--- trunk/gadget/hooke.cc	2015/07/23 19:00:38	11
+++ trunk/gadget/hooke.cc	2015/10/13 11:29:01	16
@@ -139,12 +139,12 @@
 #include "ecosystem.h"
 #include "global.h"
 
-#ifndef NO_OPENMP
+#ifdef _OPENMP
 #include "omp.h"
 #endif
 
 extern Ecosystem* EcoSystem;
-#ifndef NO_OPENMP
+#ifdef _OPENMP
 extern Ecosystem** EcoSystems;
 #endif
 
@@ -156,17 +156,9 @@
   DoubleVector z(point);
 
   minf = prevbest;
-//  for (int k=0;k<point.Size(); k++)
-//	  cout << z[k] ;
-//  cout << endl;
   for (i = 0; i < point.Size(); i++) {
-
-//	  for (int k=0;k<point.Size(); k++)
-//	  	  cout << z[k] << " " ;
-//cout << endl;
     z[param[i]] = point[param[i]] + delta[param[i]];
     ftmp = EcoSystem->SimulateAndUpdate(z);
-//    cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " << ftmp << endl;
     if (ftmp < minf) {
       minf = ftmp;
     } else {
@@ -178,7 +170,6 @@
       else
         z[param[i]] = point[param[i]];
     }
-//    cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " << ftmp  <<" - " << prevbest << endl;
   }
 
   for (i = 0; i < point.Size(); i++)
@@ -187,9 +178,13 @@
 }
 
 /* given a point, look for a better one nearby, one coord at a time */
-#ifndef NO_OPENMP
-double OptInfoHooke::bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
-//cout << "beastNearbyINI" << endl;
+#ifdef _OPENMP
+/*
+ * function bestBeraby parallelized with OpenMP
+ * · 2 threads per coord to parallelize the calculation of +delta/-delta
+ * · parallelize the calculation of the best nearby of the coord
+ */
+double OptInfoHooke::bestNearbyRepro(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
   double minf;//, ftmp;
   int i, j, k;
   DoubleVector z(point);
@@ -214,19 +209,13 @@
 	  return -1;
   }
 
+  omp_set_dynamic(0);
+  omp_set_nested(1); //permit the nested parallelization
   while ( i < nvars) {
-//	  for (int k=0;k<point.Size(); k++)
-//	  	  cout << z[k] << " " ;
-//	    cout << endl;
 	  if ((i + paral_tokens -1) >= nvars)
 		  paral_tokens = nvars - i;
-//cout << "i:" << i << endl;
-	  omp_set_dynamic(0);
-	  omp_set_nested(1);
-#pragma omp parallel for num_threads(paral_tokens) private(k)
+#pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2)
 	  for (j = 0; j < paral_tokens; ++j) {
-//		  cout << "thr_for:" << omp_get_num_threads()<<endl;
-//		  cout << "J:" << j << endl;
 		  storage[j].z = z;
 		  storage[j].delta = delta;
 		  DoubleVector v1(z);
@@ -235,42 +224,31 @@
 		  v1[k] +=  delta[k];
 		  v2[k] -=  delta[k];
 
-#pragma omp parallel sections num_threads(2)
+#pragma omp parallel sections num_threads(2) //parrallelize the +/- delta simulation for each parameter
 		  {
-//			  cout << "thr_sec:" << omp_get_num_threads()<<endl;
-  #pragma omp section
+	#pragma omp section
 			  {
 		  storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
-//		  cout << "->" << j << endl;
 			  }
-  #pragma omp section
+	#pragma omp section
 			  {
-//				  cout << "->" << j+paral_tokens << endl;
 		  storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
 			  }
 		  }
-//cout << "-----" << endl;
-//		  cout <<i+j<< "-z["<< param[i+j]<<"]:" << (storage[j].z)[param[i+j]] << " - " << storage[j].ftmp << endl;
 
 		  if (storage[j].ftmp < minf) {
 			  storage[j].iters = 1;
 			  storage[j].z[k] = v1[k];
-//			  storage[j].delta[k] += delta[k];
-//		        minf = ftmp;
 		  } else {
 			  storage[j].iters = 2;
 			  storage[j].delta[k] = 0.0 - delta[k];
 			  if (storage[j+paral_tokens].ftmp < minf) {
-//				  minf = ftmp;
 				  storage[j].ftmp = storage[j+paral_tokens].ftmp;
 				  storage[j].z[k] = v2[k];
 			  }
-//			  else
-//				  storage[j].z[k] = point[k];
 		  }
 	  }
 
-//	  cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " ;
 	  for (j = 0; j < paral_tokens; ++j) {
 		  i++;
 		  iters += storage[j].iters;
@@ -278,14 +256,11 @@
 			  minf = storage[j].ftmp;
 			  z = storage[j].z;
 			  delta = storage[j].delta;
-//			  cout << storage[j].ftmp;
 			  break;
 		  }
 	  }
   }
