Honesty about vaccination


I am not medically qualified or a statistician or epidemiologist so this is not expert opinion, but I cannot help but feel that there has been a lack of honesty in the public discourse about what vaccination against covid-19 will achieve.

First of all, I think the effort to find a vaccination will succeed and in doing so will be the greatest scientific achievement of humanity since Armstrong walked on the Moon. By next Spring I expect to see large numbers of people being vaccinated.

But I also think the following are quite likely to be true:

  • Vaccination – certainly no vaccine available in 2021 – will not give guaranteed immunity to the disease, but will probably either give some people protection and for others will lower the severity of the disease.
  • The first vaccinations in the field won’t be the best achievable and the cycle of research and improvement will go for years.
  • The vaccination may only be really effective for a short time – perhaps as short as six months.
  • All of this means that we are going to have to live with many of the restrictions on our lives for years to come – hopefully there will be a rapid de-escalation throughout 2021 but there will still be a base line that will quite likely require social distancing in many indoor venues.
  • Lowering the infection rate will mean “test and trace” finally starts to work.

The other point to note is that Covid-19 is the third novel corona virus of the 21st century – we were luckier with the other two, SARS and MERS, but unless we act to minimise the risks we should expect a fourth such disease in the next decade if not sooner.

Programmatic adventures in the hyperbolic plane


Recursion is dangerous. So dangerous that in my professional life I’m explicitly banned from using it in code. After all who would want their aeroplane to stall, their car to crash or their fridge to defrost because of a stack overflow.

But it is also, most (but not all) programmers would agree, beautiful and more than a bit mysterious. Indeed one of the biggest selling popular science/mathematics books of all time centres on recursion as a possible partial explanation of human consciousness.

In the last week or so I’ve been attempting to use recursion in R to draw a tessellating pattern of squares on the conformal representation of the hyperbolic plane – an idea from M.C. Escher via Nobel laureate Sir Roger Penrose.

The code here is where I’ve got to.

#!/usr/bin/env Rscript

library(ggplot2)
library(ggforce)

basicLengths<-seq(from<-0, to<-1.0, by<-1/250000)
projLengths<-basicLengths * (2 /(1 + basicLengths * basicLengths))

getBestMatch<-function(numIn)
{
  repVec<-basicLengths[which.min(abs(numIn - projLengths))]
  return(repVec)
}


drawConformalLine <-function(startX, startY, endX, endY)
{
  #Length and angle of line to starting conformal point
  lengthToStart <- sqrt(startX * startX + startY * startY)
  angleToStart <- atan2(startY, startX)
  #And to end Point
  lengthToEnd <- sqrt(endX * endX + endY * endY)
  angleToEnd <- atan2(endY, endX)
  #Project
  newStartLength<-lengthToStart * (2 / (1 + lengthToStart * lengthToStart))
  newEndLength<-lengthToEnd * (2 / (1 + lengthToEnd * lengthToEnd))
  newStartX<-newStartLength * cos(angleToStart)
  newStartY<-newStartLength * sin(angleToStart)
  newEndX<-newEndLength * cos(angleToEnd)
  newEndY<-newEndLength * sin(angleToEnd)

  xpoints<-seq(from<-newStartX, to<-newEndX, by<-(newEndX - newStartX)/50)
  ypoints<-seq(from<-newStartY, to<-newEndY, by<-(newEndY - newStartY)/50)
  euclidLine<-data.frame(xpoints, ypoints)
  #Conformal line (reverse projection)
  lengths<-sqrt(euclidLine$xpoints * euclidLine$xpoints + euclidLine$ypoints * euclidLine$ypoints)
  newLengths<-vector()
  for (i in 1:length(lengths))
  {
    newLengths<-c(newLengths, getBestMatch(lengths[i]))
  }
  angles<-atan2(euclidLine$ypoints, euclidLine$xpoints)
  newPoints<-data.frame(xp = cos(angles) * newLengths, yp = sin(angles) * newLengths)
  return(newPoints)
}

recursiveSquares<-function(setOfPoints, base, radius, sign)
{
  if (base + radius > 0.75){
    return (setOfPoints)
  }
  setOfPoints<-rbind(setOfPoints, recursiveSquares(setOfPoints, base + radius, radius/2, sign))
  positiveCorner = (sign * base + sign * radius)
  negativeCorner = (sign * base)
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, positiveCorner, positiveCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, negativeCorner, negativeCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, negativeCorner, negativeCorner, positiveCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, positiveCorner, positiveCorner, positiveCorner))
  return (setOfPoints)
}


