[942] | 1 | #ifndef LEMON_SIMANN_H |
---|
| 2 | #define LEMON_SIMANN_H |
---|
[918] | 3 | |
---|
[1142] | 4 | /// \ingroup experimental |
---|
| 5 | /// \file |
---|
| 6 | /// \brief Simulated annealing framework. |
---|
| 7 | /// \author Akos Ladanyi |
---|
| 8 | |
---|
[966] | 9 | #include <cstdlib> |
---|
| 10 | #include <cmath> |
---|
[1018] | 11 | #include <lemon/time_measure.h> |
---|
[966] | 12 | |
---|
[942] | 13 | namespace lemon { |
---|
[918] | 14 | |
---|
[1142] | 15 | /// \addtogroup experimental |
---|
| 16 | /// @{ |
---|
| 17 | |
---|
[942] | 18 | const double INFTY = 1e24; |
---|
[918] | 19 | |
---|
[1142] | 20 | /*! \brief Simulated annealing base class. */ |
---|
[942] | 21 | class SimAnnBase { |
---|
[918] | 22 | public: |
---|
[942] | 23 | class Controller; |
---|
| 24 | private: |
---|
[1142] | 25 | /*! Pointer to the controller. */ |
---|
[942] | 26 | Controller *controller; |
---|
| 27 | protected: |
---|
[1142] | 28 | /*! \brief Cost of the current solution. */ |
---|
[942] | 29 | double curr_cost; |
---|
[1142] | 30 | /*! \brief Cost of the best solution. */ |
---|
[1023] | 31 | double best_cost; |
---|
[1142] | 32 | /*! \brief Cost of the previous solution. */ |
---|
[942] | 33 | double prev_cost; |
---|
[1142] | 34 | /*! \brief Cost of the solution preceding the previous one. */ |
---|
[1023] | 35 | double prev_prev_cost; |
---|
[918] | 36 | |
---|
[1142] | 37 | /*! \brief Step to a neighbouring state. */ |
---|
[1150] | 38 | virtual void mutate() = 0; |
---|
[1142] | 39 | /*! \brief Reverts the last mutate(). */ |
---|
[1150] | 40 | virtual void revert() = 0; |
---|
[1142] | 41 | /*! \brief Saves the current solution as the best one. */ |
---|
[1150] | 42 | virtual void saveAsBest() = 0; |
---|
[942] | 43 | public: |
---|
[1142] | 44 | /*! \brief Constructor. */ |
---|
[942] | 45 | SimAnnBase() { |
---|
[1023] | 46 | best_cost = prev_cost = prev_prev_cost = INFTY; |
---|
[942] | 47 | } |
---|
[1142] | 48 | /*! \brief Sets the controller class to use. */ |
---|
[957] | 49 | void setController(Controller &_controller) { |
---|
| 50 | controller = &_controller; |
---|
| 51 | controller->setBase(this); |
---|
| 52 | } |
---|
[1142] | 53 | /*! \brief Returns the cost of the current solution. */ |
---|
[1018] | 54 | double getCurrCost() const { return curr_cost; } |
---|
[1142] | 55 | /*! \brief Returns the cost of the previous solution. */ |
---|
[1018] | 56 | double getPrevCost() const { return prev_cost; } |
---|
[1142] | 57 | /*! \brief Returns the cost of the best solution. */ |
---|
[1018] | 58 | double getBestCost() const { return best_cost; } |
---|
[1142] | 59 | /*! \brief Starts the annealing process. */ |
---|
[942] | 60 | void run() { |
---|
[966] | 61 | controller->init(); |
---|
[1018] | 62 | do { |
---|
[1150] | 63 | curr_cost=mutate(); |
---|
[957] | 64 | if (controller->accept()) { |
---|
[942] | 65 | controller->acceptEvent(); |
---|
| 66 | if (curr_cost < best_cost) { |
---|
| 67 | saveAsBest(); |
---|
| 68 | controller->improveEvent(); |
---|
| 69 | } |
---|
| 70 | } |
---|
| 71 | else { |
---|
| 72 | revert(); |
---|
| 73 | controller->rejectEvent(); |
---|
| 74 | } |
---|
[1018] | 75 | } while (controller->next()); |
---|
[918] | 76 | } |
---|
| 77 | |
---|
[1000] | 78 | /*! \brief A base class for controllers. */ |
---|
[942] | 79 | class Controller { |
---|
| 80 | public: |
---|
[1142] | 81 | /*! \brief Pointer to the simulated annealing base class. */ |
---|
[957] | 82 | SimAnnBase *base; |
---|
[1142] | 83 | /*! \brief Initializes the controller. */ |
---|
[966] | 84 | virtual void init() {} |
---|
[1000] | 85 | /*! \brief This is called when a neighbouring state gets accepted. */ |
---|
[942] | 86 | virtual void acceptEvent() {} |
---|
[1000] | 87 | /*! \brief This is called when the accepted neighbouring state's cost is |
---|
| 88 | * less than the best found one's. |
---|
| 89 | */ |
---|
[942] | 90 | virtual void improveEvent() {} |
---|
[1000] | 91 | /*! \brief This is called when a neighbouring state gets rejected. */ |
---|
[942] | 92 | virtual void rejectEvent() {} |
---|
[1142] | 93 | /*! \brief Sets the simulated annealing base class to use. */ |
---|
[957] | 94 | virtual void setBase(SimAnnBase *_base) { base = _base; } |
---|
[1142] | 95 | /*! \brief Decides whether to continue the annealing process or not. */ |
---|
[942] | 96 | virtual bool next() = 0; |
---|
[1142] | 97 | /*! \brief Decides whether to accept the current solution or not. */ |
---|
[957] | 98 | virtual bool accept() = 0; |
---|
[942] | 99 | }; |
---|
| 100 | }; |
---|
[918] | 101 | |
---|
[1142] | 102 | /*! \brief Simulated annealing class. */ |
---|
[942] | 103 | template <typename E> |
---|
| 104 | class SimAnn : public SimAnnBase { |
---|
| 105 | private: |
---|
[1142] | 106 | /*! \brief Pointer to the current entity. */ |
---|
[942] | 107 | E *curr_ent; |
---|
[1142] | 108 | /*! \brief Pointer to the best entity. */ |
---|
[942] | 109 | E *best_ent; |
---|
| 110 | public: |
---|
[1142] | 111 | /*! \brief Constructor. */ |
---|
[957] | 112 | SimAnn() : SimAnnBase() {} |
---|
[1142] | 113 | /*! \brief Sets the initial entity. */ |
---|
[957] | 114 | void setEntity(E &ent) { |
---|
| 115 | curr_ent = new E(ent); |
---|
| 116 | best_ent = new E(ent); |
---|
[1023] | 117 | curr_cost = curr_ent->getCost(); |
---|
[942] | 118 | } |
---|
[1142] | 119 | /*! \brief Returns the best found entity. */ |
---|
[942] | 120 | E getBestEntity() { return *best_ent; } |
---|
[1142] | 121 | /*! \brief Step to a neighbouring state. */ |
---|
[942] | 122 | void mutate() { |
---|
[1023] | 123 | prev_prev_cost = prev_cost; |
---|
[1018] | 124 | prev_cost = curr_cost; |
---|
[1023] | 125 | curr_ent->mutate(); |
---|
| 126 | curr_cost = curr_ent->getCost(); |
---|
[942] | 127 | } |
---|
[1142] | 128 | /*! \brief Reverts the last mutate(). */ |
---|
[942] | 129 | void revert() { |
---|
| 130 | curr_ent->revert(); |
---|
[1018] | 131 | curr_cost = prev_cost; |
---|
[1023] | 132 | prev_cost = prev_prev_cost; |
---|
[942] | 133 | } |
---|
[1142] | 134 | /*! \brief Saves the current solution as the best one. */ |
---|
[942] | 135 | void saveAsBest() { |
---|
[1096] | 136 | delete(best_ent); |
---|
| 137 | best_ent = new E(*curr_ent); |
---|
[942] | 138 | best_cost = curr_cost; |
---|
| 139 | } |
---|
| 140 | }; |
---|
| 141 | |
---|
[1142] | 142 | /*! \brief Skeleton of an entity class. */ |
---|
[956] | 143 | class EntitySkeleton { |
---|
[942] | 144 | public: |
---|
[1142] | 145 | /*! \brief Returns the cost of the entity. */ |
---|
[1023] | 146 | double getCost() { return 0.0; } |
---|
| 147 | /*! \brief Makes a minor change to the entity. */ |
---|
| 148 | void mutate() {} |
---|
[966] | 149 | /*! \brief Restores the entity to its previous state i.e. reverts the |
---|
[1142] | 150 | * effects of the last mutate(). |
---|
[966] | 151 | */ |
---|
[942] | 152 | void revert() {} |
---|
| 153 | }; |
---|
| 154 | |
---|
[1142] | 155 | /*! \brief A simple controller for the simulated annealing class. */ |
---|
[956] | 156 | class SimpleController : public SimAnnBase::Controller { |
---|
| 157 | public: |
---|
[1142] | 158 | /*! \brief Number of iterations. */ |
---|
| 159 | long iter; |
---|
| 160 | /*! \brief Number of iterations which did not improve the solution since |
---|
| 161 | * the last improvement. */ |
---|
| 162 | long last_impr; |
---|
| 163 | /*! \brief Maximum number of iterations. */ |
---|
| 164 | long max_iter; |
---|
| 165 | /*! \brief Maximum number of iterations which do not improve the |
---|
| 166 | * solution. */ |
---|
| 167 | long max_no_impr; |
---|
| 168 | /*! \brief Temperature. */ |
---|
| 169 | double temp; |
---|
| 170 | /*! \brief Annealing factor. */ |
---|
| 171 | double ann_fact; |
---|
| 172 | /*! \brief Constructor. |
---|
| 173 | * \param _max_iter maximum number of iterations |
---|
[1000] | 174 | * \param _max_no_impr maximum number of consecutive iterations which do |
---|
| 175 | * not yield a better solution |
---|
| 176 | * \param _temp initial temperature |
---|
| 177 | * \param _ann_fact annealing factor |
---|
| 178 | */ |
---|
| 179 | SimpleController(long _max_iter = 500000, long _max_no_impr = 20000, |
---|
[1096] | 180 | double _temp = 1000.0, double _ann_fact = 0.9999) : iter(0), last_impr(0), |
---|
[1000] | 181 | max_iter(_max_iter), max_no_impr(_max_no_impr), temp(_temp), |
---|
| 182 | ann_fact(_ann_fact) {} |
---|
[1145] | 183 | /*! \brief This is called when a neighbouring state gets accepted. */ |
---|
[956] | 184 | void acceptEvent() { |
---|
| 185 | iter++; |
---|
| 186 | } |
---|
[1142] | 187 | /*! \brief This is called when the accepted neighbouring state's cost is |
---|
| 188 | * less than the best found one's. |
---|
| 189 | */ |
---|
[956] | 190 | void improveEvent() { |
---|
| 191 | last_impr = iter; |
---|
| 192 | } |
---|
[1142] | 193 | /*! \brief This is called when a neighbouring state gets rejected. */ |
---|
[956] | 194 | void rejectEvent() { |
---|
| 195 | iter++; |
---|
| 196 | } |
---|
[1142] | 197 | /*! \brief Decides whether to continue the annealing process or not. Also |
---|
| 198 | * decreases the temperature. */ |
---|
[956] | 199 | bool next() { |
---|
[1000] | 200 | temp *= ann_fact; |
---|
[956] | 201 | bool quit = (iter > max_iter) || (iter - last_impr > max_no_impr); |
---|
| 202 | return !quit; |
---|
| 203 | } |
---|
[1142] | 204 | /*! \brief Decides whether to accept the current solution or not. */ |
---|
[957] | 205 | bool accept() { |
---|
[1018] | 206 | double cost_diff = base->getPrevCost() - base->getCurrCost(); |
---|
| 207 | if (cost_diff < 0.0) { |
---|
[1096] | 208 | bool ret = drand48() <= exp(cost_diff / temp); |
---|
| 209 | return ret; |
---|
[1018] | 210 | } |
---|
| 211 | else { |
---|
| 212 | return true; |
---|
| 213 | } |
---|
[966] | 214 | } |
---|
| 215 | }; |
---|
| 216 | |
---|
| 217 | /*! \brief A controller with preset running time for the simulated annealing |
---|
| 218 | * class. |
---|
[1145] | 219 | * |
---|
| 220 | * With this controller you can set the running time of the annealing |
---|
| 221 | * process in advance. It works the following way: the controller measures |
---|
| 222 | * a kind of divergence. The divergence is the difference of the average |
---|
| 223 | * cost of the recently found solutions the cost of the best found one. In |
---|
| 224 | * case this divergence is greater than a given threshold, then we decrease |
---|
| 225 | * the annealing factor, that is we cool the system faster. In case the |
---|
| 226 | * divergence is lower than the threshold, then we increase the temperature. |
---|
| 227 | * The threshold is a function of the elapsed time which reaches zero at the |
---|
| 228 | * desired end time. |
---|
[966] | 229 | */ |
---|
| 230 | class AdvancedController : public SimAnnBase::Controller { |
---|
| 231 | private: |
---|
[1018] | 232 | Timer timer; |
---|
[1000] | 233 | /*! \param time the elapsed time in seconds */ |
---|
[1018] | 234 | virtual double threshold(double time) { |
---|
[1096] | 235 | return (-1.0) * start_threshold / end_time * time + start_threshold; |
---|
[1018] | 236 | } |
---|
[966] | 237 | public: |
---|
[1142] | 238 | double alpha; |
---|
| 239 | double beta; |
---|
| 240 | double gamma; |
---|
[1145] | 241 | /*! \brief The time at the end of the algorithm. */ |
---|
[1142] | 242 | double end_time; |
---|
[1145] | 243 | /*! \brief The time at the start of the algorithm. */ |
---|
[1142] | 244 | double start_time; |
---|
[1145] | 245 | /*! \brief Starting threshold. */ |
---|
[1018] | 246 | double start_threshold; |
---|
[1145] | 247 | /*! \brief Average cost of recent solutions. */ |
---|
[966] | 248 | double avg_cost; |
---|
[1145] | 249 | /*! \brief Temperature. */ |
---|
[1142] | 250 | double temp; |
---|
[1145] | 251 | /*! \brief Annealing factor. */ |
---|
[1142] | 252 | double ann_fact; |
---|
[1145] | 253 | /*! \brief Initial annealing factor. */ |
---|
| 254 | double init_ann_fact; |
---|
[1018] | 255 | bool warmup; |
---|
[1142] | 256 | /*! \brief Constructor. |
---|
| 257 | * \param _end_time running time in seconds |
---|
[1000] | 258 | * \param _alpha parameter used to calculate the running average |
---|
| 259 | * \param _beta parameter used to decrease the annealing factor |
---|
| 260 | * \param _gamma parameter used to increase the temperature |
---|
[1145] | 261 | * \param _ann_fact initial annealing factor |
---|
[1000] | 262 | */ |
---|
| 263 | AdvancedController(double _end_time, double _alpha = 0.2, |
---|
[1145] | 264 | double _beta = 0.9, double _gamma = 1.6, double _ann_fact = 0.9999) : |
---|
| 265 | alpha(_alpha), beta(_beta), gamma(_gamma), end_time(_end_time), |
---|
| 266 | ann_fact(_ann_fact), init_ann_fact(_ann_fact), warmup(true) {} |
---|
[966] | 267 | void init() { |
---|
[1018] | 268 | avg_cost = base->getCurrCost(); |
---|
[966] | 269 | } |
---|
[1142] | 270 | /*! \brief This is called when a neighbouring state gets accepted. */ |
---|
[966] | 271 | void acceptEvent() { |
---|
| 272 | avg_cost = alpha * base->getCurrCost() + (1.0 - alpha) * avg_cost; |
---|
[1023] | 273 | if (warmup) { |
---|
[1096] | 274 | static int cnt = 0; |
---|
| 275 | cnt++; |
---|
| 276 | if (cnt >= 100) { |
---|
[1023] | 277 | // calculate starting threshold and starting temperature |
---|
[1096] | 278 | start_threshold = 5.0 * fabs(base->getBestCost() - avg_cost); |
---|
| 279 | temp = 10000.0; |
---|
[1023] | 280 | warmup = false; |
---|
| 281 | timer.reset(); |
---|
| 282 | } |
---|
| 283 | } |
---|
[966] | 284 | } |
---|
[1142] | 285 | /*! \brief Decides whether to continue the annealing process or not. */ |
---|
[966] | 286 | bool next() { |
---|
[1018] | 287 | if (warmup) { |
---|
| 288 | return true; |
---|
[1000] | 289 | } |
---|
| 290 | else { |
---|
[1018] | 291 | double elapsed_time = timer.getRealTime(); |
---|
| 292 | if (fabs(avg_cost - base->getBestCost()) > threshold(elapsed_time)) { |
---|
| 293 | // decrease the annealing factor |
---|
| 294 | ann_fact *= beta; |
---|
| 295 | } |
---|
| 296 | else { |
---|
| 297 | // increase the temperature |
---|
| 298 | temp *= gamma; |
---|
[1145] | 299 | // reset the annealing factor |
---|
| 300 | ann_fact = init_ann_fact; |
---|
[1018] | 301 | } |
---|
| 302 | temp *= ann_fact; |
---|
| 303 | return elapsed_time < end_time; |
---|
[1000] | 304 | } |
---|
[966] | 305 | } |
---|
[1142] | 306 | /*! \brief Decides whether to accept the current solution or not. */ |
---|
[966] | 307 | bool accept() { |
---|
[1018] | 308 | if (warmup) { |
---|
| 309 | // we accept eveything during the "warm up" phase |
---|
| 310 | return true; |
---|
| 311 | } |
---|
| 312 | else { |
---|
| 313 | double cost_diff = base->getPrevCost() - base->getCurrCost(); |
---|
| 314 | if (cost_diff < 0.0) { |
---|
| 315 | return (drand48() <= exp(cost_diff / temp)); |
---|
| 316 | } |
---|
| 317 | else { |
---|
| 318 | return true; |
---|
| 319 | } |
---|
| 320 | } |
---|
[956] | 321 | } |
---|
| 322 | }; |
---|
| 323 | |
---|
[1142] | 324 | /// @} |
---|
| 325 | |
---|
[942] | 326 | } |
---|
[918] | 327 | |
---|
| 328 | #endif |
---|