COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/akos/simann.h @ 1059:bd97feae7d90

Last change on this file since 1059:bd97feae7d90 was 1023:3268fef5d623, checked in by Akos Ladanyi, 19 years ago

Added a getCost() method to the Entity. Now prevCost() returns what its name suggests.

File size: 7.0 KB
RevLine 
[942]1#ifndef LEMON_SIMANN_H
2#define LEMON_SIMANN_H
[918]3
[966]4#include <cstdlib>
5#include <cmath>
[1018]6#include <lemon/time_measure.h>
[966]7
[942]8namespace lemon {
[918]9
[942]10  const double INFTY = 1e24;
[918]11
[942]12  class SimAnnBase {
[918]13  public:
[942]14    class Controller;
15  private:
16    Controller *controller;
17  protected:
18    double curr_cost;
[1023]19    double best_cost;
[942]20    double prev_cost;
[1023]21    double prev_prev_cost;
[918]22
[942]23    virtual void mutate() = 0;
24    virtual void revert() = 0;
25    virtual void saveAsBest() = 0;
26  public:
27    SimAnnBase() {
[1023]28      best_cost = prev_cost = prev_prev_cost = INFTY;
[942]29    }
[957]30    void setController(Controller &_controller) {
31      controller = &_controller;
32      controller->setBase(this);
33    }
[1018]34    double getCurrCost() const { return curr_cost; }
35    double getPrevCost() const { return prev_cost; }
36    double getBestCost() const { return best_cost; }
[942]37    void run() {
[966]38      controller->init();
[1018]39      do {
[942]40        mutate();
[957]41        if (controller->accept()) {
[942]42          controller->acceptEvent();
43          if (curr_cost < best_cost) {
44            saveAsBest();
45            controller->improveEvent();
46          }
47        }
48        else {
49          revert();
50          controller->rejectEvent();
51        }
[1018]52      } while (controller->next());
[918]53    }
54
[1000]55    /*! \brief A base class for controllers. */
[942]56    class Controller {
57    public:
[957]58      SimAnnBase *base;
[966]59      virtual void init() {}
[1000]60      /*! \brief This is called when a neighbouring state gets accepted. */
[942]61      virtual void acceptEvent() {}
[1000]62      /*! \brief This is called when the accepted neighbouring state's cost is
63       *  less than the best found one's.
64       */
[942]65      virtual void improveEvent() {}
[1000]66      /*! \brief This is called when a neighbouring state gets rejected. */
[942]67      virtual void rejectEvent() {}
[957]68      virtual void setBase(SimAnnBase *_base) { base = _base; }
[1000]69      /*! */
[942]70      virtual bool next() = 0;
[1000]71      /*! */
[957]72      virtual bool accept() = 0;
[942]73    };
74  };
[918]75
[942]76  template <typename E>
77  class SimAnn : public SimAnnBase {
78  private:
79    E *curr_ent;
80    E *best_ent;
81  public:
[957]82    SimAnn() : SimAnnBase() {}
83    void setEntity(E &ent) {
84      curr_ent = new E(ent);
85      best_ent = new E(ent);
[1023]86      curr_cost = curr_ent->getCost();
[942]87    }
88    E getBestEntity() { return *best_ent; }
89    void mutate() {
[1023]90      prev_prev_cost = prev_cost;
[1018]91      prev_cost = curr_cost;
[1023]92      curr_ent->mutate();
93      curr_cost = curr_ent->getCost();
[942]94    }
95    void revert() {
96      curr_ent->revert();
[1018]97      curr_cost = prev_cost;
[1023]98      prev_cost = prev_prev_cost;
[942]99    }
100    void saveAsBest() {
101      *best_ent = *curr_ent;
102      best_cost = curr_cost;
103    }
104  };
105
[956]106  class EntitySkeleton {
[942]107  public:
[1023]108    /*! \return the cost of the entity */
109    double getCost() { return 0.0; }
110    /*! \brief Makes a minor change to the entity. */
111    void mutate() {}
[966]112    /*! \brief Restores the entity to its previous state i.e. reverts the
113     *  effects of the last mutate.
114     */
[942]115    void revert() {}
116  };
117
[966]118  /*! \brief A simple controller for the simulated annealing class.
119   *  \todo Find a way to set the various parameters.
120   */
[956]121  class SimpleController : public SimAnnBase::Controller {
122  public:
123    long iter, last_impr, max_iter, max_no_impr;
[1000]124    double temp, ann_fact;
125    /*! \param _max_iter maximum number of iterations
126     *  \param _max_no_impr maximum number of consecutive iterations which do
127     *         not yield a better solution
128     *  \param _temp initial temperature
129     *  \param _ann_fact annealing factor
130     */
131    SimpleController(long _max_iter = 500000, long _max_no_impr = 20000,
132    double _temp = 1000, double _ann_fact = 0.9999) : iter(0), last_impr(0),
133    max_iter(_max_iter), max_no_impr(_max_no_impr), temp(_temp),
134    ann_fact(_ann_fact) {}
[956]135    void acceptEvent() {
136      iter++;
137    }
138    void improveEvent() {
139      last_impr = iter;
140    }
141    void rejectEvent() {
142      iter++;
143    }
144    bool next() {
[1000]145      temp *= ann_fact;
[956]146      bool quit = (iter > max_iter) || (iter - last_impr > max_no_impr);
147      return !quit;
148    }
[957]149    bool accept() {
[1018]150      double cost_diff = base->getPrevCost() - base->getCurrCost();
151      if (cost_diff < 0.0) {
152        return (drand48() <= exp(cost_diff / temp));
153      }
154      else {
155        return true;
156      }
[966]157    }
158  };
159
160  /*! \brief A controller with preset running time for the simulated annealing
161   *  class.
162   *  \todo Find a better name.
163   */
164  class AdvancedController : public SimAnnBase::Controller {
165  private:
[1018]166    Timer timer;
[1000]167    /*! \param time the elapsed time in seconds */
[1018]168    virtual double threshold(double time) {
169      // this is the function 1 / log(x) scaled and offset
170      static double xm = 5.0 / end_time;
171      static double ym = start_threshold / (1 / log(1.2) - 1 / log(5.0 + 1.2));
172      return ym * (1 / log(xm * time + 1.2) - 1 / log(5.0 + 1.2));
173    }
[966]174  public:
[1000]175    double alpha, beta, gamma;
176    double end_time, start_time;
[1018]177    double start_threshold;
[966]178    double avg_cost;
[1000]179    double temp, ann_fact;
[1018]180    bool warmup;
181    /*! \param _end_time running time in seconds
[1000]182     *  \param _alpha parameter used to calculate the running average
183     *  \param _beta parameter used to decrease the annealing factor
184     *  \param _gamma parameter used to increase the temperature
185     */
186    AdvancedController(double _end_time, double _alpha = 0.2,
187    double _beta = 0.9, double _gamma = 1.2) : alpha(_alpha), beta(_beta),
[1023]188    gamma(_gamma), end_time(_end_time), ann_fact(0.9999), warmup(true) {}
[966]189    void init() {
[1018]190      avg_cost = base->getCurrCost();
[966]191    }
192    void acceptEvent() {
193      avg_cost = alpha * base->getCurrCost() + (1.0 - alpha) * avg_cost;
[1023]194      if (warmup) {
195        static double max_cost_diff = 0.0;
196        static int incr_cnt = 0;
197        double cost_diff = base->getCurrCost() - base->getPrevCost();
198        if (cost_diff > 0.0) {
199          incr_cnt++;
200          if (cost_diff > max_cost_diff) {
201            max_cost_diff = cost_diff;
202          }
203        }
204        if (incr_cnt >= 100) {
205          // calculate starting threshold and starting temperature
206          start_threshold = fabs(base->getBestCost() - avg_cost);
207          temp = max_cost_diff / log(0.5);
208          warmup = false;
209          timer.reset();
210        }
211      }
[966]212    }
213    void improveEvent() {
214    }
215    void rejectEvent() {
216    }
217    bool next() {
[1018]218      if (warmup) {
219        return true;
[1000]220      }
221      else {
[1018]222        double elapsed_time = timer.getRealTime();
223        if (fabs(avg_cost - base->getBestCost()) > threshold(elapsed_time)) {
224          // decrease the annealing factor
225          ann_fact *= beta;
226        }
227        else {
228          // increase the temperature
229          temp *= gamma;
230        }
231        temp *= ann_fact;
232        return elapsed_time < end_time;
[1000]233      }
[966]234    }
235    bool accept() {
[1018]236      if (warmup) {
237        // we accept eveything during the "warm up" phase
238        return true;
239      }
240      else {
241        double cost_diff = base->getPrevCost() - base->getCurrCost();
242        if (cost_diff < 0.0) {
243          return (drand48() <= exp(cost_diff / temp));
244        }
245        else {
246          return true;
247        }
248      }
[956]249    }
250  };
251
[942]252}
[918]253
254#endif
Note: See TracBrowser for help on using the repository browser.