partitions<-100
distances<-seq(from=1, to=partitions, by=1)
distances<-exp(distances/partitions)/(2 * distances)
circles<-data.frame(x0=rep(0, partitions), y0=rep(0, partitions), r=distances)
t <- ggplot() + geom_circle(aes(x0=x0, y0=y0, r=1-r, color=rgb(r,1-r,1)), data=circles)
t <- t + geom_circle(aes(x0=0, y0=0, r=1), color="black", data=data.frame(x0=0, y0=0, r=1)) + coord_fixed()
distances<- 1-distances
n_distances<-(distances * 2)/(1 + distances * distances)
n_circles<-data.frame(x0=rep(0, partitions), y0=rep(0, partitions), r=n_distances)
y <- ggplot() + geom_circle(aes(x0=x0, y0=y0, r=r, color=rgb(r,1-r,1)), data=n_circles)
y <- y + geom_circle(aes(x0=0, y0=0, r=1), color="black", data=data.frame(x0=0, y0=0, r=1)) + coord_fixed()
square<-0.4
squarePoints<-data.frame()
squarePoints<-recursiveSquares(squarePoints, 0, square, 1)
nPoints<-data.frame()
nPoints<-recursiveSquares(nPoints, 0, square, -1)
squarePoints<-rbind(squarePoints, nPoints)
t<-t + geom_point(aes(color="magenta", x=squarePoints$xp, y=squarePoints$yp))
t<-t + geom_point(aes(color="magenta", x = -squarePoints$xp, y=squarePoints$yp)) 


I had several problems to solve, and recursion seemed a good fit for some of them, though in the end I’ve only been able to make partial progress.

Here are some of the problems and how I tackled them:

  1. Finding the straight lines

In the conformal representation of the hyperbolic plane ‘straight lines’ – in the sense of the shortest routes between two points – are actually represented by the segment of the circle joins the two points and that meets the boundary of the representation of the plane orthogonally (ie., at right angles).

For points that lie on the same polar angle in the circle that’s easy – the ‘straight line’ that joins them is simply the Euclidean line that passes through both of them and the origin (centre of the circular representation of the plane). The circle they are linked by in the conformal representation is of infinite radius.

But what about other, non-degenerate, cases? Perhaps there is an algorithmic method to work this out? But if there is I don’t know what it is and I couldn’t find it online either!

What I do know (again thanks to Sir Roger) is that is we substitute the projective representation for the conformal representation we can join the points by Euclidean lines. To go from the conformal to the projective representation of the hyperbolic plane by preserving the polar angle of points but magnifying the distance from the origin by \frac{2R^2}{R^2+r^2_c}, where R is the radius of the bounding circle (so R=1 in our case)and r_c is the length of the of the Euclidean line from the centre (origin) to the point.

The great advantage of the projective representation is that straight lines between points are just what we think of them as being in Euclidean geometry – in other words both the shortest length between two points and at a constant angle.

So to get the points on the line in the conformal representation the thing to do seemed to be to convert the points to the projective representation, calculate the straight line (very easy) and shrink everything back by the factor \frac{1+r^2_c}{2} .

But the problem with that is, that with the exception of the start and end points, we don’t know what r^2_c is to begin with – indeed that is more or less what we are trying to calculate.

My answer to that is to generate lookup tables (relatively quickly done using R’s vectorisation) and a function that indexes one to the other:

basicLengths<-seq(from<-0, to<-1.0, by<-1/250000)
projLengths<-basicLengths * (2 /(1 + basicLengths * basicLengths))

getBestMatch<-function(numIn)
{
  repVec<-basicLengths[which.min(abs(numIn - projLengths))]
  return(repVec)
}

2. Generate the squares

This is where the recursion comes in – I settled on the idea of four squares touching at each vertex, meaning each square as we travelled out from the centre the Euclidean lengths halve.

Here is (a simplified) version of the code:

recursiveSquares<-function(setOfPoints, base, radius)
{
  if (base + radius > 0.75){
    return (setOfPoints)
  }
  setOfPoints<-rbind(setOfPoints, recursiveSquares(setOfPoints, base + radius, radius/2, sign))
  positiveCorner = (base + radius)
  negativeCorner = (base)
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, positiveCorner, positiveCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, negativeCorner, negativeCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, negativeCorner, negativeCorner, positiveCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, positiveCorner, positiveCorner, positiveCorner))
  return (setOfPoints)
}

Those of you who have read the classic Structure and Interpretation of Computer Programs may immediately recognise a weakness here – by not using tail recursion I am putting more stress on the stack.

And it’s true that an alternative tail recursive form (this is the full code) works:

recursiveSquares<-function(setOfPoints, base, radius, sign)
{
  if (base + radius > 0.75){
    return (setOfPoints)
  }
  positiveCorner = (sign * base + sign * radius)
  negativeCorner = (sign * base)
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, positiveCorner, positiveCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, negativeCorner, negativeCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, negativeCorner, negativeCorner, positiveCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, positiveCorner, positiveCorner, positiveCorner))
  recursiveSquares(setOfPoints, base + radius, radius/2, sign)
}

But it still is liable to bust the stack – indeed the whole system quite finely balanced. Make the initial squares too small and it is easy to break the stack – though the tail recursive form does seem a little more robust.

3. Getting a tessellation/exploiting the dihedral symmetry

