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

Diff of /trunk/gadget/hooke.cc

-revision 11, Thu Jul 23 19:00:38 2015 UTC
+revision 16, Tue Oct 13 11:29:01 2015 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 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 178
+Line 170
        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 187
+Line 178
  }
  /* given a point, look for a better one nearby, one coord at a time */
- #ifndef NO_OPENMP
+ #ifdef _OPENMP
- double OptInfoHooke::bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
+ /*
- //cout << "beastNearbyINI" << endl;
+  * 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);
-Line 214
+Line 209
            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;
+ #pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2)
-           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);
-Line 235
+Line 224
                    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
                            {
                    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;
-Line 278
+Line 256
                            minf = storage[j].ftmp;
                            z = storage[j].z;
                            delta = storage[j].delta;
- //                        cout << storage[j].ftmp;
                            break;
                    }
            }
    }
+   delete[] storage;
- //  omp_set_num_threads(numThr);
    for (i = 0; i < nvars; ++i)
      point[i] = z[i];
    return minf;
-Line 314
+Line 289
    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);
-Line 337
+Line 312
      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 360
+Line 334
    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 371
+Line 345
      /* 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 386
+Line 359
      /* 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->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");
-Line 398
+Line 371
  #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");
-Line 469
+Line 441
        /* 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 482
+Line 454
        for (i = 0; i < nvars; i++)
          x[i] = trialx[i];
- #ifndef NO_OPENMP
+ #ifdef _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");
-Line 496
+Line 468
          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");
-Line 514
+Line 485
          EcoSystem->storeVariables(score, bestx);
          return;
        }
-     }
+     } // while (newf < bestf)
- #ifdef NO_OPENMP
+ #ifndef _OPENMP
      iters = EcoSystem->getFuncEval() - offset;
  #endif
      if (newf < bestf) {
-Line 552
+Line 523
    }
  }
- #ifdef GADGET_OPENMP
+ /* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
-   //TODO doc
+ #ifdef SPECULATIVE
- double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
+ double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
- //cout << "beastNearbyINI" << endl;
+   double minf;
-   double minf;//, ftmp;
    int i, j, k, ii;
    DoubleVector z(point);
    int bestId = 0;
-Line 569
+Line 539
    minf = prevbest;
    int paral_tokens, numThr, nvars = point.Size();
    numThr = omp_get_max_threads ( );
-Line 580
+Line 549
            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);
-Line 605
+Line 568
          #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)
-Line 649
+Line 601
                            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];
-Line 686
+Line 635
             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);
-Line 709
+Line 656
               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 732
+Line 678
             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 743
+Line 689
               /* 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 758
+Line 703
               /* 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 770
+Line 715
           #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");
-Line 786
+Line 730
                 for (i = 0; i < nvars; i++)
                   bestx[i] = trialx[i] * init[i];
                 EcoSystem->storeVariables(score, bestx);
- //             EcoSystem->writeOptValuesOMP();
                 return;
               }
-Line 842
+Line 785
                 /* 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 855
+Line 798
                 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 869
+Line 812
                   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");
-Line 885
+Line 827
                   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) {
-Line 916
+Line 857
                 score = bestf;
                 handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
                 EcoSystem->storeVariables(bestf, bestx);
- //             EcoSystem->writeOptValuesOMP();
                 return;
               }

 Legend:



Removed from v.11
 


changed lines


 
Added in v.16
 Legend:



Removed from v.11
 


changed lines


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

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