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

Diff of /trunk/gadget/hooke.cc

-revision 14, Thu Jul 30 12:36:46 2015 UTC
+revision 20, Fri Apr  7 09:20:55 2017 UTC
 Line 139
  #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
 Line 178
  }
  /* 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);
 Line 209
            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);
+ #pragma omp parallel for num_threads(paral_tokens*2) private(k) //parallelize the parameters (numThr)
-           omp_set_nested(1); //permit the nested parallelization
+           for (j = 0; j < (paral_tokens*2); ++j) {
- #pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2)
-           for (j = 0; j < paral_tokens; ++j) {
                    storage[j].z = z;
                    storage[j].delta = delta;
-                   DoubleVector v1(z);
+                   DoubleVector v(z);
-                   DoubleVector v2(z);
-                   k = param[i+j];
-                   v1[k] +=  delta[k];
-                   v2[k] -=  delta[k];
- #pragma omp parallel sections num_threads(2) //parrallelize the +/- delta simulation for each parameter
+                   if (j<paral_tokens) {
-                   {
+                           k = param[i+j];
-         #pragma omp section
+                           v[k] +=  delta[k];
-                           {
-                   storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
-                           }
-         #pragma omp section
-                           {
-                   storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
                            }
+                   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];;
                            }
                    }
            }
-Line 260
+Line 258
                    }
            }
    }
+   delete[] storage;
    for (i = 0; i < nvars; ++i)
      point[i] = z[i];
    return minf;
  }
- #endif
+ void OptInfoHooke::OptimiseLikelihoodREP() {
- void OptInfoHooke::OptimiseLikelihood() {
    double oldf, newf, bestf, steplength, tmp;
    int    i, offset;
-Line 289
+Line 285
    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);
-Line 334
+Line 328
    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 359
+Line 351
      /* find best new point, one coord at a time */
      for (i = 0; i < nvars; i++)
        trialx[i] = x[i];
- #ifndef NO_OPENMP
+     newf = this->bestNearbyRepro(delta, trialx, bestf, param);
-     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
      if (iters > hookeiter) {
        handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
        handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
-Line 441
+Line 427
        /* 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 454
+Line 438
        for (i = 0; i < nvars; i++)
          x[i] = trialx[i];
- #ifndef NO_OPENMP
+       newf = this->bestNearbyRepro(delta, trialx, bestf, param);
-       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
        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 463
          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];
-Line 524
+Line 718
  }
  /* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
- #ifdef GADGET_OPENMP
+ #ifdef _OPENMP
- double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
+ //#ifdef SPECULATIVE
+ double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
    double minf;
    int i, j, k, ii;
    DoubleVector z(point);
-Line 549
+Line 744
            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);
+         #pragma omp parallel for num_threads(paral_tokens*2) private(k)
-                   omp_set_nested(1);
+                   for (j = 0; j < (paral_tokens*2); ++j) {
-         #pragma omp parallel for num_threads(paral_tokens) private(k)
-                   for (j = 0; j < paral_tokens; ++j) {
                            storage[j].z = z;
                            storage[j].delta = delta;
-                           DoubleVector v1(z);
+                           DoubleVector v(z);
-                           DoubleVector v2(z);
-                           k = param[i+j];
-                           v1[k] +=  delta[k];
-                           v2[k] -=  delta[k];
- #pragma omp parallel sections num_threads(2)
+                           if (j<paral_tokens) {
-                           {
+                                   k = param[i+j];
-           #pragma omp section
+                                   v[k] +=  delta[k];
-                                   {
-                           storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
-                                   }
-           #pragma omp section
-                                   {
-                           storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
                                    }
+                           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;
                            }
-Line 605
+Line 800
                    i += paral_tokens;
            }
+           paral_tokens = numThr / 2;
          }
    delete[] storage;
-Line 678
+Line 874
             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");
-Line 703
+Line 899
               /* find best new point, one coord at a time */
               for (i = 0; i < nvars; i++)
                 trialx[i] = x[i];
-          #ifndef NO_OPENMP
+          #ifdef _OPENMP
-              newf = this->bestNearbyOMP2(delta, trialx, bestf, param);
+              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");
-Line 715
+Line 911
           #endif
               /* if too many function evaluations occur, terminate the algorithm */
-          #ifdef NO_OPENMP
+          #ifndef _OPENMP
               iters = EcoSystem->getFuncEval() - offset;
           #endif
               if (iters > hookeiter) {
-Line 785
+Line 981
                 /* 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)) {
-Line 798
+Line 994
                 for (i = 0; i < nvars; i++)
                   x[i] = trialx[i];
-          #ifndef NO_OPENMP
+          #ifdef _OPENMP
-                newf = this->bestNearbyOMP2(delta, trialx, bestf, param);
+                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");
-Line 812
+Line 1008
                   break;
                 /* if too many function evaluations occur, terminate the algorithm */
-          #ifdef NO_OPENMP
+          #ifndef _OPENMP
                 iters = EcoSystem->getFuncEval() - offset;
           #endif
                 if (iters > hookeiter) {
-Line 831
+Line 1027
                 }
               }
-          #ifdef NO_OPENMP
+          #ifndef _OPENMP
               iters = EcoSystem->getFuncEval() - offset;
           #endif
               if (newf < bestf) {
-Line 866
+Line 1062
                 delta[i] *= rho;
             }
    }
+ //#endif
  #endif

 Legend:



Removed from v.14
 


changed lines


 
Added in v.20
 Legend:



Removed from v.14
 


changed lines


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

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