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

Diff of /trunk/gadget/simann.cc

-revision 11, Thu Jul 23 19:00:38 2015 UTC
+revision 20, Fri Apr  7 09:20:55 2017 UTC
 Line 178
  #include "ecosystem.h"
  #include "global.h"
  #include "seq_optimize_template.h"
- #ifdef GADGET_OPENMP
+ #ifdef SPECULATIVE
  #include <omp.h>
  #endif
  extern Ecosystem* EcoSystem;
- #ifndef NO_OPENMP
+ #ifdef _OPENMP
  extern Ecosystem** EcoSystems;
  #endif
- //extern StochasticData* data;
+ /*sequential code replaced at seq_optimize_template.h*/
  //void OptInfoSimann::OptimiseLikelihood() {
  //
  //  //set initial values
-Line 305
+Line 304
  ////              trialx[i] = x[i];
  ////          }
  //
- ////          cout << "param:" << param[l] << "-" << trialx[param[l]] << endl;
- //
  //          //Evaluate the function with the trial point trialx and return as -trialf
  //          trialf = EcoSystem->SimulateAndUpdate(trialx);
  //          trialf = -trialf;
-Line 327
+Line 324
  //            handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score);
  //            return;
  //          }
- ////          cout << "f:" << trialf << endl;
  //          //Accept the new point if the new function value better
- ////          cout << "mustAccept:" << trialf << "|" << funcval<< endl;
  //          if ((trialf - funcval) > verysmall) {
  //            for (i = 0; i < nvars; i++)
  //              x[i] = trialx[i];
-Line 353
+Line 348
  //              nrej++;
  //            }
  //          }
- ////          cout << "goog VALUE:" << funcval << endl;
  //          // JMB added check for really silly values
  //          if (isZero(trialf)) {
  //            handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0");
-Line 442
+Line 436
  //}
- #ifdef GADGET_OPENMP
+ #ifdef _OPENMP
+ //#ifdef SPECULATIVE
  void OptInfoSimann::OptimiseLikelihoodOMP() {
    //set initial values
-Line 452
+Line 447
    double tmp, p, pp, ratio, nsdiv;
    double fopt, funcval, trialf;
-   int    a, i, j, k, l, offset, quit;
+   int    a, i, j, k, l, quit;
    int    rchange, rcheck, rnumber;  //Used to randomise the order of the parameters
+   // store the info of the different threads
    struct Storage {
            DoubleVector trialx;
            double newLikelihood;
-Line 475
+Line 471
    IntVector nacp(nvars, 0);
    EcoSystem->resetVariables();  //JMB need to reset variables in case they have been scaled
-   if (scale)
+   if (scale) {
      EcoSystem->scaleVariables();
+                 int numThr = omp_get_max_threads ( );
+                 for(i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
+                         EcoSystems[i]->scaleVariables();
+         }
    EcoSystem->getOptScaledValues(x);
    EcoSystem->getOptLowerBounds(lowerb);
    EcoSystem->getOptUpperBounds(upperb);
 Line 507
      handle.logMessage(LOGINFO, "Error starting Simulated Annealing optimisation with f(x) = infinity");
      converge = -1;
      iters = 1;
- //    EcoSystem->writeOptValuesOMP();
      return;
    }
    //the function is to be minimised so switch the sign of funcval (and trialf)
    funcval = -funcval;
-   offset = EcoSystem->getFuncEval();  //number of function evaluations done before loop
    nacc++;
    cs /= lratio;  //JMB save processing time
    nsdiv = 1.0 / ns;
-Line 532
+Line 530
    for (i=0; i<numThr; ++i)
            storage[i].trialx = trialx;
-   //Start the main loop.  Note that it terminates if
-   //(i) the algorithm succesfully optimises the function or
-   //(ii) there are too many function evaluations
-   int aux;
-  DoubleVector vns(nvars, 0);
+   int aux; //store the number of evaluations that are not useful for the algorithm
+  DoubleVector vns(nvars, 0); //vector of ns
   int ns_ = ceil(numThr/2.);
   double res;
    aux=0;
+   //Start the main loop.  Note that it terminates if
+   //(i) the algorithm succesfully optimises the function or
+   //(ii) there are too many function evaluations
    while (1) {
      for (a = 0; a < nt; ++a) {
        //Randomize the order of the parameters once in a while, to avoid
-Line 652
+Line 651
                x[i] = storage[bestId].trialx[i];
              funcval = trialf;
              nacc++;
- //            vns[param[l+bestId]]++;
-             //nacp[param[l+bestId]]+=numThr;  //JMB - not sure about this ...
- //            for (i = 0; i < numThr; ++i)
- //              nacp[param[l+i]]++;
- //            nacp[param[l+bestId]]++;
            } else {
              //Accept according to metropolis condition
              p = expRep((trialf - funcval) / t);
-Line 669
+Line 662
                funcval = trialf;
                naccmet++;
                nacp[param[(l+bestId)%nvars]]++;
- //              vns[param[l+bestId]]++;
- //              nacp[param[l+bestId]]+=numThr;
- //              for (i = 0; i < numThr; ++i)
- //                nacp[param[l+i]]++;
- //                nacp[param[l+bestId]]++;
              } else {
                //Reject point
                nrej++;
-Line 683
+Line 671
            if (isZero(trialf)) {
              handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0");
              converge = -1;
- //            EcoSystem->writeOptValuesOMP();
              delete[] storage;
              return;
            }
-Line 697
+Line 684
              if (scale) {
                for (i = 0; i < nvars; ++i)
                  scalex[i] = bestx[i] * init[i];
-               //FIXME EcoSystems[]->
                EcoSystem->storeVariables(-fopt, scalex);
              } else
                  EcoSystem->storeVariables(-fopt, bestx);
-Line 708
+Line 694
            }
          }
        }
-       //Adjust vm so that approximately half of all evaluations are accepted
- //      for (i = 0; i < nvars; ++i) {
- ////              cout << "nacp["<<i<<"]:"<<nacp[i]<<endl;
- //        ratio = nsdiv * nacp[i];
- //        nacp[i] = 0;
- //        if (ratio > uratio) {
- //          vm[i] = vm[i] * (1.0 + cs * (ratio - uratio));
- //        } else if (ratio < lratio) {
- //          vm[i] = vm[i] / (1.0 + cs * (lratio - ratio));
- //        }
- //
- //        if (vm[i] < rathersmall)
- //          vm[i] = rathersmall;
- //        if (vm[i] > (upperb[i] - lowerb[i]))
- //          vm[i] = upperb[i] - lowerb[i];
- //
- ////        cout << "vm["<<i<<"]:"<<vm[i]<<endl;
- //
- //      }
      }
      //Check termination criteria
