For some dissemination work, I want to create a nice graph to explain the exponential growth in pandemics, related to the value of R_0. Recall that R_0 corresponds to the average number of people that a contagious person can infect. Hence, with R_0=1.5, 4 people will contaminate 6 people, and those 6 will contaminate 9, etc. After n iteration, the number of contaminated people is simply R_0{}^n. As explained by Daniel Kahneman

people, certainly including myself, don’t seem to be able to think straight about exponential growth. What we see today are infections that occurred 2 or 3 weeks ago and the deaths today are people who got infected 4 or 5 weeks ago. All of this is I think beyond intuitive human comprehension

For different values of R_0 (on each row), I wanted to visualise the number of contaminated people after 3, 5 or 7 iterations, since graphs are usually the most simple way to give some intuition. The graph I had in mind was the following

(to be honest, I am quite sure I had seen it somewhere, but I cannot find where). The main challenge here is pack *optimally* k identical circles intro a unit circle: we need here the location of the points (center of the disks) and the radius. It seems to be a rather complicated mathematical problem. Nicely, on http://hydra.nat.uni-magdeburg.de/packing, it is possible to get the “best known packings of equal circles in a circle” (up to 5000, but many k‘s are missing). For instance, for k=37, we have

And interestingly, on the same website, we can get the coordinates of the centers, for example with 37 disks, so it is possible to recreate the R graph.

k = 37
base = read.table(paste("http://hydra.nat.uni-magdeburg.de/packing/cci/txt/cci",k,".txt",sep=""), header=FALSE) |

k = 37
base = read.table(paste("http://hydra.nat.uni-magdeburg.de/packing/cci/txt/cci",k,".txt",sep=""), header=FALSE)

The problem, as discussed earlier, is that some cases are not solved, yes, for instance k=2^{12}=4096: the next feasable case is 4105. To avoid that issue, one can use

T = "Error"
while(T == "Error"){
T = substr(try(base = read.table(paste("http://hydra.nat.uni-magdeburg.de/packing/cci/txt/cci",k,".txt",sep=""), header=FALSE),silent = TRUE),1,5)
k=k+1
}
k=k-1 |

T = "Error"
while(T == "Error"){
T = substr(try(base = read.table(paste("http://hydra.nat.uni-magdeburg.de/packing/cci/txt/cci",k,".txt",sep=""), header=FALSE),silent = TRUE),1,5)
k=k+1
}
k=k-1

Now we can almost plot it. The problem is that the radius of the circles is missing, here. But we can compute it

D=as.matrix(dist(x = base[,2:3]))
diag(D)=1e5
i=which(D == min(D), arr.ind = TRUE)
r = D[i[1,1],i[1,2]] |

D=as.matrix(dist(x = base[,2:3]))
diag(D)=1e5
i=which(D == min(D), arr.ind = TRUE)
r = D[i[1,1],i[1,2]]

Here the radius is

To plot it, use

plot(base$V2,base$V3,xlim=c(-1,1),ylim=c(-1,1))
n=100
theta=seq(0,pi,length=n+1)
circ= function(x,y,r,h=1){
vu=x+r*cos(theta)
vv=r*sin(theta)
cbind(c(vu,rev(vu))*h,c(y+vv,y-rev(vv))*h)
}
for(i in 1:k) polygon(circ(base[i,2],base[i,3],r/2*.95),col=colr,border=NA) |

plot(base$V2,base$V3,xlim=c(-1,1),ylim=c(-1,1))
n=100
theta=seq(0,pi,length=n+1)
circ= function(x,y,r,h=1){
vu=x+r*cos(theta)
vv=r*sin(theta)
cbind(c(vu,rev(vu))*h,c(y+vv,y-rev(vv))*h)
}
for(i in 1:k) polygon(circ(base[i,2],base[i,3],r/2*.95),col=colr,border=NA)

We can now use that code to create the graph above, with k=R_0{}^n for various values of n

And we can also use it to visualize more subtile differences, like R_0=1.1, R_0=1.3, R_0=1.5 and R_0=1.7