[mareframe] Diff of /trunk/gadget/hooke.cc

Diff of /trunk/gadget/hooke.cc

-revision 19, Wed May 25 16:36:33 2016 UTC
+revision 20, Fri Apr  7 09:20:55 2017 UTC
 Line 263
      point[i] = z[i];
    return minf;
  }
- #endif
+ void OptInfoHooke::OptimiseLikelihoodREP() {
- void OptInfoHooke::OptimiseLikelihood() {
    double oldf, newf, bestf, steplength, tmp;
    int    i, offset;
-Line 287
+Line 285
    IntVector trapped(nvars, 0);
    EcoSystem->scaleVariables();
- #ifdef _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);
-Line 332
+Line 328
    while (1) {
      if (isZero(bestf)) {
- #ifndef _OPENMP
        iters = EcoSystem->getFuncEval() - offset;
- #endif
        handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
        converge = -1;
        return;
-Line 357
+Line 351
      /* find best new point, one coord at a time */
      for (i = 0; i < nvars; i++)
        trialx[i] = x[i];
- #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");
          return;
      }
- #else
-     newf = this->bestNearby(delta, trialx, bestf, param);
- #endif
      /* if too many function evaluations occur, terminate the algorithm */
- #ifndef _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");
-Line 439
+Line 427
        /* only move forward if this is really an improvement    */
        oldf = newf;
        newf = EcoSystem->SimulateAndUpdate(trialx);
- #ifdef _OPENMP
        iters++;
- #endif
        if ((isEqual(newf, oldf)) || (newf > oldf)) {
          newf = oldf;  //JMB no improvement, so reset the value of newf
          break;
-Line 452
+Line 438
        for (i = 0; i < nvars; i++)
          x[i] = trialx[i];
- #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");
                  return;
            }
- #else
-       newf = this->bestNearby(delta, trialx, bestf, param);
- #endif
        if (isEqual(newf, bestf))
          break;
        /* if too many function evaluations occur, terminate the algorithm */
- #ifndef _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");
-Line 485
+Line 465
        }
      } // while (newf < bestf)
- #ifndef _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];
-Line 522
+Line 718
  }
  /* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
- #ifdef SPECULATIVE
+ #ifdef _OPENMP
+ //#ifdef SPECULATIVE
  double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
    double minf;
    int i, j, k, ii;
-Line 865
+Line 1062
                 delta[i] *= rho;
             }
    }
+ //#endif
  #endif

 Legend:



Removed from v.19
 


changed lines


 
Added in v.20
 Legend:



Removed from v.19
 


changed lines


 
Added in v.20
-Removed from v.19
+Added in v.20

root@forge.cesga.es	ViewVC Help
Powered by ViewVC 1.0.0