-Line 758
+Line 724
        converge = 1;
        score = EcoSystem->SimulateAndUpdate(bestx);
        handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score);
- //      EcoSystem->writeOptValuesOMP();
        delete[] storage;
        return;
      }
-Line 773
+Line 738
      for (i = 0; i < nvars; ++i)
        x[i] = bestx[i];
    }
- //  EcoSystem->writeOptValuesOMP();
  }
+ //#endif
  #endif
+ // calcule a new point
  void OptInfoSimann::newValue(int nvars, int l, IntVector& param, DoubleVector& trialx,
                  DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm)
  {
          int i, k;
- //      double old = trialx[param[l]];
- //      cout << "[" << endl;
            for (i = 0; i < nvars; ++i) {
                  if (i == param[l]) {
                    trialx[i] = x[i] + ((randomNumber(&seed) * 2.0) - 1.0) * vm[i];
-Line 800
+Line 764
                    }
                  } else
                    trialx[i] = x[i];
- //              cout <<" " << trialx[i];
            }
- //        cout << endl<< "old[" <<param[l] << "]:" << old <<"-" << trialx[param[l]]<< " - " << vm[param[l]]<<endl;
  }
+ /*################################################################################
+  * code use in the template (seq_optimize_template.h)
+  ################################################################################*/
  //Generate trialx, the trial value of x
  void newValue(int nvars, int l, IntVector& param, DoubleVector& trialx,
                  DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm, unsigned* seed)
  {
          int i, k;
- //      double old = trialx[param[l]];
- //      cout << "[" << endl;
            for (i = 0; i < nvars; ++i) {
                  if (i == param[l]) {
                    trialx[i] = x[i] + ((randomNumber(&*seed) * 2.0) - 1.0) * vm[i];
-Line 831
+Line 794
                    }
                  } else
                    trialx[i] = x[i];
- //              cout  << trialx[i]<<" ";
            }
