--- trunk/gadget/hooke.cc	2015/07/30 12:36:46	14
+++ trunk/gadget/hooke.cc	2017/04/07 09:20:55	20
@@ -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
 
@@ -178,13 +178,13 @@
 }
 
 /* given a point, look for a better one nearby, one coord at a time */
-#ifndef NO_OPENMP
+#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::bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
+double OptInfoHooke::bestNearbyRepro(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
   double minf;//, ftmp;
   int i, j, k;
   DoubleVector z(point);
@@ -209,42 +209,40 @@
 	  return -1;
   }
 
+//  omp_set_dynamic(0);
+//  omp_set_nested(1); //permit the nested parallelization
   while ( i < nvars) {
 	  if ((i + paral_tokens -1) >= nvars)
 		  paral_tokens = nvars - i;
-	  omp_set_dynamic(0);
-	  omp_set_nested(1); //permit the nested parallelization
-#pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2)
-	  for (j = 0; j < paral_tokens; ++j) {
+#pragma omp parallel for num_threads(paral_tokens*2) private(k) //parallelize the parameters (numThr)
+	  for (j = 0; j < (paral_tokens*2); ++j) {
 		  storage[j].z = z;
 		  storage[j].delta = delta;
-		  DoubleVector v1(z);
-		  DoubleVector v2(z);
-		  k = param[i+j];
-		  v1[k] +=  delta[k];
-		  v2[k] -=  delta[k];
+		  DoubleVector v(z);
 
-#pragma omp parallel sections num_threads(2) //parrallelize the +/- delta simulation for each parameter
-		  {
-	#pragma omp section
-			  {
-		  storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
-			  }
-	#pragma omp section
-			  {
-		  storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
-			  }
+		  if (j<paral_tokens) {
+			  k = param[i+j];
+			  v[k] +=  delta[k];
+		  }
+		  else {
+			  k = param[i+j-paral_tokens];
+			  v[k] -=  delta[k];
 		  }
 
-		  if (storage[j].ftmp < minf) {
+		  storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v);
+		  storage[j].z[k] = v[k];
+	  }
+	  for (j = 0; j < paral_tokens; ++j) {
+		  k = param[i+j];
+	  	  if (storage[j].ftmp < minf) {
 			  storage[j].iters = 1;
-			  storage[j].z[k] = v1[k];
+//			  storage[j].z[k] = v1[k];
 		  } else {
 			  storage[j].iters = 2;
 			  storage[j].delta[k] = 0.0 - delta[k];
 			  if (storage[j+paral_tokens].ftmp < minf) {
 				  storage[j].ftmp = storage[j+paral_tokens].ftmp;
-				  storage[j].z[k] = v2[k];
+				  storage[j].z[k] = storage[j+paral_tokens].z[k];;
 			  }
 		  }
 	  }
@@ -260,14 +258,12 @@
 		  }
 	  }
   }
-
+  delete[] storage;
   for (i = 0; i < nvars; ++i)
     point[i] = z[i];
   return minf;
 }
