// -- load packages and programs 

load "lapack"
load "UMFPACK64"
load "ffrandom"
include "getARGV.idp"
include "subprogs.idp"
include "meshprog.idp" 
include "checkboundarylabels.idp"
verbosity = 0; // to remove spurious FreeFEM prints

// -- basic parameters and variables

int k; // iteration count variable for the optimization process
int kmax = 1000; // maximum number of iterations
int countalpha, countalphaaux; // it indicates the number of iterations in the Armijo loop
int isconnected; // index for identifying if an initial guess is connected 
int countrejectedcon = 0; // it counts the number of rejected configurations when generatng one single initial guess
int indexstop = -1; // index for identifying which stopping criterion was reached
int M = 1; // M = 1 for monotone Armijo together with alphaspectral = alphainit | M > 1 in case of nonmonotone line search
int stopAlg1 = 0; // variable for algorithm 2
int startAlg2 = 0; // variable for algorithm 2
int steps, stepsinit, stepsadd; stepsinit = 2; // for set the random initial guess for the current round
real aux; // auxiliary variable 
real errorci; // error related to the centers of the balls
real alpha; // Armijo steplength
real alphaspectral; // for nonmonotone line search
real alphainit = 0.1; // initial Armijo steplength
real gamma = 1e-4; // for the Armijo loop
real tolnsd = 1e-4; // error tolerance for nsd
real tolerrorci = 1e-4; // error tolerance for errorci
real eps = 1e-20; // small value
real Tinitround, Tfinalround; // it marks the beginning of the time counting
string filename, foldername;
mesh Th; // the finite element mesh for the discretization of the Omega-domain

// -- geometric domain parameters and number of rounds
int ncircles; ncircles = getARGV("-ncircles", 1);
int nroundscalled;  nroundscalled = getARGV("-nroundscalled", 1); // nroundscalled = 3;
int nroundsevaluated; nroundsevaluated = getARGV("-nroundsevaluated", 1); // nroundsevaluated = nroundscalled; 
real squaresize = getARGV("-squaresize", 1);  // the length of the side of the square

// -- other variables 
int nev = 1; //number of computed eigen value close to sigma
int icircle;
int countlambdak, countlambdaround; // number of evaluations of the objective funcion per iteration (countlambdak) and per round (countlambdaround); 
int round, optround;
real sigma = 0; //value of the shift
real lambda, lambdatrial, lambdaold, errorlambda;
real sumcountlambdak, meancountlambdaperk, maxlastlambdas;
real lambdaround, timeround, optlambda, opttime, sumtimerounds, meantimeperround;
real sdx, sdy, nsd;

// -- all the vectors 
int[int] alllabels = 1:ncircles;
int[int] countlambdakv(kmax); 
int[int] v = 1:(ncircles-1);
real[int] diffxy(2*ncircles);
real[int] diffallsd(2*ncircles);
real[int] lastlambdas(M); lastlambdas = -tgv;
real[int] ev(nev); //to store the nev eigenvalue
real[int] optresultsAlg1(11);
real[int] optresultsAlgs1n2(11);

// -- all the matrices 
real[int,int] circles(ncircles,6);
real[int,int] circlesold(ncircles,circles.m);
real[int,int] circlestrial(ncircles,circles.m);
real[int,int] optinitguesses(ncircles,6);
real[int,int] finaltop(circles.n,circles.m);
real[int,int] allsd(ncircles,2);
real[int,int] allsdold(ncircles,2);
matrix OP, B;

// -- find best result from Algorithhm 1 

foldername = "Results" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Alg1";
filename = "summaryresults" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Alg1.txt";
ifstream file0(foldername + "/" + filename);

optlambda = tgv; sumtimerounds = 0;

for(int i = 0; i < nroundsevaluated; i++){

	file0 >> round; 
	file0 >> lambdaround; 
	file0 >> timeround; 

	sumtimerounds = sumtimerounds + timeround;

	if( lambdaround <= optlambda ){
		optround = round;
		optlambda = lambdaround;
		opttime = timeround;
	}

}

// -- set the initial guess from the txt file containing all the initial guesses

foldername = "InitGuesses" + ncircles + "bR" + nroundscalled + "Q" + squaresize;
filename = "initguesses" + ncircles + "bR" + nroundscalled + "Q" + squaresize + ".txt";

