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

Diff of /trunk/gadget/simann.cc

-revision 11, Thu Jul 23 19:00:38 2015 UTC
+revision 14, Thu Jul 30 12:36:46 2015 UTC
 Line 186
  #ifndef NO_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 452
+Line 446
    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 470
    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
+Line 506
      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 529
    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 650
                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 661
                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 670
            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 683
              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 693
            }
          }
        }
-       //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 723
        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 737
      for (i = 0; i < nvars; ++i)
        x[i] = bestx[i];
    }
- //  EcoSystem->writeOptValuesOMP();
  }
  #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 762
                    }
                  } 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 792
                    }
                  } 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 850
+Line 808
  #ifndef NO_OPENMP
          int id = omp_get_thread_num ();
          trialf = EcoSystems[id]->SimulateAndUpdate(params);
- //      cout << id << "-";
  #else
          trialf = EcoSystem->SimulateAndUpdate(params);
  #endif
- //      cout << "trailf:" << trialf<< endl;
    return -trialf;
  }
-Line 872
+Line 828
                  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 857
                  DoubleVector* bestx = seed.getBestx();
-                 //funcval = fopt;
                  for (i = 0; i < seed.getNvars(); i++)
                          x[i] = (*bestx)[i];
          }
-Line 945
+Line 895
            //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 904
                  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 928
          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 948
           bool quit = false;
            if (siman.getNT() >= siman.getNt())
            {
- //                cout << "nt:" << siman.getNT();
-                   //Check termination criteria
                    int i;
                    DoubleVector fstar = siman.getFstar();
-Line 1067
+Line 998
          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;
  }
-Line 1089
+Line 1014
          //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 1112
+Line 1037
                  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
          }
-Line 1152
+Line 1077
          //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>
          pa(s, x, simanniter);
- //      cout << "Sequential run" << endl;
- //      pa.seq_opt();
  #ifndef NO_OPENMP
+         // OpenMP parallelization
          int numThr = omp_get_max_threads ( );
-         pa.paral_opt_omp(numThr,numThr);
+         pa.paral_opt_omp(funcval,numThr,numThr);
  #else
-         pa.seq_opt();
+         // sequential code
+         pa.seq_opt(funcval);
  #endif
  }

 Legend:



Removed from v.11
 


changed lines


 
Added in v.14
 Legend:



Removed from v.11
 


changed lines


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

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