Unfortunately the tendency of the stack to blow up makes it very difficult (as far as I can see anyway) to write some elegant code that will seamlessly deliver the tessellations sought. But some simple group theory can help us get some more patterns.

The pattern above shows rotational and reflective symmetry – in other words the shape has the symmetries of square (dihedral group D_4 ).

We use this in the code through the <code>sign</code> factor in the calls to <code>recursiveSquares</code> and then in the two lines:

t<-t + geom_point(aes(color="magenta", x=squarePoints$xp, y=squarePoints$yp))
t<-t + geom_point(aes(color="magenta", x = -squarePoints$xp, y=squarePoints$yp)) 

These effectively reflect earlier results.

Representing equal distances in the conformal representation of the hyperbolic plane


(This is mainly about posting a pretty picture.)

As I have previously remarked, the last couple of times I have tried to read Roger Penrose’s The Road to Reality I have stumbled over the discussion of non-Euclidean geometry in Chapter 2 (of a 34 chapter book).

Now I feel more confident about the maths but I am also taking more time and exploring the issues – after all I’ve had the book for more than a decade so there has been no rush so far.

In Chapter 2 Penrose presents the conformal (angles are preserved) representation of a hyperbolic plane through one of Escher’s woodcuts – my picture (generated with some R code) dispenses with Escher’s fish (see this blog post from Allyson Redhunt for a further discussion) and instead just looks at the way equal distances are seen in the projection.

In this representation the black line at the edge is at infinity, but all the other circles are equidistant – I’ve applied a coloring effect so as we approach the edge they don’t just look like a solid block.

The R code is shown below.

#!/usr/bin/env Rscript
library(ggplot2)
library(ggforce)
partitions<-1000
distances<-seq(from=1, to=partitions, by=1) 
distances<-exp(distances/partitions)/(2 * distances)
circles<-data.frame(x0=rep(0, partitions), y0=rep(0, partitions), r=distances)
t <- ggplot() + geom_circle(aes(x0=x0, y0=y0, r=1-r, color=rgb(r,1-r,1)), data=circles)
t <- t + geom_circle(aes(x0=0, y0=0, r=1), color="black", data=data.frame(x0=0, y0=0, r=1)) + coord_fixed()

In praise of Roger Penrose


Roger Penrose has been awarded a share in the Nobel Prize for physics and I could not be more pleased.

It is not that I have met him or even attended a lecture by him and nor do we even see him much on TV – but I owe him a debt for his insight and also for his ability to educate and inform.

A few years ago I bought his “Fashion, Faith and Fantasy” – his stern critique of what he sees as the fashion of string theory and assorted similar ideas. I’m not going to pretend it’s an easy book, and much of it was well over my head – but it is still choc-full of insight. From it I went back to an older copy of “Cycles of Time”. This is a shorter book and was much more heavily marketed as for the “lay reader” but, actually, I found much of it much harder to get to grips with – if you don’t really understand conformal geometry (and I didn’t) it’s pretty hard to follow. But, and this is a big but, it has an absolutely brilliant explanation of entropy in its opening chapters and if, like me, you never really got to grips with this subject (because, I maintain, it was so poorly explained) as an undergraduate, this really cuts through the fog.

It mattered to me because its discussion of phase spaces gave me an insight into the stochastic nature of memory management systems and the way in which entropy can help us explain latency in computer systems. I wrote quite a bit about that in my PhD thesis and I am convinced it’s a subject that would repay a lot more investigation. Sir Roger collected an additional citation for that (not that he needed it of course).

Only last night I again made an attempt on the massive “The Road to Reality” – previously I’ve never managed to get past the introductory discussion of non-Euclidean geometry in chapter two, but I feel a bit more confident about that now, so giving it another go.

Overflow in unsigned arithmetic in C


This came up in work recently – and I got a bit confused about it, so I thought it was worth writing it down in case someone like me needs to look it up.

Let’s consider a real world example – we have a 16 bit counter that increments on every event X and which we sample periodically: at our last sampling was at value 0xFFFD and at our current sample is 0x03. How many incidents of X have occurred?

Ignoring the possibility of a ‘double rollover’ or anything similar, the answer is 6 (don’t forget the counter goes through 0). But how can we calculate this programmatically?

It turns out it’s very simple and this always works:

uint16_t numberEvents = currentSample - previousSample

In C the calculation (constrained to 16 bits) of 0x03 – 0xFFFD with unsigned arithmetic results in an overflow (we cannot calculate a negative number) but the language standard gives us a well-defined answer that works well.

The answer of -65530 is reinterpreted as modulo 0x10000 (65536) and this comes out at 6. To explain why think of a smaller number system – modulo 5 say. Then the number line (integers and integers mod 5) looks like this:

Integers: 10 9 8 7 6 5 4 3 2 1 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10
Modulo 5: 0 4 3 2 1 0 4 3 2 1 0 4 3 2 1 0 4 3 2 1 0

Under this system we can see that 1 – 4 = -3 for integers but 2 mod 5 – which if this was a simple two bit counter would be the correct number of increments that had taken place.