if( optround == 1 ){ 
	steps = stepsinit;
}
else{
	stepsadd = (optround - 1)*ncircles*circles.m;
	steps = stepsinit + stepsadd;
}

{
	ifstream file(foldername + "/" + filename);

	for(int i = 1; i <= steps; i++){ 
		file >> aux; 
	}

	for(int i = 0; i < ncircles; i++){
		for(int j = 0; j < circles.m; j++){
			file >> circles(i,j);
		}
	}

	optinitguesses = circles;
}

// -- main program

// --- create a folder and file to store results

foldername = "OptimalResult" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Algs1n2";
exec("mkdir " + foldername);

filename = "historical" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Ri" + optround + "Algs1n2.txt";
ofstream file1(foldername + "/" + filename); int npfile1 = file1.precision(12); // save info in a file .txt

// --- round initilization

Tinitround = clock();

cout << "\n Round  " << optround << ": Algorithm 2: busy \n " << endl;

file1 << " \n" << " ------------------------------------------------- Round  " << optround << " ----------------------------------------------- " << " \n";

k = 0; // iteration index -- k = 0 for the initial guess, and k > 0 otherwise
countlambdakv = 0; // reinitialize the vector countlambdakv with zero entries
errorci = tgv; // reinitialize the variable errorci with tgv

// --- find all intersection points between circles
real[int,int] allip = findallip(circles);

// --- find all intersection points on the boundary
real[int,int] bip = findbip(circles, allip);

// --- find all the inner points 
real[int,int] innerpoints = findinnerpoints(allip, bip);

// --- map all arcs on the boundary
real[int,int] allarcs = mapallarcs(circles, bip, innerpoints);

// --- map all boundaries
real[int,int] allboundaries = mapallboundaries(circles, allarcs, bip);

// --- check boundary labels
if( allboundaries.n > 1 ){ real[int,int] allboundaries = checkboundarylabels(circles, allboundaries); }

// --- mesh generation
Th = createmesh(circles, allarcs, allboundaries);

// --- fespace (load "Element_P3" for cubic interpolation functions, and load "Element_P4" for fourth-order interpolation functions)
fespace Vh(Th, P1);
Vh u1, u2, u, utrial, uold;
Vh[int] eV(nev); //to store the nev eigenvector

// --- solve Eigenvalue problem (OP = A - sigma B is the shifted matrix)
varf op (u1, u2) = int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma*u1*u2 ) + on(alllabels, u1=0);    
OP = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition
B = b(Vh, Vh, solver=CG, eps=1e-20);
EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV, tol=1e-10, maxit=0, ncv=0);
u = eV[0]; recomputeu(u[]); // 1st eigenvector
lambda = ev[0]; // 1st eigenvalue

// --- update lambda-related counters and variables
countlambdaround = 1; countlambdak = 1; countlambdakv(k) = countlambdak;
lastlambdas(k) = lambda; maxlastlambdas = lastlambdas(0); 
for(int i = 1; i < M; i++){ maxlastlambdas = max(maxlastlambdas,lastlambdas(i)); }

// --- compute shape derivatives and the respective l2-norm
for(int i = 0; i < ncircles; i++){
	sdx = - int1d(Th,i+1)((dx(u)*dx(u) + dy(u)*dy(u))*N.x);
	sdy = - int1d(Th,i+1)((dx(u)*dx(u) + dy(u)*dy(u))*N.y);
	allsd(i,0) = sdx; allsd(i,1) = sdy; 
}
allsd(0,0) = 0.0; allsd(0,1) = 0.0; allsd(1,0) = 0.0;
nsd = computevectornorml2(allsd);

// --- write quantitative results in a file
file1 << " R = " << optround << " k = " << k << " lambda = " << lambda << " nsd = " << nsd << " \n";

// -----------------------------------------------------------
// ---------------- optimization process  --------------------
// -----------------------------------------------------------
// find the minimum eigenvalue by using the method of steepest 
// descent with the armijo steplength

