createPop = function(N, D, LL, UL) {
	x = matrix(runif(N*D, LL, UL), ncol=D)
}

rosen <- function(xx) {
  d <- length(xx)
  xi <- xx[1:(d-1)]
  xnext <- xx[2:d]
	
  sum <- sum(100*(xnext-xi^2)^2 + (xi-1)^2)
	
  y <- sum
  return(y)
}



fitnessPop = function(p) {

	fp = NULL
	for (i in 1:nrow(p)) {
		fp[i] = rosen(p[i,])
	}
	return (fp)
}

select = function(p, fp, o, fo) {
	N = nrow(p)
	a = rbind(p, o)
	b = c(fp, fo)
	i = order(b)
	a = a[i,]
	a = a[1:N,]
	fa = b[i][1:N]
	return(list(a=a, fa=fa))
}

printBest = function(p, fp, g) {
	i = which.min(fp)
	fx = sprintf("% -6.2e", fp[i])
	x = sprintf("% -6.2e", p[i,])
	x = paste(x, collapse=" ")
	print(paste("G:", g, "  Best: ", fx, " x = ", x))
}

rosen2D <- function(x, y) {
  xx=cbind(x,y)
  #s = xx[,1] + xx[,2]
  s <- sum(100*(xx[,2]-xx[,1]^2)^2 + (xx[,1]-1)^2)
  return(s)
}

z = NULL

evaluateZ = function(n=50, p, fp, zoom=F, LL, UL) {

	if (is.null(z)) {
		y <<- x <<- seq(LL, UL, length.out=n)
		z <<- outer(x, y, Vectorize(rosen2D))

	}
	#par( mar=c(2,2,2,2))

	pal=colorRampPalette(c("red", "yellow", "green"), space="rgb")	
	best = p[which.min(fp), , drop=F]

	if (!zoom) {
		#d = c(0, unique(sort(z)))
		#b = d[c(1:15, round(seq(200, length(d), length.out=10)))]
		filled.contour(x,y, log(z), las=0, col=pal(20), #levels=b, 
			plot.axes={
				points(p, pch=16, col='gray60')
				points(best, pch=16, col='blue', cex=1.3)
				points(1,1, pch=1, col='yellow', cex=1.3)
				points(1,1, pch=1, col='yellow', cex=1.4)
				axis(1,cex.axis=2)
				axis(2,cex.axis=2)
			}
		)
	} else {
		d = c(0, unique(sort(z)))
		b = d[c(1,2,3,4,5, 10, round(seq(20, length(d), length.out=10)))]
		x = seq(min(p[,1]), max(p[,1]), length.out=n)
		y = seq(min(p[,2]), max(p[,2]), length.out=n)
		z <<- outer(x, y, Vectorize(rosen2D))
		filled.contour(x,y, log(z), las=0,  #levels=b, 
			plot.axes={
				points(p, pch=16, col='gray60')
				points(best, pch=16, col='blue', cex=1.3)
				points(1,1, pch=1, col='yellow', cex=1.3)
				points(1,1, pch=1, col='yellow', cex=1.4)
				axis(1,cex.axis=2)
				axis(2,cex.axis=2)
			}
		)
	}
}


DE = function(N = 100, D = 2, FF = 0.6, CR = 0.8, LL = -5, UL = 5, G = 200, plotZ=F) {
	
	pl = list()
	fl = list()

	p = createPop(N, D, LL, UL)
	o = createPop(N, D, LL, UL)
	fp = fitnessPop(p)
	g = 1
	if (plotZ) evaluateZ(n=350, p, fp, F, LL, UL)
	pl[[g]] = p
	fl[[g]] = fp
	while (g < G) {
	
		for (i in 1:N) {
			R = sample(N, 3)
			iRand = sample(D, 1)
			for (j in 1:D) {
				if ((runif(1) <= CR) | (j == iRand)) {
					o[i,j] = p[R[1], j] + FF * (p[R[2], j] - p[R[3], j])
				} else {
					o[i,j] = p[i, j]
				}			
			}
		}
		fo = fitnessPop(o)
		s = select(p, fp, o, fo)
		p = s$a
		fp = s$fa
		printBest(p, fp, g)
		g = g+1
		
		pl[[g]] = p
		fl[[g]] = fp		
		
		if (plotZ) evaluateZ(n=350, p, fp, F, LL, UL)
		
		if (min(fp) == 0) break;
	}
	
	mean = sapply(fl, mean)
	sd   = sapply(fl, sd)
	min  = sapply(fl, min)
	max  = sapply(fl, max)
	
	x = 1:length(mean)
	df = rbind(data.frame(x=x, mean=mean, min=min, max=max, sd=sd))
	return (df)

	
}

printAnalysis = function() {

	source("C:/Users/Daniel/Documents/ICMC/ssc0713/DE/de.r")
	
	DE(N = 20, D = 2, FF = 0.6, CR = 0.8, LL = -5, UL = 5, G = 200, plotZ=T)
	
	
	# Population Size
	vp = c(10, 20, 50, 100, 200)
	df = NULL
	for (i in 1:length(vp)) {
		o = DE(N = vp[i], D = 2, FF = 0.6, CR = 0.8, LL = -5, UL = 5, G = 200, plotZ=F)
		df = rbind(df, cbind.data.frame(o, vp=vp[i]))
	}
	
	require(ggplot2)
	ggplot(df, aes(x=x, y=sd, colour=as.factor(vp))) + geom_line(size=1.5) + ylim(0, 0.1)
	

	# Crossover rate
	vF = c(0.1, 0.5, 1.0, 1.5)
	df = NULL
	for (i in 1:length(vF)) {
		o = DE(N = 100, D = 2, FF = vF[i], CR = 0.8, LL = -5, UL = 5, G = 200, plotZ=F)
		df = rbind(df, cbind.data.frame(o, vp=vF[i]))
	}
	
	require(ggplot2)
	ggplot(df, aes(x=x, y=min, colour=as.factor(vp))) + geom_line(size=1.5) + ylim(0, 2)
	
	
	# test t
	df = NULL
	x = NULL
	for (i in 1:30) {
		o = DE(N = 100, D = 4, FF = 0.5, CR = 0.8, LL = -5, UL = 5, G = 100, plotZ=F)
		df = rbind(df, cbind.data.frame(o, vp=i))
		j = nrow(o)
		x = rbind(x, cbind.data.frame(mean=o$mean[j], min=o$min[j], max=o$max[j], sd=o$sd[j]))
	}
	
	df2 = NULL
	y = NULL
	for (i in 1:30) {
		o = DE(N = 100, D = 4, FF = 0.1, CR = 0.8, LL = -5, UL = 5, G = 100, plotZ=F)
		df2 = rbind(df, cbind.data.frame(o, vp=i))
		j = nrow(o)
		y = rbind(y, cbind.data.frame(mean=o$mean[j], min=o$min[j], max=o$max[j], sd=o$sd[j]))
	}
	
	require(nortest)
	ad.test(runif(1000))
	ad.test(x$mean)
	ad.test(y$mean)
	
	t.test(x$mean, y$mean)
	t.test(x$min, y$min)

	qplot(mean, data=df3, geom="density", fill=alg, alpha=I(.5))  + scale_colour_discrete(drop = FALSE) 

	wilcox.test(x$mean, y$mean)
	wilcox.test(x$min, y$min)
	
	df3 = rbind(cbind.data.frame(x, alg='x'), cbind.data.frame(y, alg='y'))
	ggplot(df3, aes(x=alg, y=min, fill=alg)) + geom_boxplot()

	
	hist(df3$mean)
	
}