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

Diff of /trunk/gadget/hooke.cc

-revision 10, Wed Jun 10 13:01:42 2015 UTC
+revision 11, Thu Jul 23 19:00:38 2015 UTC
 Line 1
  /* Nonlinear Optimization using the algorithm of Hooke and Jeeves  */
  /*  12 February 1994    author: Mark G. Johnson                    */
+ //
  /* Find a point X where the nonlinear function f(X) has a local    */
  /* minimum.  X is an n-vector and f(X) is a scalar.  In mathe-     */
  /* matical notation  f: R^n -> R^1.  The objective function f()    */
  /* is not required to be continuous.  Nor does f() need to be      */
  /* differentiable.  The program does not use or require            */
  /* derivatives of f().                                             */
+ //
  /* The software user supplies three things: a subroutine that      */
  /* computes f(X), an initial "starting guess" of the minimum point */
  /* X, and values for the algorithm convergence parameters.  Then   */
  /* the program searches for a local minimum, beginning from the    */
  /* starting guess, using the Direct Search algorithm of Hooke and  */
  /* Jeeves.                                                         */
+ //
  /* This C program is adapted from the Algol pseudocode found in    */
  /* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun-  */
  /* ications of the ACM, Vol 6. p.313 (June 1963).  It includes the */
 Line 24
  /* highly as the one by A. Kaupe, is:  R. Hooke and T. A. Jeeves,  */
  /* "Direct Search Solution of Numerical and Statistical Problems", */
  /* Journal of the ACM, Vol. 8, April 1961, pp. 212-229.            */
+ //
  /* Calling sequence:                                               */
  /*  int hooke(nvars, startpt, endpt, rho, epsilon, itermax)        */
  /*                                                                 */
 Line 59
  /*     itermax     {an integer}  A second, rarely used, halting    */
  /*         criterion.  If the algorithm uses >= itermax            */
  /*         iterations, halt.                                       */
+ //
  /* The user-supplied objective function f(x,n) should return a C   */
  /* "double".  Its  arguments are  x -- an array of doubles, and    */
  /* n -- an integer.  x is the point at which f(x) should be        */
  /* evaluated, and n is the number of coordinates of x.  That is,   */
  /* n is the number of coefficients being fitted.                   */
+ //
  /*             rho, the algorithm convergence control              */
+ //
  /*  The algorithm works by taking "steps" from one estimate of     */
  /*    a minimum, to another (hopefully better) estimate.  Taking   */
  /*    big steps gets to the minimum more quickly, at the risk of   */
 Line 95
  /*    again with a larger value of rho such as 0.85, using the     */
  /*    result of the first minimization as the starting guess to    */
  /*    begin the second minimization.                               */
+ //
  /*                        Normal use:                              */
  /*    (1) Code your function f() in the C language                 */
  /*    (2) Install your starting guess {or read it in}              */
  /*    (3) Run the program                                          */
  /*    (4) {for the skeptical}: Use the computed minimum            */
  /*        as the starting point for another run                    */
+ //
  /* Data Fitting:                                                   */
  /*  Code your function f() to be the sum of the squares of the     */
  /*  errors (differences) between the computed values and the       */
 Line 113
  /*  f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i]))     */
  /*                + (C*tan(t[i]))))^2                              */
  /*  where x[] is a 3-vector consisting of {A, B, C}.               */
+ //
  /*  The author of this software is M.G. Johnson.                   */
  /*  Permission to use, copy, modify, and distribute this software  */
  /*  for any purpose without fee is hereby granted, provided that   */
 Line 125
  /*  AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY     */
  /*  KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS    */
  /*  FITNESS FOR ANY PARTICULAR PURPOSE.                            */
+ //
  /* JMB this has been modified to work with the gadget object structure   */
  /* This means that the function has been replaced by a call to ecosystem */
  /* object, and we can use the vector objects that have been defined      */
 Line 139
  #include "ecosystem.h"
  #include "global.h"
- extern Ecosystem* EcoSystem;
+ #ifndef NO_OPENMP
+ #include "omp.h"
+ #endif
+ extern Ecosystem* EcoSystem;
+ #ifndef NO_OPENMP
+ extern Ecosystem** EcoSystems;
+ #endif
  /* given a point, look for a better one nearby, one coord at a time */
  double OptInfoHooke::bestNearby(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
-Line 150
+Line 156
    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 {
-Line 164
+Line 178
        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++)
-Line 171
+Line 186
    return minf;
  }