while( nsd > tolnsd &&  k < kmax && ( (stopAlg1 == 0 && errorci > tolerrorci) || (stopAlg1 == 1 && startAlg2 == 1) ) ){
//while( nsd > tolnsd & errorci > tolerrorci & k < kmax ){

	// ---- store results from the previous iteration
	circlesold = circles; lambdaold = lambda; uold = u; allsdold = allsd;
	countalpha = 0; // reinitialize the variable countalpha with zero value

	if( stopAlg1 == 0 ){

		// ---- update alpha
		if( k == 0 ){
			alpha = alphainit;
		}
		else{
			alpha = alphaspectral;
		}
		
		// ---- update circles
		circlestrial = computecirclestrial(alpha, circles, allsd);
		
		// ---- find all intersection points between circles
		real[int,int] allip = findallip(circlestrial);

		// ---- find all intersection points on the boundary
		real[int,int] bip = findbip(circlestrial, allip);

		// ---- find all inner points
		real[int,int] innerpoints = findinnerpoints(allip, bip);

		// ---- map all arcs on the boundary
		real[int,int] allarcs = mapallarcs(circlestrial, bip, innerpoints);

		// ---- map all boundaries
		real[int,int] allboundaries = mapallboundaries(circlestrial, allarcs, bip);

		// ---- check boundary labels
		if( allboundaries.n > 1 ){ real[int,int] allboundaries = checkboundarylabels(circlestrial, allboundaries); }

		// ---- mesh generation
		Th = createmesh(circlestrial, allarcs, allboundaries);

		// ---- solve Eigenvalue problem (OP = A - sigma B is the shifted matrix)
		varf op (u1, u2) = int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma*u1*u2 ) + on(alllabels, u1=0);    
		OP = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
		varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition
		B = b(Vh, Vh, solver=CG, eps=1e-20);
		EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV, tol=1e-10, maxit=0, ncv=0);
		utrial = eV[0]; recomputeu(utrial[]); // 1st eigenvector
		lambdatrial = ev[0]; // 1st eigenvalue

		// ---- update lambda-related counters 
		countlambdak = 1; countlambdaround = countlambdaround + 1;
		
		// ---- find the armijo steplength

		while( !( lambdatrial <= maxlastlambdas - alpha*gamma*nsd^2 ) ){

			alpha = alpha/2.0; //cout << " alpha = " << alpha << endl;
			countalpha = countalpha + 1;

			// ----- update circles
			circlestrial = computecirclestrial(alpha, circles, allsd);
			
			// ----- find all intersection points between circles
			real[int,int] allip = findallip(circlestrial);

			// ----- find all intersection points on the boundary
			real[int,int] bip = findbip(circlestrial, allip);

			// ----- find all inner points
			real[int,int] innerpoints = findinnerpoints(allip, bip);

			// ----- map all arcs on the boundary
			real[int,int] allarcs = mapallarcs(circlestrial, bip, innerpoints);

			// ----- map all boundaries
			real[int,int] allboundaries = mapallboundaries(circlestrial, allarcs, bip);

			// ----- check boundary labels
			if( allboundaries.n > 1 ){ real[int,int] allboundaries = checkboundarylabels(circlestrial, allboundaries); }

			// ----- mesh generation
			Th = createmesh(circlestrial, allarcs, allboundaries);

			// ----- solve Eigenvalue problem (OP = A - sigma B is the shifted matrix)
			varf op (u1, u2) = int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma*u1*u2 ) + on(alllabels, u1=0);    
			OP = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
			varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition
			B = b(Vh, Vh, solver=CG, eps=1e-20);
			EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV, tol=1e-10, maxit=0, ncv=0);
			utrial = eV[0]; recomputeu(utrial[]); // 1st eigenvector
			lambdatrial = ev[0]; // 1st eigenvalue

			// ----- update lambda-related counters 
			countlambdak = countlambdak + 1; countlambdaround = countlambdaround + 1;

			// check if algorithm 2 should start
			errorci = computeerrorci(eps, circles, circlestrial);
			if(errorci <= tolerrorci){ stopAlg1 = 1; }

		}

	}

	// Algorithm 2 starts here !

	if( stopAlg1 == 1 ){

		startAlg2 = 0;

		for(int ind = 0; (ind < (ncircles-1)) && (startAlg2 == 0); ind++){

			icircle = v(ind);
			alpha = alphainit;
			countalphaaux = 0; 

			for(int j = 0; j < ncircles; j++){ circlestrial(j,0:5) = circles(j,0:5); }

			if( icircle == 1 ){
				circlestrial(icircle,0) = circles(icircle,0);
				circlestrial(icircle,1) = circles(icircle,1) - alpha*allsd(icircle,1);
			}
			else{
				circlestrial(icircle,0) = circles(icircle,0) - alpha*allsd(icircle,0);
				circlestrial(icircle,1) = circles(icircle,1) - alpha*allsd(icircle,1);
			}

			// -- find all intersection points between circles
			real[int,int] allip = findallip(circlestrial);

			// -- find all intersection points on the boundary
			real[int,int] bip = findbip(circlestrial, allip);

			// -- find all inner points
			real[int,int] innerpoints = findinnerpoints(allip, bip);

			// -- map all arcs on the boundary
			real[int,int] allarcs = mapallarcs(circlestrial, bip, innerpoints);

			// -- map all boundaries
			real[int,int] allboundaries = mapallboundaries(circlestrial, allarcs, bip);

			// -- check boundary labels
			if( allboundaries.n > 1 ){ real[int,int] allboundaries = checkboundarylabels(circlestrial, allboundaries); }

			// -- mesh generation
			Th = createmesh(circlestrial, allarcs, allboundaries);

			// solve Eigenvalue problem (OP = A - sigma B is the shifted matrix)
			varf op (u1, u2) = int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma*u1*u2 ) + on(alllabels, u1=0);    
			OP = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
			varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition
			B = b(Vh, Vh, solver=CG, eps=1e-20);
			EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV, tol=1e-10, maxit=0, ncv=0);
			utrial = eV[0]; recomputeu(utrial[]); // 1st eigenvector
			lambdatrial = ev[0]; // 1st eigenvalue
			// update lambda-related counters 
			countlambdak = countlambdak + 1; countlambdaround = countlambdaround + 1;

			// compute errorci
			errorci = sqrt( (circles(icircle,0) - circlestrial(icircle,0))^2 + (circles(icircle,1) - circlestrial(icircle,1))^2 )/sqrt( (circles(icircle,0))^2 + (circles(icircle,1))^2 );	
		
			if( ( lambdatrial < lambda ) && ( countalphaaux == 0 ) ){
				lambda = lambdatrial;
				startAlg2 = 1;
			}
	
		}

		//v for next iteration
		for(int i = 0; i < (ncircles-1); i++){ v(i) = (icircle + i) % (ncircles - 1) + 1; }

	}

	// Algorithm 2 ends here !

	// ---- update variables
	if( lambdatrial <= lambda ){ 
		circlesold = circles;
		circles = circlestrial; 
		lambda = lambdatrial;
		u = utrial;
	}
	countlambdakv(k+1) = countlambdak;
	lastlambdas( (k+1) % M ) = lambda; maxlastlambdas = lastlambdas(0); 
	for(int i = 1; i < M; i++){ maxlastlambdas = max(maxlastlambdas,lastlambdas(i)); }

	// ---- compute errors
	errorlambda = abs(lambda - lambdaold);
	errorci = computeerrorci(eps, circles, circlesold);

	// ---- compute the new shape derivatives and the respective l2-norm
	for(int i = 0; i < ncircles; i++){
		sdx = - int1d(Th,i+1)((dx(u)*dx(u) + dy(u)*dy(u))*N.x);
		sdy = - int1d(Th,i+1)((dx(u)*dx(u) + dy(u)*dy(u))*N.y);
		allsd(i,0) = sdx; allsd(i,1) = sdy; 
	}
	allsd(0,0) = 0.0; allsd(0,1) = 0.0; allsd(1,0) = 0.0;
	nsd = computevectornorml2(allsd);

	// ---- compute alpha for next iteration
	
	if( M == 1 ){
		alphaspectral = alphainit;
	}
	else{
		diffxy = computediffxy(circles, circlesold);
		diffallsd = computediffallsd(allsd, allsdold);
		alphaspectral = dotproduct(diffxy,diffxy)/dotproduct(diffxy,diffallsd);
		alphaspectral = max(1e-8, min(alphaspectral,1e8));
	}
	
	// ---- verify which stopping criterion was reached, if any
	// 0 for nsd <= tolnsd | 1 for errorci <= tolerrorci | 2 for k = kmax
	if( nsd <= tolnsd ){ indexstop = 0; }
	if( errorci <= tolerrorci ){ indexstop = 1; }
	if( k == kmax ){ indexstop = 2; }
	
	// ---- write quantitative results in a file
	k = k + 1;
	file1 << " R = " << optround << " k = " << k << " lambda = " << lambda << " nsd = " << nsd << " errorci = " << errorci << " alpha = " << alpha << " calpha = " << countalpha << " \n";
	
}

