--- trunk/gadget/simann.cc	2015/07/23 19:00:38	11
+++ trunk/gadget/simann.cc	2015/07/30 12:36:46	14
@@ -186,9 +186,8 @@
 #ifndef NO_OPENMP
 extern Ecosystem** EcoSystems;
 #endif
-//extern StochasticData* data;
-
 
+/*sequential code replaced at seq_optimize_template.h*/
 //void OptInfoSimann::OptimiseLikelihood() {
 //
 //  //set initial values
@@ -305,8 +304,6 @@
 ////              trialx[i] = x[i];
 ////          }
 //
-////          cout << "param:" << param[l] << "-" << trialx[param[l]] << endl;
-//
 //          //Evaluate the function with the trial point trialx and return as -trialf
 //          trialf = EcoSystem->SimulateAndUpdate(trialx);
 //          trialf = -trialf;
@@ -327,9 +324,7 @@
 //            handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score);
 //            return;
 //          }
-////          cout << "f:" << trialf << endl;
 //          //Accept the new point if the new function value better
-////          cout << "mustAccept:" << trialf << "|" << funcval<< endl;
 //          if ((trialf - funcval) > verysmall) {
 //            for (i = 0; i < nvars; i++)
 //              x[i] = trialx[i];
@@ -353,7 +348,6 @@
 //              nrej++;
 //            }
 //          }
-////          cout << "goog VALUE:" << funcval << endl;
 //          // JMB added check for really silly values
 //          if (isZero(trialf)) {
 //            handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0");
@@ -452,9 +446,10 @@
 
   double tmp, p, pp, ratio, nsdiv;
   double fopt, funcval, trialf;
-  int    a, i, j, k, l, offset, quit;
+  int    a, i, j, k, l, quit;
   int    rchange, rcheck, rnumber;  //Used to randomise the order of the parameters
 
+  // store the info of the different threads
   struct Storage {
   	  DoubleVector trialx;
   	  double newLikelihood;
@@ -475,8 +470,12 @@
   IntVector nacp(nvars, 0);
 
   EcoSystem->resetVariables();  //JMB need to reset variables in case they have been scaled
-  if (scale)
-    EcoSystem->scaleVariables();
+  if (scale) {
+  		EcoSystem->scaleVariables();
+  		int numThr = omp_get_max_threads ( );
+  		for(i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
+  			EcoSystems[i]->scaleVariables();
+  	}
   EcoSystem->getOptScaledValues(x);
   EcoSystem->getOptLowerBounds(lowerb);
   EcoSystem->getOptUpperBounds(upperb);
@@ -507,13 +506,11 @@
     handle.logMessage(LOGINFO, "Error starting Simulated Annealing optimisation with f(x) = infinity");
     converge = -1;
     iters = 1;
-//    EcoSystem->writeOptValuesOMP();
     return;
   }
 
   //the function is to be minimised so switch the sign of funcval (and trialf)
   funcval = -funcval;
-  offset = EcoSystem->getFuncEval();  //number of function evaluations done before loop
   nacc++;
   cs /= lratio;  //JMB save processing time
   nsdiv = 1.0 / ns;
@@ -532,15 +529,16 @@
   for (i=0; i<numThr; ++i)
 	  storage[i].trialx = trialx;
 
-  //Start the main loop.  Note that it terminates if
-  //(i) the algorithm succesfully optimises the function or
-  //(ii) there are too many function evaluations
-  int aux;
 
- DoubleVector vns(nvars, 0);
+  int aux; //store the number of evaluations that are not useful for the algorithm
+
+ DoubleVector vns(nvars, 0); //vector of ns
  int ns_ = ceil(numThr/2.);
  double res;
   aux=0;
+  //Start the main loop.  Note that it terminates if
+  //(i) the algorithm succesfully optimises the function or
+  //(ii) there are too many function evaluations
   while (1) {
     for (a = 0; a < nt; ++a) {
       //Randomize the order of the parameters once in a while, to avoid
@@ -652,12 +650,6 @@
               x[i] = storage[bestId].trialx[i];
             funcval = trialf;
             nacc++;
-//            vns[param[l+bestId]]++;
-            //nacp[param[l+bestId]]+=numThr;  //JMB - not sure about this ...
-//            for (i = 0; i < numThr; ++i)
-//            	nacp[param[l+i]]++;
-//            nacp[param[l+bestId]]++;
-
           } else {
             //Accept according to metropolis condition
             p = expRep((trialf - funcval) / t);
@@ -669,11 +661,6 @@
               funcval = trialf;
               naccmet++;
               nacp[param[(l+bestId)%nvars]]++;
-//              vns[param[l+bestId]]++;
-//              nacp[param[l+bestId]]+=numThr;
-//              for (i = 0; i < numThr; ++i)
-//            	  nacp[param[l+i]]++;
-//                nacp[param[l+bestId]]++;
             } else {
               //Reject point
               nrej++;
@@ -683,7 +670,6 @@
           if (isZero(trialf)) {
             handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0");
             converge = -1;
-//            EcoSystem->writeOptValuesOMP();
             delete[] storage;
             return;
           }
@@ -697,7 +683,6 @@
             if (scale) {
               for (i = 0; i < nvars; ++i)
                 scalex[i] = bestx[i] * init[i];
-              //FIXME EcoSystems[]->
               EcoSystem->storeVariables(-fopt, scalex);
             } else
             	EcoSystem->storeVariables(-fopt, bestx);
@@ -708,26 +693,6 @@
           }
         }
       }
