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

Diff of /trunk/gadget/pso.cc

-revision 27, Mon May  8 11:58:23 2017 UTC
+revision 28, Wed Jul  5 12:08:06 2017 UTC
 Line 1
  #include "gadget.h"    //All the required standard header files are in here
  #include "optinfo.h"
  #include "mathfunc.h"
-Line 8
+Line 7
  #include "ecosystem.h"
  #include "global.h"
  #include "pso.h"
+ #include <float.h>
  #ifdef _OPENMP
  #include "omp.h"
 Line 18
  extern Ecosystem** EcoSystems;
  #endif
  //==============================================================
  //// calulate swarm size based on dimensionality
  int OptInfoPso::pso_calc_swarm_size(int dim) {
-Line 29
+Line 28
  void OptInfoPso::OptimiseLikelihood() {
    handle.logMessage(LOGINFO, "\nStarting PSO optimisation algorithm\n");
-   double oldf, newf, steplength, tmp;
+     double bestf, tmp;
    int    i, offset,rchange,rcheck,rnumber;
    int nvars = EcoSystem->numOptVariables();
    DoubleVector x(nvars);
-   DoubleVector trialx(nvars);
+ //    DoubleVector trialx(nvars);
    DoubleVector bestx(nvars);
-   struct Best {double bestf; int index;};
-   struct Best best;
-   double bestf;
    DoubleVector lowerb(nvars);
    DoubleVector upperb(nvars);
    DoubleVector init(nvars);
-Line 50
+Line 46
  /**
   * Internal variables for PSO
   */
+ //    printf("nvars:%d\n",nvars);
+     if (size == 0) size = pso_calc_swarm_size(nvars);
    DoubleMatrix pos(size,nvars,0.0); // position matrix
    DoubleMatrix vel(size,nvars,0.0); // velocity matrix
    DoubleMatrix pos_b(size,nvars,0.0); // best position matrix
-Line 66
+Line 64
    Calc_inertia_fun calc_inertia_fun; // inertia weight update function
    int psosteps=0;
-   if (this->seed!=0)
+     if (seed == 0) seed = 1234;
-      srand(this->seed);
    handle.logMessage(LOGINFO,"Starting PSO with particles ",size,"\n");
-   switch (nhood_strategy)
-         {
+ // SELECT APPROPRIATE NHOOD UPDATE FUNCTION
+     switch (nhood_strategy) {
          case PSO_NHOOD_GLOBAL :
              // comm matrix not used
              handle.logMessage(LOGINFO,"PSO will use global communication\n");
-Line 89
+Line 87
              break;
          }
    // SELECT APPROPRIATE INERTIA WEIGHT UPDATE FUNCTION
-   switch (w_strategy)
+     switch (w_strategy) {
-   {
          case PSO_W_CONST :
               handle.logMessage(LOGINFO,"PSO will use constant inertia\n");
               calc_inertia_fun = &OptInfoPso::calc_inertia_const;
-Line 109
+Line 106
    handle.logMessage(LOGINFO,"PSO c2",c2,"\n");
    handle.logMessage(LOGINFO,"PSO goal",goal,"\n");
+ //    if (scale)
    EcoSystem->scaleVariables();
- #ifdef _OPENMP
+ //#ifdef _OPENMP
-   int numThr = omp_get_max_threads ( );
+ //  int numThr = omp_get_max_threads ( );
- #pragma omp parallel for
+ //#pragma omp parallel for
-   for (i=0;i<numThr;i++)
+ //  for (i=0;i<numThr;i++)
-       EcoSystems[i]->scaleVariables();
+ //      EcoSystems[i]->scaleVariables();
- #endif
+ //#endif
    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;
+     }
+ //se escala?
+ //    if (scale)
    for (i = 0; i < nvars; i++) {
      // Scaling the bounds, because the parameters are scaled
      lowerb[i] = lowerb[i] / init[i];
-Line 130
+Line 135
        lowerb[i] = upperb[i];
        upperb[i] = tmp;
      }
-     bestx[i] = x[i];
-     trialx[i] = x[i];
-     param[i] = i;
    }
-   bestf = EcoSystem->SimulateAndUpdate(trialx);
+     bestf = EcoSystem->SimulateAndUpdate(x);
    if (bestf != bestf) { //check for NaN
      handle.logMessage(LOGINFO, "Error starting PSO optimisation with f(x) = infinity");
      converge = -1;
-Line 144
+Line 146
    }
    offset = EcoSystem->getFuncEval();  //number of function evaluations done before loop
+ //  newf = bestf;
+ //  oldf = bestf;
    iters = 0;
-Line 152
+Line 156
    //
    i=0;
    handle.logMessage(LOGINFO,"Initialising PSO ",size,"particles\n");
+     // SWARM INITIALIZATION
    //
    // Initialize the fist particle with the best known position
    //
    for (d=0; d<nvars; d++) {
              // generate two numbers within the specified range and around trialx
-             a = trialx[d] ;
+         a = x[d];
              // initialize position
              pos[i][d] = a;
              pos_b[i][d] = a;
              pos_nb[i][d] = a;
-             a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
+ //        a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand() * 1.0 / RAND_MAX);
-             b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
+         b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
              vel[i][d] = (a-b) / 2.;
    }
    fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
    fit_b[i] = fit[i]; // this is also the personal best
    if (fit[i] < bestf) {
              bestf = fit[i];
-             best.bestf=bestf;
              bestx=pos[i];
-             best.index=i;
    }
- //
  // Now other particles with a random position in the solution space
  //
- best.bestf=bestf;
+ //#pragma omp parallel for default(shared) private(i,d,a,b)
- best.index=0;
- handle.logMessage(LOGINFO,"before Best.bestf=",best.bestf);
- #pragma omp declare reduction(bestomp: struct Best  : omp_out = omp_in.bestf < omp_out.bestf ? omp_in : omp_out)
- #pragma omp parallel for default(shared) private(i,d,a,b) reduction(bestomp:best) firstprivate(bestf)
    for (i=1; i<size; i++) {
-         // for each dimension
        ostringstream spos,svel;
        spos << "Initial position of particle " << i<<" [" ;
        svel << "Initial velocity of particle " << i<<" [" ;
+         // for each dimension
          for (d=0; d<nvars; d++) {
              // generate two numbers within the specified range and around trialx
-             a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
+             a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
-             b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
+             b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
              // initialize position
              pos[i][d] = a;
              // best position is trialx
-Line 202
+Line 201
              spos << pos[i][d] << " ";
              svel << vel[i][d] << " ";
          }
- #ifdef _OPENMP
+ //#ifdef _OPENMP
-         fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]);
+ //      fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]);
-         best.bestf = fit[i];
+ //#else
- #else
          fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
- #endif
+ //#endif
          fit_b[i] = fit[i]; // this is also the personal best
          // update gbest??
  //#pragma omp critical (best1)
          if (fit[i] < bestf) {
              // update best fitness
-                 handle.logMessage(LOGINFO,"iInside Best.bestf=",best.bestf);
+             bestf = fit[i];
              // copy particle pos to gbest vector
              //memmove((void *)&bestx, (void *)&pos[i],
          //          sizeof(double) * nvars);
-             best.index=i;
+             bestx = pos[i];
-             bestf=best.bestf;
-             //bestx=pos[i];
          }
          spos<<" ]";
          svel<<" ]";
          handle.logMessage(LOGDEBUG, spos.str().c_str());
          handle.logMessage(LOGDEBUG, svel.str().c_str());
-   }
- //
- // Compare the reduction with the best value
- //
-   if (best.bestf < bestf){
-         bestx=pos[best.index];
-         bestf=best.bestf;
    }
   //
   // RUN ALGORITHM
   //
-Line 260
+Line 246
        handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of PSO steps was reached");
        handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
-       bestf = EcoSystem->SimulateAndUpdate(bestx);
+             score = EcoSystem->SimulateAndUpdate(bestx);
        for (i = 0; i < nvars; i++)
          bestx[i] = bestx[i] * init[i];
-       handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", bestf);
+             handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score);
-       EcoSystem->storeVariables(bestf, bestx);
+             EcoSystem->storeVariables(score, bestx);
        return;
      }
