1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/simann.h Tue Aug 16 20:17:43 2005 +0000
1.3 @@ -0,0 +1,348 @@
1.4 +#ifndef LEMON_SIMANN_H
1.5 +#define LEMON_SIMANN_H
1.6 +
1.7 +/// \ingroup experimental
1.8 +/// \file
1.9 +/// \brief Simulated annealing framework.
1.10 +/// \author Akos Ladanyi
1.11 +
1.12 +#include <cstdlib>
1.13 +#include <cmath>
1.14 +#include <lemon/time_measure.h>
1.15 +
1.16 +namespace lemon {
1.17 +
1.18 +/// \addtogroup experimental
1.19 +/// @{
1.20 +
1.21 + /*! \brief A base class for controllers. */
1.22 + class ControllerBase {
1.23 + friend class SimAnnBase;
1.24 + public:
1.25 + /*! \brief Pointer to the simulated annealing base class. */
1.26 + SimAnnBase *simann;
1.27 + /*! \brief Initializes the controller. */
1.28 + virtual void init() {}
1.29 + /*! \brief This is called when a neighbouring state gets accepted. */
1.30 + virtual void acceptEvent() {}
1.31 + /*! \brief This is called when the accepted neighbouring state's cost is
1.32 + * less than the best found one's.
1.33 + */
1.34 + virtual void improveEvent() {}
1.35 + /*! \brief This is called when a neighbouring state gets rejected. */
1.36 + virtual void rejectEvent() {}
1.37 + /*! \brief Decides whether to continue the annealing process or not. */
1.38 + virtual bool next() = 0;
1.39 + /*! \brief Decides whether to accept the current solution or not. */
1.40 + virtual bool accept() = 0;
1.41 + };
1.42 +
1.43 + /*! \brief Skeleton of an entity class. */
1.44 + class EntityBase {
1.45 + public:
1.46 + /*! \brief Makes a minor change to the entity.
1.47 + * \return the new cost
1.48 + */
1.49 + virtual double mutate() = 0;
1.50 + /*! \brief Restores the entity to its previous state i.e. reverts the
1.51 + * effects of the last mutate().
1.52 + */
1.53 + virtual void revert() = 0;
1.54 + /*! \brief Makes a copy of the entity. */
1.55 + virtual EntityBase* clone() = 0;
1.56 + /*! \brief Makes a major change to the entity. */
1.57 + virtual void randomize() = 0;
1.58 + };
1.59 +
1.60 + /*! \brief Simulated annealing base class. */
1.61 + class SimAnnBase {
1.62 + private:
1.63 + /*! Pointer to the controller. */
1.64 + ControllerBase *controller;
1.65 + /*! \brief Cost of the current solution. */
1.66 + double curr_cost;
1.67 + /*! \brief Cost of the best solution. */
1.68 + double best_cost;
1.69 + /*! \brief Cost of the previous solution. */
1.70 + double prev_cost;
1.71 + /*! \brief Cost of the solution preceding the previous one. */
1.72 + double prev_prev_cost;
1.73 + /*! \brief Number of iterations. */
1.74 + long iter;
1.75 + /*! \brief Number of iterations which did not improve the solution since
1.76 + * the last improvement. */
1.77 + long last_impr;
1.78 + protected:
1.79 + /*! \brief Step to a neighbouring state. */
1.80 + virtual double mutate() = 0;
1.81 + /*! \brief Reverts the last mutate(). */
1.82 + virtual void revert() = 0;
1.83 + /*! \brief Saves the current solution as the best one. */
1.84 + virtual void saveAsBest() = 0;
1.85 + /*! \brief Does initializations before each run. */
1.86 + virtual void init() {
1.87 + controller->init();
1.88 + curr_cost = prev_cost = prev_prev_cost = best_cost =
1.89 + std::numeric_limits<double>::infinity();
1.90 + iter = last_impr = 0;
1.91 + }
1.92 + public:
1.93 + /*! \brief Sets the controller class to use. */
1.94 + void setController(ControllerBase &_controller) {
1.95 + controller = &_controller;
1.96 + controller->simann = this;
1.97 + }
1.98 + /*! \brief Returns the cost of the current solution. */
1.99 + double getCurrCost() const { return curr_cost; }
1.100 + /*! \brief Returns the cost of the previous solution. */
1.101 + double getPrevCost() const { return prev_cost; }
1.102 + /*! \brief Returns the cost of the best solution. */
1.103 + double getBestCost() const { return best_cost; }
1.104 + /*! \brief Returns the number of iterations. */
1.105 + long getIter() const { return iter; }
1.106 + /*! \brief Returns the number of the last iteration when the solution was
1.107 + * improved.
1.108 + */
1.109 + long getLastImpr() const { return last_impr; }
1.110 + /*! \brief Performs one iteration. */
1.111 + bool step() {
1.112 + iter++;
1.113 + prev_prev_cost = prev_cost;
1.114 + prev_cost = curr_cost;
1.115 + curr_cost = mutate();
1.116 + if (controller->accept()) {
1.117 + controller->acceptEvent();
1.118 + last_impr = iter;
1.119 + if (curr_cost < best_cost) {
1.120 + best_cost = curr_cost;
1.121 + saveAsBest();
1.122 + controller->improveEvent();
1.123 + }
1.124 + }
1.125 + else {
1.126 + revert();
1.127 + curr_cost = prev_cost;
1.128 + prev_cost = prev_prev_cost;
1.129 + controller->rejectEvent();
1.130 + }
1.131 + return controller->next();
1.132 + }
1.133 + /*! \brief Performs a given number of iterations.
1.134 + * \param n the number of iterations
1.135 + */
1.136 + bool step(int n) {
1.137 + for(; n > 0 && step(); --n) ;
1.138 + return !n;
1.139 + }
1.140 + /*! \brief Starts the annealing process. */
1.141 + void run() {
1.142 + init();
1.143 + do { } while (step());
1.144 + }
1.145 + };
1.146 +
1.147 + /*! \brief Simulated annealing class. */
1.148 + class SimAnn : public SimAnnBase {
1.149 + private:
1.150 + /*! \brief Pointer to the current entity. */
1.151 + EntityBase *curr_ent;
1.152 + /*! \brief Pointer to the best entity. */
1.153 + EntityBase *best_ent;
1.154 + /*! \brief Does initializations before each run. */
1.155 + void init() {
1.156 + SimAnnBase::init();
1.157 + if (best_ent) delete best_ent;
1.158 + best_ent = NULL;
1.159 + curr_ent->randomize();
1.160 + }
1.161 + public:
1.162 + /*! \brief Constructor. */
1.163 + SimAnn() : curr_ent(NULL), best_ent(NULL) {}
1.164 + /*! \brief Destructor. */
1.165 + virtual ~SimAnn() {
1.166 + if (best_ent) delete best_ent;
1.167 + }
1.168 + /*! \brief Step to a neighbouring state. */
1.169 + double mutate() {
1.170 + return curr_ent->mutate();
1.171 + }
1.172 + /*! \brief Reverts the last mutate(). */
1.173 + void revert() {
1.174 + curr_ent->revert();
1.175 + }
1.176 + /*! \brief Saves the current solution as the best one. */
1.177 + void saveAsBest() {
1.178 + if (best_ent) delete best_ent;
1.179 + best_ent = curr_ent->clone();
1.180 + }
1.181 + /*! \brief Sets the current entity. */
1.182 + void setEntity(EntityBase &_ent) {
1.183 + curr_ent = &_ent;
1.184 + }
1.185 + /*! \brief Returns a copy of the best found entity. */
1.186 + EntityBase* getBestEntity() { return best_ent->clone(); }
1.187 + };
1.188 +
1.189 + /*! \brief A simple controller for the simulated annealing class. */
1.190 + class SimpleController : public ControllerBase {
1.191 + public:
1.192 + /*! \brief Maximum number of iterations. */
1.193 + long max_iter;
1.194 + /*! \brief Maximum number of iterations which do not improve the
1.195 + * solution. */
1.196 + long max_no_impr;
1.197 + /*! \brief Temperature. */
1.198 + double temp;
1.199 + /*! \brief Annealing factor. */
1.200 + double ann_fact;
1.201 + /*! \brief Constructor.
1.202 + * \param _max_iter maximum number of iterations
1.203 + * \param _max_no_impr maximum number of consecutive iterations which do
1.204 + * not yield a better solution
1.205 + * \param _temp initial temperature
1.206 + * \param _ann_fact annealing factor
1.207 + */
1.208 + SimpleController(long _max_iter = 500000, long _max_no_impr = 20000,
1.209 + double _temp = 1000.0, double _ann_fact = 0.9999) : max_iter(_max_iter),
1.210 + max_no_impr(_max_no_impr), temp(_temp), ann_fact(_ann_fact)
1.211 + {
1.212 + srand48(time(0));
1.213 + }
1.214 + /*! \brief This is called when a neighbouring state gets accepted. */
1.215 + void acceptEvent() {}
1.216 + /*! \brief This is called when the accepted neighbouring state's cost is
1.217 + * less than the best found one's.
1.218 + */
1.219 + void improveEvent() {}
1.220 + /*! \brief This is called when a neighbouring state gets rejected. */
1.221 + void rejectEvent() {}
1.222 + /*! \brief Decides whether to continue the annealing process or not. Also
1.223 + * decreases the temperature. */
1.224 + bool next() {
1.225 + temp *= ann_fact;
1.226 + bool quit = (simann->getIter() > max_iter) ||
1.227 + (simann->getIter() - simann->getLastImpr() > max_no_impr);
1.228 + return !quit;
1.229 + }
1.230 + /*! \brief Decides whether to accept the current solution or not. */
1.231 + bool accept() {
1.232 + double cost_diff = simann->getPrevCost() - simann->getCurrCost();
1.233 + return (drand48() <= exp(cost_diff / temp));
1.234 + }
1.235 + };
1.236 +
1.237 + /*! \brief A controller with preset running time for the simulated annealing
1.238 + * class.
1.239 + *
1.240 + * With this controller you can set the running time of the annealing
1.241 + * process in advance. It works the following way: the controller measures
1.242 + * a kind of divergence. The divergence is the difference of the average
1.243 + * cost of the recently found solutions the cost of the best found one. In
1.244 + * case this divergence is greater than a given threshold, then we decrease
1.245 + * the annealing factor, that is we cool the system faster. In case the
1.246 + * divergence is lower than the threshold, then we increase the temperature.
1.247 + * The threshold is a function of the elapsed time which reaches zero at the
1.248 + * desired end time.
1.249 + */
1.250 + class AdvancedController : public ControllerBase {
1.251 + private:
1.252 + Timer timer;
1.253 + /*! \param time the elapsed time in seconds */
1.254 + virtual double threshold(double time) {
1.255 + return (-1.0) * start_threshold / end_time * time + start_threshold;
1.256 + }
1.257 + public:
1.258 + double alpha;
1.259 + double beta;
1.260 + double gamma;
1.261 + /*! \brief The time at the end of the algorithm. */
1.262 + double end_time;
1.263 + /*! \brief The time at the start of the algorithm. */
1.264 + double start_time;
1.265 + /*! \brief Starting threshold. */
1.266 + double start_threshold;
1.267 + /*! \brief Average cost of recent solutions. */
1.268 + double avg_cost;
1.269 + /*! \brief Temperature. */
1.270 + double temp;
1.271 + /*! \brief Annealing factor. */
1.272 + double ann_fact;
1.273 + /*! \brief Initial annealing factor. */
1.274 + double init_ann_fact;
1.275 + bool warmup;
1.276 + /*! \brief Constructor.
1.277 + * \param _end_time running time in seconds
1.278 + * \param _alpha parameter used to calculate the running average
1.279 + * \param _beta parameter used to decrease the annealing factor
1.280 + * \param _gamma parameter used to increase the temperature
1.281 + * \param _ann_fact initial annealing factor
1.282 + */
1.283 + AdvancedController(double _end_time, double _alpha = 0.2,
1.284 + double _beta = 0.9, double _gamma = 1.6, double _ann_fact = 0.9999) :
1.285 + alpha(_alpha), beta(_beta), gamma(_gamma), end_time(_end_time),
1.286 + ann_fact(_ann_fact), init_ann_fact(_ann_fact), warmup(true)
1.287 + {
1.288 + srand48(time(0));
1.289 + }
1.290 + void init() {
1.291 + avg_cost = simann->getCurrCost();
1.292 + }
1.293 + /*! \brief This is called when a neighbouring state gets accepted. */
1.294 + void acceptEvent() {
1.295 + avg_cost = alpha * simann->getCurrCost() + (1.0 - alpha) * avg_cost;
1.296 + if (warmup) {
1.297 + static int cnt = 0;
1.298 + cnt++;
1.299 + if (cnt >= 100) {
1.300 + // calculate starting threshold and starting temperature
1.301 + start_threshold = 5.0 * fabs(simann->getBestCost() - avg_cost);
1.302 + temp = 10000.0;
1.303 + warmup = false;
1.304 + timer.reset();
1.305 + }
1.306 + }
1.307 + }
1.308 + /*! \brief Decides whether to continue the annealing process or not. */
1.309 + bool next() {
1.310 + if (warmup) {
1.311 + return true;
1.312 + }
1.313 + else {
1.314 + double elapsed_time = timer.getRealTime();
1.315 + if (fabs(avg_cost - simann->getBestCost()) > threshold(elapsed_time)) {
1.316 + // decrease the annealing factor
1.317 + ann_fact *= beta;
1.318 + }
1.319 + else {
1.320 + // increase the temperature
1.321 + temp *= gamma;
1.322 + // reset the annealing factor
1.323 + ann_fact = init_ann_fact;
1.324 + }
1.325 + temp *= ann_fact;
1.326 + return elapsed_time < end_time;
1.327 + }
1.328 + }
1.329 + /*! \brief Decides whether to accept the current solution or not. */
1.330 + bool accept() {
1.331 + if (warmup) {
1.332 + // we accept eveything during the "warm up" phase
1.333 + return true;
1.334 + }
1.335 + else {
1.336 + double cost_diff = simann->getPrevCost() - simann->getCurrCost();
1.337 + if (cost_diff < 0.0) {
1.338 + return (drand48() <= exp(cost_diff / temp));
1.339 + }
1.340 + else {
1.341 + return true;
1.342 + }
1.343 + }
1.344 + }
1.345 + };
1.346 +
1.347 +/// @}
1.348 +
1.349 +}
1.350 +
1.351 +#endif