sumcountlambdak = 0;
for(int i = 0; i <= k; i++){ sumcountlambdak = sumcountlambdak + countlambdakv(i); }
meancountlambdaperk = sumcountlambdak/(k+1);

Tfinalround = clock() - Tinitround; 

// -- round results storage
filename = "results" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Ri" + optround + "Algs1n2.txt";
{
	ofstream file2(foldername + "/" + filename); int npfile2 = file2.precision(12); // save info in a file .txt
	file2 << " " << ncircles  << " " << squaresize << " " << nroundscalled << " " << optround << " " << k << " " << countlambdaround << " " << meancountlambdaperk << " " << indexstop << " " << lambda << " " << nsd << " " << Tfinalround << " \n";
}

// -- final topology storage
finaltop = circles;
filename = "finaltop" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Ri" + optround + "Algs1n2.txt";
ofstream file3(foldername + "/" + filename); int npfile3 = file3.precision(12); 
file3 << finaltop;

foldername = "Results" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Alg1";
filename = "results" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Ri" + optround + ".txt"; 

{
	ifstream file(foldername + "/" + filename);

	for(int i = 0; i < 11; i++){
		file >> optresultsAlg1(i);
	}
}

file1 << " ------------------------ " << " \n";
file1 << " Total time of Round " << optround << " = " << Tfinalround << "s" << " \n";
file1 << " ------------------------ " << " \n";

