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(); |
|
#ifdef _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)) { |
|
#ifndef _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]; |
|
#ifdef _OPENMP |
|
354 |
newf = this->bestNearbyRepro(delta, trialx, bestf, param); |
newf = this->bestNearbyRepro(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 |
|
|
|
#ifndef _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); |
|
#ifdef _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 |
|
|
|
#ifdef _OPENMP |
|
441 |
newf = this->bestNearbyRepro(delta, trialx, bestf, param); |
newf = this->bestNearbyRepro(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 */ |
|
#ifndef _OPENMP |
|
451 |
iters = EcoSystem->getFuncEval() - offset; |
iters = EcoSystem->getFuncEval() - offset; |
|
#endif |
|
452 |
if (iters > hookeiter) { |
if (iters > hookeiter) { |
453 |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
454 |
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
465 |
} |
} |
466 |
} // while (newf < bestf) |
} // while (newf < bestf) |
467 |
|
|
|
#ifndef _OPENMP |
|
468 |
iters = EcoSystem->getFuncEval() - offset; |
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 |
|
handle.logMessage(LOGINFO, "Error in Hooke & Jeeves optimisation after", iters, "function evaluations, f(x) = 0"); |
566 |
|
converge = -1; |
567 |
|
return; |
568 |
|
} |
569 |
|
|
570 |
|
/* randomize the order of the parameters once in a while */ |
571 |
|
rchange = 0; |
572 |
|
while (rchange < nvars) { |
573 |
|
rnumber = rand() % nvars; |
574 |
|
rcheck = 1; |
575 |
|
for (i = 0; i < rchange; i++) |
576 |
|
if (param[i] == rnumber) |
577 |
|
rcheck = 0; |
578 |
|
if (rcheck) { |
579 |
|
param[rchange] = rnumber; |
580 |
|
rchange++; |
581 |
|
} |
582 |
|
} |
583 |
|
|
584 |
|
/* find best new point, one coord at a time */ |
585 |
|
for (i = 0; i < nvars; i++) |
586 |
|
trialx[i] = x[i]; |
587 |
|
newf = this->bestNearby(delta, trialx, bestf, param); |
588 |
|
/* if too many function evaluations occur, terminate the algorithm */ |
589 |
|
|
590 |
|
if (iters > hookeiter) { |
591 |
|
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
592 |
|
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
593 |
|
handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); |
594 |
|
handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations"); |
595 |
|
handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); |
596 |
|
|
597 |
|
score = EcoSystem->SimulateAndUpdate(trialx); |
598 |
|
handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); |
599 |
|
for (i = 0; i < nvars; i++) |
600 |
|
bestx[i] = trialx[i] * init[i]; |
601 |
|
EcoSystem->storeVariables(score, bestx); |
602 |
|
return; |
603 |
|
} |
604 |
|
|
605 |
|
/* if we made some improvements, pursue that direction */ |
606 |
|
while (newf < bestf) { |
607 |
|
for (i = 0; i < nvars; i++) { |
608 |
|
/* if it has been trapped but f has now gotten better (bndcheck) */ |
609 |
|
/* we assume that we are out of the trap, reset the counters */ |
610 |
|
/* and go back to the stepsize we had when we got trapped */ |
611 |
|
if ((trapped[i]) && (newf < oldf * bndcheck)) { |
612 |
|
trapped[i] = 0; |
613 |
|
lbound[i] = 0; |
614 |
|
rbounds[i] = 0; |
615 |
|
delta[i] = initialstep[i]; |
616 |
|
|
617 |
|
} else if (trialx[i] < (lowerb[i] + verysmall)) { |
618 |
|
lbound[i]++; |
619 |
|
trialx[i] = lowerb[i]; |
620 |
|
if (!trapped[i]) { |
621 |
|
initialstep[i] = delta[i]; |
622 |
|
trapped[i] = 1; |
623 |
|
} |
624 |
|
/* if it has hit the bounds 2 times then increase the stepsize */ |
625 |
|
if (lbound[i] >= 2) |
626 |
|
delta[i] /= rho; |
627 |
|
|
628 |
|
} else if (trialx[i] > (upperb[i] - verysmall)) { |
629 |
|
rbounds[i]++; |
630 |
|
trialx[i] = upperb[i]; |
631 |
|
if (!trapped[i]) { |
632 |
|
initialstep[i] = delta[i]; |
633 |
|
trapped[i] = 1; |
634 |
|
} |
635 |
|
/* if it has hit the bounds 2 times then increase the stepsize */ |
636 |
|
if (rbounds[i] >= 2) |
637 |
|
delta[i] /= rho; |
638 |
|
} |
639 |
|
} |
640 |
|
|
641 |
|
for (i = 0; i < nvars; i++) { |
642 |
|
/* firstly, arrange the sign of delta[] */ |
643 |
|
if (trialx[i] < x[i]) |
644 |
|
delta[i] = 0.0 - fabs(delta[i]); |
645 |
|
else |
646 |
|
delta[i] = fabs(delta[i]); |
647 |
|
|
648 |
|
/* now, move further in this direction */ |
649 |
|
tmp = x[i]; |
650 |
|
x[i] = trialx[i]; |
651 |
|
trialx[i] = trialx[i] + trialx[i] - tmp; |
652 |
|
} |
653 |
|
|
654 |
|
/* only move forward if this is really an improvement */ |
655 |
|
oldf = newf; |
656 |
|
newf = EcoSystem->SimulateAndUpdate(trialx); |
657 |
|
if ((isEqual(newf, oldf)) || (newf > oldf)) { |
658 |
|
newf = oldf; //JMB no improvement, so reset the value of newf |
659 |
|
break; |
660 |
|
} |
661 |
|
|
662 |
|
/* OK, it's better, so update variables and look around */ |
663 |
|
bestf = newf; |
664 |
|
for (i = 0; i < nvars; i++) |
665 |
|
x[i] = trialx[i]; |
666 |
|
|
667 |
|
if (isEqual(newf, bestf)) |
668 |
|
break; |
669 |
|
|
670 |
|
/* if too many function evaluations occur, terminate the algorithm */ |
671 |
|
if (iters > hookeiter) { |
672 |
|
handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n"); |
673 |
|
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
674 |
|
handle.logMessage(LOGINFO, "The steplength was reduced to", steplength); |
675 |
|
handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of function evaluations"); |
676 |
|
handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); |
677 |
|
|
678 |
|
score = EcoSystem->SimulateAndUpdate(trialx); |
679 |
|
handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score); |
680 |
|
for (i = 0; i < nvars; i++) |
681 |
|
bestx[i] = trialx[i] * init[i]; |
682 |
|
EcoSystem->storeVariables(score, bestx); |
683 |
|
return; |
684 |
|
} |
685 |
|
} // while (newf < bestf) |
686 |
|
|
687 |
if (newf < bestf) { |
if (newf < bestf) { |
688 |
for (i = 0; i < nvars; i++) |
for (i = 0; i < nvars; i++) |
689 |
bestx[i] = x[i] * init[i]; |
bestx[i] = x[i] * init[i]; |
718 |
} |
} |
719 |
|
|
720 |
/* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/ |
/* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/ |
721 |
#ifdef SPECULATIVE |
#ifdef _OPENMP |
722 |
|
//#ifdef SPECULATIVE |
723 |
double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { |
double OptInfoHooke::bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) { |
724 |
double minf; |
double minf; |
725 |
int i, j, k, ii; |
int i, j, k, ii; |
1062 |
delta[i] *= rho; |
delta[i] *= rho; |
1063 |
} |
} |
1064 |
} |
} |
1065 |
|
//#endif |
1066 |
#endif |
#endif |