#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