load "ff-cmaes"

int NP=10;			//Control the mesh quality
func ud=100; 		//Desired temperature in the target zone
real L=2,l=1.5,R=0.075;//L:oven length, l:oven height, R: resistors radius
real eps=1e-2,hugeval=1.e30;		//a small value governing collisions and a huge one for something else
int nbres=6;		//number of resistors (should not be changed) 
int NPres=NP*(2*pi*R);
real[int] pr(nbres+2), K(nbres+2); 
int regi=nbres, rege=nbres+1, lext=1;
K=1;	K[regi]=10;

macro Grad(u) [dx(u),dy(u)]//EOM
macro Resistor(i,cx,cy,lab)
border Res#i(t=0.,2*pi) {x=cx+R*cos(t); y=cy+R*sin(t); label=lab;}//EOM
real pena = 100;
func target = 100.;

real[int] Xmin(3*nbres),Xmax(3*nbres);
for(int i=0;i<nbres;++i)
{
	Xmin[3*i] = -hugeval;
	Xmax[3*i] = hugeval;
	Xmin[3*i+1] = (i%3)*(L/3) +  R+eps/2;
	Xmax[3*i+1] = (i%3 + 1)*(L/3) - (R+eps/2);
	Xmin[3*i+2] = R+eps + (i>2 ? 2*l/3 : 0);
	Xmax[3*i+2] = (i>2 ? l : l/3) - (R+eps);
}
cout << "Xmin: " << Xmin << endl << "Xmax: " << Xmax << endl;

int iter=0;
func real J(real[int] &X)
{
	verbosity=0;
	real[int] XX=X,diff=X;
	for(int i=0;i<X.n;++i)
	{
		XX[i] = X[i]>Xmin[i] ? X[i] : Xmin[i];
		XX[i] = XX[i]<Xmax[i] ? XX[i] : Xmax[i];
	}
	diff -= XX;
	real corr = diff.l2;
	//cout << X << endl << XX << endl << "==>" << corr << endl;
	
	border b1(t=0,1){x=L*t; y=0; label=0;}
	border b2(t=0,1){x=L;   y=l*t; label=0;}
	border b3(t=0,1){x=L*(1-t); y=l; label=0;}
	border b4(t=0,1){x=0;   y=l*(1-t); label=0;}
	Resistor(1,XX[1] ,XX[2] ,1);
	Resistor(2,XX[4] ,XX[5] ,1);
	Resistor(3,XX[7] ,XX[8] ,1);
	Resistor(4,XX[10],XX[11],1);
	Resistor(5,XX[13],XX[14],1);
	Resistor(6,XX[16],XX[17],1);
	border bi1(t=0,2){x=(0.5+t)*L/3; y=l/3; label=1;}
	border bi2(t=0,1){x=5*L/6; y=(1+t)*l/3; label=1;}
	border bi3(t=0,2){x=(2.5-t)*L/3; y=2*l/3; label=1;}
	border bi4(t=0,1){x=L/6; y=(2-t)*l/3; label=1;}
	mesh Th = buildmesh(
		 bi1(2*NP*L/3)+bi2(NP*l/3)+bi3(2*NP*L/3)+bi4(NP*l/3)
		+ Res1(NPres)+Res2(NPres)+Res3(NPres)+Res4(NPres)+Res5(NPres)+Res6(NPres)
		+ b1(NP*L)+b2(NP*l)+b3(NP*L)+b4(NP*l));
	//plot(Th,wait=1);
	
	fespace Vh(Th,P1);
	Vh u,v,reg=region;
	real[int] K(reg[].max+1),pr(reg[].max+1);
	K=1;
	K[Th(L/2,l/2).region] = 10;
	pr=0;
	for(int i=0;i<nbres;++i) pr(Th(XX[3*i+1],XX[3*i+2]).region) = XX[3*i];
	//plot(reg,value=1,fill=1,wait=1);
	
	solve Chaleur(u,v) = int2d(Th)( (dx(u)*dx(v) + dy(u)*dy(v)) * K[region]) + int1d(Th,0)(u*v) - int2d(Th)(pr[region]*v) ;
	if(iter %200 == 199)
	plot(u,fill=1,nbiso=50,wait=0,dim=2);
	real meanu = int2d(Th,Th(L/2,l/2).region)(u)/int2d(Th,Th(L/2,l/2).region)(1.);
	real res = int2d(Th,Th(L/2,l/2).region)(square(u-target));
	cout << "J: iter " << iter++ << "  -  mean=" << meanu << "  -  corr=" << corr << "  -  res=" << res << " ==> " << res+pena*corr << endl;
	cout << "pos : ";
	for(int j=0;j<X.n;++j) cout << " " << X[j];
	cout << endl;
	return res + pena*corr;
}

