--- trunk/gadget/simann.cc 2015/07/23 19:00:38 11 +++ trunk/gadget/simann.cc 2015/07/24 18:36:24 12 @@ -186,9 +186,8 @@ #ifndef NO_OPENMP extern Ecosystem** EcoSystems; #endif -//extern StochasticData* data; - +/*sequential code replaced at seq_optimize_template.h*/ //void OptInfoSimann::OptimiseLikelihood() { // // //set initial values @@ -507,7 +506,6 @@ handle.logMessage(LOGINFO, "Error starting Simulated Annealing optimisation with f(x) = infinity"); converge = -1; iters = 1; -// EcoSystem->writeOptValuesOMP(); return; } @@ -652,12 +650,6 @@ x[i] = storage[bestId].trialx[i]; funcval = trialf; 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]]++; - } else { //Accept according to metropolis condition p = expRep((trialf - funcval) / t); @@ -669,11 +661,6 @@ funcval = trialf; naccmet++; 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]]++; } else { //Reject point nrej++; @@ -683,7 +670,6 @@ if (isZero(trialf)) { handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0"); converge = -1; -// EcoSystem->writeOptValuesOMP(); delete[] storage; return; } @@ -697,7 +683,6 @@ if (scale) { for (i = 0; i < nvars; ++i) scalex[i] = bestx[i] * init[i]; - //FIXME EcoSystems[]-> EcoSystem->storeVariables(-fopt, scalex); } else EcoSystem->storeVariables(-fopt, bestx); @@ -708,26 +693,6 @@ } } } - - //Adjust vm so that approximately half of all evaluations are accepted -// for (i = 0; i < nvars; ++i) { -//// cout << "nacp["< 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["<SimulateAndUpdate(bestx); handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", score); -// EcoSystem->writeOptValuesOMP(); delete[] storage; return; } @@ -773,7 +737,6 @@ for (i = 0; i < nvars; ++i) x[i] = bestx[i]; } -// EcoSystem->writeOptValuesOMP(); } #endif @@ -812,8 +775,6 @@ DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm, unsigned* seed) { int i, k; -// double old = trialx[param[l]]; -// cout << "[" << endl; for (i = 0; i < nvars; ++i) { if (i == param[l]) { trialx[i] = x[i] + ((randomNumber(&*seed) * 2.0) - 1.0) * vm[i]; @@ -831,10 +792,7 @@ } } else trialx[i] = x[i]; -// cout << trialx[i]<<" "; } - -// cout << endl <SimulateAndUpdate(params); -// cout << id << "-"; #else trialf = EcoSystem->SimulateAndUpdate(params); #endif -// cout << "trailf:" << trialf<< endl; return -trialf; } @@ -872,12 +828,7 @@ double uratio = seed.getUratio(); double lratio = seed.getLratio(); - // cout << "uratio:" << uratio << endl; - // cout << "lratio:" << lratio << endl; - - //cout << "adjust vm" << endl; for (i = 0; i < seed.getNvars(); i++) { -// cout << "nacp[" << i << "]:" << seed.getNacp(i) << endl; ratio = nsdiv * seed.getNacp(i); seed.setNacp(i,0); if (ratio > uratio) { @@ -906,7 +857,6 @@ DoubleVector* bestx = seed.getBestx(); - //funcval = fopt; for (i = 0; i < seed.getNvars(); i++) x[i] = (*bestx)[i]; } @@ -945,7 +895,6 @@ //Accept the new point if the new function value better bool ret = true; int i; -// cout << "mustAccept:" << trialf << "|" << funcval<< endl; int aux = siman.getParam()[siman.getL()]; if ((trialf - funcval) > verysmall) { siman.incrementNacp(aux); @@ -955,32 +904,22 @@ double p, pp; //Accept according to metropolis condition -// cout << "t:" << siman.getT() << endl; - p = expRep((trialf - funcval) / siman.getT()); pp = randomNumber(siman.getSeedM()); if (pp < p) { -// cout << pp << "<" << p << endl; siman.incrementNacp(aux); siman.incrementNaccmet(); // //Accept point -// for (i = 0; i < nvars; i++) -// x[i] = trialx[i]; -// funcval = trialf; -// naccmet++; -// nacp[param[l]]++; } else { //Reject point ret = false; siman.incrementNrej(); } } -// cout << "nacp[" << aux << "]:" << siman.getNacp(aux) << endl; // JMB added check for really silly values if (isZero(trialf)) { handle.logMessage(LOGINFO, "Error in Simulated Annealing optimisation after", iters, "function evaluations, f(x) = 0"); - siman.setConverge(-1); return false; } @@ -989,16 +928,10 @@ if (siman.getL() == 0) siman.incrementNS(); if (siman.getNS() >= siman.getNs()){ - siman.setNS(0); - siman.incrementNT(); adjustVm(siman); siman.randomizeParams(); -// if (seed.getNT() >= seed.getNt()){ -// -// } - } return ret; @@ -1015,8 +948,6 @@ bool quit = false; if (siman.getNT() >= siman.getNt()) { -// cout << "nt:" << siman.getNT(); - //Check termination criteria int i; DoubleVector fstar = siman.getFstar(); @@ -1067,18 +998,12 @@ handle.logMessage(LOGINFO, "Number of rejected points", siman.getNrej()); *score = EcoSystem->SimulateAndUpdate(*bestX); handle.logMessage(LOGINFO, "\nSimulated Annealing finished with a likelihood score of", *score); -// EcoSystem->writeOptValues(); } }; // Required std::ostream &operator<<(std::ostream &os, const DoubleVector &p) { -// os << "[ AUCH"; -// for (int i = 0; i< NPARAMS; ++i) { -// std::cout << p[i] << ' '; -// } -// std::cout << "]"; os << ""; return os; } @@ -1089,12 +1014,12 @@ //set initial values double tmp, p, pp; - double fopt, funcval, trialf; + double funcval, trialf; int a, i, j, k, l, quit; int rchange, rcheck, rnumber; //Used to randomise the order of the parameters handle.logMessage(LOGINFO, - "\nStarting Simulated Annealing optimisation algorithm----\n"); + "\nStarting Simulated Annealing optimisation algorithm\n"); int nvars = EcoSystem->numOptVariables(); DoubleVector x(nvars); DoubleVector init(nvars); @@ -1152,30 +1077,23 @@ //the function is to be minimised so switch the sign of funcval (and trialf) funcval = -funcval; -// offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop -// nacc++; cs /= lratio; //JMB save processing time -// nsdiv = 1.0 / ns; - fopt = funcval; for (i = 0; i < tempcheck; i++) fstar[i] = funcval; - //unsigned seed = 1234; - Siman s(seed, seedM, seedP, nvars, nt, ns, param, &x, &lowerb, &upperb, vm, t, rt, (1.0 / ns), tempcheck, simanneps, fstar, lratio, uratio, cs, &bestx, scale, &converge, &score); ReproducibleSearch pa(s, x, simanniter); -// cout << "Sequential run" << endl; -// pa.seq_opt(); - #ifndef NO_OPENMP + // OpenMP parallelization int numThr = omp_get_max_threads ( ); - pa.paral_opt_omp(numThr,numThr); + pa.paral_opt_omp(funcval,numThr,numThr); #else - pa.seq_opt(); + // sequential code + pa.seq_opt(funcval); #endif }