source: lemon-0.x/src/work/akos/simann.h @ 1150:c20bcf71efe3

#ifndef LEMON_SIMANN_H
#define LEMON_SIMANN_H

/// \ingroup experimental
/// \file
/// \brief Simulated annealing framework.
/// \author Akos Ladanyi

#include <cstdlib>
#include <cmath>
#include <lemon/time_measure.h>

namespace lemon {

/// \addtogroup experimental
/// @{

  const double INFTY = 1e24;

  /*! \brief Simulated annealing base class. */
  class SimAnnBase {
  public:
    class Controller;
  private:
    /*! Pointer to the controller. */
    Controller *controller;
  protected:
    /*! \brief Cost of the current solution. */
    double curr_cost;
    /*! \brief Cost of the best solution. */
    double best_cost;
    /*! \brief Cost of the previous solution. */
    double prev_cost;
    /*! \brief Cost of the solution preceding the previous one. */
    double prev_prev_cost;

    /*! \brief Step to a neighbouring state. */
    virtual void mutate() = 0;
    /*! \brief Reverts the last mutate(). */
    virtual void revert() = 0;
    /*! \brief Saves the current solution as the best one. */
    virtual void saveAsBest() = 0;
  public:
    /*! \brief Constructor. */
    SimAnnBase() {
      best_cost = prev_cost = prev_prev_cost = INFTY;
    }
    /*! \brief Sets the controller class to use. */
    void setController(Controller &_controller) {
      controller = &_controller;
      controller->setBase(this);
    }
    /*! \brief Returns the cost of the current solution. */
    double getCurrCost() const { return curr_cost; }
    /*! \brief Returns the cost of the previous solution. */
    double getPrevCost() const { return prev_cost; }
    /*! \brief Returns the cost of the best solution. */
    double getBestCost() const { return best_cost; }
    /*! \brief Starts the annealing process. */
    void run() {
      controller->init();
      do {
        curr_cost=mutate();
        if (controller->accept()) {
          controller->acceptEvent();
          if (curr_cost < best_cost) {
            saveAsBest();
            controller->improveEvent();
          }
        }
        else {
          revert();
          controller->rejectEvent();
        }
      } while (controller->next());
    }

    /*! \brief A base class for controllers. */
    class Controller {
    public:
      /*! \brief Pointer to the simulated annealing base class. */
      SimAnnBase *base;
      /*! \brief Initializes the controller. */
      virtual void init() {}
      /*! \brief This is called when a neighbouring state gets accepted. */
      virtual void acceptEvent() {}
      /*! \brief This is called when the accepted neighbouring state's cost is
       *  less than the best found one's.
       */
      virtual void improveEvent() {}
      /*! \brief This is called when a neighbouring state gets rejected. */
      virtual void rejectEvent() {}
      /*! \brief Sets the simulated annealing base class to use. */
      virtual void setBase(SimAnnBase *_base) { base = _base; }
      /*! \brief Decides whether to continue the annealing process or not. */
      virtual bool next() = 0;
      /*! \brief Decides whether to accept the current solution or not. */
      virtual bool accept() = 0;
    };
  };

  /*! \brief Simulated annealing class. */
  template <typename E>
  class SimAnn : public SimAnnBase {
  private:
    /*! \brief Pointer to the current entity. */
    E *curr_ent;
    /*! \brief Pointer to the best entity. */
    E *best_ent;
  public:
    /*! \brief Constructor. */
    SimAnn() : SimAnnBase() {}
    /*! \brief Sets the initial entity. */
    void setEntity(E &ent) {
      curr_ent = new E(ent);
      best_ent = new E(ent);
      curr_cost = curr_ent->getCost();
    }
    /*! \brief Returns the best found entity. */
    E getBestEntity() { return *best_ent; }
    /*! \brief Step to a neighbouring state. */
    void mutate() {
      prev_prev_cost = prev_cost;
      prev_cost = curr_cost;
      curr_ent->mutate();
      curr_cost = curr_ent->getCost();
    }
    /*! \brief Reverts the last mutate(). */
    void revert() {
      curr_ent->revert();
      curr_cost = prev_cost;
      prev_cost = prev_prev_cost;
    }
    /*! \brief Saves the current solution as the best one. */
    void saveAsBest() {
      delete(best_ent);
      best_ent = new E(*curr_ent);
      best_cost = curr_cost;
    }
  };

  /*! \brief Skeleton of an entity class. */
  class EntitySkeleton {
  public:
    /*! \brief Returns the cost of the entity. */
    double getCost() { return 0.0; }
    /*! \brief Makes a minor change to the entity. */
    void mutate() {}
    /*! \brief Restores the entity to its previous state i.e. reverts the
     *  effects of the last mutate().
     */
    void revert() {}
  };

