## Sunday, 31 December 2017

### film review: Rogue One

Now, that’s what a prequel should be like!

As a young girl, Jyn Erso sees her mother shot and her engineer father taken by Imperial troopers. A decade later she is an embittered young criminal, on the run from everyone. The Empire wants her to put pressure on her father. The Rebellion wants her to stop the Empire getting her. Then a pilot defects from the Empire with news that could mean a new hope for the Rebellion. Jyn Erso, along with a motley crew of rebels, finds herself engaged in a desperate race against time.

This is written for, and by, Star Wars original trilogy fans. There are lots of little call outs to events in those films, and a couple of space battles in the original style. The plot explains two big questions from the original: How did Leia get hold of the Death Star plans? Why did the Death Star have such a ridiculous weakness? And it manages this without introducing too many new questions, although it does introduce a new series of gloomy muddy planets. Nowadays it is a bit harder to fully engage with an armed insurrection quite as whole-heartedly as in the more innocent days of 1977, though.

One great piece of continuity is Peter Cushing reprising his role as Grand Moff Tarkin from the 40-year old original. Not a bad trick given that Peter Cushing died in 1994. That’s what CGI can do nowadays: take one actor, and plaster the face of another over the top. The face worked brilliantly, but the eyes seemed a bit glassy. They played the same trick at the end with Leia, but that was more uncanny valley territory: a grizzled old face in the gloom is easier to fake that a young unlined face in full light, it seems. Mom Mothma also recurs, recast with a good look-alike. Darth Vader was a bit easier to recast behind that black mask, but still has his great James Earl Jones voice.

Speaking of Leia and Mom Mothma, does this film, with its female lead, pass the Bechdel test? Opinion is divided: the young Jyn does briefly talk with her mother at the beginning, but they are talking about her hiding from the men coming for them. Later, Jyn addresses the Rebellion council headed by Mom Mothma. Are they talking to each other, or is Jyn just addressing the assembled, overwhelmingly male, crowd? Even if it can be argued a technical pass, it does needs to be argued. The crew Jyn flies with is all male. About the only other women are another council member, and a couple of fighter pilots, who get a few lines. Would it have killed them to have had another woman in the crew? Or even to have had Jyn’s mother be the abducted engineer? Ah well. Baby steps.

Technically, and plot-wise, the film is great. Not as good as the original (but what is?), but much superior to the prequel trilogy. (We probably should have watched it in the cinema to better appreciate the space battles and mushroom clouds, though.) The plot of Rogue One finishes just before that of Star Wars starts. It’s a very satisfying prequel that adds to the canon.

For all my film reviews, see my main website.

## Saturday, 30 December 2017

### books, books, everywhere you looks

10 Famous Book Hoarders

 Karl Lagerfeld: 300,000 books

We come in somewhere between number 3 and 4 on the list. More than 10k, but less than 20k (so far).

For all my social networking posts, see my Google+ page

## Thursday, 28 December 2017

### fleshly material-semiotic presences, anyone?

Continuing strange recommendations by Amazon include a book who’s description includes the phrase "dogs are fleshly material-semiotic presences in the body of technoscience". Don’t judge a book by its blurb?

For all my social networking posts, see my Google+ page

## Tuesday, 26 December 2017

### against a dark background

The Christmas jigsaw event is complete for another year.  This year’s task: the 1000-piece Highland Glen.  The puzzle was bought from a charity shop, and we discovered it is actually the 999-piece Highland Glen.  Can you spot the missing piece?  (This is an unretouched photo, but we did do something to the jigsaw before snapping it.)

If you need a clue, we slid an appropriately coloured piece of paper under the puzzle.

If you need a further clue, look at the post title.

If that doesn’t help, here’s a shot without that piece of paper.

## Monday, 25 December 2017

### Christmas done

So that’s all the Christmas Day traditions done for another year:
1. exchanging presents
2. eating turkey
3. watching The Queen
4. watching the Doctor Who Christmas special (with two regenerations for the price of one! excellent!)

### sequestering carbon, one Christmas at a time V

His‘n’hers Christmas presents:

Something to read whilst unable to move from overeating turkey.

## Thursday, 21 December 2017

### Cel-ee-brAAAAte!

Even the House Dalek is getting into the spirit of the season!

## Wednesday, 20 December 2017

### tree by night

Christmas comes but once a year.  But when it comes, it’s accompanied by trees, and lights, and twinkly stuff.

## Friday, 15 December 2017

### sequestering carbon, several books at a time LXXVIII

