Log In | Get Help   
Home My Page Projects Code Snippets Project Openings Mareframe
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files
[mareframe] Diff of /trunk/gadget/hooke.cc
[mareframe] / trunk / gadget / hooke.cc Repository:
ViewVC logotype

Diff of /trunk/gadget/hooke.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1, Mon Feb 10 17:09:07 2014 UTC revision 14, Thu Jul 30 12:36:46 2015 UTC
# Line 1  Line 1 
1  /* Nonlinear Optimization using the algorithm of Hooke and Jeeves  */  /* Nonlinear Optimization using the algorithm of Hooke and Jeeves  */
2  /*  12 February 1994    author: Mark G. Johnson                    */  /*  12 February 1994    author: Mark G. Johnson                    */
3    //
4  /* Find a point X where the nonlinear function f(X) has a local    */  /* Find a point X where the nonlinear function f(X) has a local    */
5  /* minimum.  X is an n-vector and f(X) is a scalar.  In mathe-     */  /* minimum.  X is an n-vector and f(X) is a scalar.  In mathe-     */
6  /* matical notation  f: R^n -> R^1.  The objective function f()    */  /* matical notation  f: R^n -> R^1.  The objective function f()    */
7  /* is not required to be continuous.  Nor does f() need to be      */  /* is not required to be continuous.  Nor does f() need to be      */
8  /* differentiable.  The program does not use or require            */  /* differentiable.  The program does not use or require            */
9  /* derivatives of f().                                             */  /* derivatives of f().                                             */
10    //
11  /* The software user supplies three things: a subroutine that      */  /* The software user supplies three things: a subroutine that      */
12  /* computes f(X), an initial "starting guess" of the minimum point */  /* computes f(X), an initial "starting guess" of the minimum point */
13  /* X, and values for the algorithm convergence parameters.  Then   */  /* X, and values for the algorithm convergence parameters.  Then   */
14  /* the program searches for a local minimum, beginning from the    */  /* the program searches for a local minimum, beginning from the    */
15  /* starting guess, using the Direct Search algorithm of Hooke and  */  /* starting guess, using the Direct Search algorithm of Hooke and  */
16  /* Jeeves.                                                         */  /* Jeeves.                                                         */
17    //
18  /* This C program is adapted from the Algol pseudocode found in    */  /* This C program is adapted from the Algol pseudocode found in    */
19  /* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun-  */  /* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun-  */
20  /* ications of the ACM, Vol 6. p.313 (June 1963).  It includes the */  /* ications of the ACM, Vol 6. p.313 (June 1963).  It includes the */
# Line 24  Line 24 
24  /* highly as the one by A. Kaupe, is:  R. Hooke and T. A. Jeeves,  */  /* highly as the one by A. Kaupe, is:  R. Hooke and T. A. Jeeves,  */
25  /* "Direct Search Solution of Numerical and Statistical Problems", */  /* "Direct Search Solution of Numerical and Statistical Problems", */
26  /* Journal of the ACM, Vol. 8, April 1961, pp. 212-229.            */  /* Journal of the ACM, Vol. 8, April 1961, pp. 212-229.            */
27    //
28  /* Calling sequence:                                               */  /* Calling sequence:                                               */
29  /*  int hooke(nvars, startpt, endpt, rho, epsilon, itermax)        */  /*  int hooke(nvars, startpt, endpt, rho, epsilon, itermax)        */
30  /*                                                                 */  /*                                                                 */
# Line 59  Line 59 
59  /*     itermax     {an integer}  A second, rarely used, halting    */  /*     itermax     {an integer}  A second, rarely used, halting    */
60  /*         criterion.  If the algorithm uses >= itermax            */  /*         criterion.  If the algorithm uses >= itermax            */
61  /*         iterations, halt.                                       */  /*         iterations, halt.                                       */
62    //
63  /* The user-supplied objective function f(x,n) should return a C   */  /* The user-supplied objective function f(x,n) should return a C   */
64  /* "double".  Its  arguments are  x -- an array of doubles, and    */  /* "double".  Its  arguments are  x -- an array of doubles, and    */
65  /* n -- an integer.  x is the point at which f(x) should be        */  /* n -- an integer.  x is the point at which f(x) should be        */
66  /* evaluated, and n is the number of coordinates of x.  That is,   */  /* evaluated, and n is the number of coordinates of x.  That is,   */
67  /* n is the number of coefficients being fitted.                   */  /* n is the number of coefficients being fitted.                   */
68    //
69  /*             rho, the algorithm convergence control              */  /*             rho, the algorithm convergence control              */
70    //
71  /*  The algorithm works by taking "steps" from one estimate of     */  /*  The algorithm works by taking "steps" from one estimate of     */
72  /*    a minimum, to another (hopefully better) estimate.  Taking   */  /*    a minimum, to another (hopefully better) estimate.  Taking   */
73  /*    big steps gets to the minimum more quickly, at the risk of   */  /*    big steps gets to the minimum more quickly, at the risk of   */
# Line 95  Line 95 
95  /*    again with a larger value of rho such as 0.85, using the     */  /*    again with a larger value of rho such as 0.85, using the     */
96  /*    result of the first minimization as the starting guess to    */  /*    result of the first minimization as the starting guess to    */
97  /*    begin the second minimization.                               */  /*    begin the second minimization.                               */
98    //
99  /*                        Normal use:                              */  /*                        Normal use:                              */
100  /*    (1) Code your function f() in the C language                 */  /*    (1) Code your function f() in the C language                 */
101  /*    (2) Install your starting guess {or read it in}              */  /*    (2) Install your starting guess {or read it in}              */
102  /*    (3) Run the program                                          */  /*    (3) Run the program                                          */
103  /*    (4) {for the skeptical}: Use the computed minimum            */  /*    (4) {for the skeptical}: Use the computed minimum            */
104  /*        as the starting point for another run                    */  /*        as the starting point for another run                    */
105    //
106  /* Data Fitting:                                                   */  /* Data Fitting:                                                   */
107  /*  Code your function f() to be the sum of the squares of the     */  /*  Code your function f() to be the sum of the squares of the     */
108  /*  errors (differences) between the computed values and the       */  /*  errors (differences) between the computed values and the       */
# Line 113  Line 113 
113  /*  f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i]))     */  /*  f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i]))     */
114  /*                + (C*tan(t[i]))))^2                              */  /*                + (C*tan(t[i]))))^2                              */
115  /*  where x[] is a 3-vector consisting of {A, B, C}.               */  /*  where x[] is a 3-vector consisting of {A, B, C}.               */
116    //
117  /*  The author of this software is M.G. Johnson.                   */  /*  The author of this software is M.G. Johnson.                   */
118  /*  Permission to use, copy, modify, and distribute this software  */  /*  Permission to use, copy, modify, and distribute this software  */
119  /*  for any purpose without fee is hereby granted, provided that   */  /*  for any purpose without fee is hereby granted, provided that   */
# Line 125  Line 125 
125  /*  AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY     */  /*  AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY     */
126  /*  KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS    */  /*  KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS    */
127  /*  FITNESS FOR ANY PARTICULAR PURPOSE.                            */  /*  FITNESS FOR ANY PARTICULAR PURPOSE.                            */
128    //
129  /* JMB this has been modified to work with the gadget object structure   */  /* JMB this has been modified to work with the gadget object structure   */
130  /* This means that the function has been replaced by a call to ecosystem */  /* This means that the function has been replaced by a call to ecosystem */
131  /* object, and we can use the vector objects that have been defined      */  /* object, and we can use the vector objects that have been defined      */
# Line 139  Line 139 
139  #include "ecosystem.h"  #include "ecosystem.h"
140  #include "global.h"  #include "global.h"
141    
142  extern Ecosystem* EcoSystem;  #ifndef NO_OPENMP
143    #include "omp.h"
144    #endif
145    
146    extern Ecosystem* EcoSystem;
147    #ifndef NO_OPENMP
148    extern Ecosystem** EcoSystems;
149    #endif
150    
151  /* given a point, look for a better one nearby, one coord at a time */  /* given a point, look for a better one nearby, one coord at a time */
152  double OptInfoHooke::bestNearby(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {  double OptInfoHooke::bestNearby(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
# Line 171  Line 177 
177    return minf;    return minf;
178  }  }
179    
180    /* given a point, look for a better one nearby, one coord at a time */
181    #ifndef NO_OPENMP
182    /*
183     * function bestBeraby parallelized with OpenMP
184     * · 2 threads per coord to parallelize the calculation of +delta/-delta
185     * · parallelize the calculation of the best nearby of the coord
186     */
187    double OptInfoHooke::bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
188      double minf;//, ftmp;
189      int i, j, k;
190      DoubleVector z(point);
191    
192      struct Storage {
193              DoubleVector z;
194              DoubleVector delta;
195              double ftmp;
196              int iters;
197      };
198    
199      minf = prevbest;
200      i = 0;
201    
202      int paral_tokens, numThr, nvars = point.Size();
203      numThr = omp_get_max_threads ( );
204    
205      Storage* storage = new Storage[numThr];
206      if ((numThr % 2) == 0)
207              paral_tokens = numThr / 2;
208      else {
209              return -1;
210      }
211    
212      while ( i < nvars) {
213              if ((i + paral_tokens -1) >= nvars)
214                      paral_tokens = nvars - i;
215              omp_set_dynamic(0);
216              omp_set_nested(1); //permit the nested parallelization
217    #pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2)
218              for (j = 0; j < paral_tokens; ++j) {
219                      storage[j].z = z;
220                      storage[j].delta = delta;
221                      DoubleVector v1(z);
222                      DoubleVector v2(z);
223                      k = param[i+j];
224                      v1[k] +=  delta[k];
225                      v2[k] -=  delta[k];
226    
227    #pragma omp parallel sections num_threads(2) //parrallelize the +/- delta simulation for each parameter
228                      {
229            #pragma omp section
230                              {
231                      storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
232                              }
233            #pragma omp section
234                              {
235                      storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
236                              }
237                      }
238    
239                      if (storage[j].ftmp < minf) {
240                              storage[j].iters = 1;
241                              storage[j].z[k] = v1[k];
242                      } else {
243                              storage[j].iters = 2;
244                              storage[j].delta[k] = 0.0 - delta[k];
245                              if (storage[j+paral_tokens].ftmp < minf) {
246                                      storage[j].ftmp = storage[j+paral_tokens].ftmp;
247                                      storage[j].z[k] = v2[k];
248                              }
249                      }
250              }
251    
252              for (j = 0; j < paral_tokens; ++j) {
253                      i++;
254                      iters += storage[j].iters;
255                      if (storage[j].ftmp < minf) {
256                              minf = storage[j].ftmp;
257                              z = storage[j].z;
258                              delta = storage[j].delta;
259                              break;
260                      }
261              }
262      }
263    
264      for (i = 0; i < nvars; ++i)
265        point[i] = z[i];
266      return minf;
267    }
268    #endif
269    
270  void OptInfoHooke::OptimiseLikelihood() {  void OptInfoHooke::OptimiseLikelihood() {
271    
272    double oldf, newf, bestf, steplength, tmp;    double oldf, newf, bestf, steplength, tmp;
# Line 193  Line 289 
289    IntVector trapped(nvars, 0);    IntVector trapped(nvars, 0);
290    
291    EcoSystem->scaleVariables();    EcoSystem->scaleVariables();
292    #ifndef NO_OPENMP
293      int numThr = omp_get_max_threads ( );
294      for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
295              EcoSystems[i]->scaleVariables();
296    #endif
297    EcoSystem->getOptScaledValues(x);    EcoSystem->getOptScaledValues(x);
298    EcoSystem->getOptLowerBounds(lowerb);    EcoSystem->getOptLowerBounds(lowerb);
299    EcoSystem->getOptUpperBounds(upperb);    EcoSystem->getOptUpperBounds(upperb);
# Line 229  Line 330 
330    if (isZero(steplength))    if (isZero(steplength))
331      steplength = rho;      steplength = rho;
332    
333      iters = 0;
334    
335    while (1) {    while (1) {
336      if (isZero(bestf)) {      if (isZero(bestf)) {
337    #ifdef NO_OPENMP
338        iters = EcoSystem->getFuncEval() - offset;        iters = EcoSystem->getFuncEval() - offset;
339    #endif
340        handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");        handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
341        converge = -1;        converge = -1;
342        return;        return;
# Line 254  Line 359 
359      /* find best new point, one coord at a time */      /* find best new point, one coord at a time */
360      for (i = 0; i < nvars; i++)      for (i = 0; i < nvars; i++)
361        trialx[i] = x[i];        trialx[i] = x[i];
362    #ifndef NO_OPENMP
363        newf = this->bestNearbyOMP(delta, trialx, bestf, param);
364        if (newf == -1) {
365            handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
366            handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
367            return;
368        }
369    #else
370        newf = this->bestNearby(delta, trialx, bestf, param);
371    #endif
372        /* if too many function evaluations occur, terminate the algorithm */
373    
374    #ifdef NO_OPENMP
375        iters = EcoSystem->getFuncEval() - offset;
376    #endif
377        if (iters > hookeiter) {
378          handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
379          handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
380          handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
381          handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations");
382          handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
383    
384          score = EcoSystem->SimulateAndUpdate(trialx);
385          handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
386          for (i = 0; i < nvars; i++)
387            bestx[i] = trialx[i] * init[i];
388          EcoSystem->storeVariables(score, bestx);
389          return;
390        }
391    
392        /* if we made some improvements, pursue that direction */
393        while (newf < bestf) {
394          for (i = 0; i < nvars; i++) {
395            /* if it has been trapped but f has now gotten better (bndcheck) */
396            /* we assume that we are out of the trap, reset the counters     */
397            /* and go back to the stepsize we had when we got trapped        */
398            if ((trapped[i]) && (newf < oldf * bndcheck)) {
399              trapped[i] = 0;
400              lbound[i] = 0;
401              rbounds[i] = 0;
402              delta[i] = initialstep[i];
403    
404            } else if (trialx[i] < (lowerb[i] + verysmall)) {
405              lbound[i]++;
406              trialx[i] = lowerb[i];
407              if (!trapped[i]) {
408                initialstep[i] = delta[i];
409                trapped[i] = 1;
410              }
411              /* if it has hit the bounds 2 times then increase the stepsize */
412              if (lbound[i] >= 2)
413                delta[i] /= rho;
414    
415            } else if (trialx[i] > (upperb[i] - verysmall)) {
416              rbounds[i]++;
417              trialx[i] = upperb[i];
418              if (!trapped[i]) {
419                initialstep[i] = delta[i];
420                trapped[i] = 1;
421              }
422              /* if it has hit the bounds 2 times then increase the stepsize */
423              if (rbounds[i] >= 2)
424                delta[i] /= rho;
425            }
426          }
427    
428          for (i = 0; i < nvars; i++) {
429            /* firstly, arrange the sign of delta[] */
430            if (trialx[i] < x[i])
431              delta[i] = 0.0 - fabs(delta[i]);
432            else
433              delta[i] = fabs(delta[i]);
434    
435            /* now, move further in this direction  */
436            tmp = x[i];
437            x[i] = trialx[i];
438            trialx[i] = trialx[i] + trialx[i] - tmp;
439          }
440    
441          /* only move forward if this is really an improvement    */
442          oldf = newf;
443          newf = EcoSystem->SimulateAndUpdate(trialx);
444    #ifndef NO_OPENMP
445          iters++;
446    #endif
447          if ((isEqual(newf, oldf)) || (newf > oldf)) {
448            newf = oldf;  //JMB no improvement, so reset the value of newf
449            break;
450          }
451    
452          /* OK, it's better, so update variables and look around  */
453          bestf = newf;
454          for (i = 0; i < nvars; i++)
455            x[i] = trialx[i];
456    
457    #ifndef NO_OPENMP
458          newf = this->bestNearbyOMP(delta, trialx, bestf, param);
459          if (newf == -1) {
460                    handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
461                    handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
462                    return;
463              }
464    #else
465      newf = this->bestNearby(delta, trialx, bestf, param);      newf = this->bestNearby(delta, trialx, bestf, param);
466    #endif
467          if (isEqual(newf, bestf))
468            break;
469    
470      /* if too many function evaluations occur, terminate the algorithm */      /* if too many function evaluations occur, terminate the algorithm */
471    #ifdef NO_OPENMP
472      iters = EcoSystem->getFuncEval() - offset;      iters = EcoSystem->getFuncEval() - offset;
473    #endif
474          if (iters > hookeiter) {
475            handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
476            handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
477            handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
478            handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations");
479            handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
480    
481            score = EcoSystem->SimulateAndUpdate(trialx);
482            handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
483            for (i = 0; i < nvars; i++)
484              bestx[i] = trialx[i] * init[i];
485            EcoSystem->storeVariables(score, bestx);
486            return;
487          }
488        }
489    
490    #ifdef NO_OPENMP
491        iters = EcoSystem->getFuncEval() - offset;
492    #endif
493        if (newf < bestf) {
494          for (i = 0; i < nvars; i++)
495            bestx[i] = x[i] * init[i];
496          bestf = newf;
497          handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
498          handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point");
499          EcoSystem->storeVariables(bestf, bestx);
500          EcoSystem->writeBestValues();
501    
502        } else
503          handle.logMessage(LOGINFO, "Checking convergence criteria after", iters, "function evaluations ...");
504    
505        /* if the step length is less than hookeeps, terminate the algorithm */
506        if (steplength < hookeeps) {
507          handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
508          handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
509          handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
510          handle.logMessage(LOGINFO, "The optimisation stopped because an optimum was found for this run");
511    
512          converge = 1;
513          score = bestf;
514          handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
515          EcoSystem->storeVariables(bestf, bestx);
516          return;
517        }
518    
519        steplength *= rho;
520        handle.logMessage(LOGINFO, "Reducing the steplength to", steplength);
521        for (i = 0; i < nvars; i++)
522          delta[i] *= rho;
523      }
524    }
525    
526    /* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
527    #ifdef GADGET_OPENMP
528    double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
529      double minf;
530      int i, j, k, ii;
531      DoubleVector z(point);
532      int bestId = 0;
533      struct Storage {
534              DoubleVector z;
535              DoubleVector delta;
536              double ftmp;
537              int iters;
538      };
539    
540      minf = prevbest;
541    
542      int paral_tokens, numThr, nvars = point.Size();
543      numThr = omp_get_max_threads ( );
544    
545      Storage* storage = new Storage[numThr];
546      if ((numThr % 2) == 0)
547              paral_tokens = numThr / 2;
548      else {
549              return -1;
550      }
551    
552      for (ii=0; ii< paral_tokens; ii++) {
553              i = 0;
554              while ( i < nvars) {
555                      if ((i + paral_tokens -1) >= nvars)
556                              paral_tokens = nvars - i;
557                      omp_set_dynamic(0);
558                      omp_set_nested(1);
559            #pragma omp parallel for num_threads(paral_tokens) private(k)
560                      for (j = 0; j < paral_tokens; ++j) {
561                              storage[j].z = z;
562                              storage[j].delta = delta;
563                              DoubleVector v1(z);
564                              DoubleVector v2(z);
565                              k = param[i+j];
566                              v1[k] +=  delta[k];
567                              v2[k] -=  delta[k];
568    
569    #pragma omp parallel sections num_threads(2)
570                              {
571              #pragma omp section
572                                      {
573                              storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
574                                      }
575              #pragma omp section
576                                      {
577                              storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
578                                      }
579                              }
580                              if (storage[j].ftmp < minf) {
581                                      storage[j].iters = 1;
582                                      storage[j].z[k] = v1[k];
583                              } else {
584                                      storage[j].iters = 2;
585                                      storage[j].delta[k] = 0.0 - delta[k];
586                                      if (storage[j+paral_tokens].ftmp < minf) {
587                                              storage[j].ftmp = storage[j+paral_tokens].ftmp;
588                                              storage[j].z[k] = v2[k];
589                                      }
590                                      else iters += 2;
591                              }
592                      }
593    
594                      bestId = 0;
595                      for (j = 1; j < paral_tokens; ++j) {
596                              if (storage[j].ftmp < storage[bestId].ftmp)
597                                      bestId = j;
598                      }
599                      if (storage[bestId].ftmp < minf) {
600                              iters += storage[bestId].iters;
601                              minf = storage[bestId].ftmp;
602                              z = storage[bestId].z;
603                              delta = storage[bestId].delta;
604                      }
605    
606                      i += paral_tokens;
607              }
608            }
609    
610      delete[] storage;
611      for (i = 0; i < nvars; ++i)
612        point[i] = z[i];
613    
614      return minf;
615    }
616    
617      void OptInfoHooke::OptimiseLikelihoodOMP() {
618              double oldf, newf, bestf, steplength, tmp;
619               int    i, offset;
620               int    rchange, rcheck, rnumber;  //Used to randomise the order of the parameters
621    
622               handle.logMessage(LOGINFO, "\nStarting Hooke & Jeeves optimisation algorithm\n");
623               int nvars = EcoSystem->numOptVariables();
624               DoubleVector x(nvars);
625               DoubleVector trialx(nvars);
626               DoubleVector bestx(nvars);
627               DoubleVector lowerb(nvars);
628               DoubleVector upperb(nvars);
629               DoubleVector init(nvars);
630               DoubleVector initialstep(nvars, rho);
631               DoubleVector delta(nvars);
632               IntVector param(nvars, 0);
633               IntVector lbound(nvars, 0);
634               IntVector rbounds(nvars, 0);
635               IntVector trapped(nvars, 0);
636    
637               EcoSystem->scaleVariables();
638               int numThr = omp_get_max_threads ( );
639               for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
640                      EcoSystems[i]->scaleVariables();
641               EcoSystem->getOptScaledValues(x);
642               EcoSystem->getOptLowerBounds(lowerb);
643               EcoSystem->getOptUpperBounds(upperb);
644               EcoSystem->getOptInitialValues(init);
645    
646               for (i = 0; i < nvars; i++) {
647                 // Scaling the bounds, because the parameters are scaled
648                 lowerb[i] = lowerb[i] / init[i];
649                 upperb[i] = upperb[i] / init[i];
650                 if (lowerb[i] > upperb[i]) {
651                   tmp = lowerb[i];
652                   lowerb[i] = upperb[i];
653                   upperb[i] = tmp;
654                 }
655    
656                 bestx[i] = x[i];
657                 trialx[i] = x[i];
658                 param[i] = i;
659                 delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign
660               }
661    
662               bestf = EcoSystem->SimulateAndUpdate(trialx);
663               if (bestf != bestf) { //check for NaN
664                 handle.logMessage(LOGINFO, "Error starting Hooke & Jeeves optimisation with f(x) = infinity");
665                 converge = -1;
666                 iters = 1;
667                 return;
668               }
669    
670               offset = EcoSystem->getFuncEval();  //number of function evaluations done before loop
671               newf = bestf;
672               oldf = bestf;
673               steplength = lambda;
674               if (isZero(steplength))
675                 steplength = rho;
676    
677               iters = 0;
678    
679               while (1) {
680                 if (isZero(bestf)) {
681             #ifdef NO_OPENMP
682                   iters = EcoSystem->getFuncEval() - offset;
683             #endif
684                   handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
685                   converge = -1;
686                   return;
687                 }
688    
689                 /* randomize the order of the parameters once in a while */
690                 rchange = 0;
691                 while (rchange < nvars) {
692                   rnumber = rand() % nvars;
693                   rcheck = 1;
694                   for (i = 0; i < rchange; i++)
695                     if (param[i] == rnumber)
696                       rcheck = 0;
697                   if (rcheck) {
698                     param[rchange] = rnumber;
699                     rchange++;
700                   }
701                 }
702    
703                 /* find best new point, one coord at a time */
704                 for (i = 0; i < nvars; i++)
705                   trialx[i] = x[i];
706             #ifndef NO_OPENMP
707                 newf = this->bestNearbyOMP2(delta, trialx, bestf, param);
708                 if (newf == -1) {
709                    handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
710                    handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
711                    return;
712                 }
713             #else
714                 newf = this->bestNearby(delta, trialx, bestf, param);
715             #endif
716                 /* if too many function evaluations occur, terminate the algorithm */
717    
718             #ifdef NO_OPENMP
719                 iters = EcoSystem->getFuncEval() - offset;
720             #endif
721      if (iters > hookeiter) {      if (iters > hookeiter) {
722        handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");        handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
723        handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");        handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
# Line 325  Line 785 
785        /* only move forward if this is really an improvement    */        /* only move forward if this is really an improvement    */
786        oldf = newf;        oldf = newf;
787        newf = EcoSystem->SimulateAndUpdate(trialx);        newf = EcoSystem->SimulateAndUpdate(trialx);
788             #ifndef NO_OPENMP
789                   iters++;
790             #endif
791        if ((isEqual(newf, oldf)) || (newf > oldf)) {        if ((isEqual(newf, oldf)) || (newf > oldf)) {
792          newf = oldf;  //JMB no improvement, so reset the value of newf          newf = oldf;  //JMB no improvement, so reset the value of newf
793          break;          break;
# Line 334  Line 797 
797        bestf = newf;        bestf = newf;
798        for (i = 0; i < nvars; i++)        for (i = 0; i < nvars; i++)
799          x[i] = trialx[i];          x[i] = trialx[i];
800    
801             #ifndef NO_OPENMP
802                   newf = this->bestNearbyOMP2(delta, trialx, bestf, param);
803                   if (newf == -1) {
804                            handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
805                            handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n");
806                            return;
807                       }
808             #else
809        newf = this->bestNearby(delta, trialx, bestf, param);        newf = this->bestNearby(delta, trialx, bestf, param);
810             #endif
811        if (isEqual(newf, bestf))        if (isEqual(newf, bestf))
812          break;          break;
813    
814        /* if too many function evaluations occur, terminate the algorithm */        /* if too many function evaluations occur, terminate the algorithm */
815             #ifdef NO_OPENMP
816        iters = EcoSystem->getFuncEval() - offset;        iters = EcoSystem->getFuncEval() - offset;
817             #endif
818        if (iters > hookeiter) {        if (iters > hookeiter) {
819          handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");          handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
820          handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");          handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
# Line 356  Line 831 
831        }        }
832      }      }
833    
834             #ifdef NO_OPENMP
835      iters = EcoSystem->getFuncEval() - offset;      iters = EcoSystem->getFuncEval() - offset;
836             #endif
837      if (newf < bestf) {      if (newf < bestf) {
838        for (i = 0; i < nvars; i++)        for (i = 0; i < nvars; i++)
839          bestx[i] = x[i] * init[i];          bestx[i] = x[i] * init[i];
# Line 389  Line 866 
866        delta[i] *= rho;        delta[i] *= rho;
867    }    }
868  }  }
869    #endif

Legend:
Removed from v.1  
changed lines
  Added in v.14

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

Powered By FusionForge