From 7ae7468428248d0db936ebaf5b8fb9c7fd9fc1a5 Mon Sep 17 00:00:00 2001 From: LPTK Date: Fri, 19 Jul 2013 13:38:20 +0200 Subject: [PATCH] Cleaning, doc & fixes for trikiSA --- .../coolingSchedule/moTrikiCoolingSchedule.h | 751 ++++++++++-------- 1 file changed, 419 insertions(+), 332 deletions(-) diff --git a/mo/src/coolingSchedule/moTrikiCoolingSchedule.h b/mo/src/coolingSchedule/moTrikiCoolingSchedule.h index 259f1747c..bbb60ee2b 100644 --- a/mo/src/coolingSchedule/moTrikiCoolingSchedule.h +++ b/mo/src/coolingSchedule/moTrikiCoolingSchedule.h @@ -1,363 +1,450 @@ /* - - Copyright (C) DOLPHIN Project-Team, INRIA Futurs, 2006-2008 - (C) OPAC Team, LIFL, 2002-2008 - Sébastien Cahon, Jean-Charles Boisson (Jean-Charles.Boisson@lifl.fr) +(c) Thales group, 2010 - This software is governed by the CeCILL license under French law and - abiding by the rules of distribution of free software. You can use, - modify and/ or redistribute the software under the terms of the CeCILL - license as circulated by CEA, CNRS and INRIA at the following URL - "http://www.cecill.info". + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; + version 2 of the License. - As a counterpart to the access to the source code and rights to copy, - modify and redistribute granted by the license, users are provided only - with a limited warranty and the software's author, the holder of the - economic rights, and the successive licensors have only limited liability. + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. - In this respect, the user's attention is drawn to the risks associated - with loading, using, modifying and/or developing or reproducing the - software by the user in light of its specific status of free software, - that may mean that it is complicated to manipulate, and that also - therefore means that it is reserved for developers and experienced - professionals having in-depth computer knowledge. Users are therefore - encouraged to load and test the software's suitability as regards their - requirements in conditions enabling the security of their systems and/or - data to be ensured and, more generally, to use and operate it in the - same conditions as regards security. - The fact that you are presently reading this means that you have had - knowledge of the CeCILL license and that you accept its terms. + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - ParadisEO WebSite : http://paradiseo.gforge.inria.fr - Contact: paradiseo-help@lists.gforge.inria.fr - */ +Contact: http://eodev.sourceforge.net + +Authors: +Lionel Parreaux + +*/ #ifndef _moTrikiCoolingSchedule_h #define _moTrikiCoolingSchedule_h #include +#include + #include #include -#include #include -#include +#include - -/* -#include - -#include -#include -#include -#include +/** + * Cooling Schedule, adapted from E.Triki, Y.Collette, P.Siarry (2004) + * This CS is based on an initial estimation of the standard deviation + * and an expected decrease in cost between each Markov chain, + * possibly re-estimating the prameters + * + * + * A detailed explanation follows: + * + * + * Initialization + * + * A random walk of n steps should be performed to estimate the std dev + * of the fitness. The init temp is set to this value. + * + * + * Algorithm + * + * The CS is based on Markov chains, during which the temp is constant. + * A Markov chain ends when a given number of solutions 'max_accepted' + * has been reached, or when a given number 'max_generated' of solutions + * have been generated. + * + * After each chain, the average cost of the solutions of the chain is + * expected to have decreased of delta (compared to the previous chain) + * with delta initialized to = stddev/mu2. + * If it's the case (ie: avgCost/(prevAvgCost-delta) < xi, where xi == 1+epsilon) + * then we say we're at equilibrium and we apply a normal temperature + * decrease, of ratio + * alpha = 1-_temp*delta/variance + * + * Negative temperatures (when alpha < 0) happen when the initial std + * dev was misestimated or, according to the article, when "the minimum + * cost is greater than the current average cost minus delta"(sic)(??). + * In case of a neg temp, we increment a counter 'negative_temp' and we + * reinitialize the algorithm until the temperature is no more negative, + * unless the counter of neg temp reaches a constant K2, in which case + * the behavior has been chosen to be this of a greedy algorithm (SA + * with nil temperature). + * + * Reinitializing the algo consists in: + * - setting the temperature "decrease" factor alpha to a + * constant lambda1 > 1 in order to in fact increase it + * - setting delta to sigma/mu1 (sigma being the std dev of the + * current Markov chain) + * + * Note that when not reinitializing, the expected decrease in cost + * 'delta' is never supposed to change. + * + * If the eq is not reached after the current chain, we increment a + * counter 'equilibrium_not_reached', and when it reaches K1 we + * reinitialize the algorithm and reset the counter to 0 (resetting + * to 0 was an added behavior and not part of the article; but without + * it the algo got trapped). + * + * + * Termination + * + * Currently, the algo terminates when the current average cost stops + * changing for 'theta' chains, or when the current std dev becomes + * null (added behavior; indeed, when the std dev is null it is no more + * possible to compute alpha), or when there is no accepted solution + * for the current "chain" (added, cf in this case we can't compute a + * std dev or an average). + * In practice, the algorithm never seems to terminate by "freezing" + * (first case), obviously because we need an implementation of + * approximate double comparison instead of exact comparison. + * + * */ - -#include -using namespace std; - - -//! -/*! - */ -template< class EOT, class Neighbor > //, class Neighborhood > +template< class EOT > class moTrikiCoolingSchedule: public moCoolingSchedule< EOT > { public: - //typedef typename Neighbor::EOT EOT ; - typedef moNeighborhood Neighborhood ; + + /** + * Constructor for the cooling schedule + * @param _initTemp the temperature at which the CS begins; a recommended value is _stdDevEstimation + * @param _stdDevEstimation an estimation of the standard deviation of the fitness. Typically, a random walk of n steps is performed to estimate the std dev of the fitness + * @param _max_accepted maximum number of solutions to accept before ending the Markov chain and reducing the temperature; depends on the pb/neighborhood + * @param _max_generated maximum number of solutions to generate before ending the Markov chain and reducing the temperature; depends on the pb/neighborhood + * @param _mu2 target decrease in cost factor, mu2 typically belongs to [1; 20] + * @param _mu1 target decrease in cost factor when reinitializing, in [2; 20] + * @param _lambda1 the increase in temperature (reheating factor), typically in [1.5; 4] + * @param _lambda2 lambda2 in [0.5; 0.99] + * @param _xi typically belongs to [1; 1.1] + * @param _theta typically set to 10 + * @param _K1 in [1; 4], the number of chains without reaching equilibrium before we raise the temperature + * @param _K2 maximul number of consecutive negative temperatures before switching to a greedy algorithm + */ + moTrikiCoolingSchedule ( + double _initTemp, + double _stdDevEstimation, + int _max_accepted = 50, + int _max_generated = 100, + double _mu2 = 2.5, + double _mu1 = 10, + double _lambda1 = 2, + double _lambda2 = .7, + double _xi = 1.05, + int _theta = 10, + int _K1 = 10, + int _K2 = 5 + ) + : initTemp(_initTemp), + initStdDev(_stdDevEstimation), + mu2(_mu2), + K1(_K1), + K2(_K2), + lambda1(_lambda1), + lambda2(_lambda2), + mu1(_mu1), + xi(_xi), + max_accepted(_max_accepted), + max_generated(_max_generated), + theta(_theta), + statIsInitialized(false) + { + chainStat.temperature = initTemp; + } + + /** + * Initialization + * @param _solution initial solution + */ + double init(EOT & _solution) { - //! Constructor - /*! - */ + chainStat.temperature = initTemp; + + accepted = generated = 0; + + negative_temp = equilibrium_not_reached = frozen = 0; + + chainStat.delta = initStdDev/mu2; + + reinitializing = false; + + return initTemp; + } - moTrikiCoolingSchedule (Neighborhood& _neighborhood, moEval& _eval, double _initTemp) - : initTemp(_initTemp), - mu2(10), // mu2 typically belongs to [1; 20] - K1(2), // K1 in [1; 4], the number of chains without reaching equilibrium before we raise the temperature - K2(5), // ??? - lambda1(2), // the increase in temperature, typically in [1.5; 4] - lambda2(.7), // lambda2 in [0.5; 0.99] - mu1(10), // target decrease in cost factor, in [2; 20] - xi(1.05), // xi typically belongs to [1; 1.1] - max_accepted(50), // depends on pb/neighborhood - max_generated(100), // depends on pb/neighborhood - theta(10), // theta is typically set to 10 - statIsInitialized(false), - outf("out.data") - { } - - moTrikiCoolingSchedule ( - Neighborhood& _neighborhood, moEval& _eval, double _initTemp, - double _max_accepted, - double _max_generated - ) - : initTemp(_initTemp), - mu2(10), // mu2 typically belongs to [1; 20] - K1(2), // K1 in [1; 4], the number of chains without reaching equilibrium before we raise the temperature - K2(5), // ??? - lambda1(2), // the increase in temperature, typically in [1.5; 4] - lambda2(.7), // lambda2 in [0.5; 0.99] - mu1(10), // target decrease in cost factor, in [2; 20] - xi(1.05), // xi typically belongs to [1; 1.1] - max_accepted(_max_accepted), // depends on pb/neighborhood - max_generated(_max_generated), // depends on pb/neighborhood - theta(10), // theta is typically set to 10 - statIsInitialized(false), - outf("out.data") - { } - - /** - * Initial temperature - * @param _solution initial solution - */ - double init(EOT & _solution) { - - accepted = generated = costs_sum = 0; - - negative_temp = equilibrium_not_reached = frozen = 0; - - reinitializing = false; - terminated = false; - statIsInitialized = false; - - /// - cout << "INIT T=" << initTemp << endl; - /// - - //outf.open("out"); - //outf << "ok"; - //outf.close(); - - - return initTemp; - } + /** + * update the temperature by a factor + * @param _temp current temperature to update + * @param _acceptedMove true when the move is accepted, false otherwise + */ + void update(double& _temp, bool _acceptedMove, EOT & _solution) { + + /* + * In the following code, things were added or modified from + * the original (incomplete) version of the algorithm + * described in [2004, Triki et al.] + * Each added/modified behavior is labelled + * with a "// ADDED!" comment. + */ + + chainStat.temperature = _temp; + chainStat.stoppingReason = NULL; + chainStat.chainEndingReason = NULL; + chainStat.equilibriumNotReached = false; + chainStat.negativeTemp = false; + chainStat.generatedSolutions = generated; + chainStat.acceptedSolutions = accepted; + + generated++; + + if (_acceptedMove) + { + accepted++; + if (statIsInitialized) + momentStat(_solution); + else momentStat.init(_solution), statIsInitialized = true; + } + + if (accepted > max_accepted || generated > max_generated) { + + chainStat.chainEndingReason = accepted > max_accepted ? chainEndingReasons[0]: chainEndingReasons[1]; + + double avgFitness = momentStat.value().first; + double prevAvgFitness = chainStat.avgFitness; + + double alpha = 0; + + if (accepted == 0) // ADDED! Otherwise the computed std dev is null; we're probably at equilibrium + { + chainStat.stoppingReason = stoppingReasons[0]; + + // Note: we could also not stop and just become greedy (temperature set to 0) + } + else + { + double avgFitness = momentStat.value().first; + double variance = momentStat.value().second; + chainStat.stdDev = sqrt(variance); + double sigma = chainStat.stdDev; + + accepted = generated = 0; + statIsInitialized = false; + + if (negative_temp < K2) + { + if (!reinitializing) + { + if (avgFitness/(prevAvgFitness-chainStat.delta) > xi) + equilibrium_not_reached++, chainStat.equilibriumNotReached = true; + else equilibrium_not_reached = 0; + } + if (equilibrium_not_reached > K1) + { + reinitializing = true; + + alpha = lambda1; + chainStat.delta = sigma/mu1; + equilibrium_not_reached = 0; // ADDED! Otherwise the algo gets trapped here! + } + else if (_temp*chainStat.delta/(sigma*sigma) >= 1) + { + negative_temp++; + reinitializing = true; + chainStat.negativeTemp = true; + + if (negative_temp < K2) + { + alpha = lambda1; + chainStat.delta = sigma/mu1; + } else + alpha = lambda2; + } + else + { + reinitializing = false; + alpha = 1-_temp*chainStat.delta/variance; + + if (sigma == 0) // ADDED! When std dev is null, the solution is probably at eq, and the algo can't go on anyways + chainStat.stoppingReason = stoppingReasons[1]; + } + } + else + { /* Note: the paper doesn't specify a value for alpha in this case. + We've chosen to let it set to 0, which means the algorithm becomes greedy. */ + alpha = 0; + } + } + + _temp *= alpha; + + chainStat.currentFitness = _solution.fitness(); + chainStat.alpha = alpha; + chainStat.avgFitness = avgFitness; + + + // TODO use a relative-epsilon comparison to approximate equality + if (avgFitness == prevAvgFitness) + frozen++; + else frozen = 0; + + if (frozen >= theta) + chainStat.stoppingReason = stoppingReasons[2]; + + } + + } + + /* + * operator() Determines if the cooling schedule shall end or continue + * @param temperature the current temperature + */ + bool operator() (double temperature) + { + + return frozen < theta + && !chainStat.stoppingReason ; // ADDED! because 'frozen' isn't a sufficient terminating criterion (yet?) + + } + + /* + * Definition of getter functions useful for monitoring the algorithm + * using an eoGetterUpdater. + */ +#define __triki_makeGetter(name, type) type name() { return chainStat.name; } + + __triki_makeGetter(stdDev, double) + __triki_makeGetter(avgFitness, double) + __triki_makeGetter(temperature, double) + __triki_makeGetter(currentFitness, double) + __triki_makeGetter(alpha, double) + __triki_makeGetter(delta, double) + + __triki_makeGetter(generatedSolutions, int) + __triki_makeGetter(acceptedSolutions, int) - /** - * update the temperature by a factor - * @param _temp current temperature to update - * @param _acceptedMove true when the move is accepted, false otherwise - */ - void update(double& _temp, bool _acceptedMove, EOT & _solution) { - - //cout << _temp << " g " << generated << endl; - - generated++; - - if (_acceptedMove) - { - accepted++; - //costs_sum += _solution.fitness(); - //varStat(_solution); - if (statIsInitialized) - momStat(_solution); - else momStat.init(_solution), statIsInitialized = true; - - //cout << _solution.fitness() << " avgCost=" << momStat.value().first << endl; - } - - - if (accepted > max_accepted || generated > max_generated) { - - if (accepted == 0) // ADDED! Otherwise the computed std dev is null; we're probably at equilibrium - { - /// - cout << "Stopping: no accepted solution" << endl; - /// - - terminated = true; - return; - } - - /// - cout << (accepted > max_accepted? "MAXACC ": "MAXGEN "); - /// - - //double avgCost = costs_sum/(double)accepted; - //double stdDev = sqrt(varStat.value()); // WARNING: IT'S NO MORE THE AVG COST, NOW IT'S THE STD DEV! - //double variance = varStat.value(); - double avgCost = momStat.value().first; - double variance = momStat.value().second; - double stdDev = sqrt(variance); - double sigma = stdDev; - double delta = sigma/mu2; - - - //outf << avgCost << endl; - //outf << _temp << endl; - outf << prevAvgCost-delta << endl; - - - accepted = generated = costs_sum = 0; - //varStat.init(_solution);//TODON use next chain's first sol - //momStat.init(_solution);//TODONE use next chain's first sol - statIsInitialized = false; - - /// - cout << "T=" << _temp << " avgCost=" << avgCost << " stdDev=" << stdDev << " currCost=" << _solution.fitness() << endl; - /// - - double alpha; - double oldprevAvgCost = prevAvgCost; - - /// - cout << "negTemp: " << negative_temp << " / " << K2 << endl; - /// - - if (negative_temp < K2) - { - if (!reinitializing) - { - /// - if (avgCost/(prevAvgCost-delta) > xi) cout << "/!\\ eq not reached!" << endl; - /// - - if (avgCost/(prevAvgCost-delta) > xi) - equilibrium_not_reached++; - else equilibrium_not_reached = 0; - } - if (equilibrium_not_reached > K1) - { - /// - cout << "/!\\ Reinitializing (eq not reached)" << endl; - /// - - reinitializing = true; - alpha = lambda1; - delta = sigma/mu1; - equilibrium_not_reached = 0; // ADDED! Otherwise the algo gets trapped here! - } - else if (_temp*delta/(sigma*sigma) >= 1) - { - /// - cout << "/!\\ neg temp!" << endl; - /// - - negative_temp++; - reinitializing = true; - if (negative_temp < K2) - { - alpha = lambda1; - delta = sigma/mu1; - } else - alpha = lambda2; - } - - // First interpretation of the pseudocode indentation: (seems obviously false because it makes the above code unreachable) - /* - } - else - { - cout << "ccc" << endl; - reinitializing = false; - prevAvgCost = avgCost; - alpha = 1-_temp*delta/(sigma*sigma); - } - */ - - // Second interpretation of the pseudocode indentation: - else - { - /// - cout << "[normal decrease]" << endl; - /// - - reinitializing = false; - prevAvgCost = avgCost; - //alpha = 1-_temp*delta/(sigma*sigma); - alpha = 1-_temp*delta/variance; - - //alpha = (sigma==0? 1: 1-_temp*delta/(sigma*sigma)); // ADDED! but removed - - if (sigma == 0) // ADDED! When std dev is null, the solution is probably at eq, and the algo can't go on anyways - terminated = true, cout << "Stopping: null std dev" << endl; - } - } - // FIXME: else what? alpha=? - - + __triki_makeGetter(equilibriumNotReached, bool) + __triki_makeGetter(negativeTemp, bool) + __triki_makeGetter(terminated, bool) + + __triki_makeGetter(stoppingReason, const char *) + __triki_makeGetter(chainEndingReason, const char *) - /// - cout << "*=" << alpha << endl; - /// - - _temp *= alpha; - - - // Never seems to be used - if (avgCost == oldprevAvgCost) // use a neighborhood to approximate double equality? - frozen++; - else frozen = 0; - - - - //exit(0); - //cin.get(); - - - } - - } +#undef __triki_makeGetter + - //! Function which proceeds to the cooling - /*! - */ - bool operator() (double temperature) - { - /// - if (terminated) cout << "TERMINATED" << endl; - /// - - return frozen < theta - && !terminated ; // ADDED! because 'frozen' doesn't terminate anything - } + const bool markovChainEnded() const + { + return chainStat.chainEndingReason != NULL; + } + + struct MarkovChainStats; + const MarkovChainStats& markovChainStats() const + { + return chainStat; + } + + struct MarkovChainStats { + + MarkovChainStats() + : stdDev(0), avgFitness(0), temperature(-1), currentFitness(-1), alpha(0), delta(0), + equilibriumNotReached(false), negativeTemp(false) + { } + + double + stdDev, + avgFitness, + temperature, + currentFitness, + alpha, + delta + ; + int generatedSolutions, acceptedSolutions; + bool equilibriumNotReached, negativeTemp; + const char * stoppingReason; // if NULL, the algo has not stopped + const char * chainEndingReason; // if NULL, the chain has not ended + + void print(std::ostream& os = std::cout, bool onlyWhenChainEnds = true) const + { + if (chainEndingReason != NULL || !onlyWhenChainEnds) + { + os << "T=" << temperature << " avgFitness=" << avgFitness << " stdDev=" << stdDev + << " currentFitness=" << currentFitness << " expected decrease in cost=" << delta + << std::endl; + if (chainEndingReason) + os << "T*=" << alpha << " chain ended, because " << chainEndingReason; + if (equilibriumNotReached) + os << " /!\\ equilibrium not reached"; + if (negativeTemp) + os << " /!\\ negative temperature"; + os << std::endl; + if (stoppingReason) + os << "Terminated, because " << stoppingReason << std::endl; + } + } + + }; + private: - //moNeighborhoodStat nhStat; - //moStdFitnessNeighborStat stdDevStat; - const double - // parameters of the algorithm - //currentTemp, - initTemp, - //ratio, - //threshold, - mu2, // mu2 typically belongs to [1; 20] - K1, // K1 in [1; 4], the number of chains without reaching equilibrium before we raise the temperature - K2, - lambda1, // the increase in temperature, typically in [1.5; 4] - lambda2, // lambda2 in [0.5; 0.99] - mu1, // target decrease in cost factor, in [2; 20] - xi // xi typically belongs to [1; 1.1] - // private variables - ; - double - stdDev, - prevAvgCost, - expectedDecreaseInCost, // delta - costs_sum - ; - const int - max_accepted, - max_generated, - theta // theta is typically set to 10 - ; - int - accepted, - generated, - equilibrium_not_reached, - negative_temp, - frozen - ; - bool reinitializing, terminated; - - //moFitnessVarianceStat varStat; - moFitnessMomentsStat momStat; - bool statIsInitialized; - - ofstream outf; - + + // parameters of the algorithm + + const double + initTemp, + initStdDev, + mu2, + K1, + K2, + lambda1, + lambda2, + mu1, + xi + ; + const int + max_accepted, + max_generated, + theta + ; + + // Variables of the algorithm + + MarkovChainStats chainStat; + + int + accepted, + generated, + equilibrium_not_reached, + negative_temp, + frozen + ; + + bool statIsInitialized, reinitializing; + + moFitnessMomentsStat momentStat; + + + // Possible reasons why the algorithm has stopped + static const char * stoppingReasons[]; + + // Possible reasons why the previous Markov chain has ended + static const char * chainEndingReasons[]; + }; + +/* + * Definition of the static members of the class + */ +template< class Neighbor > +const char * moTrikiCoolingSchedule::stoppingReasons[] = {"no accepted solution", "null std dev" , "frozen >= theta"}; +template< class Neighbor > +const char * moTrikiCoolingSchedule::chainEndingReasons[] = {"MAX ACCepted solutions", "MAX GENerated solutions"}; + + #endif + +