-
-      //Adjust vm so that approximately half of all evaluations are accepted
-//      for (i = 0; i < nvars; ++i) {
-////    	  cout << "nacp["<<i<<"]:"<<nacp[i]<<endl;
-//        ratio = nsdiv * nacp[i];
-//        nacp[i] = 0;
-//        if (ratio > uratio) {
-//          vm[i] = vm[i] * (1.0 + cs * (ratio - uratio));
-//        } else if (ratio < lratio) {
-//          vm[i] = vm[i] / (1.0 + cs * (lratio - ratio));
-//        }
-//
-//        if (vm[i] < rathersmall)
-//          vm[i] = rathersmall;
-//        if (vm[i] > (upperb[i] - lowerb[i]))
-//          vm[i] = upperb[i] - lowerb[i];
-//
-////        cout << "vm["<<i<<"]:"<<vm[i]<<endl;
-//
-//      }
     }
 
     //Check termination criteria
@@ -758,7 +723,6 @@
       converge = 1;
       score = EcoSystem->SimulateAndUpdate(bestx);
       handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score);
-//      EcoSystem->writeOptValuesOMP();
       delete[] storage;
       return;
     }
@@ -773,16 +737,14 @@
     for (i = 0; i < nvars; ++i)
       x[i] = bestx[i];
   }
-//  EcoSystem->writeOptValuesOMP();
 }
 #endif
 
+// calcule a new point
 void OptInfoSimann::newValue(int nvars, int l, IntVector& param, DoubleVector& trialx,
 		DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm)
 {
 	int i, k;
-//	double old = trialx[param[l]];
-//	cout << "[" << endl;
 	  for (i = 0; i < nvars; ++i) {
 		if (i == param[l]) {
 		  trialx[i] = x[i] + ((randomNumber(&seed) * 2.0) - 1.0) * vm[i];
@@ -800,20 +762,19 @@
 		  }
 		} else
 		  trialx[i] = x[i];
-//		cout <<" " << trialx[i];
 	  }
-
-//	  cout << endl<< "old[" <<param[l] << "]:" << old <<"-" << trialx[param[l]]<< " - " << vm[param[l]]<<endl;
 }
 
+/*################################################################################
+ * code use in the template (seq_optimize_template.h)
+ ################################################################################*/
+
 
 //Generate trialx, the trial value of x
 void newValue(int nvars, int l, IntVector& param, DoubleVector& trialx,
 		DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm, unsigned* seed)
 {
 	int i, k;
-//	double old = trialx[param[l]];
-//	cout << "[" << endl;
 	  for (i = 0; i < nvars; ++i) {
 		if (i == param[l]) {
 		  trialx[i] = x[i] + ((randomNumber(&*seed) * 2.0) - 1.0) * vm[i];
@@ -831,10 +792,7 @@
 		  }
 		} else
 		  trialx[i] = x[i];
-//		cout  << trialx[i]<<" ";
 	  }
-
-//	  cout << endl <<l<< "-old[" <<param[l] << "]:" << old <<"-" << trialx[param[l]]<< " - " << vm[param[l]]<<endl;
 }
 
 void buildNewParams_f(Siman& seed, DoubleVector& params) {
@@ -850,11 +808,9 @@
 #ifndef NO_OPENMP
 	int id = omp_get_thread_num ();
 	trialf = EcoSystems[id]->SimulateAndUpdate(params);
-//	cout << id << "-";
 #else
 	trialf = EcoSystem->SimulateAndUpdate(params);
 #endif
-//	cout << "trailf:" << trialf<< endl;
   return -trialf;
 }
 
@@ -872,12 +828,7 @@
 		double uratio = seed.getUratio();
 		double lratio = seed.getLratio();
 
