ladanyi@942: #ifndef LEMON_SIMANN_H
ladanyi@942: #define LEMON_SIMANN_H
ladanyi@918: 
ladanyi@1142: /// \ingroup experimental
ladanyi@1142: /// \file
ladanyi@1142: /// \brief Simulated annealing framework.
ladanyi@1142: /// \author Akos Ladanyi
ladanyi@1142: 
ladanyi@966: #include <cstdlib>
ladanyi@966: #include <cmath>
ladanyi@1018: #include <lemon/time_measure.h>
ladanyi@966: 
ladanyi@942: namespace lemon {
ladanyi@918: 
ladanyi@1142: /// \addtogroup experimental
ladanyi@1142: /// @{
ladanyi@1142: 
ladanyi@942:   const double INFTY = 1e24;
ladanyi@918: 
ladanyi@1142:   /*! \brief Simulated annealing base class. */
ladanyi@942:   class SimAnnBase {
ladanyi@918:   public:
ladanyi@942:     class Controller;
ladanyi@942:   private:
ladanyi@1142:     /*! Pointer to the controller. */
ladanyi@942:     Controller *controller;
ladanyi@942:   protected:
ladanyi@1142:     /*! \brief Cost of the current solution. */
ladanyi@942:     double curr_cost;
ladanyi@1142:     /*! \brief Cost of the best solution. */
ladanyi@1023:     double best_cost;
ladanyi@1142:     /*! \brief Cost of the previous solution. */
ladanyi@942:     double prev_cost;
ladanyi@1142:     /*! \brief Cost of the solution preceding the previous one. */
ladanyi@1023:     double prev_prev_cost;
ladanyi@918: 
ladanyi@1142:     /*! \brief Step to a neighbouring state. */
alpar@1150:     virtual void mutate() = 0;
ladanyi@1142:     /*! \brief Reverts the last mutate(). */
alpar@1150:     virtual void revert() = 0;
ladanyi@1142:     /*! \brief Saves the current solution as the best one. */
alpar@1150:     virtual void saveAsBest() = 0;
ladanyi@942:   public:
ladanyi@1142:     /*! \brief Constructor. */
ladanyi@942:     SimAnnBase() {
ladanyi@1023:       best_cost = prev_cost = prev_prev_cost = INFTY;
ladanyi@942:     }
ladanyi@1142:     /*! \brief Sets the controller class to use. */
ladanyi@957:     void setController(Controller &_controller) {
ladanyi@957:       controller = &_controller;
ladanyi@957:       controller->setBase(this);
ladanyi@957:     }
ladanyi@1142:     /*! \brief Returns the cost of the current solution. */
ladanyi@1018:     double getCurrCost() const { return curr_cost; }
ladanyi@1142:     /*! \brief Returns the cost of the previous solution. */
ladanyi@1018:     double getPrevCost() const { return prev_cost; }
ladanyi@1142:     /*! \brief Returns the cost of the best solution. */
ladanyi@1018:     double getBestCost() const { return best_cost; }
ladanyi@1142:     /*! \brief Starts the annealing process. */
ladanyi@942:     void run() {
ladanyi@966:       controller->init();
ladanyi@1018:       do {
alpar@1150:         curr_cost=mutate();
ladanyi@957:         if (controller->accept()) {
ladanyi@942:           controller->acceptEvent();
ladanyi@942:           if (curr_cost < best_cost) {
ladanyi@942:             saveAsBest();
ladanyi@942:             controller->improveEvent();
ladanyi@942:           }
ladanyi@942:         }
ladanyi@942:         else {
ladanyi@942:           revert();
ladanyi@942:           controller->rejectEvent();
ladanyi@942:         }
ladanyi@1018:       } while (controller->next());
ladanyi@918:     }
ladanyi@918: 
ladanyi@1000:     /*! \brief A base class for controllers. */
ladanyi@942:     class Controller {
ladanyi@942:     public:
ladanyi@1142:       /*! \brief Pointer to the simulated annealing base class. */
ladanyi@957:       SimAnnBase *base;
ladanyi@1142:       /*! \brief Initializes the controller. */
ladanyi@966:       virtual void init() {}
ladanyi@1000:       /*! \brief This is called when a neighbouring state gets accepted. */
ladanyi@942:       virtual void acceptEvent() {}
ladanyi@1000:       /*! \brief This is called when the accepted neighbouring state's cost is
ladanyi@1000:        *  less than the best found one's.
ladanyi@1000:        */
ladanyi@942:       virtual void improveEvent() {}
ladanyi@1000:       /*! \brief This is called when a neighbouring state gets rejected. */
ladanyi@942:       virtual void rejectEvent() {}
ladanyi@1142:       /*! \brief Sets the simulated annealing base class to use. */
ladanyi@957:       virtual void setBase(SimAnnBase *_base) { base = _base; }
ladanyi@1142:       /*! \brief Decides whether to continue the annealing process or not. */
ladanyi@942:       virtual bool next() = 0;
ladanyi@1142:       /*! \brief Decides whether to accept the current solution or not. */
ladanyi@957:       virtual bool accept() = 0;
ladanyi@942:     };
ladanyi@942:   };
ladanyi@918: 
ladanyi@1142:   /*! \brief Simulated annealing class. */
ladanyi@942:   template <typename E>
ladanyi@942:   class SimAnn : public SimAnnBase {
ladanyi@942:   private:
ladanyi@1142:     /*! \brief Pointer to the current entity. */
ladanyi@942:     E *curr_ent;
ladanyi@1142:     /*! \brief Pointer to the best entity. */
ladanyi@942:     E *best_ent;
ladanyi@942:   public:
ladanyi@1142:     /*! \brief Constructor. */
ladanyi@957:     SimAnn() : SimAnnBase() {}
ladanyi@1142:     /*! \brief Sets the initial entity. */
ladanyi@957:     void setEntity(E &ent) {
ladanyi@957:       curr_ent = new E(ent);
ladanyi@957:       best_ent = new E(ent);
ladanyi@1023:       curr_cost = curr_ent->getCost();
ladanyi@942:     }
ladanyi@1142:     /*! \brief Returns the best found entity. */
ladanyi@942:     E getBestEntity() { return *best_ent; }
ladanyi@1142:     /*! \brief Step to a neighbouring state. */
ladanyi@942:     void mutate() {
ladanyi@1023:       prev_prev_cost = prev_cost;
ladanyi@1018:       prev_cost = curr_cost;
ladanyi@1023:       curr_ent->mutate();
ladanyi@1023:       curr_cost = curr_ent->getCost();
ladanyi@942:     }
ladanyi@1142:     /*! \brief Reverts the last mutate(). */
ladanyi@942:     void revert() {
ladanyi@942:       curr_ent->revert();
ladanyi@1018:       curr_cost = prev_cost;
ladanyi@1023:       prev_cost = prev_prev_cost;
ladanyi@942:     }
ladanyi@1142:     /*! \brief Saves the current solution as the best one. */
ladanyi@942:     void saveAsBest() {
ladanyi@1096:       delete(best_ent);
ladanyi@1096:       best_ent = new E(*curr_ent);
ladanyi@942:       best_cost = curr_cost;
ladanyi@942:     }
ladanyi@942:   };
ladanyi@942: 
ladanyi@1142:   /*! \brief Skeleton of an entity class. */
ladanyi@956:   class EntitySkeleton {
ladanyi@942:   public:
ladanyi@1142:     /*! \brief Returns the cost of the entity. */
ladanyi@1023:     double getCost() { return 0.0; }
ladanyi@1023:     /*! \brief Makes a minor change to the entity. */
ladanyi@1023:     void mutate() {}
ladanyi@966:     /*! \brief Restores the entity to its previous state i.e. reverts the
ladanyi@1142:      *  effects of the last mutate().
ladanyi@966:      */
ladanyi@942:     void revert() {}
ladanyi@942:   };
ladanyi@942: 
ladanyi@1142:   /*! \brief A simple controller for the simulated annealing class. */
ladanyi@956:   class SimpleController : public SimAnnBase::Controller {
ladanyi@956:   public:
ladanyi@1142:     /*! \brief Number of iterations. */
ladanyi@1142:     long iter;
ladanyi@1142:     /*! \brief Number of iterations which did not improve the solution since
ladanyi@1142:      *  the last improvement. */
ladanyi@1142:     long last_impr;
ladanyi@1142:     /*! \brief Maximum number of iterations. */
ladanyi@1142:     long max_iter;
ladanyi@1142:     /*! \brief Maximum number of iterations which do not improve the
ladanyi@1142:      *  solution. */
ladanyi@1142:     long max_no_impr;
ladanyi@1142:     /*! \brief Temperature. */
ladanyi@1142:     double temp;
ladanyi@1142:     /*! \brief Annealing factor. */
ladanyi@1142:     double ann_fact;
ladanyi@1142:     /*! \brief Constructor.
ladanyi@1142:      *  \param _max_iter maximum number of iterations
ladanyi@1000:      *  \param _max_no_impr maximum number of consecutive iterations which do
ladanyi@1000:      *         not yield a better solution
ladanyi@1000:      *  \param _temp initial temperature
ladanyi@1000:      *  \param _ann_fact annealing factor
ladanyi@1000:      */
ladanyi@1000:     SimpleController(long _max_iter = 500000, long _max_no_impr = 20000,
ladanyi@1096:     double _temp = 1000.0, double _ann_fact = 0.9999) : iter(0), last_impr(0),
ladanyi@1000:     max_iter(_max_iter), max_no_impr(_max_no_impr), temp(_temp),
ladanyi@1000:     ann_fact(_ann_fact) {}
ladanyi@1145:     /*! \brief This is called when a neighbouring state gets accepted. */
ladanyi@956:     void acceptEvent() {
ladanyi@956:       iter++;
ladanyi@956:     }
ladanyi@1142:     /*! \brief This is called when the accepted neighbouring state's cost is
ladanyi@1142:      *  less than the best found one's.
ladanyi@1142:      */
ladanyi@956:     void improveEvent() {
ladanyi@956:       last_impr = iter;
ladanyi@956:     }
ladanyi@1142:     /*! \brief This is called when a neighbouring state gets rejected. */
ladanyi@956:     void rejectEvent() {
ladanyi@956:       iter++;
ladanyi@956:     }
ladanyi@1142:     /*! \brief Decides whether to continue the annealing process or not. Also
ladanyi@1142:      *  decreases the temperature. */
ladanyi@956:     bool next() {
ladanyi@1000:       temp *= ann_fact;
ladanyi@956:       bool quit = (iter > max_iter) || (iter - last_impr > max_no_impr);
ladanyi@956:       return !quit;
ladanyi@956:     }
ladanyi@1142:     /*! \brief Decides whether to accept the current solution or not. */
ladanyi@957:     bool accept() {
ladanyi@1018:       double cost_diff = base->getPrevCost() - base->getCurrCost();
ladanyi@1018:       if (cost_diff < 0.0) {
ladanyi@1096:         bool ret = drand48() <= exp(cost_diff / temp);
ladanyi@1096:         return ret;
ladanyi@1018:       }
ladanyi@1018:       else {
ladanyi@1018:         return true;
ladanyi@1018:       }
ladanyi@966:     }
ladanyi@966:   };
ladanyi@966: 
ladanyi@966:   /*! \brief A controller with preset running time for the simulated annealing
ladanyi@966:    *  class.
ladanyi@1145:    *
ladanyi@1145:    *  With this controller you can set the running time of the annealing
ladanyi@1145:    *  process in advance. It works the following way: the controller measures
ladanyi@1145:    *  a kind of divergence. The divergence is the difference of the average
ladanyi@1145:    *  cost of the recently found solutions the cost of the best found one. In
ladanyi@1145:    *  case this divergence is greater than a given threshold, then we decrease
ladanyi@1145:    *  the annealing factor, that is we cool the system faster. In case the
ladanyi@1145:    *  divergence is lower than the threshold, then we increase the temperature.
ladanyi@1145:    *  The threshold is a function of the elapsed time which reaches zero at the
ladanyi@1145:    *  desired end time.
ladanyi@966:    */
ladanyi@966:   class AdvancedController : public SimAnnBase::Controller {
ladanyi@966:   private:
ladanyi@1018:     Timer timer;
ladanyi@1000:     /*! \param time the elapsed time in seconds */
ladanyi@1018:     virtual double threshold(double time) {
ladanyi@1096:       return (-1.0) * start_threshold / end_time * time + start_threshold;
ladanyi@1018:     }
ladanyi@966:   public:
ladanyi@1142:     double alpha;
ladanyi@1142:     double beta;
ladanyi@1142:     double gamma;
ladanyi@1145:     /*! \brief The time at the end of the algorithm. */
ladanyi@1142:     double end_time;
ladanyi@1145:     /*! \brief The time at the start of the algorithm. */
ladanyi@1142:     double start_time;
ladanyi@1145:     /*! \brief Starting threshold. */
ladanyi@1018:     double start_threshold;
ladanyi@1145:     /*! \brief Average cost of recent solutions. */
ladanyi@966:     double avg_cost;
ladanyi@1145:     /*! \brief Temperature. */
ladanyi@1142:     double temp;
ladanyi@1145:     /*! \brief Annealing factor. */
ladanyi@1142:     double ann_fact;
ladanyi@1145:     /*! \brief Initial annealing factor. */
ladanyi@1145:     double init_ann_fact;
ladanyi@1018:     bool warmup;
ladanyi@1142:     /*! \brief Constructor.
ladanyi@1142:      *  \param _end_time running time in seconds
ladanyi@1000:      *  \param _alpha parameter used to calculate the running average
ladanyi@1000:      *  \param _beta parameter used to decrease the annealing factor
ladanyi@1000:      *  \param _gamma parameter used to increase the temperature
ladanyi@1145:      *  \param _ann_fact initial annealing factor
ladanyi@1000:      */
ladanyi@1000:     AdvancedController(double _end_time, double _alpha = 0.2,
ladanyi@1145:     double _beta = 0.9, double _gamma = 1.6, double _ann_fact = 0.9999) :
ladanyi@1145:     alpha(_alpha), beta(_beta), gamma(_gamma), end_time(_end_time),
ladanyi@1145:     ann_fact(_ann_fact), init_ann_fact(_ann_fact), warmup(true) {}
ladanyi@966:     void init() {
ladanyi@1018:       avg_cost = base->getCurrCost();
ladanyi@966:     }
ladanyi@1142:     /*! \brief This is called when a neighbouring state gets accepted. */
ladanyi@966:     void acceptEvent() {
ladanyi@966:       avg_cost = alpha * base->getCurrCost() + (1.0 - alpha) * avg_cost;
ladanyi@1023:       if (warmup) {
ladanyi@1096:         static int cnt = 0;
ladanyi@1096:         cnt++;
ladanyi@1096:         if (cnt >= 100) {
ladanyi@1023:           // calculate starting threshold and starting temperature
ladanyi@1096:           start_threshold = 5.0 * fabs(base->getBestCost() - avg_cost);
ladanyi@1096:           temp = 10000.0;
ladanyi@1023:           warmup = false;
ladanyi@1023:           timer.reset();
ladanyi@1023:         }
ladanyi@1023:       }
ladanyi@966:     }
ladanyi@1142:     /*! \brief Decides whether to continue the annealing process or not. */
ladanyi@966:     bool next() {
ladanyi@1018:       if (warmup) {
ladanyi@1018:         return true;
ladanyi@1000:       }
ladanyi@1000:       else {
ladanyi@1018:         double elapsed_time = timer.getRealTime();
ladanyi@1018:         if (fabs(avg_cost - base->getBestCost()) > threshold(elapsed_time)) {
ladanyi@1018:           // decrease the annealing factor
ladanyi@1018:           ann_fact *= beta;
ladanyi@1018:         }
ladanyi@1018:         else {
ladanyi@1018:           // increase the temperature
ladanyi@1018:           temp *= gamma;
ladanyi@1145:           // reset the annealing factor
ladanyi@1145:           ann_fact = init_ann_fact;
ladanyi@1018:         }
ladanyi@1018:         temp *= ann_fact;
ladanyi@1018:         return elapsed_time < end_time;
ladanyi@1000:       }
ladanyi@966:     }
ladanyi@1142:     /*! \brief Decides whether to accept the current solution or not. */
ladanyi@966:     bool accept() {
ladanyi@1018:       if (warmup) {
ladanyi@1018:         // we accept eveything during the "warm up" phase
ladanyi@1018:         return true;
ladanyi@1018:       }
ladanyi@1018:       else {
ladanyi@1018:         double cost_diff = base->getPrevCost() - base->getCurrCost();
ladanyi@1018:         if (cost_diff < 0.0) {
ladanyi@1018:           return (drand48() <= exp(cost_diff / temp));
ladanyi@1018:         }
ladanyi@1018:         else {
ladanyi@1018:           return true;
ladanyi@1018:         }
ladanyi@1018:       }
ladanyi@956:     }
ladanyi@956:   };
ladanyi@956: 
ladanyi@1142: /// @}
ladanyi@1142: 
ladanyi@942: }
ladanyi@918: 
ladanyi@918: #endif