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

// -- basic parameters

string filename, foldername; // name of the file
int isconnected; // index 0 or 1 for checking if there is a disjoint ball in a set of balls
int orowindex, crowindex, iaux;
real countrejectedcon; countrejectedcon = 0;
real meancountrejectedcon;

// -- geometric domain parameters and number of rounds
int ncircles; ncircles = getARGV("-ncircles", 1);
real[int,int] circles(ncircles,6);
int nroundi, nroundscalled; nroundscalled = getARGV("-nroundscalled", 1);
int[int,int] infoinitguesses(nroundscalled,5); infoinitguesses = 0;
real[int,int] initguesses(ncircles*nroundscalled,circles.m);

// -- random real numbers generation
real[int] randv(5000000); 
int randvi; randvi = 0;
randinit(100); //srandomdev(); // set a true random seed 
for(int i = 0; i < randv.n; i++){
	randv(i) = randreal1();
}

//-- square domain for generation of initial guesses
// --> square domain = [oep,oep] x [cep,cep] with its baricenter fixed at the origin
real squaresize, oep, cep; // oep/cep --> open/close end-point of the interval
squaresize = getARGV("-squaresize", 1); 
oep = -squaresize/2.0; cep = squaresize/2.0;

// -- create a folder for results storage
foldername = "InitGuesses" + ncircles + "bR" + nroundscalled + "Q" + squaresize;
exec("mkdir " + foldername);

// -- main program

for(int nroundi = 1; nroundi <= nroundscalled; nroundi++){ // loop for nroundscalled

	infoinitguesses(nroundi-1,0) = nroundi;
	infoinitguesses(nroundi-1,1) = squaresize;

	orowindex = (nroundi-1)*ncircles; crowindex = nroundi*ncircles - 1;

	//---------------------- random initial guess ------------------------------

	isconnected = 0;
	
	while( isconnected == 0 ){

		// 1st circle
		circles(0,0) = 0.0; circles(0,1) = 0.0; 
		circles(0,2) = 1.0; circles(0,3) = 0; 
		circles(0,4) = squaresize; circles(0,5) = nroundi;

		// 2nd circle
		if( randvi != 0 ){ randvi = randvi + 1; }
		infoinitguesses(nroundi-1,2) = randvi;
		circles(1,0) = 0.0; circles(1,1) = randv(randvi);
		circles(1,2) = 1.0; circles(1,3) = 1; 
		circles(1,4) = squaresize; circles(1,5) = nroundi;

		// the remaining (ncircles - 2) circles
		for(int i = 2; i < ncircles; i++){

			// create a new circle
			randvi = randvi + 1; 
			circles(i,0) = oep + 2*cep*randv(randvi);
			randvi = randvi + 1;
			circles(i,1) = oep + 2*cep*randv(randvi);
			circles(i,2) = 1.0; 
			circles(i,3) = i;
			circles(i,4) = squaresize;
			circles(i,5) = nroundi;

		}

		// -- check if the domain is connected
		isconnected = connected(circles);
		if( isconnected == 0 ){
			countrejectedcon = countrejectedcon + 1;
		}
		infoinitguesses(nroundi-1,4) = countrejectedcon;

	}

	infoinitguesses(nroundi-1,3) = infoinitguesses(nroundi-1,2) + (ncircles - 2)*2;

	// initial guess storage
	for(int i = orowindex; i <= crowindex; i++){
		for(int j = 0; j < circles.m; j++){
			iaux = i % ncircles;
			initguesses(i,j) = circles(iaux,j);
		}
	}

}

meancountrejectedcon = countrejectedcon/nroundscalled;

cout << "\n Random initial guesses generated. \n " << endl;

// -- create files to store the random initial guesses and their respective information
filename = "infoinitguesses" + ncircles + "bR" + nroundscalled + "Q" + squaresize + ".txt";
ofstream fileinfoinitguesses(foldername + "/" + filename); int npfileinfoinitguesses = fileinfoinitguesses.precision(12);
fileinfoinitguesses << infoinitguesses;

filename = "initguesses" + ncircles + "bR" + nroundscalled + "Q" + squaresize + ".txt";
ofstream fileinitguesses(foldername + "/" + filename); int npfileinitguesses = fileinitguesses.precision(12);
fileinitguesses << initguesses;

// -- create a folder and file to store results
foldername = "Results" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Alg1"; 
exec("mkdir " + foldername);
filename = "summaryresults" + ncircles + "bR" + nroundscalled + "Q" + squaresize + "Alg1.txt";
ofstream filesummaryresults(foldername + "/" + filename); int npfilesummaryresults = filesummaryresults.precision(12);