  /*! \brief A simple controller for the simulated annealing class. */
  class SimpleController : public SimAnnBase::Controller {
  public:
    /*! \brief Number of iterations. */
    long iter;
    /*! \brief Number of iterations which did not improve the solution since
     *  the last improvement. */
    long last_impr;
    /*! \brief Maximum number of iterations. */
    long max_iter;
    /*! \brief Maximum number of iterations which do not improve the
     *  solution. */
    long max_no_impr;
    /*! \brief Temperature. */
    double temp;
    /*! \brief Annealing factor. */
    double ann_fact;
    /*! \brief Constructor.
     *  \param _max_iter maximum number of iterations
     *  \param _max_no_impr maximum number of consecutive iterations which do
     *         not yield a better solution
     *  \param _temp initial temperature
     *  \param _ann_fact annealing factor
     */
    SimpleController(long _max_iter = 500000, long _max_no_impr = 20000,
    double _temp = 1000.0, double _ann_fact = 0.9999) : iter(0), last_impr(0),
    max_iter(_max_iter), max_no_impr(_max_no_impr), temp(_temp),
    ann_fact(_ann_fact) {}
    /*! \brief This is called when a neighbouring state gets accepted. */
    void acceptEvent() {
      iter++;
    }
    /*! \brief This is called when the accepted neighbouring state's cost is
     *  less than the best found one's.
     */
    void improveEvent() {
      last_impr = iter;
    }
    /*! \brief This is called when a neighbouring state gets rejected. */
    void rejectEvent() {
      iter++;
    }
    /*! \brief Decides whether to continue the annealing process or not. Also
     *  decreases the temperature. */
    bool next() {
      temp *= ann_fact;
      bool quit = (iter > max_iter) || (iter - last_impr > max_no_impr);
      return !quit;
    }
    /*! \brief Decides whether to accept the current solution or not. */
    bool accept() {
      double cost_diff = base->getPrevCost() - base->getCurrCost();
      if (cost_diff < 0.0) {
        bool ret = drand48() <= exp(cost_diff / temp);
        return ret;
      }
      else {
        return true;
      }
    }
  };

  /*! \brief A controller with preset running time for the simulated annealing
   *  class.
   *
   *  With this controller you can set the running time of the annealing
   *  process in advance. It works the following way: the controller measures
   *  a kind of divergence. The divergence is the difference of the average
   *  cost of the recently found solutions the cost of the best found one. In
   *  case this divergence is greater than a given threshold, then we decrease
   *  the annealing factor, that is we cool the system faster. In case the
   *  divergence is lower than the threshold, then we increase the temperature.
   *  The threshold is a function of the elapsed time which reaches zero at the
   *  desired end time.
   */
  class AdvancedController : public SimAnnBase::Controller {
  private:
    Timer timer;
    /*! \param time the elapsed time in seconds */
    virtual double threshold(double time) {
      return (-1.0) * start_threshold / end_time * time + start_threshold;
    }
  public:
    double alpha;
    double beta;
    double gamma;
    /*! \brief The time at the end of the algorithm. */
    double end_time;
    /*! \brief The time at the start of the algorithm. */
    double start_time;
    /*! \brief Starting threshold. */
    double start_threshold;
    /*! \brief Average cost of recent solutions. */
    double avg_cost;
    /*! \brief Temperature. */
    double temp;
    /*! \brief Annealing factor. */
    double ann_fact;
    /*! \brief Initial annealing factor. */
    double init_ann_fact;
    bool warmup;
    /*! \brief Constructor.
     *  \param _end_time running time in seconds
     *  \param _alpha parameter used to calculate the running average
     *  \param _beta parameter used to decrease the annealing factor
     *  \param _gamma parameter used to increase the temperature
     *  \param _ann_fact initial annealing factor
     */
    AdvancedController(double _end_time, double _alpha = 0.2,
    double _beta = 0.9, double _gamma = 1.6, double _ann_fact = 0.9999) :
    alpha(_alpha), beta(_beta), gamma(_gamma), end_time(_end_time),
    ann_fact(_ann_fact), init_ann_fact(_ann_fact), warmup(true) {}
    void init() {
      avg_cost = base->getCurrCost();
    }
    /*! \brief This is called when a neighbouring state gets accepted. */
    void acceptEvent() {
      avg_cost = alpha * base->getCurrCost() + (1.0 - alpha) * avg_cost;
      if (warmup) {
        static int cnt = 0;
        cnt++;
        if (cnt >= 100) {
          // calculate starting threshold and starting temperature
          start_threshold = 5.0 * fabs(base->getBestCost() - avg_cost);
          temp = 10000.0;
          warmup = false;
          timer.reset();
        }
      }
    }
    /*! \brief Decides whether to continue the annealing process or not. */
    bool next() {
      if (warmup) {
        return true;
      }
      else {
        double elapsed_time = timer.getRealTime();
        if (fabs(avg_cost - base->getBestCost()) > threshold(elapsed_time)) {
          // decrease the annealing factor
          ann_fact *= beta;
        }
        else {
          // increase the temperature
          temp *= gamma;
          // reset the annealing factor
          ann_fact = init_ann_fact;
        }
        temp *= ann_fact;
        return elapsed_time < end_time;
      }
    }
    /*! \brief Decides whether to accept the current solution or not. */
    bool accept() {
      if (warmup) {
        // we accept eveything during the "warm up" phase
        return true;
      }
      else {
        double cost_diff = base->getPrevCost() - base->getCurrCost();
        if (cost_diff < 0.0) {
          return (drand48() <= exp(cost_diff / temp));
        }
        else {
          return true;
        }
      }
    }
  };

/// @}

}

#endif