+         // update inertia weight
-     if (w_strategy)
+         if (w_strategy) w = (this->*calc_inertia_fun)(psosteps);
-       w = (this->*calc_inertia_fun)(psosteps);
          // check optimization goal
      if (bestf <= goal) {
        handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n");
        handle.logMessage(LOGINFO, "Goal achieved!!! @ step ", psosteps);
-       bestf = EcoSystem->SimulateAndUpdate(bestx);
+             score = EcoSystem->SimulateAndUpdate(bestx);
-       handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", bestf);
+             handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score);
        break;
        }
-Line 286
+Line 271
      improved = 0;
          // update all particles
-     best.bestf=bestf;
+ //#pragma omp parallel for private(d,rho1,rho2) default(shared)
- #pragma omp parallel for default(shared) private(i,d,rho1,rho2)  reduction(bestomp:best) firstprivate(bestf)
      for (i=0; i<size; i++) {
-           // for each dimension
+ //      handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n");
-       handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n");
        ostringstream spos,svel;
        spos << "New position of particle " << i<<" [" ;
        svel << "New velocity of particle " << i<<" [" ;
+             // for each dimension
        for (d=0; d<nvars; d++) {
          // calculate stochastic coefficients
-         rho1 = c1 * (rand()*1.0)/RAND_MAX;
+                 rho1 = c1 * (rand_r(&seed) * 1.0) / RAND_MAX;
-         rho2 = c2 * (rand()*1.0)/RAND_MAX;
+                 rho2 = c2 * (rand_r(&seed) * 1.0) / RAND_MAX;
          // update velocity
-         vel[i][d] = w * vel[i][d] +     \
+                 vel[i][d] = w * vel[i][d] + rho1 * (pos_b[i][d] - pos[i][d]) + rho2 * (pos_nb[i][d] - pos[i][d]);
-             rho1 * (pos_b[i][d] - pos[i][d]) +  \
-             rho2 * (pos_nb[i][d] - pos[i][d]);
                  // update position
          pos[i][d] += vel[i][d];
          spos << pos[i][d] <<" ";
-Line 311
+Line 293
              if (pos[i][d] < lowerb[d]) {
                  pos[i][d] = lowerb[d];
                  vel[i][d] = 0;
-             } else if (pos[i][d] > upperb[d]) {
+                     }
+                     else if (pos[i][d] > upperb[d]) {
                  pos[i][d] = upperb[d];
                  vel[i][d] = 0;
              }
-         } else {
+                 }
+                 else {
                      // enforce periodic boundary conditions
              if (pos[i][d] < lowerb[d]) {
+                         pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d], upperb[d] - lowerb[d]);
-                 pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d],
-                                           upperb[d] - lowerb[d]);
                  vel[i][d] = 0;