cout << " \n " << endl;
cout << " ---------------------------------------------------------------- " << endl;
cout << " -----      Optimal result obtained from Algorithhm 1       ----- " << endl;
cout << " ---------------------------------------------------------------- " << endl;
cout << " Optimal round = " << optround  << endl;
cout << " Optimal lambda = " << optresultsAlg1(8) << endl;
cout << " Total round time (only Alg 1) = " << optresultsAlg1(10) << "s" << endl;
cout << " ---------------------------------------------------------------- " << endl;
cout << " -----   Optimal result obtained from Algorithhms 1 and 2   ----- " << endl;
cout << " ---------------------------------------------------------------- " << endl;
cout << " Optimal lambda = " << lambda << endl;
cout << " Total round time (Algs 1 and 2) = " << Tfinalround << "s" << endl;
cout << " Total round time (only Alg 2)   = " << Tfinalround - optresultsAlg1(10) << "s" << endl;
cout << " ---------------------------------------------------------------- " << endl;
cout << " \n " << endl;

foldername = "InitGuesses" + ncircles + "bR" + nroundscalled + "Q" + squaresize;
filename = "initguesses" + ncircles + "bR" + nroundscalled + "Q" + squaresize + ".txt";

if( optround == 1 ){ 
	steps = stepsinit;
}
else{
	stepsadd = (optround - 1)*ncircles*optinitguesses.m;
	steps = stepsinit + stepsadd;
}

{
	ifstream file(foldername + "/" + filename);

	for(int i = 1; i <= steps; i++){ 
		file >> aux; 
	}

	for(int i = 0; i < ncircles; i++){
		for(int j = 0; j < optinitguesses.m; j++){
			file >> optinitguesses(i,j);
		}
	}
}

foldername = "OptimalResult" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Algs1n2";
filename = "results" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Ri" + optround + "Algs1n2.txt"; 

{
	ifstream file(foldername + "/" + filename);

	for(int i = 0; i < 11; i++){
		file >> optresultsAlgs1n2(i);
	}
}

foldername = "OptimalResult" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Algs1n2";
filename = "optimalresult" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Algs1n2.txt";