-	//	cout << "uratio:" << uratio << endl;
-	//	cout << "lratio:" << lratio << endl;
-
-	//cout << "adjust vm" << endl;
 		for (i = 0; i < seed.getNvars(); i++) {
-//			cout << "nacp[" << i << "]:" << seed.getNacp(i) << endl;
 			ratio = nsdiv * seed.getNacp(i);
 			seed.setNacp(i,0);
 			if (ratio > uratio) {
@@ -906,7 +857,6 @@
 
 		DoubleVector* bestx = seed.getBestx();
 
-		//funcval = fopt;
 		for (i = 0; i < seed.getNvars(); i++)
 			x[i] = (*bestx)[i];
 	}
@@ -945,7 +895,6 @@
 	  //Accept the new point if the new function value better
 	  bool ret = true;
 	  int i;
-//	  cout << "mustAccept:" << trialf << "|" << funcval<< endl;
 	  int aux = siman.getParam()[siman.getL()];
 	if ((trialf - funcval) > verysmall) {
 		siman.incrementNacp(aux);
@@ -955,32 +904,22 @@
 		double p, pp;
 		//Accept according to metropolis condition
 
-//		cout << "t:" << siman.getT() << endl;
-
 		p = expRep((trialf - funcval) / siman.getT());
 		pp = randomNumber(siman.getSeedM());
 		if (pp < p) {
-//			cout << pp << "<" << p << endl;
 			siman.incrementNacp(aux);
 			siman.incrementNaccmet();
 //			//Accept point
-//			for (i = 0; i < nvars; i++)
-//				x[i] = trialx[i];
-//			funcval = trialf;
-//			naccmet++;
-//			nacp[param[l]]++;
 		} else {
 			//Reject point
 			ret = false;
 			siman.incrementNrej();
 		}
 	}
-//	cout << "nacp[" << aux << "]:" << siman.getNacp(aux) << endl;
 	// JMB added check for really silly values
 	if (isZero(trialf)) {
 		handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after",
 				iters, "function evaluations, f(x) = 0");
-
 		siman.setConverge(-1);
 		return false;
 	}
@@ -989,16 +928,10 @@
 	if (siman.getL() == 0)
 		siman.incrementNS();
 	if (siman.getNS() >= siman.getNs()){
-
 		siman.setNS(0);
-
 		siman.incrementNT();
 		adjustVm(siman);
 		siman.randomizeParams();
-//		if (seed.getNT() >= seed.getNt()){
-//
-//		}
-
 	}
 
 	return ret;
@@ -1015,8 +948,6 @@
 	 bool quit = false;
 	  if (siman.getNT() >= siman.getNt())
 	  {
-//		  cout << "nt:" << siman.getNT();
-		  //Check termination criteria
 		  int i;
 		  DoubleVector fstar = siman.getFstar();
 
@@ -1067,18 +998,12 @@
 	handle.logMessage(LOGINFO, "Number of rejected points", siman.getNrej());
 	*score = EcoSystem->SimulateAndUpdate(*bestX);
 	handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", *score);
-//	EcoSystem->writeOptValues();
   }
 };
 
 // Required
 std::ostream &operator<<(std::ostream &os, const DoubleVector &p)
 {
-//  os <<  "[ AUCH";
-//  for (int i = 0; i< NPARAMS; ++i) {
-//    std::cout << p[i] << ' ';
-//  }
-//  std::cout << "]";
 	os << "";
   return os;
 }
@@ -1089,12 +1014,12 @@
 	//set initial values
 
 	double tmp, p, pp;
-	double fopt, funcval, trialf;
+	double funcval, trialf;
 	int a, i, j, k, l, quit;
 	int rchange, rcheck, rnumber; //Used to randomise the order of the parameters
 
 	handle.logMessage(LOGINFO,
-			"\nStarting Simulated Annealing optimisation algorithm----\n");
+			"\nStarting Simulated Annealing optimisation algorithm\n");
 	int nvars = EcoSystem->numOptVariables();
 	DoubleVector x(nvars);
 	DoubleVector init(nvars);
@@ -1112,7 +1037,7 @@
 		EcoSystem->scaleVariables();
 #ifndef NO_OPENMP
 		int numThr = omp_get_max_threads ( );
-		for(i = 0; i < numThr; i++)
+		for(i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
 			EcoSystems[i]->scaleVariables();
 #endif
 	}
@@ -1152,30 +1077,23 @@
 
 	//the function is to be minimised so switch the sign of funcval (and trialf)
 	funcval = -funcval;
-//	offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop
-//	nacc++;
 	cs /= lratio;  //JMB save processing time
-//	nsdiv = 1.0 / ns;
-	fopt = funcval;
 	for (i = 0; i < tempcheck; i++)
 		fstar[i] = funcval;
 
-	//unsigned seed = 1234;
-
 	Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns),
 			tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score);
 
 	ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_f, buildNewParams_f>
 	pa(s, x, simanniter);
 
-//	cout << "Sequential run" << endl;
-//	pa.seq_opt();
-
 #ifndef NO_OPENMP
+	// OpenMP parallelization
 	int numThr = omp_get_max_threads ( );
-	pa.paral_opt_omp(numThr,numThr);
+	pa.paral_opt_omp(funcval,numThr,numThr);
 #else
-	pa.seq_opt();
+	// sequential code
+	pa.seq_opt(funcval);
 #endif
 
 }