+                     }
-             } else if (pos[i][d] > upperb[d]) {
+                     else if (pos[i][d] > upperb[d]) {
+                         pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d], upperb[d] - lowerb[d]);
-                 pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d],
-                                           upperb[d] - lowerb[d]);
                  vel[i][d] = 0;
              }
            }
+             } //end for nvars
-          } //end for d
         spos<< "]" << endl;
         svel<< "]" << endl;
         handle.logMessage(LOGDEBUG,spos.str().c_str());
         handle.logMessage(LOGDEBUG,svel.str().c_str());
              // update particle fitness
-        fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]);
+ //       fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]);
+             fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
              // update personal best position?
-        best.bestf = fit[i];
         if (fit[i] < fit_b[i]) {
             fit_b[i] = fit[i];
-                 // copy contents of pos[i] to pos_b[i]
-                 //memmove((void *)&pos_b[i], (void *)&pos[i],
-         //              sizeof(double) * nvars);
             pos_b[i]=pos[i];
         }
              // update gbest??
-Line 358
+Line 334
          if (fit[i] < bestf) {
             improved = 1;
                  // update best fitness
-            best.index=i;
             bestf=fit[i];
-                 // copy particle pos to gbest vector
-                 //memmove((void *)&bestx, (void *)&pos[i],
-         //              sizeof(double) * nvars);
-         //   for (d = 0; d < nvars; d++)
-         //     bestx[d] = pos[i][d] * init[d];
-         //   iters = EcoSystem->getFuncEval() - offset;
-         //   handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
-         //   handle.logMessage(LOGINFO, "The likelihood score is", best.bestf, "at the point");
-         //   EcoSystem->storeVariables(best.bestf, bestx);
-         //   EcoSystem->writeBestValues();
-         //   bestx=pos[i];
-         }
-       } // end i
-       if (best.bestf < bestf){
-            bestf=best.bestf;
-            bestx=pos[best.index];
             for (d = 0; d < nvars; d++)
-              bestx[d] = bestx[d] * init[d];
+                     bestx[d] = pos[i][d] * init[d];
-            iters = EcoSystem->getFuncEval() - offset; //To change. It is incorrect when is parallel.
+                 iters = EcoSystem->getFuncEval() - offset;
             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();
-            bestx=pos[best.index];
+                 bestx = pos[i];
          }