-
-//  omp_set_num_threads(numThr);
-
+  delete[] storage;
   for (i = 0; i < nvars; ++i)
     point[i] = z[i];
   return minf;
@@ -314,9 +289,9 @@
   IntVector trapped(nvars, 0);
 
   EcoSystem->scaleVariables();
-#ifndef NO_OPENMP
+#ifdef _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
   EcoSystem->getOptScaledValues(x);
@@ -337,7 +312,6 @@
     bestx[i] = x[i];
     trialx[i] = x[i];
     param[i] = i;
-    //FIXME rand_r
     delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign
   }
 
@@ -360,7 +334,7 @@
 
   while (1) {
     if (isZero(bestf)) {
-#ifdef NO_OPENMP
+#ifndef _OPENMP
       iters = EcoSystem->getFuncEval() - offset;
 #endif
       handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
@@ -371,7 +345,6 @@
     /* randomize the order of the parameters once in a while */
     rchange = 0;
     while (rchange < nvars) {
-    	//FIXME rand_r
       rnumber = rand() % nvars;
       rcheck = 1;
       for (i = 0; i < rchange; i++)
@@ -386,8 +359,8 @@
     /* find best new point, one coord at a time */
     for (i = 0; i < nvars; i++)
       trialx[i] = x[i];
-#ifndef NO_OPENMP
-    newf = this->bestNearbyOMP(delta, trialx, bestf, param);
+#ifdef _OPENMP
+    newf = this->bestNearbyRepro(delta, trialx, bestf, param);
     if (newf == -1) {
     	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
     	handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
@@ -398,10 +371,9 @@
 #endif
     /* if too many function evaluations occur, terminate the algorithm */
 
-#ifdef NO_OPENMP
+#ifndef _OPENMP
     iters = EcoSystem->getFuncEval() - offset;
 #endif
-//    cout << "newf:" << newf << "-" << iters<<endl;
     if (iters > hookeiter) {
       handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
       handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
@@ -469,7 +441,7 @@
       /* only move forward if this is really an improvement    */
       oldf = newf;
       newf = EcoSystem->SimulateAndUpdate(trialx);
-#ifndef NO_OPENMP
+#ifdef _OPENMP
       iters++;
 #endif
       if ((isEqual(newf, oldf)) || (newf > oldf)) {
@@ -482,8 +454,8 @@
       for (i = 0; i < nvars; i++)
         x[i] = trialx[i];
 
-#ifndef NO_OPENMP
-      newf = this->bestNearbyOMP(delta, trialx, bestf, param);
+#ifdef _OPENMP
+      newf = this->bestNearbyRepro(delta, trialx, bestf, param);
       if (newf == -1) {
           	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
           	handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
@@ -496,10 +468,9 @@
         break;
 
       /* if too many function evaluations occur, terminate the algorithm */
-#ifdef NO_OPENMP
+#ifndef _OPENMP
       iters = EcoSystem->getFuncEval() - offset;
 #endif
-//      cout << "newf:" << newf << "-" << iters<<endl;
       if (iters > hookeiter) {
         handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
         handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
@@ -514,9 +485,9 @@
         EcoSystem->storeVariables(score, bestx);
         return;
       }
-    }
+    } // while (newf < bestf)
 
-#ifdef NO_OPENMP
+#ifndef _OPENMP
     iters = EcoSystem->getFuncEval() - offset;
 #endif
     if (newf < bestf) {
@@ -552,11 +523,10 @@
   }
 }
 
-#ifdef GADGET_OPENMP
-  //TODO doc
-double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
-//cout << "beastNearbyINI" << endl;
-  double minf;//, ftmp;
+/* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
+#ifdef SPECULATIVE
+double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
+  double minf;
   int i, j, k, ii;
   DoubleVector z(point);
   int bestId = 0;
@@ -569,7 +539,6 @@
 
   minf = prevbest;
 
-
   int paral_tokens, numThr, nvars = point.Size();
   numThr = omp_get_max_threads ( );
 
@@ -580,21 +549,15 @@
 	  return -1;
   }
 
+  omp_set_dynamic(0);
+  omp_set_nested(1); //permit the nested parallelization
   for (ii=0; ii< paral_tokens; ii++) {
 	  i = 0;
 	  while ( i < nvars) {
-	//	  for (int k=0;k<point.Size(); k++)
-	//	  	  cout << z[k] << " " ;
-	//	    cout << endl;
 		  if ((i + paral_tokens -1) >= nvars)
 			  paral_tokens = nvars - i;
-	//cout << "i:" << i << endl;
-		  omp_set_dynamic(0);
-		  omp_set_nested(1);
 	#pragma omp parallel for num_threads(paral_tokens) private(k)
 		  for (j = 0; j < paral_tokens; ++j) {
-	//		  cout << "thr_for:" << omp_get_num_threads()<<endl;
-	//		  cout << "J:" << j << endl;
 			  storage[j].z = z;
 			  storage[j].delta = delta;
 			  DoubleVector v1(z);
@@ -603,42 +566,31 @@
 			  v1[k] +=  delta[k];
 			  v2[k] -=  delta[k];
 
-	#pragma omp parallel sections num_threads(2)
+#pragma omp parallel sections num_threads(2)
 			  {
-	//			  cout << "thr_sec:" << omp_get_num_threads()<<endl;
 	  #pragma omp section
 				  {
 			  storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
-	//		  cout << "->" << j << endl;
 				  }
 	  #pragma omp section
 				  {
-	//				  cout << "->" << j+paral_tokens << endl;
 			  storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
 				  }
 			  }
-	//cout << "-----" << endl;
-	//		  cout <<i+j<< "-z["<< param[i+j]<<"]:" << (storage[j].z)[param[i+j]] << " - " << storage[j].ftmp << endl;
-
 			  if (storage[j].ftmp < minf) {
 				  storage[j].iters = 1;
 				  storage[j].z[k] = v1[k];
-	//			  storage[j].delta[k] += delta[k];
-	//		        minf = ftmp;
 			  } else {
 				  storage[j].iters = 2;
 				  storage[j].delta[k] = 0.0 - delta[k];
 				  if (storage[j+paral_tokens].ftmp < minf) {
-	//				  minf = ftmp;
 					  storage[j].ftmp = storage[j+paral_tokens].ftmp;
 					  storage[j].z[k] = v2[k];
 				  }
 				  else iters += 2;
-	//				  storage[j].z[k] = point[k];
 			  }
 		  }
 
-	//	  cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " ;
 		  bestId = 0;
 		  for (j = 1; j < paral_tokens; ++j) {
 			  if (storage[j].ftmp < storage[bestId].ftmp)
@@ -649,15 +601,12 @@
 			  minf = storage[bestId].ftmp;
 			  z = storage[bestId].z;
 			  delta = storage[bestId].delta;
-	//			  cout << storage[j].ftmp;
-
 		  }
 
 		  i += paral_tokens;
 	  }
 	}
 
-//  omp_set_num_threads(numThr);
   delete[] storage;
   for (i = 0; i < nvars; ++i)
     point[i] = z[i];
@@ -686,11 +635,9 @@
 	   IntVector trapped(nvars, 0);
 
 	   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
 	   EcoSystem->getOptScaledValues(x);
 	   EcoSystem->getOptLowerBounds(lowerb);
 	   EcoSystem->getOptUpperBounds(upperb);
@@ -709,7 +656,6 @@
 	     bestx[i] = x[i];
 	     trialx[i] = x[i];
 	     param[i] = i;
-	     //FIXME rand_r
 	     delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign
 	   }
 
@@ -732,7 +678,7 @@
 
 	   while (1) {
 	     if (isZero(bestf)) {
-	 #ifdef NO_OPENMP
+	 #ifndef _OPENMP
 	       iters = EcoSystem->getFuncEval() - offset;
 	 #endif
 	       handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
@@ -743,7 +689,6 @@
 	     /* randomize the order of the parameters once in a while */
 	     rchange = 0;
 	     while (rchange < nvars) {
-	     	//FIXME rand_r
 	       rnumber = rand() % nvars;
 	       rcheck = 1;
 	       for (i = 0; i < rchange; i++)
@@ -758,8 +703,8 @@
 	     /* find best new point, one coord at a time */
 	     for (i = 0; i < nvars; i++)
 	       trialx[i] = x[i];
-	 #ifndef NO_OPENMP
-	     newf = this->bestNearbyOMP2(delta, trialx, bestf, param);
+	 #ifdef _OPENMP
+	     newf = this->bestNearbySpec(delta, trialx, bestf, param);
 	     if (newf == -1) {
 	     	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
 	     	handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
@@ -770,10 +715,9 @@
 	 #endif
 	     /* if too many function evaluations occur, terminate the algorithm */
 
-	 #ifdef NO_OPENMP
+	 #ifndef _OPENMP
 	     iters = EcoSystem->getFuncEval() - offset;
 	 #endif
-	 //    cout << "newf:" << newf << "-" << iters<<endl;
 	     if (iters > hookeiter) {
 	       handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
 	       handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
@@ -786,7 +730,6 @@
 	       for (i = 0; i < nvars; i++)
 	         bestx[i] = trialx[i] * init[i];
 	       EcoSystem->storeVariables(score, bestx);
-//	       EcoSystem->writeOptValuesOMP();
 	       return;
 	     }
 
@@ -842,7 +785,7 @@
 	       /* only move forward if this is really an improvement    */
 	       oldf = newf;
 	       newf = EcoSystem->SimulateAndUpdate(trialx);
-	 #ifndef NO_OPENMP
+	 #ifdef _OPENMP
 	       iters++;
 	 #endif
 	       if ((isEqual(newf, oldf)) || (newf > oldf)) {
@@ -855,8 +798,8 @@
 	       for (i = 0; i < nvars; i++)
 	         x[i] = trialx[i];
 
-	 #ifndef NO_OPENMP
-	       newf = this->bestNearbyOMP2(delta, trialx, bestf, param);
+	 #ifdef _OPENMP
+	       newf = this->bestNearbySpec(delta, trialx, bestf, param);
 	       if (newf == -1) {
 	           	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
 	           	handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
@@ -869,10 +812,9 @@
 	         break;
 
 	       /* if too many function evaluations occur, terminate the algorithm */
-	 #ifdef NO_OPENMP
+	 #ifndef _OPENMP
 	       iters = EcoSystem->getFuncEval() - offset;
 	 #endif
-	 //      cout << "newf:" << newf << "-" << iters<<endl;
 	       if (iters > hookeiter) {
 	         handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
 	         handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
@@ -885,12 +827,11 @@
 	         for (i = 0; i < nvars; i++)
 	           bestx[i] = trialx[i] * init[i];
 	         EcoSystem->storeVariables(score, bestx);
-//	         EcoSystem->writeOptValuesOMP();
 	         return;
 	       }
 	     }
 
-	 #ifdef NO_OPENMP
+	 #ifndef _OPENMP
 	     iters = EcoSystem->getFuncEval() - offset;
 	 #endif
 	     if (newf < bestf) {
@@ -916,7 +857,6 @@
 	       score = bestf;
 	       handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
 	       EcoSystem->storeVariables(bestf, bestx);
-//	       EcoSystem->writeOptValuesOMP();
 	       return;
 	     }