suppressMessages(library(GA))
library(data.tree)

dynVisc = 0.0000181 
density = 1.225 
roughness = 0.01

tree <- CreateRegularTree(4, 2)
tree$Set(flowRate = 0.1, filterFun = isLeaf)
tree$Do(function(x) x$flowRate <- Aggregate(x, "flowRate", sum))

bitnum=7
nvars=tree$totalCount

string = sample(rep(0:1,each=(bitnum*nvars)/2))

fitness = function(string) {
        x=split(string, ceiling(seq_along(string)/bitnum))
	y=t(as.data.frame(x))
	ductDia=vector()
	for(i in 1:nvars) ductDia[i] <- 150+(2*(((binary2decimal(y[i,]))-1)))
	tree$Set(ductDia = ductDia)
	tree$Do(function(x) x$velocity = x$flowRate/(pi*((((x$ductDia)/1000)/2)^2)))
	tree$Do(function(x) x$Re = (density * x$velocity * ((x$ductDia)/1000))/dynVisc)
	tree$Do(function(x) x$ff = 1/((-1.8*log10((6.9/x$Re)+(((roughness/x$ductDia)/3.71)^1.11)))^2))
	tree$Do(function(x) x$Pa =  x$ff*(1/(x$ductDia/1000))*(1/2)*density*((x$velocity)^2) )
	Pas = tree$Get("Pa")
	fitness = 1/(((sum((Pas-1)^2))/(length(Pas)-1))^(1/2))
	}	

GA <- ga(type = "binary", nBits = bitnum*nvars, maxFitness=20000, pmutation = 0.6, fitness = fitness,popSize = 100000, maxiter = 10, run = 1, monitor=TRUE)

#summary(GA)
#
string = as.vector(GA@solution)
x=split(string, ceiling(seq_along(string)/bitnum))
y=t(as.data.frame(x))
ductDia2=vector()
for(i in 1:nvars) ductDia2[i] <- 150+(2*((binary2decimal(y[i,]))-1))

tree$Set(ductDia2 = ductDia2)
tree$Do(function(x) x$velocity2 = x$flowRate/(pi*((((x$ductDia2)/1000)/2)^2)))
tree$Do(function(x) x$Re2 = (density * x$velocity2 * ((x$ductDia2)/1000))/dynVisc)
tree$Do(function(x) x$ff2 = 1/((-1.8*log10((6.9/x$Re2)+(((roughness/x$ductDia2)/3.71)^1.11)))^2))
tree$Do(function(x) x$Pa2 =  x$ff2*(1/(x$ductDia2/1000))*(1/2)*density*((x$velocity2)^2) )

print(tree, "flowRate", "ductDia2", "velocity2", "Pa2")