+         } // end i
        psosteps++;
-Line 398
+Line 358
      EcoSystem->storeVariables(bestf, bestx);
      EcoSystem->writeBestValues();
  }
+ #ifdef _OPENMP
  void OptInfoPso::OptimiseLikelihoodREP() {
        OptimiseLikelihood();
  }
  void OptInfoPso::OptimiseLikelihoodOMP() {
-       OptimiseLikelihood();
+     handle.logMessage(LOGINFO, "\nStarting PSO optimisation algorithm\n");
+     double tmp;
+     int i, offset, rchange, rcheck, rnumber;
+     int nvars = EcoSystem->numOptVariables();
+     DoubleVector x(nvars);
+ //    DoubleVector trialx(nvars);
+     DoubleVector bestx(nvars);
+     DoubleVector scalex(nvars);
+     DoubleVector lowerb(nvars);
+     DoubleVector upperb(nvars);
+     DoubleVector init(nvars);
+ //DoubleVector initialstep(nvars, rho);
+     IntVector param(nvars, 0);
+     IntVector lbound(nvars, 0);
+     IntVector rbounds(nvars, 0);
+     IntVector trapped(nvars, 0);
+     struct Best { double bestf; int index;};
+     struct Best best;
+     /**
+      * Internal variables for PSO
+      */
+ //    printf("nvars:%d\n",nvars);
+     if (size == 0) size = pso_calc_swarm_size(nvars);
+     DoubleMatrix pos(size, nvars, 0.0); // position matrix
+     DoubleMatrix vel(size, nvars, 0.0); // velocity matrix
+     DoubleMatrix pos_b(size, nvars, 0.0); // best position matrix
+     DoubleVector fit(size); // particle fitness vector
+     DoubleVector fit_b(size); // best fitness vector
+     DoubleMatrix pos_nb(size, nvars, 0.0); // what is the best informed position for each particle
+     IntMatrix comm(size, size, 0); // communications:who informs who rows : those who inform cols : those who are informed
+     int improved; // whether solution->error was improved during the last iteration
+     int d, step;
+     double a, b;
+     double rho1, rho2; //random numbers (coefficients)
+     double w; //current omega
+     Inform_fun inform_fun;
+     Calc_inertia_fun calc_inertia_fun; // inertia weight update function
+     int psosteps = 0; //!= iters?
+ //    if (seed == 0) seed = 1234;
+     handle.logMessage(LOGINFO, "Starting PSO with particles ", size, "\n");
+ // SELECT APPROPRIATE NHOOD UPDATE FUNCTION
+     switch (nhood_strategy) {
+         case PSO_NHOOD_GLOBAL:
+             // comm matrix not used
+             handle.logMessage(LOGINFO, "PSO will use global communication\n");
+             inform_fun = &OptInfoPso::inform_global;
+             break;
+         case PSO_NHOOD_RING:
+             handle.logMessage(LOGINFO, "PSO will use ring communication\n");
+             init_comm_ring(comm);
+             inform_fun = &OptInfoPso::inform_ring;
+             break;
+         case PSO_NHOOD_RANDOM:
+             handle.logMessage(LOGINFO, "PSO will use random communication\n");
+             init_comm_random(comm);
+             inform_fun = &OptInfoPso::inform_random;
+             break;
+     }
+ // SELECT APPROPRIATE INERTIA WEIGHT UPDATE FUNCTION
+     switch (w_strategy) {
+         case PSO_W_CONST:
+             handle.logMessage(LOGINFO, "PSO will use constant inertia\n");
+             calc_inertia_fun = &OptInfoPso::calc_inertia_const;
+             break;
+         case PSO_W_LIN_DEC:
+             handle.logMessage(LOGINFO, "PSO will use linear decreasing inertia\n");
+             calc_inertia_fun = &OptInfoPso::calc_inertia_lin_dec;
+             break;
+     }
+ //
+ //initialize omega using standard value
+ //
+     w = PSO_INERTIA;
+     handle.logMessage(LOGINFO, "PSO initial inertia ", w, "\n");
+     handle.logMessage(LOGINFO, "PSO c1", c1, "\n");
+     handle.logMessage(LOGINFO, "PSO c2", c2, "\n");
+     handle.logMessage(LOGINFO, "PSO goal", goal, "\n");
+     EcoSystem->scaleVariables();
+     int numThr = omp_get_max_threads ( );
+ //    if (scale) {
+ #pragma omp parallel for
+         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) {
+         bestx[i] = x[i];
+         param[i] = i;
+       }
+ //se escala?
+ //    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;
+             }
+         }
+     best.bestf = EcoSystem->SimulateAndUpdate(x);
+     best.index=0;
+     if (best.bestf != best.bestf) { //check for NaN
+         handle.logMessage(LOGINFO, "Error starting PSO optimisation with f(x) = infinity");
+         converge = -1;
+         iters = 1;
+         return;
+     }
+     offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop
+ //  newf = bestf;
+ //  oldf = bestf;
+     iters = 0;
+ //
+ // Initialize the particles
+ //
+     i = 0;
+     handle.logMessage(LOGINFO, "Initialising PSO ", size, "particles\n");
+     // SWARM INITIALIZATION
+ //
+ // Initialize the fist particle with the best known position
+ //
+ //    best.bestf=bestf;
+     for (d = 0; d < nvars; d++) {
+         // generate two numbers within the specified range and around trialx
+         a = x[d];
+         // initialize position
+         pos[i][d] = a;
+         pos_b[i][d] = a;
+         pos_nb[i][d] = a;
+ //        a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand() * 1.0 / RAND_MAX);
+         b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
+         vel[i][d] = (a - b) / 2.;
+     }
+     fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
+     fit_b[i] = fit[i]; // this is also the personal best
+     if (fit[i] < best.bestf) {
+         best.bestf = fit[i];
+ //        best.bestf=bestf;
+         bestx = pos[i];
+         best.index=i;
    }
