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