-#endif
-
-void OptInfoHooke::OptimiseLikelihood() {
+void OptInfoHooke::OptimiseLikelihoodREP() {
 
   double oldf, newf, bestf, steplength, tmp;
   int    i, offset;
@@ -289,11 +285,9 @@
   IntVector trapped(nvars, 0);
 
   EcoSystem->scaleVariables();
-#ifndef NO_OPENMP
   int numThr = omp_get_max_threads ( );
   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);
@@ -334,9 +328,7 @@
 
   while (1) {
     if (isZero(bestf)) {
-#ifdef NO_OPENMP
       iters = EcoSystem->getFuncEval() - offset;
-#endif
       handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
       converge = -1;
       return;
@@ -359,21 +351,15 @@
     /* 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);
+    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");
     	return;
     }
-#else
-    newf = this->bestNearby(delta, trialx, bestf, param);
-#endif
     /* if too many function evaluations occur, terminate the algorithm */
 
-#ifdef NO_OPENMP
     iters = EcoSystem->getFuncEval() - offset;
-#endif
     if (iters > hookeiter) {
       handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
       handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
@@ -441,9 +427,7 @@
       /* only move forward if this is really an improvement    */
       oldf = newf;
       newf = EcoSystem->SimulateAndUpdate(trialx);
-#ifndef NO_OPENMP
       iters++;
-#endif
       if ((isEqual(newf, oldf)) || (newf > oldf)) {
         newf = oldf;  //JMB no improvement, so reset the value of newf
         break;
@@ -454,23 +438,17 @@
       for (i = 0; i < nvars; i++)
         x[i] = trialx[i];
 
-#ifndef NO_OPENMP
-      newf = this->bestNearbyOMP(delta, trialx, bestf, param);
+      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");
           	return;
           }
-#else
-      newf = this->bestNearby(delta, trialx, bestf, param);
-#endif
       if (isEqual(newf, bestf))
         break;
 
       /* if too many function evaluations occur, terminate the algorithm */
-#ifdef NO_OPENMP
       iters = EcoSystem->getFuncEval() - offset;
-#endif
       if (iters > hookeiter) {
         handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
         handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
@@ -485,11 +463,227 @@
         EcoSystem->storeVariables(score, bestx);
         return;
       }
-    }
+    } // while (newf < bestf)
 
-#ifdef NO_OPENMP
     iters = EcoSystem->getFuncEval() - offset;
+    if (newf < bestf) {
+      for (i = 0; i < nvars; i++)
+        bestx[i] = x[i] * init[i];
+      bestf = newf;
+      handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
+      handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point");
+      EcoSystem->storeVariables(bestf, bestx);
+      EcoSystem->writeBestValues();
+
+    } else
+      handle.logMessage(LOGINFO, "Checking convergence criteria after", iters, "function evaluations ...");
+
+    /* if the step length is less than hookeeps, terminate the algorithm */
+    if (steplength < hookeeps) {
+      handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
+      handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
+      handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
+      handle.logMessage(LOGINFO, "The optimisation stopped because an optimum was found for this run");
+
+      converge = 1;
+      score = bestf;
+      handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
+      EcoSystem->storeVariables(bestf, bestx);
+      return;
+    }
+
+    steplength *= rho;
+    handle.logMessage(LOGINFO, "Reducing the steplength to", steplength);
+    for (i = 0; i < nvars; i++)
+      delta[i] *= rho;
+  }
+}
 #endif
