Skip to content

Instantly share code, notes, and snippets.

@CnrLwlss
Last active December 4, 2016 15:47
Show Gist options
  • Save CnrLwlss/4571810 to your computer and use it in GitHub Desktop.
Save CnrLwlss/4571810 to your computer and use it in GitHub Desktop.
Functions for packing N circles into a rectangle of width W and height H, together with a function for plotting solution and some example code fitting 13 circles into a square. http://cnr.lwlss.net/CircleObjectivesR/
###########
# Optimum circle packing in R
###########
# Function closure, sets up list of possible combinations
# and returns objective function
createObj<-function(N,W,H){
print("Generating combinations...")
# Generate all possible circle pairings
cpairs=data.frame(t(combn(1:N,2)))
colnames(cpairs)=c("A","B")
print("Combinations complete.")
return({
newObj<-function(z){
# This vector is split up into triplets of x,y,r
# z[1],z[2],z[3] -> x1,y1,r1
x=z[1:N*3-2]
y=z[1:N*3-1]
r=z[1:N*3-0]
res=0
# Some linear inequality constraints to be satisfied
c1=x+r-W # <0
c2=y+r-H # <0
c3=r-x # <= 0
c4=r-y # <=0
# Many nonlinear inequality constraints
c5=r[cpairs$A]+r[cpairs$B]-sqrt((x[cpairs$A]-x[cpairs$B])^2+(y[cpairs$A]-y[cpairs$B])^2) # <=0
c1=c1[c1>0]
c2=c2[c2>0]
c3=c3[c3>0]
c4=c4[c4>0]
c5=c5[c5>0]
constraints=c(c1,c2,c3,c4,c5)
if(length(constraints)>0){
constraints=constraints^2
res=res-sum(constraints)
}
# Actual objective function (fraction of area covered)
res=res+(pi*sum(r^2))/(W*H)
return(-1*res)
}
})
}
# Function closure, sets up list of possible combinations
# and returns objective function
createHardObj<-function(N,W,H){
print("Generating combinations...")
# Generate all possible circle pairings
cpairs=data.frame(t(combn(1:N,2)))
colnames(cpairs)=c("A","B")
print("Combinations complete.")
return({
newHardObj<-function(z){
# This vector is split up into triplets of x,y,r
# z[1],z[2],z[3] -> x1,y1,r1
x=z[1:N*3-2]
y=z[1:N*3-1]
r=z[1:N*3-0]
res=0
# Some nonlinear inequality constraints to be satisfied
if(FALSE%in%(x+r<W)) res = Inf
if(FALSE%in%(y+r<H)) res = Inf
if(FALSE%in%(x-r>=0)) res = Inf
if(FALSE%in%(y-r>=0)) res = Inf
# Many nonlinear inequality constraints
if(FALSE%in%(((x[cpairs$A]-x[cpairs$B])^2+(y[cpairs$A]-y[cpairs$B])^2)>=(r[cpairs$A]+r[cpairs$B])^2)) res = Inf
# Actual objective function (fraction of area covered)
res=res+(pi*sum(r^2))/(W*H)
return(-1*res)
}
})
}
# Draw circles (centres and radii speficied in df)
# within a bounding rectangle (width * height)
drawCirclesRect<-function(df,width,height,colour="black",pdfName=NULL,numbers=FALSE){
if((is.null(pdfName))){
graphics.off()
dev.new(width=7, height=7*height/width)
}else{
pdf(pdfName,width=7,height=7*height/width)
}
maxdim=max(width,height)
width=width/maxdim
height=height/maxdim
df=df/maxdim
# Change to coordinate system based on centre of rectangle
df$x=df$x-width/2
df$y=df$y-height/2
# Draws circles on a square
par(plt=c(0,1,0,1))
plot(NULL,xlim=c(-width/2,width/2),ylim=c(-height/2,height/2),axes=FALSE,xaxt='n',yaxt='n',ann=FALSE)
symbols(df$x,df$y,circles=df$r,fg=colour,bg=NA,add=TRUE,inches=FALSE)
rect(-width/2,-height/2,width/2,height/2)
if(numbers) text(df$x,df$y,1:N)
if(!is.null(pdfName)) dev.off()
return((pi*sum(df$r^2))/(width*height))
}
drawCircles<-function(df,width,height,colour="black",pdfName=NULL,numbers=FALSE){
maxdim=max(width,height)
width=width/maxdim
height=height/maxdim
df=df/maxdim
# Change to coordinate system based on centre of rectangle
df$x=df$x-width/2
df$y=df$y-height/2
# Draws circles on a square
opc<-par(plt=c(0,1,0,1))
plot(NULL,xlim=c(-width/2,width/2),ylim=c(-height/2,height/2),axes=FALSE,xaxt='n',yaxt='n',ann=FALSE)
symbols(df$x,df$y,circles=df$r,fg=colour,bg=NA,add=TRUE,inches=FALSE)
rect(-width/2,-height/2,width/2,height/2)
if(numbers) text(df$x,df$y,1:N)
par(opc)
return((pi*sum(df$r^2))/(width*height))
}
# Generate sensible initial guess (random circle coordinates and radii)
genGuess<-function(N,W,H){
z=rep(0,3*N)
z[1:N*3-2]=runif(N,0,W)
z[1:N*3-1]=runif(N,0,H)
z[1:N*3-0]=runif(N,0.01*min(W,H),0.1*min(W,H))
return(z)
}
# Draw circles from vector z
plotCircles<-function(z,W,H,colour="black",pdfName=NULL,numbers=FALSE){
N=length(z)/3
x=z[1:N*3-2]; y=z[1:N*3-1]; r=abs(z[1:N*3-0])
df=data.frame(x=x,y=y,r=r)
drawCircles(df,W,H,colour,pdfName,numbers)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Demonstration
# Number of circles and dimensions of bounding rectangle
N=13; W=10; H=10
# Create objective function
obj=createObj(N,W,H)
# Lowest acceptable values for dimensions or radii (x,y,r) is zero
low=rep(0,3*N)
# x can go as high as W, y as high as H
# The biggest allowable circle radius is the smaller of W and H
up=rep(c(W,H,min(W,H)),N)
# Initial guess (random circle coordinates and radii)
z=genGuess(N,W,H)
# Optimium circle packing
out=optim(par=z,fn=obj,method="L-BFGS-B",lower=low,upper=up,control=list(maxit=2000))
print(c("Fractional area:",-1*out$value))
plotCircles(out$par,W,H)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment