139 |
#include "ecosystem.h" |
#include "ecosystem.h" |
140 |
#include "global.h" |
#include "global.h" |
141 |
|
|
142 |
#ifndef NO_OPENMP |
#ifdef _OPENMP |
143 |
#include "omp.h" |
#include "omp.h" |
144 |
#endif |
#endif |
145 |
|
|
146 |
extern Ecosystem* EcoSystem; |
extern Ecosystem* EcoSystem; |
147 |
#ifndef NO_OPENMP |
#ifdef _OPENMP |
148 |
extern Ecosystem** EcoSystems; |
extern Ecosystem** EcoSystems; |
149 |
#endif |
#endif |
150 |
|
|
178 |
} |
} |
179 |
|
|
180 |
/* 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 */ |
181 |
#ifndef NO_OPENMP |
#ifdef _OPENMP |
182 |
/* |
/* |
183 |
* function bestBeraby parallelized with OpenMP |
* function bestBeraby parallelized with OpenMP |
184 |
* · 2 threads per coord to parallelize the calculation of +delta/-delta |
* · 2 threads per coord to parallelize the calculation of +delta/-delta |
185 |
* · parallelize the calculation of the best nearby of the coord |
* · parallelize the calculation of the best nearby of the coord |
186 |
*/ |
*/ |
187 |
double OptInfoHooke::bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { |
double OptInfoHooke::bestNearbyRepro(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { |
188 |
double minf;//, ftmp; |
double minf;//, ftmp; |
189 |
int i, j, k; |
int i, j, k; |
190 |
DoubleVector z(point); |
DoubleVector z(point); |
209 |
return -1; |
return -1; |
210 |
} |
} |
211 |
|
|
212 |
|
// omp_set_dynamic(0); |
213 |
|
// omp_set_nested(1); //permit the nested parallelization |
214 |
while ( i < nvars) { |
while ( i < nvars) { |
215 |
if ((i + paral_tokens -1) >= nvars) |
if ((i + paral_tokens -1) >= nvars) |
216 |
paral_tokens = nvars - i; |
paral_tokens = nvars - i; |
217 |
omp_set_dynamic(0); |
#pragma omp parallel for num_threads(paral_tokens*2) private(k) //parallelize the parameters (numThr) |
218 |
omp_set_nested(1); //permit the nested parallelization |
for (j = 0; j < (paral_tokens*2); ++j) { |
|
#pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2) |
|
|
for (j = 0; j < paral_tokens; ++j) { |
|
219 |
storage[j].z = z; |
storage[j].z = z; |
220 |
storage[j].delta = delta; |
storage[j].delta = delta; |
221 |
DoubleVector v1(z); |
DoubleVector v(z); |
|
DoubleVector v2(z); |
|
|
k = param[i+j]; |
|
|
v1[k] += delta[k]; |
|
|
v2[k] -= delta[k]; |
|
222 |
|
|
223 |
#pragma omp parallel sections num_threads(2) //parrallelize the +/- delta simulation for each parameter |
if (j<paral_tokens) { |
224 |
{ |
k = param[i+j]; |
225 |
#pragma omp section |
v[k] += delta[k]; |
|
{ |
|
|
storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1); |
|
|
} |
|
|
#pragma omp section |
|
|
{ |
|
|
storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2); |
|
226 |
} |
} |
227 |
|
else { |
228 |
|
k = param[i+j-paral_tokens]; |
229 |
|
v[k] -= delta[k]; |
230 |
} |
} |
231 |
|
|
232 |
|
storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v); |
233 |
|
storage[j].z[k] = v[k]; |
234 |
|
} |
235 |
|
for (j = 0; j < paral_tokens; ++j) { |
236 |
|
k = param[i+j]; |
237 |
if (storage[j].ftmp < minf) { |
if (storage[j].ftmp < minf) { |
238 |
storage[j].iters = 1; |
storage[j].iters = 1; |
239 |
storage[j].z[k] = v1[k]; |
// storage[j].z[k] = v1[k]; |
240 |
} else { |
} else { |
241 |
storage[j].iters = 2; |
storage[j].iters = 2; |
242 |
storage[j].delta[k] = 0.0 - delta[k]; |
storage[j].delta[k] = 0.0 - delta[k]; |
243 |
if (storage[j+paral_tokens].ftmp < minf) { |
if (storage[j+paral_tokens].ftmp < minf) { |
244 |
storage[j].ftmp = storage[j+paral_tokens].ftmp; |
storage[j].ftmp = storage[j+paral_tokens].ftmp; |
245 |
storage[j].z[k] = v2[k]; |
storage[j].z[k] = storage[j+paral_tokens].z[k];; |
246 |
} |
} |
247 |
} |
} |
248 |
} |
} |
258 |
} |
} |
259 |
} |
} |
260 |
} |
} |
261 |
|
delete[] storage; |
262 |
for (i = 0; i < nvars; ++i) |
for (i = 0; i < nvars; ++i) |
263 |
point[i] = z[i]; |
point[i] = z[i]; |
264 |
return minf; |
return minf; |
265 |
} |
} |
266 |
#endif |
void OptInfoHooke::OptimiseLikelihoodREP() { |
|
|
|
|
void OptInfoHooke::OptimiseLikelihood() { |
|
267 |
|
|
268 |
double oldf, newf, bestf, steplength, tmp; |
double oldf, newf, bestf, steplength, tmp; |
269 |
int i, offset; |
int i, offset; |
285 |
IntVector trapped(nvars, 0); |
IntVector trapped(nvars, 0); |
286 |
|
|
287 |
EcoSystem->scaleVariables(); |
EcoSystem->scaleVariables(); |
|
#ifndef NO_OPENMP |
|
288 |
int numThr = omp_get_max_threads ( ); |
int numThr = omp_get_max_threads ( ); |
289 |
for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread |
for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread |
290 |
EcoSystems[i]->scaleVariables(); |
EcoSystems[i]->scaleVariables(); |
|
#endif |
|
291 |
EcoSystem->getOptScaledValues(x); |
EcoSystem->getOptScaledValues(x); |
292 |
EcoSystem->getOptLowerBounds(lowerb); |
EcoSystem->getOptLowerBounds(lowerb); |
293 |
EcoSystem->getOptUpperBounds(upperb); |
EcoSystem->getOptUpperBounds(upperb); |
328 |
|
|
329 |
while (1) { |
while (1) { |
330 |
if (isZero(bestf)) { |
if (isZero(bestf)) { |
|
#ifdef NO_OPENMP |
|
331 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
|
#endif |
|
332 |
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"); |
333 |
converge = -1; |
converge = -1; |
334 |
return; |
return; |
351 |
/* find best new point, one coord at a time */ |
/* find best new point, one coord at a time */ |
352 |
for (i = 0; i < nvars; i++) |
for (i = 0; i < nvars; i++) |
353 |
trialx[i] = x[i]; |
trialx[i] = x[i]; |
354 |
#ifndef NO_OPENMP |
newf = this->bestNearbyRepro(delta, trialx, bestf, param); |
|
newf = this->bestNearbyOMP(delta, trialx, bestf, param); |
|
355 |
if (newf == -1) { |
if (newf == -1) { |
356 |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
357 |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
358 |
return; |
return; |
359 |
} |
} |
|
#else |
|
|
newf = this->bestNearby(delta, trialx, bestf, param); |
|
|
#endif |
|
360 |
/* if too many function evaluations occur, terminate the algorithm */ |
/* if too many function evaluations occur, terminate the algorithm */ |
361 |
|
|
|
#ifdef NO_OPENMP |
|
362 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
|
#endif |
|
363 |
if (iters > hookeiter) { |
if (iters > hookeiter) { |
364 |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
365 |
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
427 |
/* only move forward if this is really an improvement */ |
/* only move forward if this is really an improvement */ |
428 |
oldf = newf; |
oldf = newf; |
429 |
newf = EcoSystem->SimulateAndUpdate(trialx); |
newf = EcoSystem->SimulateAndUpdate(trialx); |
|
#ifndef NO_OPENMP |
|
430 |
iters++; |
iters++; |
|
#endif |
|
431 |
if ((isEqual(newf, oldf)) || (newf > oldf)) { |
if ((isEqual(newf, oldf)) || (newf > oldf)) { |
432 |
newf = oldf; //JMB no improvement, so reset the value of newf |
newf = oldf; //JMB no improvement, so reset the value of newf |
433 |
break; |
break; |
438 |
for (i = 0; i < nvars; i++) |
for (i = 0; i < nvars; i++) |
439 |
x[i] = trialx[i]; |
x[i] = trialx[i]; |
440 |
|
|
441 |
#ifndef NO_OPENMP |
newf = this->bestNearbyRepro(delta, trialx, bestf, param); |
|
newf = this->bestNearbyOMP(delta, trialx, bestf, param); |
|
442 |
if (newf == -1) { |
if (newf == -1) { |
443 |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
444 |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
445 |
return; |
return; |
446 |
} |
} |
|
#else |
|
|
newf = this->bestNearby(delta, trialx, bestf, param); |
|
|
#endif |
|
447 |
if (isEqual(newf, bestf)) |
if (isEqual(newf, bestf)) |
448 |
break; |
break; |
449 |
|
|
450 |
/* if too many function evaluations occur, terminate the algorithm */ |
/* if too many function evaluations occur, terminate the algorithm */ |
|
#ifdef NO_OPENMP |
|
451 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
452 |
|
if (iters > hookeiter) { |
453 |
|
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
454 |
|
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
455 |
|
handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); |
456 |
|
handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations"); |
457 |
|
handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); |
458 |
|
|
459 |
|
score = EcoSystem->SimulateAndUpdate(trialx); |
460 |
|
handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); |
461 |
|
for (i = 0; i < nvars; i++) |
462 |
|
bestx[i] = trialx[i] * init[i]; |
463 |
|
EcoSystem->storeVariables(score, bestx); |
464 |
|
return; |
465 |
|
} |
466 |
|
} // while (newf < bestf) |
467 |
|
|
468 |
|
iters = EcoSystem->getFuncEval() - offset; |
469 |
|
if (newf < bestf) { |
470 |
|
for (i = 0; i < nvars; i++) |
471 |
|
bestx[i] = x[i] * init[i]; |
472 |
|
bestf = newf; |
473 |
|
handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations"); |
474 |
|
handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point"); |
475 |
|
EcoSystem->storeVariables(bestf, bestx); |
476 |
|
EcoSystem->writeBestValues(); |
477 |
|
|
478 |
|
} else |
479 |
|
handle.logMessage(LOGINFO, "Checking convergence criteria after", iters, "function evaluations ..."); |
480 |
|
|
481 |
|
/* if the step length is less than hookeeps, terminate the algorithm */ |
482 |
|
if (steplength < hookeeps) { |
483 |
|
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
484 |
|
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
485 |
|
handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); |
486 |
|
handle.logMessage(LOGINFO, "The optimisation stopped because an optimum was found for this run"); |
487 |
|
|
488 |
|
converge = 1; |
489 |
|
score = bestf; |
490 |
|
handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); |
491 |
|
EcoSystem->storeVariables(bestf, bestx); |
492 |
|
return; |
493 |
|
} |
494 |
|
|
495 |
|
steplength *= rho; |
496 |
|
handle.logMessage(LOGINFO, "Reducing the steplength to", steplength); |
497 |
|
for (i = 0; i < nvars; i++) |
498 |
|
delta[i] *= rho; |
499 |
|
} |
500 |
|
} |
501 |
#endif |
#endif |
502 |
|
|
503 |
|
void OptInfoHooke::OptimiseLikelihood() { |
504 |
|
|
505 |
|
double oldf, newf, bestf, steplength, tmp; |
506 |
|
int i, offset; |
507 |
|
int rchange, rcheck, rnumber; //Used to randomise the order of the parameters |
508 |
|
|
509 |
|
handle.logMessage(LOGINFO, "\nStarting Hooke & Jeeves optimisation algorithm\n"); |
510 |
|
int nvars = EcoSystem->numOptVariables(); |
511 |
|
DoubleVector x(nvars); |
512 |
|
DoubleVector trialx(nvars); |
513 |
|
DoubleVector bestx(nvars); |
514 |
|
DoubleVector lowerb(nvars); |
515 |
|
DoubleVector upperb(nvars); |
516 |
|
DoubleVector init(nvars); |
517 |
|
DoubleVector initialstep(nvars, rho); |
518 |
|
DoubleVector delta(nvars); |
519 |
|
IntVector param(nvars, 0); |
520 |
|
IntVector lbound(nvars, 0); |
521 |
|
IntVector rbounds(nvars, 0); |
522 |
|
IntVector trapped(nvars, 0); |
523 |
|
|
524 |
|
EcoSystem->scaleVariables(); |
525 |
|
EcoSystem->getOptScaledValues(x); |
526 |
|
EcoSystem->getOptLowerBounds(lowerb); |
527 |
|
EcoSystem->getOptUpperBounds(upperb); |
528 |
|
EcoSystem->getOptInitialValues(init); |
529 |
|
|
530 |
|
for (i = 0; i < nvars; i++) { |
531 |
|
// Scaling the bounds, because the parameters are scaled |
532 |
|
lowerb[i] = lowerb[i] / init[i]; |
533 |
|
upperb[i] = upperb[i] / init[i]; |
534 |
|
if (lowerb[i] > upperb[i]) { |
535 |
|
tmp = lowerb[i]; |
536 |
|
lowerb[i] = upperb[i]; |
537 |
|
upperb[i] = tmp; |
538 |
|
} |
539 |
|
|
540 |
|
bestx[i] = x[i]; |
541 |
|
trialx[i] = x[i]; |
542 |
|
param[i] = i; |
543 |
|
delta[i] = ((2 * (rand() % 2)) - 1) * rho; //JMB - randomise the sign |
544 |
|
} |
545 |
|
|
546 |
|
bestf = EcoSystem->SimulateAndUpdate(trialx); |
547 |
|
if (bestf != bestf) { //check for NaN |
548 |
|
handle.logMessage(LOGINFO, "Error starting Hooke & Jeeves optimisation with f(x) = infinity"); |
549 |
|
converge = -1; |
550 |
|
iters = 1; |
551 |
|
return; |
552 |
|
} |
553 |
|
|
554 |
|
offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop |
555 |
|
newf = bestf; |
556 |
|
oldf = bestf; |
557 |
|
steplength = lambda; |
558 |
|
if (isZero(steplength)) |
559 |
|
steplength = rho; |
560 |
|
|
561 |
|
iters = 0; |
562 |
|
|
563 |
|
while (1) { |
564 |
|
if (isZero(bestf)) { |
565 |
|
iters = EcoSystem->getFuncEval() - offset; |
566 |
|
handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0"); |
567 |
|
converge = -1; |
568 |
|
return; |
569 |
|
} |
570 |
|
|
571 |
|
/* randomize the order of the parameters once in a while */ |
572 |
|
rchange = 0; |
573 |
|
while (rchange < nvars) { |
574 |
|
rnumber = rand() % nvars; |
575 |
|
rcheck = 1; |
576 |
|
for (i = 0; i < rchange; i++) |
577 |
|
if (param[i] == rnumber) |
578 |
|
rcheck = 0; |
579 |
|
if (rcheck) { |
580 |
|
param[rchange] = rnumber; |
581 |
|
rchange++; |
582 |
|
} |
583 |
|
} |
584 |
|
|
585 |
|
/* find best new point, one coord at a time */ |
586 |
|
for (i = 0; i < nvars; i++) |
587 |
|
trialx[i] = x[i]; |
588 |
|
newf = this->bestNearby(delta, trialx, bestf, param); |
589 |
|
/* if too many function evaluations occur, terminate the algorithm */ |
590 |
|
|
591 |
|
iters = EcoSystem->getFuncEval() - offset; |
592 |
if (iters > hookeiter) { |
if (iters > hookeiter) { |
593 |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
594 |
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
603 |
EcoSystem->storeVariables(score, bestx); |
EcoSystem->storeVariables(score, bestx); |
604 |
return; |
return; |
605 |
} |
} |
606 |
|
|
607 |
|
/* if we made some improvements, pursue that direction */ |
608 |
|
while (newf < bestf) { |
609 |
|
for (i = 0; i < nvars; i++) { |
610 |
|
/* if it has been trapped but f has now gotten better (bndcheck) */ |
611 |
|
/* we assume that we are out of the trap, reset the counters */ |
612 |
|
/* and go back to the stepsize we had when we got trapped */ |
613 |
|
if ((trapped[i]) && (newf < oldf * bndcheck)) { |
614 |
|
trapped[i] = 0; |
615 |
|
lbound[i] = 0; |
616 |
|
rbounds[i] = 0; |
617 |
|
delta[i] = initialstep[i]; |
618 |
|
|
619 |
|
} else if (trialx[i] < (lowerb[i] + verysmall)) { |
620 |
|
lbound[i]++; |
621 |
|
trialx[i] = lowerb[i]; |
622 |
|
if (!trapped[i]) { |
623 |
|
initialstep[i] = delta[i]; |
624 |
|
trapped[i] = 1; |
625 |
|
} |
626 |
|
/* if it has hit the bounds 2 times then increase the stepsize */ |
627 |
|
if (lbound[i] >= 2) |
628 |
|
delta[i] /= rho; |
629 |
|
|
630 |
|
} else if (trialx[i] > (upperb[i] - verysmall)) { |
631 |
|
rbounds[i]++; |
632 |
|
trialx[i] = upperb[i]; |
633 |
|
if (!trapped[i]) { |
634 |
|
initialstep[i] = delta[i]; |
635 |
|
trapped[i] = 1; |
636 |
|
} |
637 |
|
/* if it has hit the bounds 2 times then increase the stepsize */ |
638 |
|
if (rbounds[i] >= 2) |
639 |
|
delta[i] /= rho; |
640 |
|
} |
641 |
|
} |
642 |
|
|
643 |
|
for (i = 0; i < nvars; i++) { |
644 |
|
/* firstly, arrange the sign of delta[] */ |
645 |
|
if (trialx[i] < x[i]) |
646 |
|
delta[i] = 0.0 - fabs(delta[i]); |
647 |
|
else |
648 |
|
delta[i] = fabs(delta[i]); |
649 |
|
|
650 |
|
/* now, move further in this direction */ |
651 |
|
tmp = x[i]; |
652 |
|
x[i] = trialx[i]; |
653 |
|
trialx[i] = trialx[i] + trialx[i] - tmp; |
654 |
|
} |
655 |
|
|
656 |
|
/* only move forward if this is really an improvement */ |
657 |
|
oldf = newf; |
658 |
|
newf = EcoSystem->SimulateAndUpdate(trialx); |
659 |
|
if ((isEqual(newf, oldf)) || (newf > oldf)) { |
660 |
|
newf = oldf; //JMB no improvement, so reset the value of newf |
661 |
|
break; |
662 |
} |
} |
663 |
|
|
664 |
#ifdef NO_OPENMP |
/* OK, it's better, so update variables and look around */ |
665 |
|
bestf = newf; |
666 |
|
for (i = 0; i < nvars; i++) |
667 |
|
x[i] = trialx[i]; |
668 |
|
newf = this->bestNearby(delta, trialx, bestf, param); |
669 |
|
if (isEqual(newf, bestf)) |
670 |
|
break; |
671 |
|
|
672 |
|
/* if too many function evaluations occur, terminate the algorithm */ |
673 |
|
iters = EcoSystem->getFuncEval() - offset; |
674 |
|
if (iters > hookeiter) { |
675 |
|
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
676 |
|
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
677 |
|
handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); |
678 |
|
handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations"); |
679 |
|
handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); |
680 |
|
|
681 |
|
score = EcoSystem->SimulateAndUpdate(trialx); |
682 |
|
handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); |
683 |
|
for (i = 0; i < nvars; i++) |
684 |
|
bestx[i] = trialx[i] * init[i]; |
685 |
|
EcoSystem->storeVariables(score, bestx); |
686 |
|
return; |
687 |
|
} |
688 |
|
} // while (newf < bestf) |
689 |
|
|
690 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
|
#endif |
|
691 |
if (newf < bestf) { |
if (newf < bestf) { |
692 |
for (i = 0; i < nvars; i++) |
for (i = 0; i < nvars; i++) |
693 |
bestx[i] = x[i] * init[i]; |
bestx[i] = x[i] * init[i]; |
722 |
} |
} |
723 |
|
|
724 |
/* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/ |
/* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/ |
725 |
#ifdef GADGET_OPENMP |
#ifdef _OPENMP |
726 |
double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { |
//#ifdef SPECULATIVE |
727 |
|
double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { |
728 |
double minf; |
double minf; |
729 |
int i, j, k, ii; |
int i, j, k, ii; |
730 |
DoubleVector z(point); |
DoubleVector z(point); |
748 |
return -1; |
return -1; |
749 |
} |
} |
750 |
|
|
751 |
|
// omp_set_dynamic(0); |
752 |
|
// omp_set_nested(1); //permit the nested parallelization |
753 |
for (ii=0; ii< paral_tokens; ii++) { |
for (ii=0; ii< paral_tokens; ii++) { |
754 |
i = 0; |
i = 0; |
755 |
while ( i < nvars) { |
while ( i < nvars) { |
756 |
if ((i + paral_tokens -1) >= nvars) |
if ((i + paral_tokens -1) >= nvars) |
757 |
paral_tokens = nvars - i; |
paral_tokens = nvars - i; |
758 |
omp_set_dynamic(0); |
#pragma omp parallel for num_threads(paral_tokens*2) private(k) |
759 |
omp_set_nested(1); |
for (j = 0; j < (paral_tokens*2); ++j) { |
|
#pragma omp parallel for num_threads(paral_tokens) private(k) |
|
|
for (j = 0; j < paral_tokens; ++j) { |
|
760 |
storage[j].z = z; |
storage[j].z = z; |
761 |
storage[j].delta = delta; |
storage[j].delta = delta; |
762 |
DoubleVector v1(z); |
DoubleVector v(z); |
|
DoubleVector v2(z); |
|
|
k = param[i+j]; |
|
|
v1[k] += delta[k]; |
|
|
v2[k] -= delta[k]; |
|
763 |
|
|
764 |
#pragma omp parallel sections num_threads(2) |
if (j<paral_tokens) { |
765 |
{ |
k = param[i+j]; |
766 |
#pragma omp section |
v[k] += delta[k]; |
|
{ |
|
|
storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1); |
|
|
} |
|
|
#pragma omp section |
|
|
{ |
|
|
storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2); |
|
767 |
} |
} |
768 |
|
else { |
769 |
|
k = param[i+j-paral_tokens]; |
770 |
|
v[k] -= delta[k]; |
771 |
|
} |
772 |
|
|
773 |
|
storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v); |
774 |
|
storage[j].z[k] = v[k]; |
775 |
} |
} |
776 |
|
|
777 |
|
for (j = 0; j < paral_tokens; ++j) { |
778 |
|
k = param[i+j]; |
779 |
if (storage[j].ftmp < minf) { |
if (storage[j].ftmp < minf) { |
780 |
storage[j].iters = 1; |
storage[j].iters = 1; |
781 |
storage[j].z[k] = v1[k]; |
// storage[j].z[k] = v1[k]; |
782 |
} else { |
} else { |
783 |
storage[j].iters = 2; |
storage[j].iters = 2; |
784 |
storage[j].delta[k] = 0.0 - delta[k]; |
storage[j].delta[k] = 0.0 - delta[k]; |
785 |
if (storage[j+paral_tokens].ftmp < minf) { |
if (storage[j+paral_tokens].ftmp < minf) { |
786 |
storage[j].ftmp = storage[j+paral_tokens].ftmp; |
storage[j].ftmp = storage[j+paral_tokens].ftmp; |
787 |
storage[j].z[k] = v2[k]; |
storage[j].z[k] = storage[j+paral_tokens].z[k]; |
788 |
} |
} |
789 |
else iters += 2; |
else iters += 2; |
790 |
} |
} |
804 |
|
|
805 |
i += paral_tokens; |
i += paral_tokens; |
806 |
} |
} |
807 |
|
paral_tokens = numThr / 2; |
808 |
} |
} |
809 |
|
|
810 |
delete[] storage; |
delete[] storage; |
878 |
|
|
879 |
while (1) { |
while (1) { |
880 |
if (isZero(bestf)) { |
if (isZero(bestf)) { |
881 |
#ifdef NO_OPENMP |
#ifndef _OPENMP |
882 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
883 |
#endif |
#endif |
884 |
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"); |
903 |
/* find best new point, one coord at a time */ |
/* find best new point, one coord at a time */ |
904 |
for (i = 0; i < nvars; i++) |
for (i = 0; i < nvars; i++) |
905 |
trialx[i] = x[i]; |
trialx[i] = x[i]; |
906 |
#ifndef NO_OPENMP |
#ifdef _OPENMP |
907 |
newf = this->bestNearbyOMP2(delta, trialx, bestf, param); |
newf = this->bestNearbySpec(delta, trialx, bestf, param); |
908 |
if (newf == -1) { |
if (newf == -1) { |
909 |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
910 |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
915 |
#endif |
#endif |
916 |
/* if too many function evaluations occur, terminate the algorithm */ |
/* if too many function evaluations occur, terminate the algorithm */ |
917 |
|
|
918 |
#ifdef NO_OPENMP |
#ifndef _OPENMP |
919 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
920 |
#endif |
#endif |
921 |
if (iters > hookeiter) { |
if (iters > hookeiter) { |
985 |
/* only move forward if this is really an improvement */ |
/* only move forward if this is really an improvement */ |
986 |
oldf = newf; |
oldf = newf; |
987 |
newf = EcoSystem->SimulateAndUpdate(trialx); |
newf = EcoSystem->SimulateAndUpdate(trialx); |
988 |
#ifndef NO_OPENMP |
#ifdef _OPENMP |
989 |
iters++; |
iters++; |
990 |
#endif |
#endif |
991 |
if ((isEqual(newf, oldf)) || (newf > oldf)) { |
if ((isEqual(newf, oldf)) || (newf > oldf)) { |
998 |
for (i = 0; i < nvars; i++) |
for (i = 0; i < nvars; i++) |
999 |
x[i] = trialx[i]; |
x[i] = trialx[i]; |
1000 |
|
|
1001 |
#ifndef NO_OPENMP |
#ifdef _OPENMP |
1002 |
newf = this->bestNearbyOMP2(delta, trialx, bestf, param); |
newf = this->bestNearbySpec(delta, trialx, bestf, param); |
1003 |
if (newf == -1) { |
if (newf == -1) { |
1004 |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
1005 |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
handle.logMessage(LOGINFO, "\nThe number of threads must be a multiple of 2\n"); |
1012 |
break; |
break; |
1013 |
|
|
1014 |
/* if too many function evaluations occur, terminate the algorithm */ |
/* if too many function evaluations occur, terminate the algorithm */ |
1015 |
#ifdef NO_OPENMP |
#ifndef _OPENMP |
1016 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
1017 |
#endif |
#endif |
1018 |
if (iters > hookeiter) { |
if (iters > hookeiter) { |
1031 |
} |
} |
1032 |
} |
} |
1033 |
|
|
1034 |
#ifdef NO_OPENMP |
#ifndef _OPENMP |
1035 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
1036 |
#endif |
#endif |
1037 |
if (newf < bestf) { |
if (newf < bestf) { |
1066 |
delta[i] *= rho; |
delta[i] *= rho; |
1067 |
} |
} |
1068 |
} |
} |
1069 |
|
//#endif |
1070 |
#endif |
#endif |