real[int] grad(3*nbres);
bool plotgrad=1;
func real[int] dJ(real[int] &X)
{
	verbosity=0;
	real[int] XX=X,diff=X;
	for(int i=0;i<X.n;++i)
	{
		XX[i] = X[i]>Xmin[i] ? X[i] : Xmin[i];
		XX[i] = XX[i]<Xmax[i] ? XX[i] : Xmax[i];
	}
	diff -= XX;
	real corr = diff.l2;
	//cout << X << endl << XX << endl << "==>" << corr << endl;
	
	border b1(t=0,1){x=L*t; y=0; label=7;}
	border b2(t=0,1){x=L;   y=l*t; label=7;}
	border b3(t=0,1){x=L*(1-t); y=l; label=7;}
	border b4(t=0,1){x=0;   y=l*(1-t); label=7;}
	Resistor(1,XX[1] ,XX[2] ,0);
	Resistor(2,XX[4] ,XX[5] ,1);
	Resistor(3,XX[7] ,XX[8] ,2);
	Resistor(4,XX[10],XX[11],3);
	Resistor(5,XX[13],XX[14],4);
	Resistor(6,XX[16],XX[17],5);
	border bi1(t=0,2){x=(0.5+t)*L/3; y=l/3; label=6;}
	border bi2(t=0,1){x=5*L/6; y=(1+t)*l/3; label=6;}
	border bi3(t=0,2){x=(2.5-t)*L/3; y=2*l/3; label=6;}
	border bi4(t=0,1){x=L/6; y=(2-t)*l/3; label=6;}
	mesh Th = buildmesh(
		 bi1(2*NP*L/3)+bi2(NP*l/3)+bi3(2*NP*L/3)+bi4(NP*l/3)
		+ Res1(NPres)+Res2(NPres)+Res3(NPres)+Res4(NPres)+Res5(NPres)+Res6(NPres)
		+ b1(NP*L)+b2(NP*l)+b3(NP*L)+b4(NP*l));
	//plot(Th,wait=1);
	
	fespace Vh(Th,P1);
	Vh reg=region;
	real[int] K(reg[].max+1),pr(reg[].max+1);
	K=1;
	K[Th(L/2,l/2).region] = 10;
	pr=0;
	for(int i=0;i<nbres;++i) pr(Th(XX[3*i+1],XX[3*i+2]).region) = XX[3*i];
	//plot(reg,value=1,fill=1,wait=1);
	Vh u,v;
	solve Chaleur(u,v) = int2d(Th)( (dx(u)*dx(v) + dy(u)*dy(v)) * K[region]) + int1d(Th,7)(u*v) - int2d(Th)(pr[region]*v) ;
	int regt = Th(L/2,l/2).region;
	for(int i=0;i<nbres;++i)
	{
		Vh up,ux,uy;
		int thereg = Th(XX[3*i+1],XX[3*i+2]).region;
		//cout << "pr:" << pr[thereg] << endl;
		solve dpchaleur(up,v) = int2d(Th)( (dx(up)*dx(v) + dy(up)*dy(v)) * K[region]) + int1d(Th,7)(up*v) - int2d(Th,thereg)(v);
		grad[3*i] = int2d(Th,regt)(2*up*(u-target)) + 2*pena*diff[3*i];
		solve dxchaleur(ux,v) = int2d(Th)( (dx(ux)*dx(v) + dy(ux)*dy(v)) * K[region]) + int1d(Th,7)(ux*v) - int1d(Th,i)(pr[thereg]*v*N.x);
		grad[3*i+1] = int2d(Th,regt)(2*ux*(u-target)) + 2*pena*diff[3*i+1];
		solve dychaleur(uy,v) = int2d(Th)( (dx(uy)*dx(v) + dy(uy)*dy(v)) * K[region]) + int1d(Th,7)(uy*v) - int1d(Th,i)(pr[thereg]*v*N.y);
		grad[3*i+2] = int2d(Th,regt)(2*uy*(u-target)) + 2*pena*diff[3*i+2];
		if(plotgrad)
		{
			plot(up,nbiso=50,wait=0,cmm="Derivative w.r.t. pr["+i+"]");
			plot(ux,nbiso=50,wait=0,cmm="Derivative w.r.t. x["+i+"]");
			plot(uy,nbiso=50,wait=0,cmm="Derivative w.r.t. y["+i+"]");
		}
	}
	cout << "dJ : " << grad << endl;
	return grad;
}


real pst=5000;
real isdxy=3*L/24.;
real[int] isd  =[1000,isdxy,isdxy,1000,isdxy,isdxy,1000,isdxy,isdxy,1000,isdxy,5*l/6,1000,isdxy,isdxy,1000,isdxy,isdxy];
real[int] start=[pst,L/6,l/6,pst,L/2,l/6,pst,5*L/6,l/6,pst,L/6,5*l/6,pst,L/2,5*l/6,pst,5*L/6,5*l/6];

//dJ(start);
cout << cmaes(J,start,initialStdDevs=isd,stopTolFun=1e-1,stopMaxIter=200,stopMaxFunEval=1000);
//BFGS(J,dJ,start,nbiterline=10,nbiter=3000,eps=1.e-8);
iter = 1000;

// FFCS - regression test reference value
real regtest=J(start);
