Annotation of /trunk/gadget/hooke.cc

Parent Directory | Revision Log

---

Revision 16 - (view) (download)

1 :	agomez	1	/* Nonlinear Optimization using the algorithm of Hooke and Jeeves */
2 :			/* 12 February 1994 author: Mark G. Johnson */
3 :	ulcessvp	11	//
4 :	agomez	1	/* 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- */
6 :			/* matical notation f: R^n -> R^1. The objective function f() */
7 :			/* is not required to be continuous. Nor does f() need to be */
8 :			/* differentiable. The program does not use or require */
9 :			/* derivatives of f(). */
10 :	ulcessvp	11	//
11 :	agomez	1	/* The software user supplies three things: a subroutine that */
12 :			/* computes f(X), an initial "starting guess" of the minimum point */
13 :			/* X, and values for the algorithm convergence parameters. Then */
14 :			/* the program searches for a local minimum, beginning from the */
15 :			/* starting guess, using the Direct Search algorithm of Hooke and */
16 :			/* Jeeves. */
17 :	ulcessvp	11	//
18 :	agomez	1	/* This C program is adapted from the Algol pseudocode found in */
19 :			/* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun- */
20 :			/* ications of the ACM, Vol 6. p.313 (June 1963). It includes the */
21 :			/* improvements suggested by Bell and Pike (CACM v.9, p. 684, Sept */
22 :			/* 1966) and those of Tomlin and Smith, "Remark on Algorithm 178" */
23 :			/* (CACM v.12). The original paper, which I don't recommend as */
24 :			/* highly as the one by A. Kaupe, is: R. Hooke and T. A. Jeeves, */
25 :			/* "Direct Search Solution of Numerical and Statistical Problems", */
26 :			/* Journal of the ACM, Vol. 8, April 1961, pp. 212-229. */
27 :	ulcessvp	11	//
28 :	agomez	1	/* Calling sequence: */
29 :			/* int hooke(nvars, startpt, endpt, rho, epsilon, itermax) */
30 :			/* */
31 :			/* nvars {an integer} */
32 :			/* This is the number of dimensions in the domain of f(). */
33 :			/* It is the number of coordinates of the starting point */
34 :			/* (and the minimum point.) */
35 :			/* startpt {an array of doubles} */
36 :			/* This is the user-supplied guess at the minimum. */
37 :			/* endpt {an array of doubles} */
38 :			/* This is the calculated location of the local minimum */
39 :			/* rho {a double} */
40 :			/* This is a user-supplied convergence parameter (more */
41 :			/* detail below), which should be set to a value between */
42 :			/* 0.0 and 1.0. Larger values of rho give greater */
43 :			/* probability of convergence on highly nonlinear */
44 :			/* functions, at a cost of more function evaluations. */
45 :			/* Smaller values of rho reduces the number of evaluations */
46 :			/* (and the program running time), but increases the risk */
47 :			/* of nonconvergence. See below. */
48 :			/* epsilon {a double} */
49 :			/* This is the criterion for halting the search for a */
50 :			/* minimum. When the algorithm begins to make less and */
51 :			/* less progress on each iteration, it checks the halting */
52 :			/* criterion: if the stepsize is below epsilon, terminate */
53 :			/* the iteration and return the current best estimate of */
54 :			/* the minimum. Larger values of epsilon (such as 1.0e-4) */
55 :			/* give quicker running time, but a less accurate estimate */
56 :			/* of the minimum. Smaller values of epsilon (such as */
57 :			/* 1.0e-7) give longer running time, but a more accurate */
58 :			/* estimate of the minimum. */
59 :			/* itermax {an integer} A second, rarely used, halting */
60 :			/* criterion. If the algorithm uses >= itermax */
61 :			/* iterations, halt. */
62 :	ulcessvp	11	//
63 :	agomez	1	/* The user-supplied objective function f(x,n) should return a C */
64 :			/* "double". Its arguments are x -- an array of doubles, and */
65 :			/* 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, */
67 :			/* n is the number of coefficients being fitted. */
68 :	ulcessvp	11	//
69 :	agomez	1	/* rho, the algorithm convergence control */
70 :	ulcessvp	11	//
71 :	agomez	1	/* The algorithm works by taking "steps" from one estimate of */
72 :			/* a minimum, to another (hopefully better) estimate. Taking */
73 :			/* big steps gets to the minimum more quickly, at the risk of */
74 :			/* "stepping right over" an excellent point. The stepsize is */
75 :			/* controlled by a user supplied parameter called rho. At each */
76 :			/* iteration, the stepsize is multiplied by rho (0 < rho < 1), */
77 :			/* so the stepsize is successively reduced. */
78 :			/* Small values of rho correspond to big stepsize changes, */
79 :			/* which make the algorithm run more quickly. However, there */
80 :			/* is a chance (especially with highly nonlinear functions) */
81 :			/* that these big changes will accidentally overlook a */
82 :			/* promising search vector, leading to nonconvergence. */
83 :			/* Large values of rho correspond to small stepsize changes, */
84 :			/* which force the algorithm to carefully examine nearby points */
85 :			/* instead of optimistically forging ahead. This improves the */
86 :			/* probability of convergence. */
87 :			/* The stepsize is reduced until it is equal to (or smaller */
88 :			/* than) epsilon. So the number of iterations performed by */
89 :			/* Hooke-Jeeves is determined by rho and epsilon: */
90 :			/* rho**(number_of_iterations) = epsilon */
91 :			/* In general it is a good idea to set rho to an aggressively */
92 :			/* small value like 0.5 (hoping for fast convergence). Then, */
93 :			/* if the user suspects that the reported minimum is incorrect */
94 :			/* (or perhaps not accurate enough), the program can be run */
95 :			/* again with a larger value of rho such as 0.85, using the */
96 :			/* result of the first minimization as the starting guess to */
97 :			/* begin the second minimization. */
98 :	ulcessvp	11	//
99 :	agomez	1	/* Normal use: */
100 :			/* (1) Code your function f() in the C language */
101 :			/* (2) Install your starting guess {or read it in} */
102 :			/* (3) Run the program */
103 :			/* (4) {for the skeptical}: Use the computed minimum */
104 :			/* as the starting point for another run */
105 :	ulcessvp	11	//
106 :	agomez	1	/* Data Fitting: */
107 :			/* Code your function f() to be the sum of the squares of the */
108 :			/* errors (differences) between the computed values and the */
109 :			/* measured values. Then minimize f() using Hooke-Jeeves. */
110 :			/* EXAMPLE: you have 20 datapoints (ti, yi) and you want to */
111 :			/* find A,B,C such that (A*t*t) + (B*exp(t)) + (C*tan(t)) */
112 :			/* fits the data as closely as possible. Then f() is just */
113 :			/* f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i])) */
114 :			/* + (C*tan(t[i]))))^2 */
115 :			/* where x[] is a 3-vector consisting of {A, B, C}. */
116 :	ulcessvp	11	//
117 :	agomez	1	/* The author of this software is M.G. Johnson. */
118 :			/* Permission to use, copy, modify, and distribute this software */
119 :			/* for any purpose without fee is hereby granted, provided that */
120 :			/* this entire notice is included in all copies of any software */
121 :			/* which is or includes a copy or modification of this software */
122 :			/* and in all copies of the supporting documentation for such */
123 :			/* software. THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT */
124 :			/* ANY EXPRESS OR IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE */
125 :			/* AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY */
126 :			/* KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS */
127 :			/* FITNESS FOR ANY PARTICULAR PURPOSE. */
128 :	ulcessvp	11	//
129 :	agomez	1	/* 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 */
131 :			/* object, and we can use the vector objects that have been defined */
132 :
133 :			#include "gadget.h"
134 :			#include "optinfo.h"
135 :			#include "mathfunc.h"
136 :			#include "doublevector.h"
137 :			#include "intvector.h"
138 :			#include "errorhandler.h"
139 :			#include "ecosystem.h"
140 :			#include "global.h"
141 :
142 :	ulcessvp	16	#ifdef _OPENMP
143 :	ulcessvp	11	#include "omp.h"
144 :			#endif
145 :
146 :	agomez	1	extern Ecosystem* EcoSystem;
147 :	ulcessvp	16	#ifdef _OPENMP
148 :	ulcessvp	11	extern Ecosystem** EcoSystems;
149 :			#endif
150 :	agomez	1
151 :			/* 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) {
153 :
154 :			double minf, ftmp;
155 :			int i;
156 :			DoubleVector z(point);
157 :
158 :			minf = prevbest;
159 :			for (i = 0; i < point.Size(); i++) {
160 :			z[param[i]] = point[param[i]] + delta[param[i]];
161 :			ftmp = EcoSystem->SimulateAndUpdate(z);
162 :			if (ftmp < minf) {
163 :			minf = ftmp;
164 :			} else {
165 :			delta[param[i]] = 0.0 - delta[param[i]];
166 :			z[param[i]] = point[param[i]] + delta[param[i]];
167 :			ftmp = EcoSystem->SimulateAndUpdate(z);
168 :			if (ftmp < minf)
169 :			minf = ftmp;
170 :			else
171 :			z[param[i]] = point[param[i]];
172 :			}
173 :			}
174 :
175 :			for (i = 0; i < point.Size(); i++)
176 :			point[i] = z[i];
177 :			return minf;
178 :			}
179 :
180 :	ulcessvp	11	/* given a point, look for a better one nearby, one coord at a time */
181 :	ulcessvp	16	#ifdef _OPENMP
182 :	ulcessvp	12	/*
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 :	ulcessvp	15	double OptInfoHooke::bestNearbyRepro(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
188 :	ulcessvp	11	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 :	ulcessvp	15	omp_set_dynamic(0);
213 :			omp_set_nested(1); //permit the nested parallelization
214 :	ulcessvp	11	while ( i < nvars) {
215 :			if ((i + paral_tokens -1) >= nvars)
216 :			paral_tokens = nvars - i;
217 :	ulcessvp	14	#pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2)
218 :	ulcessvp	11	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 :	ulcessvp	14	#pragma omp parallel sections num_threads(2) //parrallelize the +/- delta simulation for each parameter
228 :	ulcessvp	11	{
229 :	ulcessvp	12	#pragma omp section
230 :	ulcessvp	11	{
231 :			storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
232 :			}
233 :	ulcessvp	12	#pragma omp section
234 :	ulcessvp	11	{
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 :	ulcessvp	15	delete[] storage;
264 :	ulcessvp	11	for (i = 0; i < nvars; ++i)
265 :			point[i] = z[i];
266 :			return minf;
267 :			}
268 :			#endif
269 :
270 :	agomez	1	void OptInfoHooke::OptimiseLikelihood() {
271 :
272 :			double oldf, newf, bestf, steplength, tmp;
273 :			int i, offset;
274 :			int rchange, rcheck, rnumber; //Used to randomise the order of the parameters
275 :
276 :			handle.logMessage(LOGINFO, "\nStarting Hooke & Jeeves optimisation algorithm\n");
277 :			int nvars = EcoSystem->numOptVariables();
278 :			DoubleVector x(nvars);
279 :			DoubleVector trialx(nvars);
280 :			DoubleVector bestx(nvars);
281 :			DoubleVector lowerb(nvars);
282 :			DoubleVector upperb(nvars);
283 :			DoubleVector init(nvars);
284 :			DoubleVector initialstep(nvars, rho);
285 :			DoubleVector delta(nvars);
286 :			IntVector param(nvars, 0);
287 :			IntVector lbound(nvars, 0);
288 :			IntVector rbounds(nvars, 0);
289 :			IntVector trapped(nvars, 0);
290 :
291 :			EcoSystem->scaleVariables();
292 :	ulcessvp	16	#ifdef _OPENMP
293 :	ulcessvp	11	int numThr = omp_get_max_threads ( );
294 :	ulcessvp	14	for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
295 :	ulcessvp	11	EcoSystems[i]->scaleVariables();
296 :			#endif
297 :	agomez	1	EcoSystem->getOptScaledValues(x);
298 :			EcoSystem->getOptLowerBounds(lowerb);
299 :			EcoSystem->getOptUpperBounds(upperb);
300 :			EcoSystem->getOptInitialValues(init);
301 :
302 :			for (i = 0; i < nvars; i++) {
303 :			// Scaling the bounds, because the parameters are scaled
304 :			lowerb[i] = lowerb[i] / init[i];
305 :			upperb[i] = upperb[i] / init[i];
306 :			if (lowerb[i] > upperb[i]) {
307 :			tmp = lowerb[i];
308 :			lowerb[i] = upperb[i];
309 :			upperb[i] = tmp;
310 :			}
311 :
312 :			bestx[i] = x[i];
313 :			trialx[i] = x[i];
314 :			param[i] = i;
315 :			delta[i] = ((2 * (rand() % 2)) - 1) * rho; //JMB - randomise the sign
316 :			}
317 :
318 :			bestf = EcoSystem->SimulateAndUpdate(trialx);
319 :			if (bestf != bestf) { //check for NaN
320 :			handle.logMessage(LOGINFO, "Error starting Hooke & Jeeves optimisation with f(x) = infinity");
321 :			converge = -1;
322 :			iters = 1;
323 :			return;
324 :			}
325 :
326 :			offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop
327 :			newf = bestf;
328 :			oldf = bestf;
329 :			steplength = lambda;
330 :			if (isZero(steplength))
331 :			steplength = rho;
332 :
333 :	ulcessvp	11	iters = 0;
334 :
335 :	agomez	1	while (1) {
336 :			if (isZero(bestf)) {
337 :	ulcessvp	16	#ifndef _OPENMP
338 :	agomez	1	iters = EcoSystem->getFuncEval() - offset;
339 :	ulcessvp	11	#endif
340 :	agomez	1	handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0");
341 :			converge = -1;
342 :			return;
343 :			}
344 :
345 :			/* randomize the order of the parameters once in a while */
346 :			rchange = 0;
347 :			while (rchange < nvars) {
348 :			rnumber = rand() % nvars;
349 :			rcheck = 1;
350 :			for (i = 0; i < rchange; i++)
351 :			if (param[i] == rnumber)
352 :			rcheck = 0;
353 :			if (rcheck) {
354 :			param[rchange] = rnumber;
355 :			rchange++;
356 :			}
357 :			}
358 :
359 :			/* find best new point, one coord at a time */
360 :			for (i = 0; i < nvars; i++)
361 :			trialx[i] = x[i];
362 :	ulcessvp	16	#ifdef _OPENMP
363 :	ulcessvp	15	newf = this->bestNearbyRepro(delta, trialx, bestf, param);
364 :	ulcessvp	11	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 :	agomez	1	newf = this->bestNearby(delta, trialx, bestf, param);
371 :	ulcessvp	11	#endif
372 :			/* if too many function evaluations occur, terminate the algorithm */
373 :	agomez	1
374 :	ulcessvp	16	#ifndef _OPENMP
375 :	agomez	1	iters = EcoSystem->getFuncEval() - offset;
376 :	ulcessvp	11	#endif
377 :	agomez	1	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 :	ulcessvp	16	#ifdef _OPENMP
445 :	ulcessvp	11	iters++;
446 :			#endif
447 :	agomez	1	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 :	ulcessvp	11
457 :	ulcessvp	16	#ifdef _OPENMP
458 :	ulcessvp	15	newf = this->bestNearbyRepro(delta, trialx, bestf, param);
459 :	ulcessvp	11	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 :	agomez	1	newf = this->bestNearby(delta, trialx, bestf, param);
466 :	ulcessvp	11	#endif
467 :	agomez	1	if (isEqual(newf, bestf))
468 :			break;
469 :
470 :			/* if too many function evaluations occur, terminate the algorithm */
471 :	ulcessvp	16	#ifndef _OPENMP
472 :	agomez	1	iters = EcoSystem->getFuncEval() - offset;
473 :	ulcessvp	11	#endif
474 :	agomez	1	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 :	ulcessvp	16	} // while (newf < bestf)
489 :	agomez	1
490 :	ulcessvp	16	#ifndef _OPENMP
491 :	agomez	1	iters = EcoSystem->getFuncEval() - offset;
492 :	ulcessvp	11	#endif
493 :	agomez	1	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 :	ulcessvp	11
526 :	ulcessvp	12	/* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
527 :	ulcessvp	15	#ifdef SPECULATIVE
528 :			double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
529 :	ulcessvp	12	double minf;
530 :	ulcessvp	11	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 :	ulcessvp	15	omp_set_dynamic(0);
553 :			omp_set_nested(1); //permit the nested parallelization
554 :	ulcessvp	11	for (ii=0; ii< paral_tokens; ii++) {
555 :			i = 0;
556 :			while ( i < nvars) {
557 :			if ((i + paral_tokens -1) >= nvars)
558 :			paral_tokens = nvars - i;
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 :	ulcessvp	12	#pragma omp parallel sections num_threads(2)
570 :	ulcessvp	11	{
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 :	ulcessvp	14	for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
640 :	ulcessvp	11	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 :	ulcessvp	16	#ifndef _OPENMP
682 :	ulcessvp	11	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 :	ulcessvp	16	#ifdef _OPENMP
707 :	ulcessvp	15	newf = this->bestNearbySpec(delta, trialx, bestf, param);
708 :	ulcessvp	11	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 :	ulcessvp	16	#ifndef _OPENMP
719 :	ulcessvp	11	iters = EcoSystem->getFuncEval() - offset;
720 :			#endif
721 :			if (iters > hookeiter) {
722 :			handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
723 :			handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
724 :			handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
725 :			handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations");
726 :			handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
727 :
728 :			score = EcoSystem->SimulateAndUpdate(trialx);
729 :			handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
730 :			for (i = 0; i < nvars; i++)
731 :			bestx[i] = trialx[i] * init[i];
732 :			EcoSystem->storeVariables(score, bestx);
733 :			return;
734 :			}
735 :
736 :			/* if we made some improvements, pursue that direction */
737 :			while (newf < bestf) {
738 :			for (i = 0; i < nvars; i++) {
739 :			/* if it has been trapped but f has now gotten better (bndcheck) */
740 :			/* we assume that we are out of the trap, reset the counters */
741 :			/* and go back to the stepsize we had when we got trapped */
742 :			if ((trapped[i]) && (newf < oldf * bndcheck)) {
743 :			trapped[i] = 0;
744 :			lbound[i] = 0;
745 :			rbounds[i] = 0;
746 :			delta[i] = initialstep[i];
747 :
748 :			} else if (trialx[i] < (lowerb[i] + verysmall)) {
749 :			lbound[i]++;
750 :			trialx[i] = lowerb[i];
751 :			if (!trapped[i]) {
752 :			initialstep[i] = delta[i];
753 :			trapped[i] = 1;
754 :			}
755 :			/* if it has hit the bounds 2 times then increase the stepsize */
756 :			if (lbound[i] >= 2)
757 :			delta[i] /= rho;
758 :
759 :			} else if (trialx[i] > (upperb[i] - verysmall)) {
760 :			rbounds[i]++;
761 :			trialx[i] = upperb[i];
762 :			if (!trapped[i]) {
763 :			initialstep[i] = delta[i];
764 :			trapped[i] = 1;
765 :			}
766 :			/* if it has hit the bounds 2 times then increase the stepsize */
767 :			if (rbounds[i] >= 2)
768 :			delta[i] /= rho;
769 :			}
770 :			}
771 :
772 :			for (i = 0; i < nvars; i++) {
773 :			/* firstly, arrange the sign of delta[] */
774 :			if (trialx[i] < x[i])
775 :			delta[i] = 0.0 - fabs(delta[i]);
776 :			else
777 :			delta[i] = fabs(delta[i]);
778 :
779 :			/* now, move further in this direction */
780 :			tmp = x[i];
781 :			x[i] = trialx[i];
782 :			trialx[i] = trialx[i] + trialx[i] - tmp;
783 :			}
784 :
785 :			/* only move forward if this is really an improvement */
786 :			oldf = newf;
787 :			newf = EcoSystem->SimulateAndUpdate(trialx);
788 :	ulcessvp	16	#ifdef _OPENMP
789 :	ulcessvp	11	iters++;
790 :			#endif
791 :			if ((isEqual(newf, oldf)) || (newf > oldf)) {
792 :			newf = oldf; //JMB no improvement, so reset the value of newf
793 :			break;
794 :			}
795 :
796 :			/* OK, it's better, so update variables and look around */
797 :			bestf = newf;
798 :			for (i = 0; i < nvars; i++)
799 :			x[i] = trialx[i];
800 :
801 :	ulcessvp	16	#ifdef _OPENMP
802 :	ulcessvp	15	newf = this->bestNearbySpec(delta, trialx, bestf, param);
803 :	ulcessvp	11	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);
810 :			#endif
811 :			if (isEqual(newf, bestf))
812 :			break;
813 :
814 :			/* if too many function evaluations occur, terminate the algorithm */
815 :	ulcessvp	16	#ifndef _OPENMP
816 :	ulcessvp	11	iters = EcoSystem->getFuncEval() - offset;
817 :			#endif
818 :			if (iters > hookeiter) {
819 :			handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
820 :			handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
821 :			handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
822 :			handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations");
823 :			handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
824 :
825 :			score = EcoSystem->SimulateAndUpdate(trialx);
826 :			handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
827 :			for (i = 0; i < nvars; i++)
828 :			bestx[i] = trialx[i] * init[i];
829 :			EcoSystem->storeVariables(score, bestx);
830 :			return;
831 :			}
832 :			}
833 :
834 :	ulcessvp	16	#ifndef _OPENMP
835 :	ulcessvp	11	iters = EcoSystem->getFuncEval() - offset;
836 :			#endif
837 :			if (newf < bestf) {
838 :			for (i = 0; i < nvars; i++)
839 :			bestx[i] = x[i] * init[i];
840 :			bestf = newf;
841 :			handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
842 :			handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point");
843 :			EcoSystem->storeVariables(bestf, bestx);
844 :			EcoSystem->writeBestValues();
845 :
846 :			} else
847 :			handle.logMessage(LOGINFO, "Checking convergence criteria after", iters, "function evaluations ...");
848 :
849 :			/* if the step length is less than hookeeps, terminate the algorithm */
850 :			if (steplength < hookeeps) {
851 :			handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
852 :			handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
853 :			handle.logMessage(LOGINFO, "The steplength was reduced to", steplength);
854 :			handle.logMessage(LOGINFO, "The optimisation stopped because an optimum was found for this run");
855 :
856 :			converge = 1;
857 :			score = bestf;
858 :			handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
859 :			EcoSystem->storeVariables(bestf, bestx);
860 :			return;
861 :			}
862 :
863 :			steplength *= rho;
864 :			handle.logMessage(LOGINFO, "Reducing the steplength to", steplength);
865 :			for (i = 0; i < nvars; i++)
866 :			delta[i] *= rho;
867 :			}
868 :			}
869 :			#endif

---

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