+ // Now other particles with a random position in the solution space
+ //
+ //    handle.logMessage(LOGINFO, "ini 1 best.bestf=",best.bestf, "\n");
+ #pragma omp declare reduction(bestomp: Best : omp_out = omp_in.bestf < omp_out.bestf ? omp_in : omp_out) initializer(omp_priv=omp_orig)
+ #pragma omp parallel for default(shared) private(i,d,a,b) reduction(bestomp:best)  //firstprivate(bestf)
+     for (i = 1; i < size; i++) {
+         ostringstream spos, svel;
+         spos << "Initial position of particle " << i << " [";
+         svel << "Initial velocity of particle " << i << " [";
+         // for each dimension
+         for (d = 0; d < nvars; d++) {
+             // generate two numbers within the specified range and around trialx
+             a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
+             b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
+             // initialize position
+             pos[i][d] = a;
+             // best position is trialx
+             pos_b[i][d] = a;
+             pos_nb[i][d] = a;
+             // initialize velocity
+             vel[i][d] = (a - b) / 2.;
+             spos << pos[i][d] << " ";
+             svel << vel[i][d] << " ";
+         }
+ //#ifdef _OPENMP
+       fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]);
+ //      best.bestf = fit[i];
+ //#else
+ //        fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
+ //#endif
+         fit_b[i] = fit[i]; // this is also the personal best
+         // update gbest??
+ //#pragma omp critical (best1)
+         if (fit[i] < best.bestf) {
+             // update best fitness
+ //            bestf = fit[i];
+ // copy particle pos to gbest vector
+ //memmove((void *)&bestx, (void *)&pos[i],
+ //          sizeof(double) * nvars);
+ //            bestx = pos[i];
+             best.index=i;
+ //            bestf=best.bestf;
+             best.bestf = fit[i];
+         }
+         spos << " ]";
+         svel << " ]";
+         handle.logMessage(LOGDEBUG, spos.str().c_str());
+         handle.logMessage(LOGDEBUG, svel.str().c_str());
+     }
+ //
+ // Compare the reduction with the best value
+ //
+ //      if (best.bestf < bestf){
+ //            bestx=pos[best.index];
+ //            bestf=best.bestf;
+ //      }
+ //
+ // RUN ALGORITHM
+ //
+ //    handle.logMessage(LOGINFO, "ini 2 best.bestf=",best.bestf, "\n");
+     psosteps = 0;
+     while (1) {
+         handle.logMessage(LOGINFO, "PSO optimisation after", psosteps * size, "\n");
+         if (isZero(best.bestf)) {
+             iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset;
+             handle.logMessage(LOGINFO, "Error in PSO optimisation after", iters, "function evaluations, f(x) = 0");
+             converge = -1;
+             return;
+         }
+         /* if too many function evaluations occur, terminate the algorithm */
+         iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset;
+         if (psosteps > psoiter) {
+             handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n");
+             handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
+             handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of PSO steps was reached");
+             handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
+             score = EcoSystem->SimulateAndUpdate(bestx);
+             for (i = 0; i < nvars; i++)
+                 bestx[i] = bestx[i] * init[i];
+             handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score);
+             EcoSystem->storeVariables(score, bestx);
+             return;
+         }
+         // update inertia weight
+         if (w_strategy) w = (this->*calc_inertia_fun)(psosteps);
+         // check optimization goal
+         if (best.bestf <= goal) {
+             handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n");
+             handle.logMessage(LOGINFO, "Goal achieved!!! @ step ", psosteps);
+             score = EcoSystem->SimulateAndUpdate(bestx);
+             handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score);
+             break;
+         }
+         // update pos_nb matrix (find best of neighborhood for all particles)
+         (this->*inform_fun)(comm, pos_nb, pos_b, fit_b, bestx, improved);
+         // the value of improved was just used; reset it
+         improved = 0;
+         // update all particles
+ //        best.bestf=bestf;
+ //#pragma omp declare reduction(bestomp2: Best : omp_out = omp_in.bestf < omp_out.bestf ? omp_in : omp_out) initializer(omp_priv=omp_orig)
+ #pragma omp parallel for default(shared) private(i,d,rho1,rho2)  reduction(bestomp:best) //firstprivate(best)
+         for (i = 0; i < size; i++) {
+ //      handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n");
+             ostringstream spos, svel;
+             spos << "New position of particle " << i << " [";
+             svel << "New velocity of particle " << i << " [";
+             // for each dimension
+             for (d = 0; d < nvars; d++) {
+                 // calculate stochastic coefficients
+                 rho1 = c1 * (rand_r(&seed) * 1.0) / RAND_MAX;
+                 rho2 = c2 * (rand_r(&seed) * 1.0) / RAND_MAX;
+                 // update velocity
+                 vel[i][d] = w * vel[i][d] + rho1 * (pos_b[i][d] - pos[i][d]) + rho2 * (pos_nb[i][d] - pos[i][d]);
+                 // update position
+                 pos[i][d] += vel[i][d];
+                 spos << pos[i][d] << " ";
+                 svel << vel[i][d] << " ";
+                 // clamp position within bounds?
+                 if (clamp_pos) {
+                     if (pos[i][d] < lowerb[d]) {
+                         pos[i][d] = lowerb[d];
+                         vel[i][d] = 0;
+                     }
+                     else if (pos[i][d] > upperb[d]) {
+                         pos[i][d] = upperb[d];
+                         vel[i][d] = 0;
+                     }
+                 }
+                 else {
+                     // enforce periodic boundary conditions
+                     if (pos[i][d] < lowerb[d]) {
+                         pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d], upperb[d] - lowerb[d]);
+                         vel[i][d] = 0;
+                     }
+                     else if (pos[i][d] > upperb[d]) {
+                         pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d], upperb[d] - lowerb[d]);
+                         vel[i][d] = 0;
+                     }
+                 }
+             } //end for nvars
+             spos << "]" << endl;
+             svel << "]" << endl;
+             handle.logMessage(LOGDEBUG, spos.str().c_str());
+             handle.logMessage(LOGDEBUG, svel.str().c_str());
+             // update particle fitness
+             fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]);
+ //            fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
+             // update personal best position?
+ //            best.bestf = fit[i];
+             if (fit[i] < fit_b[i]) {
+                 fit_b[i] = fit[i];
+                 pos_b[i] = pos[i];
+             }
+             // update gbest??
+             handle.logMessage(LOGDEBUG, "Particle ", i);
+             handle.logMessage(LOGDEBUG, "The likelihood score is", best.bestf, "at the point");
+             handle.logMessage(LOGDEBUG, "Its likelihood score is", fit[i], "at the point");
+ //#pragma omp critical (best2)
+             if (fit[i] < best.bestf) {
+                 improved = 1;
+                 // update best fitness
+                 best.index=i;
+                 best.bestf = fit[i];
+                 for (d = 0; d < nvars; d++)
+                     bestx[d] = pos[i][d] * init[d];
+                 iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset;
+                 handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
+                 handle.logMessage(LOGINFO, "The likelihood score is", best.bestf, "at the point");
+                 EcoSystem->storeVariables(best.bestf, bestx);
+                 EcoSystem->writeBestValues();
+                 bestx = pos[i];
+             }
+         } // end i
+ //        if (best.bestf < bestf){
+ //           bestf=best.bestf;
+ //           bestx=pos[best.index];
+ //           for (d = 0; d < nvars; d++)
+ //             bestx[d] = bestx[d] * init[d];
+ //           iters = EcoSystem->getFuncEval() - offset; //To change. It is incorrect when is parallel.
+ //           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();
+ //           bestx=pos[best.index];
+ //        }
+ //        handle.logMessage(LOGINFO, "best.bestf=",best.bestf, "\n");
+         psosteps += numThr;
+     } // while (newf < bestf)
+     iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset;
+     handle.logMessage(LOGINFO, "Existing PSO after ", iters, "function evaluations ...");
+     for (d = 0; d < nvars; d++)
+         bestx[d] = bestx[d] * init[d];
+     EcoSystem->storeVariables(best.bestf, bestx);
+     EcoSystem->writeBestValues();
+ }
+ #endif
  /**
   * \brief calculate constant  inertia weight equal to w_max
   */