+
+void OptInfoHooke::OptimiseLikelihood() {
+
+  double oldf, newf, bestf, steplength, tmp;
+  int    i, offset;
+  int    rchange, rcheck, rnumber;  //Used to randomise the order of the parameters
+
+  handle.logMessage(LOGINFO, "\nStarting Hooke & Jeeves optimisation algorithm\n");
+  int nvars = EcoSystem->numOptVariables();
+  DoubleVector x(nvars);
+  DoubleVector trialx(nvars);
+  DoubleVector bestx(nvars);
+  DoubleVector lowerb(nvars);
+  DoubleVector upperb(nvars);
+  DoubleVector init(nvars);
+  DoubleVector initialstep(nvars, rho);
+  DoubleVector delta(nvars);
+  IntVector param(nvars, 0);
+  IntVector lbound(nvars, 0);
+  IntVector rbounds(nvars, 0);
+  IntVector trapped(nvars, 0);
+
+  EcoSystem->scaleVariables();
+  EcoSystem->getOptScaledValues(x);
+  EcoSystem->getOptLowerBounds(lowerb);
+  EcoSystem->getOptUpperBounds(upperb);
+  EcoSystem->getOptInitialValues(init);
+
+  for (i = 0; i < nvars; i++) {
+    // Scaling the bounds, because the parameters are scaled
+    lowerb[i] = lowerb[i] / init[i];
+    upperb[i] = upperb[i] / init[i];
+    if (lowerb[i] > upperb[i]) {
+      tmp = lowerb[i];
+      lowerb[i] = upperb[i];
+      upperb[i] = tmp;
+    }
+
+    bestx[i] = x[i];
+    trialx[i] = x[i];
+    param[i] = i;
+    delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign
+  }
+
+  bestf = EcoSystem->SimulateAndUpdate(trialx);
+  if (bestf != bestf) { //check for NaN
+    handle.logMessage(LOGINFO, "Error starting Hooke & Jeeves optimisation with f(x) = infinity");
+    converge = -1;
+    iters = 1;
+    return;
+  }
+
+  offset = EcoSystem->getFuncEval();  //number of function evaluations done before loop
+  newf = bestf;
+  oldf = bestf;
+  steplength = lambda;
+  if (isZero(steplength))
+    steplength = rho;
+
+  iters = 0;
+
+  while (1) {
+    if (isZero(bestf)) {
+      handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
+      converge = -1;
+      return;
+    }
+
+    /* randomize the order of the parameters once in a while */
+    rchange = 0;
+    while (rchange < nvars) {
+      rnumber = rand() % nvars;
+      rcheck = 1;
+      for (i = 0; i < rchange; i++)
+        if (param[i] == rnumber)
+          rcheck = 0;
+      if (rcheck) {
+        param[rchange] = rnumber;
+        rchange++;
+      }
+    }
+
+    /* find best new point, one coord at a time */
+    for (i = 0; i < nvars; i++)
+      trialx[i] = x[i];
+    newf = this->bestNearby(delta, trialx, bestf, param);
+    /* if too many function evaluations occur, terminate the algorithm */
+
+    if (iters > hookeiter) {
+      handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
+      handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
+      handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
+      handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations");
+      handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
+
+      score = EcoSystem->SimulateAndUpdate(trialx);
+      handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
+      for (i = 0; i < nvars; i++)
+        bestx[i] = trialx[i] * init[i];
+      EcoSystem->storeVariables(score, bestx);
+      return;
+    }
+
+    /* if we made some improvements, pursue that direction */
+    while (newf < bestf) {
+      for (i = 0; i < nvars; i++) {
+        /* if it has been trapped but f has now gotten better (bndcheck) */
+        /* we assume that we are out of the trap, reset the counters     */
+        /* and go back to the stepsize we had when we got trapped        */
+        if ((trapped[i]) && (newf < oldf * bndcheck)) {
+          trapped[i] = 0;
+          lbound[i] = 0;
+          rbounds[i] = 0;
+          delta[i] = initialstep[i];
+
+        } else if (trialx[i] < (lowerb[i] + verysmall)) {
+          lbound[i]++;
+          trialx[i] = lowerb[i];
+          if (!trapped[i]) {
+            initialstep[i] = delta[i];
+            trapped[i] = 1;
+          }
+          /* if it has hit the bounds 2 times then increase the stepsize */
+          if (lbound[i] >= 2)
+            delta[i] /= rho;
+
+        } else if (trialx[i] > (upperb[i] - verysmall)) {
+          rbounds[i]++;
+          trialx[i] = upperb[i];
+          if (!trapped[i]) {
+            initialstep[i] = delta[i];
+            trapped[i] = 1;
+          }
+          /* if it has hit the bounds 2 times then increase the stepsize */
+          if (rbounds[i] >= 2)
+            delta[i] /= rho;
+        }
+      }
+
+      for (i = 0; i < nvars; i++) {
+        /* firstly, arrange the sign of delta[] */
+        if (trialx[i] < x[i])
+          delta[i] = 0.0 - fabs(delta[i]);
+        else
+          delta[i] = fabs(delta[i]);
+
+        /* now, move further in this direction  */
+        tmp = x[i];
+        x[i] = trialx[i];
+        trialx[i] = trialx[i] + trialx[i] - tmp;
+      }
+
+      /* only move forward if this is really an improvement    */
+      oldf = newf;
+      newf = EcoSystem->SimulateAndUpdate(trialx);
+      if ((isEqual(newf, oldf)) || (newf > oldf)) {
+        newf = oldf;  //JMB no improvement, so reset the value of newf
+        break;
+      }
+
+      /* OK, it's better, so update variables and look around  */
+      bestf = newf;
+      for (i = 0; i < nvars; i++)
+        x[i] = trialx[i];
+
+      if (isEqual(newf, bestf))
+        break;
+
+      /* if too many function evaluations occur, terminate the algorithm */
+      if (iters > hookeiter) {
+        handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
+        handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
+        handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
+        handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations");
+        handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
+
+        score = EcoSystem->SimulateAndUpdate(trialx);
+        handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
+        for (i = 0; i < nvars; i++)
+          bestx[i] = trialx[i] * init[i];
+        EcoSystem->storeVariables(score, bestx);
+        return;
+      }
+    } // while (newf < bestf)
+
     if (newf < bestf) {
       for (i = 0; i < nvars; i++)
         bestx[i] = x[i] * init[i];
@@ -524,8 +718,9 @@
 }
 
 /* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
-#ifdef GADGET_OPENMP
-double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
+#ifdef _OPENMP
+//#ifdef SPECULATIVE
+double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
   double minf;
   int i, j, k, ii;
   DoubleVector z(point);
@@ -549,43 +744,43 @@
 	  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) {
 		  if ((i + paral_tokens -1) >= nvars)
 			  paral_tokens = nvars - i;
-		  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) {
+	#pragma omp parallel for num_threads(paral_tokens*2) private(k)
+		  for (j = 0; j < (paral_tokens*2); ++j) {
 			  storage[j].z = z;
 			  storage[j].delta = delta;
-			  DoubleVector v1(z);
-			  DoubleVector v2(z);
-			  k = param[i+j];
-			  v1[k] +=  delta[k];
-			  v2[k] -=  delta[k];
+			  DoubleVector v(z);
 
-#pragma omp parallel sections num_threads(2)
-			  {
-	  #pragma omp section
-				  {
-			  storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
-				  }
-	  #pragma omp section
-				  {
-			  storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
-				  }
+			  if (j<paral_tokens) {
+				  k = param[i+j];
+				  v[k] +=  delta[k];
 			  }
+			  else {
+				  k = param[i+j-paral_tokens];
+				  v[k] -=  delta[k];
+			  }
+
+			  storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v);
+			  storage[j].z[k] = v[k];
+		  }
+
+		  for (j = 0; j < paral_tokens; ++j) {
+			  k = param[i+j];
 			  if (storage[j].ftmp < minf) {
 				  storage[j].iters = 1;
-				  storage[j].z[k] = v1[k];
+//				  storage[j].z[k] = v1[k];
 			  } else {
 				  storage[j].iters = 2;
 				  storage[j].delta[k] = 0.0 - delta[k];
 				  if (storage[j+paral_tokens].ftmp < minf) {
 					  storage[j].ftmp = storage[j+paral_tokens].ftmp;
-					  storage[j].z[k] = v2[k];
+					  storage[j].z[k] = storage[j+paral_tokens].z[k];
 				  }
 				  else iters += 2;
 			  }
@@ -605,6 +800,7 @@
 
 		  i += paral_tokens;
 	  }
+	  paral_tokens = numThr / 2;
 	}
 
   delete[] storage;
@@ -678,7 +874,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");
@@ -703,8 +899,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");
@@ -715,7 +911,7 @@
 	 #endif
 	     /* if too many function evaluations occur, terminate the algorithm */
 
-	 #ifdef NO_OPENMP
+	 #ifndef _OPENMP
 	     iters = EcoSystem->getFuncEval() - offset;
 	 #endif
 	     if (iters > hookeiter) {
@@ -785,7 +981,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)) {
@@ -798,8 +994,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");
@@ -812,7 +1008,7 @@
 	         break;
 
 	       /* if too many function evaluations occur, terminate the algorithm */
-	 #ifdef NO_OPENMP
+	 #ifndef _OPENMP
 	       iters = EcoSystem->getFuncEval() - offset;
 	 #endif
 	       if (iters > hookeiter) {
@@ -831,7 +1027,7 @@
 	       }
 	     }
 
-	 #ifdef NO_OPENMP
+	 #ifndef _OPENMP
 	     iters = EcoSystem->getFuncEval() - offset;
 	 #endif
 	     if (newf < bestf) {
@@ -866,4 +1062,5 @@
 	       delta[i] *= rho;
 	   }
   }
+//#endif
 #endif