library(DEoptim)
library(rgenoud)
library(snow)
source("CirclePacking.R")

# Number of circles and dimensions of bounding rectangle 
N=100; W=10; H=10
root="N64_100"
arrnum=64
arrdim=ceiling(sqrt(arrnum))

# Create objective function
obj=createObj(N,W,H)

# Lowest acceptable values for dimensions (x,y) is zero
# Minimum radius (to prevent circles being optimised out of existance)
low=rep(c(0,0,0.01*min(W,H)),N)
# x can go as high as W, y as high as H
# Radius of biggest possible circle contained in box is the smaller of W/2 and H/2
up=rep(c(W,H,min(W,H)/2),N)

# Generate initial population of guesses
NumPart=10*3*N
pop=matrix(0,nrow=NumPart,ncol=3*N)
for(i in 1:NumPart) pop[i,]=genGuess(N,W,H)

# L-BFGS-B test
results=c()
times=c()
pdf(paste("L-BFGS-B",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	z=genGuess(N,W,H); #plotCircles(z,W,H,numbers=TRUE) # Starting guess
	tim=system.time({
		out=optim(par=z,fn=obj,method="L-BFGS-B",lower=low,upper=up,control=list(maxit=300)) # Optimise
		z=as.numeric(out$par);  results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()

# parallel genoud Differential optimisation test
doms=matrix(c(low,up),nrow=3*N,ncol=2,byrow=FALSE) # Format bounds
results=c()
times=c()
cl <- makeCluster(12, type = "SOCK") # Initiate cluster
# Unfortunately, snow clusters don't seem to respect local R variables, so we pass them as global variables to each node, like this...
clusterExport(cl,list("N","W","H"))
pdf(paste("Parallelrgenoud",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	tim=system.time({
		z=genGuess(N,W,H); #plotCircles(z,W,H,numbers=TRUE) # Starting guess
		out=genoud(obj,3*N,starting.values=z,boundary.enforcement=2,Domains=doms, cluster=cl, max.generations=10,print.level=0,control=list(maxit=300))# Optimise
		z=out$par; results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
stopCluster(cl)
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()

# genoud Differential optimisation test
doms=matrix(c(low,up),nrow=3*N,ncol=2,byrow=FALSE) # Format bounds
results=c()
times=c()
pdf(paste("rgenoud",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	tim=system.time({
		z=genGuess(N,W,H); #plotCircles(z,W,H,numbers=TRUE) # Starting guess
		out=genoud(obj,3*N,starting.values=z,boundary.enforcement=2,Domains=doms, max.generations=10,print.level=0,control=list(maxit=300))# Optimise
		z=out$par; results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()

# DEoptim Differential optimisation test
results=c()
times=c()
pdf(paste("DEoptim",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	pop=matrix(0,nrow=NumPart,ncol=3*N)
	for(i in 1:NumPart) pop[i,]=genGuess(N,W,H); # Starting guesses
	tim=system.time({
		out=DEoptim(obj, low, up, DEoptim.control(trace=0,itermax=2000,NP=NumPart,initialpop=pop)) # Optimise
		z=as.numeric(out$optim$bestmem);  results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()