- //        cout << endl <<l<< "-old[" <<param[l] << "]:" << old <<"-" << trialx[param[l]]<< " - " << vm[param[l]]<<endl;
  }
  void buildNewParams_f(Siman& seed, DoubleVector& params) {
-Line 845
+Line 805
  }
  /// Represents the function that computes how good the parameters are
- double evaluate_f(const DoubleVector& params) {
+ #ifdef _OPENMP
+ double evaluate_par_f(const DoubleVector& params) {
          double trialf;
- #ifndef NO_OPENMP
          int id = omp_get_thread_num ();
          trialf = EcoSystems[id]->SimulateAndUpdate(params);
- //      cout << id << "-";
+   return -trialf;
- #else
+ }
-         trialf = EcoSystem->SimulateAndUpdate(params);
  #endif
- //      cout << "trailf:" << trialf<< endl;
+ double evaluate_f(const DoubleVector& params) {
+         double trialf;
+         trialf = EcoSystem->SimulateAndUpdate(params);
    return -trialf;
  }
-Line 872
+Line 833
                  double uratio = seed.getUratio();
                  double lratio = seed.getLratio();
-         //      cout << "uratio:" << uratio << endl;
-         //      cout << "lratio:" << lratio << endl;
-         //cout << "adjust vm" << endl;
                  for (i = 0; i < seed.getNvars(); i++) {
- //                      cout << "nacp[" << i << "]:" << seed.getNacp(i) << endl;
                          ratio = nsdiv * seed.getNacp(i);
                          seed.setNacp(i,0);
                          if (ratio > uratio) {
-Line 906
+Line 862
                  DoubleVector* bestx = seed.getBestx();
-                 //funcval = fopt;
                  for (i = 0; i < seed.getNvars(); i++)
                          x[i] = (*bestx)[i];
          }
-Line 945
+Line 900
            //Accept the new point if the new function value better
            bool ret = true;
            int i;
- //        cout << "mustAccept:" << trialf << "|" << funcval<< endl;
            int aux = siman.getParam()[siman.getL()];
          if ((trialf - funcval) > verysmall) {
                  siman.incrementNacp(aux);
-Line 955
+Line 909
                  double p, pp;
                  //Accept according to metropolis condition
- //              cout << "t:" << siman.getT() << endl;
                  p = expRep((trialf - funcval) / siman.getT());
                  pp = randomNumber(siman.getSeedM());
                  if (pp < p) {
- //                      cout << pp << "<" << p << endl;
                          siman.incrementNacp(aux);
                          siman.incrementNaccmet();
  //                      //Accept point
- //                      for (i = 0; i < nvars; i++)
- //                              x[i] = trialx[i];
- //                      funcval = trialf;
- //                      naccmet++;
- //                      nacp[param[l]]++;
                  } else {
                          //Reject point
                          ret = false;
                          siman.incrementNrej();
                  }
          }
- //      cout << "nacp[" << aux << "]:" << siman.getNacp(aux) << endl;
          // JMB added check for really silly values
          if (isZero(trialf)) {
                  handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after",
                                  iters, "function evaluations, f(x) = 0");
                  siman.setConverge(-1);
                  return false;
          }
-Line 989
+Line 933
          if (siman.getL() == 0)
                  siman.incrementNS();
          if (siman.getNS() >= siman.getNs()){
                  siman.setNS(0);
                  siman.incrementNT();
                  adjustVm(siman);
                  siman.randomizeParams();
- //              if (seed.getNT() >= seed.getNt()){
- //
- //              }
          }
          return ret;
-Line 1015
+Line 953
           bool quit = false;
            if (siman.getNT() >= siman.getNt())
            {
- //                cout << "nt:" << siman.getNT();
-                   //Check termination criteria
                    int i;
                    DoubleVector fstar = siman.getFstar();
-Line 1067
+Line 1003
          handle.logMessage(LOGINFO, "Number of rejected points", siman.getNrej());
          *score = EcoSystem->SimulateAndUpdate(*bestX);
          handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", *score);
- //      EcoSystem->writeOptValues();
    }
  };
  // Required
  std::ostream &operator<<(std::ostream &os, const DoubleVector &p)
  {
- //  os <<  "[ AUCH";
- //  for (int i = 0; i< NPARAMS; ++i) {
- //    std::cout << p[i] << ' ';
- //  }
- //  std::cout << "]";
          os << "";
    return os;
  }
