186 |
#ifndef NO_OPENMP |
#ifndef NO_OPENMP |
187 |
extern Ecosystem** EcoSystems; |
extern Ecosystem** EcoSystems; |
188 |
#endif |
#endif |
|
//extern StochasticData* data; |
|
|
|
|
189 |
|
|
190 |
|
/*sequential code replaced at seq_optimize_template.h*/ |
191 |
//void OptInfoSimann::OptimiseLikelihood() { |
//void OptInfoSimann::OptimiseLikelihood() { |
192 |
// |
// |
193 |
// //set initial values |
// //set initial values |
506 |
handle.logMessage(LOGINFO, "Error starting Simulated Annealing optimisation with f(x) = infinity"); |
handle.logMessage(LOGINFO, "Error starting Simulated Annealing optimisation with f(x) = infinity"); |
507 |
converge = -1; |
converge = -1; |
508 |
iters = 1; |
iters = 1; |
|
// EcoSystem->writeOptValuesOMP(); |
|
509 |
return; |
return; |
510 |
} |
} |
511 |
|
|
650 |
x[i] = storage[bestId].trialx[i]; |
x[i] = storage[bestId].trialx[i]; |
651 |
funcval = trialf; |
funcval = trialf; |
652 |
nacc++; |
nacc++; |
|
// vns[param[l+bestId]]++; |
|
|
//nacp[param[l+bestId]]+=numThr; //JMB - not sure about this ... |
|
|
// for (i = 0; i < numThr; ++i) |
|
|
// nacp[param[l+i]]++; |
|
|
// nacp[param[l+bestId]]++; |
|
|
|
|
653 |
} else { |
} else { |
654 |
//Accept according to metropolis condition |
//Accept according to metropolis condition |
655 |
p = expRep((trialf - funcval) / t); |
p = expRep((trialf - funcval) / t); |
661 |
funcval = trialf; |
funcval = trialf; |
662 |
naccmet++; |
naccmet++; |
663 |
nacp[param[(l+bestId)%nvars]]++; |
nacp[param[(l+bestId)%nvars]]++; |
|
// vns[param[l+bestId]]++; |
|
|
// nacp[param[l+bestId]]+=numThr; |
|
|
// for (i = 0; i < numThr; ++i) |
|
|
// nacp[param[l+i]]++; |
|
|
// nacp[param[l+bestId]]++; |
|
664 |
} else { |
} else { |
665 |
//Reject point |
//Reject point |
666 |
nrej++; |
nrej++; |
670 |
if (isZero(trialf)) { |
if (isZero(trialf)) { |
671 |
handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0"); |
handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0"); |
672 |
converge = -1; |
converge = -1; |
|
// EcoSystem->writeOptValuesOMP(); |
|
673 |
delete[] storage; |
delete[] storage; |
674 |
return; |
return; |
675 |
} |
} |
683 |
if (scale) { |
if (scale) { |
684 |
for (i = 0; i < nvars; ++i) |
for (i = 0; i < nvars; ++i) |
685 |
scalex[i] = bestx[i] * init[i]; |
scalex[i] = bestx[i] * init[i]; |
|
//FIXME EcoSystems[]-> |
|
686 |
EcoSystem->storeVariables(-fopt, scalex); |
EcoSystem->storeVariables(-fopt, scalex); |
687 |
} else |
} else |
688 |
EcoSystem->storeVariables(-fopt, bestx); |
EcoSystem->storeVariables(-fopt, bestx); |
693 |
} |
} |
694 |
} |
} |
695 |
} |
} |
|
|
|
|
//Adjust vm so that approximately half of all evaluations are accepted |
|
|
// for (i = 0; i < nvars; ++i) { |
|
|
//// cout << "nacp["<<i<<"]:"<<nacp[i]<<endl; |
|
|
// ratio = nsdiv * nacp[i]; |
|
|
// nacp[i] = 0; |
|
|
// if (ratio > uratio) { |
|
|
// vm[i] = vm[i] * (1.0 + cs * (ratio - uratio)); |
|
|
// } else if (ratio < lratio) { |
|
|
// vm[i] = vm[i] / (1.0 + cs * (lratio - ratio)); |
|
|
// } |
|
|
// |
|
|
// if (vm[i] < rathersmall) |
|
|
// vm[i] = rathersmall; |
|
|
// if (vm[i] > (upperb[i] - lowerb[i])) |
|
|
// vm[i] = upperb[i] - lowerb[i]; |
|
|
// |
|
|
//// cout << "vm["<<i<<"]:"<<vm[i]<<endl; |
|
|
// |
|
|
// } |
|
696 |
} |
} |
697 |
|
|
698 |
//Check termination criteria |
//Check termination criteria |
723 |
converge = 1; |
converge = 1; |
724 |
score = EcoSystem->SimulateAndUpdate(bestx); |
score = EcoSystem->SimulateAndUpdate(bestx); |
725 |
handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score); |
handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score); |
|
// EcoSystem->writeOptValuesOMP(); |
|
726 |
delete[] storage; |
delete[] storage; |
727 |
return; |
return; |
728 |
} |
} |
737 |
for (i = 0; i < nvars; ++i) |
for (i = 0; i < nvars; ++i) |
738 |
x[i] = bestx[i]; |
x[i] = bestx[i]; |
739 |
} |
} |
|
// EcoSystem->writeOptValuesOMP(); |
|
740 |
} |
} |
741 |
#endif |
#endif |
742 |
|
|
775 |
DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm, unsigned* seed) |
DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm, unsigned* seed) |
776 |
{ |
{ |
777 |
int i, k; |
int i, k; |
|
// double old = trialx[param[l]]; |
|
|
// cout << "[" << endl; |
|
778 |
for (i = 0; i < nvars; ++i) { |
for (i = 0; i < nvars; ++i) { |
779 |
if (i == param[l]) { |
if (i == param[l]) { |
780 |
trialx[i] = x[i] + ((randomNumber(&*seed) * 2.0) - 1.0) * vm[i]; |
trialx[i] = x[i] + ((randomNumber(&*seed) * 2.0) - 1.0) * vm[i]; |
792 |
} |
} |
793 |
} else |
} else |
794 |
trialx[i] = x[i]; |
trialx[i] = x[i]; |
|
// cout << trialx[i]<<" "; |
|
795 |
} |
} |
|
|
|
|
// cout << endl <<l<< "-old[" <<param[l] << "]:" << old <<"-" << trialx[param[l]]<< " - " << vm[param[l]]<<endl; |
|
796 |
} |
} |
797 |
|
|
798 |
void buildNewParams_f(Siman& seed, DoubleVector& params) { |
void buildNewParams_f(Siman& seed, DoubleVector& params) { |
808 |
#ifndef NO_OPENMP |
#ifndef NO_OPENMP |
809 |
int id = omp_get_thread_num (); |
int id = omp_get_thread_num (); |
810 |
trialf = EcoSystems[id]->SimulateAndUpdate(params); |
trialf = EcoSystems[id]->SimulateAndUpdate(params); |
|
// cout << id << "-"; |
|
811 |
#else |
#else |
812 |
trialf = EcoSystem->SimulateAndUpdate(params); |
trialf = EcoSystem->SimulateAndUpdate(params); |
813 |
#endif |
#endif |
|
// cout << "trailf:" << trialf<< endl; |
|
814 |
return -trialf; |
return -trialf; |
815 |
} |
} |
816 |
|
|
828 |
double uratio = seed.getUratio(); |
double uratio = seed.getUratio(); |
829 |
double lratio = seed.getLratio(); |
double lratio = seed.getLratio(); |
830 |
|
|
|
// cout << "uratio:" << uratio << endl; |
|
|
// cout << "lratio:" << lratio << endl; |
|
|
|
|
|
//cout << "adjust vm" << endl; |
|
831 |
for (i = 0; i < seed.getNvars(); i++) { |
for (i = 0; i < seed.getNvars(); i++) { |
|
// cout << "nacp[" << i << "]:" << seed.getNacp(i) << endl; |
|
832 |
ratio = nsdiv * seed.getNacp(i); |
ratio = nsdiv * seed.getNacp(i); |
833 |
seed.setNacp(i,0); |
seed.setNacp(i,0); |
834 |
if (ratio > uratio) { |
if (ratio > uratio) { |
857 |
|
|
858 |
DoubleVector* bestx = seed.getBestx(); |
DoubleVector* bestx = seed.getBestx(); |
859 |
|
|
|
//funcval = fopt; |
|
860 |
for (i = 0; i < seed.getNvars(); i++) |
for (i = 0; i < seed.getNvars(); i++) |
861 |
x[i] = (*bestx)[i]; |
x[i] = (*bestx)[i]; |
862 |
} |
} |
895 |
//Accept the new point if the new function value better |
//Accept the new point if the new function value better |
896 |
bool ret = true; |
bool ret = true; |
897 |
int i; |
int i; |
|
// cout << "mustAccept:" << trialf << "|" << funcval<< endl; |
|
898 |
int aux = siman.getParam()[siman.getL()]; |
int aux = siman.getParam()[siman.getL()]; |
899 |
if ((trialf - funcval) > verysmall) { |
if ((trialf - funcval) > verysmall) { |
900 |
siman.incrementNacp(aux); |
siman.incrementNacp(aux); |
904 |
double p, pp; |
double p, pp; |
905 |
//Accept according to metropolis condition |
//Accept according to metropolis condition |
906 |
|
|
|
// cout << "t:" << siman.getT() << endl; |
|
|
|
|
907 |
p = expRep((trialf - funcval) / siman.getT()); |
p = expRep((trialf - funcval) / siman.getT()); |
908 |
pp = randomNumber(siman.getSeedM()); |
pp = randomNumber(siman.getSeedM()); |
909 |
if (pp < p) { |
if (pp < p) { |
|
// cout << pp << "<" << p << endl; |
|
910 |
siman.incrementNacp(aux); |
siman.incrementNacp(aux); |
911 |
siman.incrementNaccmet(); |
siman.incrementNaccmet(); |
912 |
// //Accept point |
// //Accept point |
|
// for (i = 0; i < nvars; i++) |
|
|
// x[i] = trialx[i]; |
|
|
// funcval = trialf; |
|
|
// naccmet++; |
|
|
// nacp[param[l]]++; |
|
913 |
} else { |
} else { |
914 |
//Reject point |
//Reject point |
915 |
ret = false; |
ret = false; |
916 |
siman.incrementNrej(); |
siman.incrementNrej(); |
917 |
} |
} |
918 |
} |
} |
|
// cout << "nacp[" << aux << "]:" << siman.getNacp(aux) << endl; |
|
919 |
// JMB added check for really silly values |
// JMB added check for really silly values |
920 |
if (isZero(trialf)) { |
if (isZero(trialf)) { |
921 |
handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", |
handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", |
922 |
iters, "function evaluations, f(x) = 0"); |
iters, "function evaluations, f(x) = 0"); |
|
|
|
923 |
siman.setConverge(-1); |
siman.setConverge(-1); |
924 |
return false; |
return false; |
925 |
} |
} |
928 |
if (siman.getL() == 0) |
if (siman.getL() == 0) |
929 |
siman.incrementNS(); |
siman.incrementNS(); |
930 |
if (siman.getNS() >= siman.getNs()){ |
if (siman.getNS() >= siman.getNs()){ |
|
|
|
931 |
siman.setNS(0); |
siman.setNS(0); |
|
|
|
932 |
siman.incrementNT(); |
siman.incrementNT(); |
933 |
adjustVm(siman); |
adjustVm(siman); |
934 |
siman.randomizeParams(); |
siman.randomizeParams(); |
|
// if (seed.getNT() >= seed.getNt()){ |
|
|
// |
|
|
// } |
|
|
|
|
935 |
} |
} |
936 |
|
|
937 |
return ret; |
return ret; |
948 |
bool quit = false; |
bool quit = false; |
949 |
if (siman.getNT() >= siman.getNt()) |
if (siman.getNT() >= siman.getNt()) |
950 |
{ |
{ |
|
// cout << "nt:" << siman.getNT(); |
|
|
//Check termination criteria |
|
951 |
int i; |
int i; |
952 |
DoubleVector fstar = siman.getFstar(); |
DoubleVector fstar = siman.getFstar(); |
953 |
|
|
998 |
handle.logMessage(LOGINFO, "Number of rejected points", siman.getNrej()); |
handle.logMessage(LOGINFO, "Number of rejected points", siman.getNrej()); |
999 |
*score = EcoSystem->SimulateAndUpdate(*bestX); |
*score = EcoSystem->SimulateAndUpdate(*bestX); |
1000 |
handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", *score); |
handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", *score); |
|
// EcoSystem->writeOptValues(); |
|
1001 |
} |
} |
1002 |
}; |
}; |
1003 |
|
|
1004 |
// Required |
// Required |
1005 |
std::ostream &operator<<(std::ostream &os, const DoubleVector &p) |
std::ostream &operator<<(std::ostream &os, const DoubleVector &p) |
1006 |
{ |
{ |
|
// os << "[ AUCH"; |
|
|
// for (int i = 0; i< NPARAMS; ++i) { |
|
|
// std::cout << p[i] << ' '; |
|
|
// } |
|
|
// std::cout << "]"; |
|
1007 |
os << ""; |
os << ""; |
1008 |
return os; |
return os; |
1009 |
} |
} |
1014 |
//set initial values |
//set initial values |
1015 |
|
|
1016 |
double tmp, p, pp; |
double tmp, p, pp; |
1017 |
double fopt, funcval, trialf; |
double funcval, trialf; |
1018 |
int a, i, j, k, l, quit; |
int a, i, j, k, l, quit; |
1019 |
int rchange, rcheck, rnumber; //Used to randomise the order of the parameters |
int rchange, rcheck, rnumber; //Used to randomise the order of the parameters |
1020 |
|
|
1021 |
handle.logMessage(LOGINFO, |
handle.logMessage(LOGINFO, |
1022 |
"\nStarting Simulated Annealing optimisation algorithm----\n"); |
"\nStarting Simulated Annealing optimisation algorithm\n"); |
1023 |
int nvars = EcoSystem->numOptVariables(); |
int nvars = EcoSystem->numOptVariables(); |
1024 |
DoubleVector x(nvars); |
DoubleVector x(nvars); |
1025 |
DoubleVector init(nvars); |
DoubleVector init(nvars); |
1077 |
|
|
1078 |
//the function is to be minimised so switch the sign of funcval (and trialf) |
//the function is to be minimised so switch the sign of funcval (and trialf) |
1079 |
funcval = -funcval; |
funcval = -funcval; |
|
// offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop |
|
|
// nacc++; |
|
1080 |
cs /= lratio; //JMB save processing time |
cs /= lratio; //JMB save processing time |
|
// nsdiv = 1.0 / ns; |
|
|
fopt = funcval; |
|
1081 |
for (i = 0; i < tempcheck; i++) |
for (i = 0; i < tempcheck; i++) |
1082 |
fstar[i] = funcval; |
fstar[i] = funcval; |
1083 |
|
|
|
//unsigned seed = 1234; |
|
|
|
|
1084 |
Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns), |
Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns), |
1085 |
tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score); |
tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score); |
1086 |
|
|
1087 |
ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_f, buildNewParams_f> |
ReproducibleSearch<Siman, DoubleVector, ControlClass, evaluate_f, buildNewParams_f> |
1088 |
pa(s, x, simanniter); |
pa(s, x, simanniter); |
1089 |
|
|
|
// cout << "Sequential run" << endl; |
|
|
// pa.seq_opt(); |
|
|
|
|
1090 |
#ifndef NO_OPENMP |
#ifndef NO_OPENMP |
1091 |
|
// OpenMP parallelization |
1092 |
int numThr = omp_get_max_threads ( ); |
int numThr = omp_get_max_threads ( ); |
1093 |
pa.paral_opt_omp(numThr,numThr); |
pa.paral_opt_omp(funcval,numThr,numThr); |
1094 |
#else |
#else |
1095 |
pa.seq_opt(); |
// sequential code |
1096 |
|
pa.seq_opt(funcval); |
1097 |
#endif |
#endif |
1098 |
|
|
1099 |
} |
} |