--- trunk/gadget/hooke.cc 2014/02/10 17:09:07 1 +++ trunk/gadget/hooke.cc 2015/07/23 19:00:38 11 @@ -1,20 +1,20 @@ /* Nonlinear Optimization using the algorithm of Hooke and Jeeves */ /* 12 February 1994 author: Mark G. Johnson */ - +// /* Find a point X where the nonlinear function f(X) has a local */ /* minimum. X is an n-vector and f(X) is a scalar. In mathe- */ /* matical notation f: R^n -> R^1. The objective function f() */ /* is not required to be continuous. Nor does f() need to be */ /* differentiable. The program does not use or require */ /* derivatives of f(). */ - +// /* The software user supplies three things: a subroutine that */ /* computes f(X), an initial "starting guess" of the minimum point */ /* X, and values for the algorithm convergence parameters. Then */ /* the program searches for a local minimum, beginning from the */ /* starting guess, using the Direct Search algorithm of Hooke and */ /* Jeeves. */ - +// /* This C program is adapted from the Algol pseudocode found in */ /* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun- */ /* ications of the ACM, Vol 6. p.313 (June 1963). It includes the */ @@ -24,7 +24,7 @@ /* highly as the one by A. Kaupe, is: R. Hooke and T. A. Jeeves, */ /* "Direct Search Solution of Numerical and Statistical Problems", */ /* Journal of the ACM, Vol. 8, April 1961, pp. 212-229. */ - +// /* Calling sequence: */ /* int hooke(nvars, startpt, endpt, rho, epsilon, itermax) */ /* */ @@ -59,15 +59,15 @@ /* itermax {an integer} A second, rarely used, halting */ /* criterion. If the algorithm uses >= itermax */ /* iterations, halt. */ - +// /* The user-supplied objective function f(x,n) should return a C */ /* "double". Its arguments are x -- an array of doubles, and */ /* n -- an integer. x is the point at which f(x) should be */ /* evaluated, and n is the number of coordinates of x. That is, */ /* n is the number of coefficients being fitted. */ - +// /* rho, the algorithm convergence control */ - +// /* The algorithm works by taking "steps" from one estimate of */ /* a minimum, to another (hopefully better) estimate. Taking */ /* big steps gets to the minimum more quickly, at the risk of */ @@ -95,14 +95,14 @@ /* again with a larger value of rho such as 0.85, using the */ /* result of the first minimization as the starting guess to */ /* begin the second minimization. */ - +// /* Normal use: */ /* (1) Code your function f() in the C language */ /* (2) Install your starting guess {or read it in} */ /* (3) Run the program */ /* (4) {for the skeptical}: Use the computed minimum */ /* as the starting point for another run */ - +// /* Data Fitting: */ /* Code your function f() to be the sum of the squares of the */ /* errors (differences) between the computed values and the */ @@ -113,7 +113,7 @@ /* f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i])) */ /* + (C*tan(t[i]))))^2 */ /* where x[] is a 3-vector consisting of {A, B, C}. */ - +// /* The author of this software is M.G. Johnson. */ /* Permission to use, copy, modify, and distribute this software */ /* for any purpose without fee is hereby granted, provided that */ @@ -125,7 +125,7 @@ /* AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY */ /* KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS */ /* FITNESS FOR ANY PARTICULAR PURPOSE. */ - +// /* JMB this has been modified to work with the gadget object structure */ /* This means that the function has been replaced by a call to ecosystem */ /* object, and we can use the vector objects that have been defined */ @@ -139,8 +139,14 @@ #include "ecosystem.h" #include "global.h" -extern Ecosystem* EcoSystem; +#ifndef NO_OPENMP +#include "omp.h" +#endif +extern Ecosystem* EcoSystem; +#ifndef NO_OPENMP +extern Ecosystem** EcoSystems; +#endif /* given a point, look for a better one nearby, one coord at a time */ double OptInfoHooke::bestNearby(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { @@ -150,9 +156,17 @@ DoubleVector z(point); minf = prevbest; +// for (int k=0;kSimulateAndUpdate(z); +// cout << i <<"-z["<< param[i]<<"]:" <= nvars) + paral_tokens = nvars - i; +//cout << "i:" << i << endl; + omp_set_dynamic(0); + omp_set_nested(1); +#pragma omp parallel for num_threads(paral_tokens) private(k) + for (j = 0; j < paral_tokens; ++j) { +// cout << "thr_for:" << omp_get_num_threads()<SimulateAndUpdate(v1); +// cout << "->" << j << endl; + } + #pragma omp section + { +// cout << "->" << j+paral_tokens << endl; + storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2); + } + } +//cout << "-----" << endl; +// cout <scaleVariables(); +#ifndef NO_OPENMP + int numThr = omp_get_max_threads ( ); + for (i = 0; i < numThr; i++) + EcoSystems[i]->scaleVariables(); +#endif EcoSystem->getOptScaledValues(x); EcoSystem->getOptLowerBounds(lowerb); EcoSystem->getOptUpperBounds(upperb); @@ -211,6 +337,7 @@ bestx[i] = x[i]; trialx[i] = x[i]; param[i] = i; + //FIXME rand_r delta[i] = ((2 * (rand() % 2)) - 1) * rho; //JMB - randomise the sign } @@ -229,9 +356,13 @@ if (isZero(steplength)) steplength = rho; + iters = 0; + while (1) { if (isZero(bestf)) { +#ifdef NO_OPENMP iters = EcoSystem->getFuncEval() - offset; +#endif handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0"); converge = -1; return; @@ -240,6 +371,7 @@ /* randomize the order of the parameters once in a while */ rchange = 0; while (rchange < nvars) { + //FIXME rand_r rnumber = rand() % nvars; rcheck = 1; for (i = 0; i < rchange; i++) @@ -254,10 +386,22 @@ /* find best new point, one coord at a time */ for (i = 0; i < nvars; i++) trialx[i] = x[i]; +#ifndef NO_OPENMP + newf = this->bestNearbyOMP(delta, trialx, bestf, param); + if (newf == -1) { + handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); + handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); + return; + } +#else newf = this->bestNearby(delta, trialx, bestf, param); - +#endif /* if too many function evaluations occur, terminate the algorithm */ + +#ifdef NO_OPENMP iters = EcoSystem->getFuncEval() - offset; +#endif +// cout << "newf:" << newf << "-" << iters< hookeiter) { handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); @@ -325,6 +469,9 @@ /* only move forward if this is really an improvement */ oldf = newf; newf = EcoSystem->SimulateAndUpdate(trialx); +#ifndef NO_OPENMP + iters++; +#endif if ((isEqual(newf, oldf)) || (newf > oldf)) { newf = oldf; //JMB no improvement, so reset the value of newf break; @@ -334,12 +481,25 @@ bestf = newf; for (i = 0; i < nvars; i++) x[i] = trialx[i]; + +#ifndef NO_OPENMP + newf = this->bestNearbyOMP(delta, trialx, bestf, param); + if (newf == -1) { + handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); + handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); + return; + } +#else newf = this->bestNearby(delta, trialx, bestf, param); +#endif if (isEqual(newf, bestf)) break; /* if too many function evaluations occur, terminate the algorithm */ +#ifdef NO_OPENMP iters = EcoSystem->getFuncEval() - offset; +#endif +// cout << "newf:" << newf << "-" << iters< hookeiter) { handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); @@ -356,7 +516,9 @@ } } +#ifdef NO_OPENMP iters = EcoSystem->getFuncEval() - offset; +#endif if (newf < bestf) { for (i = 0; i < nvars; i++) bestx[i] = x[i] * init[i]; @@ -389,3 +551,379 @@ delta[i] *= rho; } } + +#ifdef GADGET_OPENMP + //TODO doc +double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { +//cout << "beastNearbyINI" << endl; + double minf;//, ftmp; + int i, j, k, ii; + DoubleVector z(point); + int bestId = 0; + struct Storage { + DoubleVector z; + DoubleVector delta; + double ftmp; + int iters; + }; + + minf = prevbest; + + + int paral_tokens, numThr, nvars = point.Size(); + numThr = omp_get_max_threads ( ); + + Storage* storage = new Storage[numThr]; + if ((numThr % 2) == 0) + paral_tokens = numThr / 2; + else { + return -1; + } + + for (ii=0; ii< paral_tokens; ii++) { + i = 0; + while ( i < nvars) { + // for (int k=0;k= nvars) + paral_tokens = nvars - i; + //cout << "i:" << i << endl; + omp_set_dynamic(0); + omp_set_nested(1); + #pragma omp parallel for num_threads(paral_tokens) private(k) + for (j = 0; j < paral_tokens; ++j) { + // cout << "thr_for:" << omp_get_num_threads()<SimulateAndUpdate(v1); + // cout << "->" << j << endl; + } + #pragma omp section + { + // cout << "->" << j+paral_tokens << endl; + storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2); + } + } + //cout << "-----" << endl; + // cout <numOptVariables(); + DoubleVector x(nvars); + DoubleVector trialx(nvars); + DoubleVector bestx(nvars); + DoubleVector lowerb(nvars); + DoubleVector upperb(nvars); + DoubleVector init(nvars); + DoubleVector initialstep(nvars, rho); + DoubleVector delta(nvars); + IntVector param(nvars, 0); + IntVector lbound(nvars, 0); + IntVector rbounds(nvars, 0); + IntVector trapped(nvars, 0); + + EcoSystem->scaleVariables(); + #ifndef NO_OPENMP + int numThr = omp_get_max_threads ( ); + 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++) { + // 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; + } + + bestx[i] = x[i]; + trialx[i] = x[i]; + param[i] = i; + //FIXME rand_r + delta[i] = ((2 * (rand() % 2)) - 1) * rho; //JMB - randomise the sign + } + + bestf = EcoSystem->SimulateAndUpdate(trialx); + if (bestf != bestf) { //check for NaN + handle.logMessage(LOGINFO, "Error starting Hooke & Jeeves optimisation with f(x) = infinity"); + converge = -1; + iters = 1; + return; + } + + offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop + newf = bestf; + oldf = bestf; + steplength = lambda; + if (isZero(steplength)) + steplength = rho; + + iters = 0; + + while (1) { + if (isZero(bestf)) { + #ifdef NO_OPENMP + iters = EcoSystem->getFuncEval() - offset; + #endif + handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0"); + converge = -1; + return; + } + + /* randomize the order of the parameters once in a while */ + rchange = 0; + while (rchange < nvars) { + //FIXME rand_r + rnumber = rand() % nvars; + rcheck = 1; + for (i = 0; i < rchange; i++) + if (param[i] == rnumber) + rcheck = 0; + if (rcheck) { + param[rchange] = rnumber; + rchange++; + } + } + + /* find best new point, one coord at a time */ + for (i = 0; i < nvars; i++) + trialx[i] = x[i]; + #ifndef NO_OPENMP + newf = this->bestNearbyOMP2(delta, trialx, bestf, param); + if (newf == -1) { + handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); + handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); + return; + } + #else + newf = this->bestNearby(delta, trialx, bestf, param); + #endif + /* if too many function evaluations occur, terminate the algorithm */ + + #ifdef NO_OPENMP + iters = EcoSystem->getFuncEval() - offset; + #endif + // cout << "newf:" << newf << "-" << iters< hookeiter) { + handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); + handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); + handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); + handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations"); + handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); + + score = EcoSystem->SimulateAndUpdate(trialx); + handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); + for (i = 0; i < nvars; i++) + bestx[i] = trialx[i] * init[i]; + EcoSystem->storeVariables(score, bestx); +// EcoSystem->writeOptValuesOMP(); + return; + } + + /* if we made some improvements, pursue that direction */ + while (newf < bestf) { + for (i = 0; i < nvars; i++) { + /* if it has been trapped but f has now gotten better (bndcheck) */ + /* we assume that we are out of the trap, reset the counters */ + /* and go back to the stepsize we had when we got trapped */ + if ((trapped[i]) && (newf < oldf * bndcheck)) { + trapped[i] = 0; + lbound[i] = 0; + rbounds[i] = 0; + delta[i] = initialstep[i]; + + } else if (trialx[i] < (lowerb[i] + verysmall)) { + lbound[i]++; + trialx[i] = lowerb[i]; + if (!trapped[i]) { + initialstep[i] = delta[i]; + trapped[i] = 1; + } + /* if it has hit the bounds 2 times then increase the stepsize */ + if (lbound[i] >= 2) + delta[i] /= rho; + + } else if (trialx[i] > (upperb[i] - verysmall)) { + rbounds[i]++; + trialx[i] = upperb[i]; + if (!trapped[i]) { + initialstep[i] = delta[i]; + trapped[i] = 1; + } + /* if it has hit the bounds 2 times then increase the stepsize */ + if (rbounds[i] >= 2) + delta[i] /= rho; + } + } + + for (i = 0; i < nvars; i++) { + /* firstly, arrange the sign of delta[] */ + if (trialx[i] < x[i]) + delta[i] = 0.0 - fabs(delta[i]); + else + delta[i] = fabs(delta[i]); + + /* now, move further in this direction */ + tmp = x[i]; + x[i] = trialx[i]; + trialx[i] = trialx[i] + trialx[i] - tmp; + } + + /* only move forward if this is really an improvement */ + oldf = newf; + newf = EcoSystem->SimulateAndUpdate(trialx); + #ifndef NO_OPENMP + iters++; + #endif + if ((isEqual(newf, oldf)) || (newf > oldf)) { + newf = oldf; //JMB no improvement, so reset the value of newf + break; + } + + /* OK, it's better, so update variables and look around */ + bestf = newf; + for (i = 0; i < nvars; i++) + x[i] = trialx[i]; + + #ifndef NO_OPENMP + newf = this->bestNearbyOMP2(delta, trialx, bestf, param); + if (newf == -1) { + handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); + handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); + return; + } + #else + newf = this->bestNearby(delta, trialx, bestf, param); + #endif + if (isEqual(newf, bestf)) + break; + + /* if too many function evaluations occur, terminate the algorithm */ + #ifdef NO_OPENMP + iters = EcoSystem->getFuncEval() - offset; + #endif + // cout << "newf:" << newf << "-" << iters< hookeiter) { + handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); + handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); + handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); + handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations"); + handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); + + score = EcoSystem->SimulateAndUpdate(trialx); + handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); + for (i = 0; i < nvars; i++) + bestx[i] = trialx[i] * init[i]; + EcoSystem->storeVariables(score, bestx); +// EcoSystem->writeOptValuesOMP(); + return; + } + } + + #ifdef NO_OPENMP + iters = EcoSystem->getFuncEval() - offset; + #endif + if (newf < bestf) { + for (i = 0; i < nvars; i++) + bestx[i] = x[i] * init[i]; + bestf = newf; + 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(); + + } else + handle.logMessage(LOGINFO, "Checking convergence criteria after", iters, "function evaluations ..."); + + /* if the step length is less than hookeeps, terminate the algorithm */ + if (steplength < hookeeps) { + handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); + handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); + handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); + handle.logMessage(LOGINFO, "The optimisation stopped because an optimum was found for this run"); + + converge = 1; + score = bestf; + handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); + EcoSystem->storeVariables(bestf, bestx); +// EcoSystem->writeOptValuesOMP(); + return; + } + + steplength *= rho; + handle.logMessage(LOGINFO, "Reducing the steplength to", steplength); + for (i = 0; i < nvars; i++) + delta[i] *= rho; + } + } +#endif