-Line 423
+Line 757
      int dec_stage = 3 * psoiter / 4;
      if (step <= dec_stage)
-         return w_min + (w_max - w_min) *        \
+         return w_min + (w_max - w_min) * (dec_stage - step) / dec_stage;
-             (dec_stage - step) / dec_stage;
      else
          return w_min;
  }
  /**
   * \brief NEIGHBORHOOD (COMM) MATRIX STRATEGIES. global neighborhood
   */
- void OptInfoPso::inform_global(IntMatrix& comm, DoubleMatrix& pos_nb,
+ void OptInfoPso::inform_global(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) {
-                    DoubleMatrix& pos_b, DoubleVector& fit_b,
-                    DoubleVector& gbest,int improved)
- {
      int i;
      // all particles have the same attractor (gbest)
-Line 448
+Line 776
          pos_nb[i]=gbest;
  }
  /**
   * \brief  general inform function :: according to the connectivity matrix COMM, it copies the best position (from pos_b) of the informers of each particle to the pos_nb matrix
   */
- void OptInfoPso::inform(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, int improved)
+ void OptInfoPso::inform(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, int improved) {
- {
      int i, j;
      int b_n; // best neighbor in terms of fitness
-Line 491
+Line 817
              comm[i][1] = 1;
              // look left
              comm[i][size-1] = 1;
-         } else if (i==size-1) {
+         }
+         else if (i == size - 1) {
              // look right
              comm[i][0] = 1;
              // look left
              comm[i][i-1] = 1;
-         } else {
+         }
+         else {
              // look right
              comm[i][i+1] = 1;
              // look left
-Line 508
+Line 836
      handle.logMessage(LOGDEBUG, "\nPSO ring communication\n");
      for (i=0;i<size;i++){
         string message="[";
-        for (int j=0;j<size;j++)
+         for (int j = 0; j < size; j++) {
-        {
            std::ostringstream ss;
            ss<< comm[i][j];
            message=message+" "+ss.str();
-Line 520
+Line 847
  }
+ void OptInfoPso::inform_ring(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) {
- void OptInfoPso::inform_ring(IntMatrix& comm, DoubleMatrix& pos_nb,
-                  DoubleMatrix& pos_b, DoubleVector& fit_b,
-                  DoubleVector& gbest, int improved)
- {
      // update pos_nb matrix
      inform(comm, pos_nb, pos_b, fit_b, improved);
-Line 550
+Line 871
          for (k=0; k<nhood_size; k++) {
              // generate a random index
              //j = gsl_rng_uniform_int(settings->rng, settings->size);
-             j=rand()*size/RAND_MAX;
+             j = rand_r(&this->seed) * size / RAND_MAX;
              // particle i informs particle j
              comm[i] [j] = 1;
          }
      }
  }
+ void OptInfoPso::inform_random(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) {
- void OptInfoPso::inform_random(IntMatrix& comm, DoubleMatrix& pos_nb,
-                    DoubleMatrix& pos_b, DoubleVector& fit_b,
-                    DoubleVector& gbest, int improved)
- {
      // regenerate connectivity??
-     if (!improved)
+     if (!improved) init_comm_random(comm);
-         init_comm_random(comm);
      inform(comm, pos_nb, pos_b, fit_b, improved);
  }

 Legend:



Removed from v.27
 


changed lines


 
Added in v.28
 Legend:



Removed from v.27
 


changed lines


 
Added in v.28
-Removed from v.27
+Added in v.28

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