#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