- void OptInfoSimann::OptimiseLikelihood() {
+ #ifdef _OPENMP
+ void OptInfoSimann::OptimiseLikelihoodREP() {
          //set initial values
          double tmp, p, pp;
-         double fopt, funcval, trialf;
+         double funcval, trialf;
          int a, i, j, k, l, quit;
          int rchange, rcheck, rnumber; //Used to randomise the order of the parameters
          handle.logMessage(LOGINFO,
-                         "\nStarting Simulated Annealing optimisation algorithm----\n");
+                         "\nStarting Simulated Annealing optimisation algorithm\n");
          int nvars = EcoSystem->numOptVariables();
          DoubleVector x(nvars);
          DoubleVector init(nvars);
-Line 1110
+Line 1041
          EcoSystem->resetVariables(); //JMB need to reset variables in case they have been scaled
          if (scale) {
                  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);
-Line 1152
+Line 1082
          //the function is to be minimised so switch the sign of funcval (and trialf)
          funcval = -funcval;
- //      offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop
- //      nacc++;
          cs /= lratio;  //JMB save processing time
- //      nsdiv = 1.0 / ns;
-         fopt = funcval;
          for (i = 0; i < tempcheck; i++)
                  fstar[i] = funcval;
-         //unsigned seed = 1234;
          Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns),
                          tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score);
-         ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_f, buildNewParams_f>
+         ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_par_f, buildNewParams_f>
          pa(s, x, simanniter);
- //      cout << "Sequential run" << endl;
+         // OpenMP parallelization
- //      pa.seq_opt();
- #ifndef NO_OPENMP
          int numThr = omp_get_max_threads ( );
-         pa.paral_opt_omp(numThr,numThr);
+         pa.paral_opt_omp(funcval,numThr,numThr);
- #else
+         iters = pa.iterations();
-         pa.seq_opt();
+ }
  #endif
+ void OptInfoSimann::OptimiseLikelihood() {
+         //set initial values
+         double tmp, p, pp;
+         double funcval, trialf;
+         int a, i, j, k, l, quit;
+         int rchange, rcheck, rnumber; //Used to randomise the order of the parameters
+         handle.logMessage(LOGINFO,
+                         "\nStarting Simulated Annealing optimisation algorithm\n");
+         int nvars = EcoSystem->numOptVariables();
+         DoubleVector x(nvars);
+         DoubleVector init(nvars);
+         DoubleVector trialx(nvars, 0.0);
+         DoubleVector bestx(nvars);
+         DoubleVector scalex(nvars);
+         DoubleVector lowerb(nvars);
+         DoubleVector upperb(nvars);
+         DoubleVector fstar(tempcheck);
+         DoubleVector vm(nvars, vminit);
+         IntVector param(nvars, 0);
+         EcoSystem->resetVariables(); //JMB need to reset variables in case they have been scaled
+         if (scale) {
+                 EcoSystem->scaleVariables();
+         }
+         EcoSystem->getOptScaledValues(x);
+         EcoSystem->getOptLowerBounds(lowerb);
+         EcoSystem->getOptUpperBounds(upperb);
+         EcoSystem->getOptInitialValues(init);
+         for (i = 0; i < nvars; i++) {
+                 bestx[i] = x[i];
+                 param[i] = i;
  }
+         if (scale) {
+                 for (i = 0; i < nvars; i++) {
+                         scalex[i] = x[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;
+                         }
+                 }
+         }
+         //funcval is the function value at x
+         funcval = EcoSystem->SimulateAndUpdate(x);
+         if (funcval != funcval) { //check for NaN
+                 handle.logMessage(LOGINFO,
+                                 "Error starting Simulated Annealing optimisation with f(x) = infinity");
+                 converge = -1;
+                 iters = 1;
+                 return;
+         }
+         //the function is to be minimised so switch the sign of funcval (and trialf)
+         funcval = -funcval;
+         cs /= lratio;  //JMB save processing time
+         for (i = 0; i < tempcheck; i++)
+                 fstar[i] = funcval;
+         Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns),
+                         tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score);
+         ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_f, buildNewParams_f>
+         pa(s, x, simanniter);
+         // sequential code
+         pa.seq_opt(funcval);
+         iters = pa.iterations();
+ }

 Legend:



Removed from v.11
 


changed lines


 
Added in v.20
 Legend:



Removed from v.11
 


changed lines


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

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