+ /* 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;
+   double minf;//, ftmp;
+   int i, j, k;
+   DoubleVector z(point);
+   struct Storage {
+           DoubleVector z;
+           DoubleVector delta;
+           double ftmp;
+           int iters;
+   };
+   minf = prevbest;
+   i = 0;
+   int paral_tokens, numThr, nvars = point.Size();
+   numThr = omp_get_max_threads ( );
+   Storage* storage = new Storage[numThr];
+   if ((numThr % 2) == 0)
+           paral_tokens = numThr / 2;
+   else {
+           return -1;
+   }
+   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);
+                   DoubleVector v2(z);
+                   k = param[i+j];
+                   v1[k] +=  delta[k];
+                   v2[k] -=  delta[k];
+ #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
+ //                                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;
+                   if (storage[j].ftmp < minf) {
+                           minf = storage[j].ftmp;
+                           z = storage[j].z;
+                           delta = storage[j].delta;
+ //                        cout << storage[j].ftmp;
+                           break;
+                   }
+           }
+   }
+ //  omp_set_num_threads(numThr);
+   for (i = 0; i < nvars; ++i)
+     point[i] = z[i];
+   return minf;
+ }
+ #endif
  void OptInfoHooke::OptimiseLikelihood() {
    double oldf, newf, bestf, steplength, tmp;
-Line 193
+Line 314
    IntVector trapped(nvars, 0);
    EcoSystem->scaleVariables();
+ #ifndef NO_OPENMP
+   int numThr = omp_get_max_threads ( );
+   for (i = 0; i < numThr; i++)
+           EcoSystems[i]->scaleVariables();
+ #endif
    EcoSystem->getOptScaledValues(x);
    EcoSystem->getOptLowerBounds(lowerb);
    EcoSystem->getOptUpperBounds(upperb);
-Line 211
+Line 337
      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
    }
-Line 229
+Line 356
    if (isZero(steplength))
      steplength = rho;
+   iters = 0;
    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;
-Line 240
+Line 371
      /* 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++)
-Line 254
+Line 386
      /* 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);
+     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
+ //    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");
+       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);
+ #ifndef NO_OPENMP
+       iters++;
+ #endif
+       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];
+ #ifndef NO_OPENMP
+       newf = this->bestNearbyOMP(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
+ //      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");
-Line 272
+Line 514
        EcoSystem->storeVariables(score, bestx);
        return;
      }
+     }
+ #ifdef NO_OPENMP
+     iters = EcoSystem->getFuncEval() - offset;
+ #endif
+     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;
+   }
+ }
+ #ifdef GADGET_OPENMP
+   //TODO doc
+ double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
+ //cout << "beastNearbyINI" << endl;
+   double minf;//, ftmp;
+   int i, j, k, ii;
+   DoubleVector z(point);
+   int bestId = 0;
+   struct Storage {
+           DoubleVector z;
+           DoubleVector delta;
+           double ftmp;
+           int iters;
+   };
+   minf = prevbest;
+   int paral_tokens, numThr, nvars = point.Size();
+   numThr = omp_get_max_threads ( );
+   Storage* storage = new Storage[numThr];
+   if ((numThr % 2) == 0)
+           paral_tokens = numThr / 2;
+   else {
+           return -1;
+   }
+   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);
+                           DoubleVector v2(z);
+                           k = param[i+j];
+                           v1[k] +=  delta[k];
+                           v2[k] -=  delta[k];
+         #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)
+                                   bestId = j;
+                   }
+                   if (storage[bestId].ftmp < minf) {
+                           iters += storage[bestId].iters;
+                           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];
+   return minf;
+ }
+   void OptInfoHooke::OptimiseLikelihoodOMP() {
+           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();
+          #ifndef NO_OPENMP
+            int numThr = omp_get_max_threads ( );
+            for (i = 0; i < numThr; i++)
+                   EcoSystems[i]->scaleVariables();
+          #endif
+            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;
+              //FIXME rand_r
+              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)) {
+          #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;
+              }
+              /* 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++)
+                  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];
+          #ifndef NO_OPENMP
+              newf = this->bestNearbyOMP2(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
+          //    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");
+                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);
+ //             EcoSystem->writeOptValuesOMP();
+                return;
+              }
      /* if we made some improvements, pursue that direction */
      while (newf < bestf) {
-Line 325
+Line 842
        /* 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;
-Line 334
+Line 854
        bestf = newf;
        for (i = 0; i < nvars; i++)
          x[i] = trialx[i];
+          #ifndef NO_OPENMP
+                newf = this->bestNearbyOMP2(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
+          //      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");
-Line 352
+Line 885
          for (i = 0; i < nvars; i++)
            bestx[i] = trialx[i] * init[i];
          EcoSystem->storeVariables(score, bestx);
+ //               EcoSystem->writeOptValuesOMP();
          return;
        }
      }
+          #ifdef NO_OPENMP
      iters = EcoSystem->getFuncEval() - offset;
+          #endif
      if (newf < bestf) {
        for (i = 0; i < nvars; i++)
          bestx[i] = x[i] * init[i];
-Line 380
+Line 916
        score = bestf;
        handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
        EcoSystem->storeVariables(bestf, bestx);
+ //             EcoSystem->writeOptValuesOMP();
        return;
      }
-Line 389
+Line 926
        delta[i] *= rho;
    }
  }
+ #endif

 Legend:



Removed from v.10
 


changed lines


 
Added in v.11
 Legend:



Removed from v.10
 


changed lines


 
Added in v.11
-Removed from v.10
+Added in v.11

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