{
	ofstream fileoptAlgs1n2(foldername + "/" + filename); int npfileoptAlgs1n2 = fileoptAlgs1n2.precision(12);

	fileoptAlgs1n2 << "\n---------------------------------------------------------------";
	fileoptAlgs1n2 << "\n------- Optimal result obtained from Algorithms 1 and 2 -------";
	fileoptAlgs1n2 << "\n---------------------------------------------------------------";

	fileoptAlgs1n2 << "\n\nNumber of circles = " << optresultsAlgs1n2(0);
	fileoptAlgs1n2 << "\nSquare side length = " << optresultsAlgs1n2(1);

	fileoptAlgs1n2 << "\n\nOptimal Round = " << optresultsAlgs1n2(3);
	fileoptAlgs1n2 << "\nOptimal lambda = " << optresultsAlgs1n2(8);
	fileoptAlgs1n2 << "\nGradients' Euclidean norm = " << optresultsAlgs1n2(9);
	fileoptAlgs1n2 << "\nTotal number of iterations = " << optresultsAlgs1n2(4);
	fileoptAlgs1n2 << "\nTotal number of evaluations of the objective function = " << optresultsAlgs1n2(5);
	fileoptAlgs1n2 << "\nAverage number of evaluations of the objective function per iteration = " << optresultsAlgs1n2(6);
	fileoptAlgs1n2 << "\nStopping criterion fulfilled* = " << optresultsAlgs1n2(7);
	fileoptAlgs1n2 << "\n    *[0 if nsd <= tolnsd | 1 if errorci <= tolerrorci | 2 if k == kmax]";
	fileoptAlgs1n2 << "\nRound time = " << optresultsAlgs1n2(10) << "s";

	fileoptAlgs1n2 << "\n\nInitial guess corresponding to the Round " << optround << " in MetaPost format:";
	fileoptAlgs1n2 << "\n-----------------------------------------------------------------------------------";
	fileoptAlgs1n2 << "\n\\begin{figure} " << "\n";
	fileoptAlgs1n2 << "\\centering " << "\n";
	fileoptAlgs1n2 << "\\begin{mpost} " << "\n";
	fileoptAlgs1n2 << "input mpcolornames;" << "\n";
	fileoptAlgs1n2 << "u := 1.0cm;" << "\n";
	fileoptAlgs1n2 << "numeric m;" << "\n";
	fileoptAlgs1n2 << "pair c[];" << "\n";
	fileoptAlgs1n2 << "m := " << ncircles << ";" << "\n";
	for(int i=0; i < ncircles; i++){
		fileoptAlgs1n2 << "c[" << i+1 << "] := (" << optinitguesses(i,0) << "u," << optinitguesses(i,1) << "u);" << " \n";
	}
	fileoptAlgs1n2 << "for i=1 upto m: fill fullcircle scaled 2.0u shifted c[i] withcolor 0.9white; endfor" << "\n";
	fileoptAlgs1n2 << "for i=1 upto m: draw fullcircle scaled 2.0u shifted c[i] withcolor RoyalBlue withpen pencircle scaled 0.7; endfor" << "\n";
	fileoptAlgs1n2 << "\\end{mpost} " << "\n";
	fileoptAlgs1n2 << "\\label{fig:initguess" << ncircles << "bRi" << optround << "Q" << squaresize << "}" << " \n";
	fileoptAlgs1n2 << "\\end{figure} ";
	fileoptAlgs1n2 << "\n-----------------------------------------------------------------------------------\n";

	fileoptAlgs1n2 << "\nOptimal topology obtained from Round " << optround << " in MetaPost format:";
	fileoptAlgs1n2 << "\n-----------------------------------------------------------------------------------";
	fileoptAlgs1n2 << "\n\\begin{figure} " << "\n";
	fileoptAlgs1n2 << "\\centering " << "\n";
	fileoptAlgs1n2 << "\\begin{mpost} " << "\n";
	fileoptAlgs1n2 << "input mpcolornames;" << "\n";
	fileoptAlgs1n2 << "u := 1.0cm;" << "\n";
	fileoptAlgs1n2 << "numeric m;" << "\n";
	fileoptAlgs1n2 << "pair c[];" << "\n";
	fileoptAlgs1n2 << "m := " << ncircles << ";" << "\n";
	for(int i=0; i < ncircles; i++){
		fileoptAlgs1n2 << "c[" << i+1 << "] := (" <<  finaltop(i,0) << "u," << finaltop(i,1) << "u);" << " \n";
	}
	fileoptAlgs1n2 << "for i=1 upto m: fill fullcircle scaled 2.0u shifted c[i] withcolor 0.9white; endfor" << "\n";
	fileoptAlgs1n2 << "for i=1 upto m: draw fullcircle scaled 2.0u shifted c[i] withcolor RoyalBlue withpen pencircle scaled 0.7; endfor" << "\n";
	fileoptAlgs1n2 << "\\end{mpost} " << "\n";
	fileoptAlgs1n2 << "\\label{fig:finaltop" << ncircles << "bRi" << optround << "Q" << squaresize << "Algs1n2}" << " \n";
	fileoptAlgs1n2 << "\\end{figure} ";
	fileoptAlgs1n2 << "\n-----------------------------------------------------------------------------------\n";

}