The latest batch.  The next batch will appear Monday week, being those that have been smuggled into the house over the last month or so.

This batch includes Inspired by Nature, and Introduction to Coalgebra, a review copy, which means I actually have to read it!  The Professional Git is for my other half, who, in his SysAdmin role, has decided we need our own Git repository on our internal network.

## Tuesday, 12 December 2017

### seeing is believing - not

Mind.  Blown.

The Remarkable “Curvature Blindness” Illusion

[h/t Danny Yee's blog]

For all my social networking posts, see my Google+ page

## Monday, 11 December 2017

### mid-winter mid-afternoon

 4:02pm, GMT
The camera never captures the true deep red intensity of these sunsets.

## Sunday, 10 December 2017

### imprints of birds

We can’t see any birds at the moment, but we have strong evidence of their existence.

 footprints in the snow

## Saturday, 9 December 2017

### Generators have state

Python’s generator expressions (discussed in two previous posts) are very useful for stream programming. But sometimes, you want some state to be preserved between calls. That is where we need full Python generators.

A generator expression is actually a special case of a generator that can be written inline. The generator expression (<expr> for i in <iter>) is shorthand for the generator:
def some_gen():
for i in <iter>:
yield <expr>
For example, (n*n for n in count(1)) can be equivalently written as:
from itertools import *

def squares():
for n in count(1):
yield n*n

print_for(squares())

1, 4, 9, 16, 25, 36, 49, 64, 81, 100, ...

The yield statement acts somewhat like a return in providing its value. But the crucial difference is what happens on the next call. An ordinary function starts off again from the top; a generator starts again directly after the previous yield. This becomes important if the generator includes some internal state: this state is maintained between calls. That is, we can have memory of previous state on the next call. This is particularly useful for recurrence relations, where the current value is expressed in terms of previous values (kept as remembered state).

## Running total

The accumulate() generator provides a running total of its iterator argument. We can write this as an explicit sum:
$$T_N = \sum_{i=1}^N x_i$$ We can instead write this sum as a recurrence relation: $$T_0 = 0 ; T_n = T_{n-1} + x_n$$ Expanding this out explicitly we get $$T_0 = 0 ; T_1 = T_0 + x_1 = 0 + x_1 = x_1 ; T_2 = T_1 + x_2 = x_1 + x_2 ; \ldots$$ $T_0$ is the initial state, which is observed directly. The recurrence term $T_n$ tells us what needs to be remembered of the previous state(s), and how to use that to generate the current value.
def running_total(gen):
tot = 0
for x in gen:
tot += x
yield tot

print_for(repeat(3))
print_for(running_total(repeat(3)))
print_for(running_total(count(1)))

3, 3, 3, 3, 3, 3, 3, 3, 3, 3, ...
3, 6, 9, 12, 15, 18, 21, 24, 27, 30, ...
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...


## Factorial

Similarly we can write the factorial as an explicit product:$$N! = F_N = \prod_{i=1}^N i$$ And we can write it as a recurrence relation: $$F_1 = 1; F_n = n F_{n-1}$$ $F_1$ is the initial state, which should be the first output of the generator. We can yield this directly on the first call, then (on the next call) go into a loop yielding the following values.
def fact():
f = 1
yield f
for n in count(2):
f *= n
yield f

print_for(fact())

1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, ...

We can remove the special case, and instead calculate the next value after the yield in the loop. The subsequent call then picks up at that calculation.
def fact():
f = 1
for n in count(2):
yield f
f *= n

print_for(fact())

1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, ...


## Fibonacci numbers

The perennially popular Fibonacci numbers are naturally defined using a recurrence relation, involving two previous states:$$F_1 = F_2 = 1 ; F_n = F_{n-1} + F_{n-2}$$ This demonstrates how we can store more state than just the result of the previous yield.
def fib(start=(1,1)):
a,b = start
while True:
yield a
a,b = b,a+b

print_for(fib(),20)
print_for(fib((0,3)),20)

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, ...
0, 3, 3, 6, 9, 15, 24, 39, 63, 102, 165, 267, 432, 699, 1131, 1830, 2961, 4791, 7752,
12543, ...


## Logistic map

Difference equations, discrete-time analogues of differential equations, are a form of recurrence relation: the value of the state at the next time step is defined in terms of its value at previous timesteps.

The logistic map is a famous difference equation, exhibiting a range of periodic and chaotic behaviours, depending on the value of its paramenter $r$. $$x_0 \in (0,1) ; x_{n+1} = r x_n(1-x_n)$$
def logistic_map(r=4):
x = 0.1
while True:
yield x
x = r * x * (1-x)

print_for(logistic_map())
print_for(logistic_map(2))

0.1, 0.36000000000000004, 0.9216, 0.28901376000000006, 0.8219392261226498,
0.5854205387341974, 0.970813326249438, 0.11333924730376121, 0.4019738492975123,
0.9615634951138128, ...
0.1, 0.18000000000000002, 0.2952, 0.41611392, 0.4859262511644672, 0.49960385918742867,
0.49999968614491325, 0.49999999999980305, 0.5, 0.5, ...

We can plot these values to see the chaotic ($r=4$) and periodic ($r=2$) behaviours over time. (The %matplotlib inline “magic” allows the plot to display in a Jupyter notebook.)
%matplotlib inline
import matplotlib.pyplot as py

py.plot(list(islice(logistic_map(),200)))
py.plot(list(islice(logistic_map(3.5),200)))

We can also produce the more usual plot, with the parameter $r$ running along the $x$-axis, highlighting the various areas of periodicity and chaos.
from numpy import arange

start = 2.8
stop = 4
step = (stop-start)*0.002

skip = int(1/step) # no of initial vals to skip (converged to attractor)
npts = 150 # no of values to plot per value of lambda

for r in arange(start,stop,step):
yl = logistic_map(r)
list(islice(yl,skip)) # consume first few items
py.scatter(list(islice(repeat(r,npts),npts)), list(islice(yl,npts)), marker='.', s=1)

py.xlim(start,stop)
py.ylim(0,1)

(Here I have used s=1 to get a small point, rather than use a comma marker to get a pixel, because there is currently a bug in the use of pixels in scatter plots.)

## Faster $\pi$

Many series for generating $\pi$ converge very slowly.  One that converges extremely quickly is:$$\pi = \sum_{i=0}^\infty \frac{(i!)^2 \, 2^{i+1}}{(2i+1)!}$$We could use generators for each component of the term to code this us as:
import operator

def pi_term():
fact = accumulate(count(1), operator.mul)
factsq = (i*i for i in fact)

twoi1 = (2**(i+1) for i in count(1))

def fact2i1():
i,f = 3,6
while True:
yield f
f = f * (i+1) * (i+2)
i += 2

yield 2   # the i=0 term (needed because of 0! = 1 issues)
for i in map(lambda x,y,z: x*y/z, factsq, twoi1, fact2i1()):
yield i

print_for(accumulate(pi_term()),40)

2, 2.6666666666666665, 2.933333333333333, 3.0476190476190474, 3.098412698412698,
3.121500721500721, 3.132156732156732, 3.1371295371295367, 3.1394696806461506,
3.140578169680336, 3.141106021601377, 3.1413584725201353, 3.1414796489611394,
3.1415379931734746, 3.1415661593449467, 3.1415797881375944, 3.1415863960370602,
3.141589605588229, 3.1415911669915006, 3.1415919276751456, 3.1415922987403384,
3.141592479958223, 3.1415925685536337, 3.1415926119088344, 3.1415926331440347,
3.1415926435534467, 3.141592648659951, 3.14159265116678, 3.141592652398205,
3.1415926530034817, 3.1415926533011587, 3.141592653447635, 3.141592653519746,
3.1415926535552634, 3.141592653572765, 3.1415926535813923, 3.141592653585647,
3.141592653587746, 3.141592653588782, 3.141592653589293, ...

However, these terms get large, and take a long time to calculate.
import timeit

%time [ i for i in islice(pi_term(),10000) if i < 0 ];

Wall time: 10.9 s

We can instead write the terms as a recurrence relation.

Let $$F_{k-1} = \frac{((k-1)!)^2 \, 2^{k}}{(2k-1)!}$$ Then $$\begin{eqnarray} F_{k} &=& \frac{((k)!)^2 \, 2^{k+1}}{(2k+1)!} \ &=& \frac{((k-1)! k)^2 \, 2 \times 2^{k}}{(2k-1)!(2k)(2k+1)} \ &=& \frac{((k-1)!)^2 \, 2^{k} 2k^2}{(2k-1)! 2k(2k+1)} \ &=& F_{k-1} \frac{k}{2k+1} \end{eqnarray}$$ So $$F_0 = 2 ; F_{n} = F_{n-1} \frac{n}{2n+1}$$
import math

def pi_term_rec():
pt = 2
for n in count(1):
yield pt
pt = pt * n/(2*n+1)

print_for(pi_term_rec(),30)
print_for(accumulate(pi_term_rec()),30)
print(math.pi)

We can see how quickly the terms shrink.
2, 0.6666666666666666, 0.26666666666666666, 0.1142857142857143, 0.0507936507936508,
0.02308802308802309, 0.010656010656010658, 0.004972804972804974, 0.002340143516614105,
0.0011084890341856288, 0.0005278519210407756, 0.0002524509187586318,
0.00012117644100414327, 5.8344212335328244e-05, 2.816617147222743e-05,
1.3628792647851983e-05, 6.607899465625204e-06, 3.209551169017956e-06,
1.561403271414141e-06, 7.606836450479149e-07, 3.710651927063e-07,
1.8121788481005348e-07, 8.85954103515817e-08, 4.335520081034849e-08,
2.123520039690538e-08, 1.0409411959267343e-08, 5.1065039800179424e-09,
2.5068292265542627e-09, 1.2314248832196377e-09, 6.052766375147372e-10, ...
2, 2.6666666666666665, 2.933333333333333, 3.0476190476190474, 3.098412698412698,
3.121500721500721, 3.132156732156732, 3.1371295371295367, 3.1394696806461506,
3.140578169680336, 3.141106021601377, 3.1413584725201353, 3.1414796489611394,
3.1415379931734746, 3.1415661593449467, 3.1415797881375944, 3.1415863960370602,
3.141589605588229, 3.1415911669915006, 3.1415919276751456, 3.1415922987403384,
3.141592479958223, 3.1415925685536337, 3.1415926119088344, 3.1415926331440347,
3.1415926435534467, 3.141592648659951, 3.14159265116678, 3.141592652398205,
3.1415926530034817, ...
3.141592653589793

Not only is this code simpler, it is also much faster:
%time [ i for i in islice(pi_term_rec(),10000) if i < 0 ]

Wall time: 3.98 ms


### Sieving for primes

No discussion of generators would be complete without including a prime number generator. The standard algorithm is quite straightforward:
• Generate consecutive whole numbers, and test each for divisibility.  If the current number isn’t divisible by anything, yield a new prime.
• Only test for divisibility by primes, and only up to the square root of the number being tested; this requires keeping a record of the primes found so far.
• Optimisation: treat 2 as a special case, and generate and test only the odd numbers.
def primessqrt():
primessofar = []
yield 2
for n in count(3,2): # check odd numbers only, starting with 3
sqrtn = int(math.sqrt(n))
testprimes = takewhile(lambda i: i<=sqrtn, primessofar)
isprime = all(n % p for p in testprimes) # n % p == 0 if n is divisible
if isprime:  # if new prime, add to list and yield, else continue
yield n
primessofar.append(n)

print_while(primessqrt(), 200)

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
181, 191, 193, 197, 199, ...

This algorithm is often referred to as the Sieve of Eratosthenes, but the true sieve is somewhat different in its operation.  The sieve doesn’t require any square roots or divisions: it uses addition only, moving a marker through the numbers, striking out the multiples.  Using generators we can do this lazily, using a dictionary to store which marks have struck out a particular value (for example, when we get to 15, the dictionary entry will show that it has been struck out by 3).
• Generate consecutive whole numbers, and check whether each one’s dictionary entry has any markers.
• If there are no markers, yield a new prime $p$; start a new marker to strike out multiples of $p$ (start it at $p^2$, for the same reason the previous algorithm only needs to test values up to $\sqrt n$).
• If there are markers (such as for 15), the value isn’t prime; move each marker on to the next value it strikes out (so for 15, move the 3 marker on 6 places to strike out 21: each marker is moved on by twice its value, since we are optimising by not considering even values)
• Delete the current dictionary entry (so that the dictionary grows only as the number of primes found so far, rather than as the number of values checked, speeding up dictionary access).
from collections import defaultdict
def primesieve():
yield 2
sieve = defaultdict(set) # dict of n:{divisor} elems
for n in count(3,2): # check odd numbers only
if sieve[n] : # there are divisors, so not prime
for d in sieve[n]: # move the sieve markers on
else: #set empty, so prime
yield n
# remove current dict item
del sieve[n]

print_for(primesieve(), 100)

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277,
281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
503, 509, 521, 523, 541, ...

And this “true” sieve is faster than the divisor version.
%time [ i for i in islice(primessqrt(),100000) if i < 0 ]
%time [ i for i in islice(primesieve(),100000) if i < 0 ]

Wall time: 4.69 s
Wall time: 735 ms

So full Python generators are even more fun than generator expressions!

And in fact, there’s yet more fun to be had with generators, because it is possible to send values to them on each call, too; but that’s another story for another time.