Djinni
2.2
|
00001 #ifndef ANNEALER_H 00002 #define ANNEALER_H 00003 00004 #include "../Twister/Twister.h" 00005 #include "Penalties.h" 00006 #include <iosfwd> 00007 #include <string> 00008 #include <sstream> 00009 #include <iterator> 00010 #include <algorithm> 00011 #include <cmath> 00012 #include <iostream> 00013 #include <sys/types.h> 00014 00015 #if defined(__GNUC__) && (__GNUC__ >= 4) 00016 #include <tr1/memory> 00017 #define SHARED_PTR std::tr1::shared_ptr 00018 #else 00019 #include <boost/shared_ptr.hpp> 00020 #define SHARED_PTR boost::shared_ptr 00021 #endif 00022 00024 00029 template <class PenaltyFunc, class SolutionType> 00030 class Annealer { 00031 public: 00034 typedef typename PenaltyFunc::ReturnType PenaltyType; 00038 typedef unsigned long u_int32_t; 00039 00052 Annealer(PenaltyFunc& pfunc, SolutionType& sol, double multT, 00053 double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter): 00054 _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)), 00055 _bestIter(0), _iterations(0), _terminalBestIter(0), 00056 _pfunc(pfunc) 00057 { 00058 setParameters(multT, accept, tBI, minIter, maxIter); 00059 } 00060 00071 Annealer(SolutionType& solution, double multT, 00072 double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter) : 00073 _best(new SolutionType(solution)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)), 00074 _bestIter(0), _iterations(0), _terminalBestIter(0) 00075 { 00076 setParameters(multT, accept, tBI, minIter, maxIter); 00077 } 00078 00084 Annealer(PenaltyFunc& pfunc, SolutionType& sol) : 00085 _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)), 00086 _bestIter(0), _iterations(0), _terminalBestIter(0), 00087 _pfunc(pfunc) 00088 { 00089 } 00090 00094 virtual ~Annealer() {} 00095 00097 PenaltyFunc& getPenaltyFunc() { return _pfunc; } 00098 00100 const SolutionType& best() const { return *_best; } 00101 00108 const SolutionType& current() const { return *_current; } 00109 00111 void solve() 00112 { 00113 *_current = *_best; 00114 _best->setP(1000000); 00115 initializeParam(); 00116 tuneTemperature(); 00117 _iterations = 0; 00118 while ((_iterations <= _minIterations) || 00119 (_bestIter < _terminalBestIter)) 00120 { 00121 ++_iterations; 00122 for (u_int32_t count = 0; count < _maxIterations ; ++count) 00123 { 00124 _current->generateNeighbor(*_neighbor); 00125 testNeigh(); 00126 if(_current->getP() < _best->getP()) 00127 { 00128 (*_best) = (*_current); 00129 _bestIter = 1; 00130 } 00131 else if(_current->getP() == _best->getP() && _current->getF() < _best->getF()) 00132 { 00133 (*_best) = (*_current); 00134 _bestIter = 1; 00135 } 00136 } 00137 ++_bestIter; 00138 updateParam(); 00139 } 00140 } 00141 00145 std::string solution() const 00146 { 00147 std::stringstream ss; 00148 ss << (*_best); 00149 std::string result = ss.str(); 00150 return result; 00151 } 00152 00160 void setParameters(double multT, double accept, u_int32_t tBI, 00161 u_int32_t minIterations, u_int32_t maxIterations) 00162 { 00163 _multiplierT = multT; 00164 _acceptProb = accept; 00165 _terminalBestIter = tBI; 00166 _minIterations = minIterations; 00167 _maxIterations = maxIterations; 00168 } 00169 00173 void setSolutionParameters(const char* foo) 00174 { 00175 _best = SHARED_PTR<SolutionType>(new SolutionType(foo)); 00176 _current = SHARED_PTR<SolutionType>(new SolutionType(*_best)); 00177 _neighbor = SHARED_PTR<SolutionType>(new SolutionType(*_best)); 00178 } 00179 00186 std::ostream& dump(std::ostream& os) const 00187 { 00188 os << "Best solution: " << (_best->getF()) << "\t" 00189 << (_best->getP()) << std::endl; 00190 os << "Best iter: " << _bestIter << std::endl; 00191 os << "Iterations: " << _iterations << std::endl; 00192 os << "Count limit: " << _maxIterations << std::endl; 00193 os << "Minimum iterations: " << _minIterations << std::endl; 00194 os << "Sample size: " << _sampleSize << std::endl; 00195 os << "Multiplier: " << _multiplierT << std::endl; 00196 os << "Accept probability: " << _acceptProb << std::endl; 00197 os << "Terminal best iter: " << _terminalBestIter << std::endl; 00198 os << "Pressure: " << _lambda << std::endl; 00199 return os; 00200 } 00201 00205 double cost() const { return _best->getF(); } 00206 00210 double penalty() const { return _best->getP(); } 00214 u_int32_t bestIter() const { return _bestIter; } 00215 00219 u_int32_t iterations() const { return _iterations; } 00220 00224 u_int32_t maxIterations() const { return _maxIterations; } 00225 00229 u_int32_t minIterations() const { return _minIterations; } 00233 double multiplier() const { return _multiplierT; } 00234 00238 double probability() const { return _acceptProb; } 00242 PenaltyType getLambda() const { return _lambda; } 00243 00244 protected: 00246 void initializeParam() 00247 { 00248 double lambda1 = 0, sum = 0; 00249 for (u_int32_t j = 0 ; j < _sampleSize - 1 ; j += 2) 00250 { 00251 _current->randomize(); 00252 _current->compute(); 00253 _current->generateNeighbor(*_neighbor); 00254 double u = (_current->getF() + lambda1 * _current->getP()) - 00255 (_neighbor->getF() + lambda1 * _neighbor->getP()); 00256 sum += u > 0 ? u : (-1 * u); 00257 } 00258 00259 sum /= _sampleSize; 00260 _currentT = ((-1 * sum) / log(_acceptProb)); 00261 } 00262 00264 void tuneTemperature() 00265 { 00266 int acceptedWorse, uphill; 00267 do 00268 { 00269 acceptedWorse = uphill = 0; 00270 for (u_int32_t count = 0 ; count < _maxIterations ; count++) 00271 { 00272 _current->generateNeighbor(*_neighbor); 00273 double delta, u; 00274 delta = (_neighbor->getF() + _lambda * _neighbor->getP()) - 00275 (_current->getF() + _lambda*_current->getP()); 00276 if (delta < 0) 00277 _current.swap(_neighbor); 00278 else 00279 { 00280 uphill++; 00281 u = 0 - delta / _currentT; 00282 if (Twister::generateDouble() < exp(u)) 00283 { 00284 _current.swap(_neighbor); 00285 acceptedWorse++; 00286 } 00287 } 00288 if ( _current->getP() < _best->getP() ) 00289 *_best = *_current; 00290 else if ( (_current->getP() == _best->getP() ) && 00291 (_current->getF() < _best->getF()) ) 00292 *_best = *_current; 00293 } 00294 if ((static_cast<double>(acceptedWorse) / 00295 static_cast<double>(uphill)) < _acceptProb) 00296 _currentT = 1.5 * _currentT; 00297 } while ((static_cast<double>(acceptedWorse) / 00298 static_cast<double>(uphill)) < _acceptProb); 00299 } 00300 00302 void testNeigh() 00303 { 00304 double delta, u; 00305 delta = (_neighbor->getF() + _lambda * _neighbor->getP()) - 00306 (_current->getF() + _lambda * _current->getP()); 00307 if (delta < 0) 00308 _current.swap(_neighbor); 00309 else 00310 { 00311 u = 0 - delta / _currentT; 00312 if (Twister::generateDouble() < exp(u)) 00313 _current.swap(_neighbor); 00314 } 00315 } 00316 00318 void updateParam() 00319 { 00320 _currentT = _multiplierT * _currentT; 00321 _lambda = _pfunc(_iterations); 00322 } 00323 00324 SHARED_PTR<SolutionType> _best, _current, _neighbor; 00325 00326 u_int32_t _bestIter, _iterations, _maxIterations, _minIterations, 00327 _terminalBestIter; 00328 static const int _sampleSize = 10000; 00329 double _multiplierT, _acceptProb, _currentT; 00330 PenaltyFunc _pfunc; 00331 PenaltyType _lambda; 00332 }; 00333 00334 00335 00337 00342 template <class SolutionType> 00343 class Annealer<Compressed, SolutionType> { 00344 public: 00345 00346 typedef Compressed PenaltyFunc; 00347 typedef typename PenaltyFunc::ReturnType PenaltyType; 00348 typedef unsigned long u_int32_t; 00349 00362 Annealer(PenaltyFunc& pfunc, SolutionType& sol, double multT, 00363 double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter): 00364 _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)), 00365 _bestIter(0), _iterations(0), _terminalBestIter(0), 00366 _pfunc(pfunc) 00367 { 00368 setParameters(multT, accept, tBI, minIter, maxIter); 00369 } 00370 00381 Annealer(SolutionType& solution, double multT, 00382 double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter) : 00383 _best(new SolutionType(solution)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)), 00384 _bestIter(0), _iterations(0), _terminalBestIter(0) 00385 { 00386 setParameters(multT, accept, tBI, minIter, maxIter); 00387 } 00388 00394 Annealer(PenaltyFunc& pfunc, SolutionType& sol) : 00395 _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)), 00396 _bestIter(0), _iterations(0), _terminalBestIter(0), 00397 _pfunc(pfunc) 00398 { 00399 } 00400 00404 virtual ~Annealer() {} 00405 00407 PenaltyFunc& getPenaltyFunc() { return _pfunc; } 00408 00410 const SolutionType& best() const { return *_best; } 00411 00418 const SolutionType& current() const { return *_current; } 00419 00421 void solve() 00422 { 00423 *_current = *_best; 00424 _best->setP(1000000); 00425 initializeParam(); 00426 tuneTemperature(); 00427 _iterations = 0; 00428 while ((_iterations <= _minIterations) || 00429 (_bestIter < _terminalBestIter)) 00430 { 00431 ++_iterations; 00432 for (u_int32_t count = 0; count < _maxIterations ; ++count) 00433 { 00434 _current->generateNeighbor(*_neighbor); 00435 testNeigh(); 00436 if(_current->getP() < _best->getP()) 00437 { 00438 (*_best) = (*_current); 00439 _bestIter = 1; 00440 } 00441 else if(_current->getP() == _best->getP() && _current->getF() < _best->getF()) 00442 { 00443 (*_best) = (*_current); 00444 _bestIter = 1; 00445 } 00446 } 00447 ++_bestIter; 00448 updateParam(); 00449 } 00450 } 00451 00455 std::string solution() const 00456 { 00457 std::stringstream ss; 00458 ss << (*_best); 00459 std::string result = ss.str(); 00460 return result; 00461 } 00462 00470 void setParameters(double multT, double accept, u_int32_t tBI, 00471 u_int32_t minIterations, u_int32_t maxIterations) 00472 { 00473 _multiplierT = multT; 00474 _acceptProb = accept; 00475 _terminalBestIter = tBI; 00476 _minIterations = minIterations; 00477 _maxIterations = maxIterations; 00478 } 00479 00483 void setSolutionParameters(const char* foo) 00484 { 00485 _best = SHARED_PTR<SolutionType>(new SolutionType(foo)); 00486 _current = SHARED_PTR<SolutionType>(new SolutionType(*_best)); 00487 _neighbor = SHARED_PTR<SolutionType>(new SolutionType(*_best)); 00488 } 00489 00496 std::ostream& dump(std::ostream& os) const 00497 { 00498 os << "Best solution: " << (_best->getF()) << "\t" 00499 << (_best->getP()) << std::endl; 00500 os << "Best iter: " << _bestIter << std::endl; 00501 os << "Iterations: " << _iterations << std::endl; 00502 os << "Count limit: " << _maxIterations << std::endl; 00503 os << "Minimum iterations: " << _minIterations << std::endl; 00504 os << "Sample size: " << _sampleSize << std::endl; 00505 os << "Multiplier: " << _multiplierT << std::endl; 00506 os << "Accept probability: " << _acceptProb << std::endl; 00507 os << "Terminal best iter: " << _terminalBestIter << std::endl; 00508 os << "Pressure: " << _lambda << std::endl; 00509 return os; 00510 } 00511 00515 double cost() const { return _best->getF(); } 00516 00520 double penalty() const { return _best->getP(); } 00524 u_int32_t bestIter() const { return _bestIter; } 00525 00529 u_int32_t iterations() const { return _iterations; } 00530 00534 u_int32_t maxIterations() const { return _maxIterations; } 00535 00539 u_int32_t minIterations() const { return _minIterations; } 00543 double multiplier() const { return _multiplierT; } 00544 00548 double probability() const { return _acceptProb; } 00552 PenaltyType getLambda() const { return _lambda; } 00553 00554 protected: 00556 void initializeParam() 00557 { 00558 double lambda1 = 0, sum = 0; 00559 double lambda0 = 0, cap = 0; 00560 double solutionScalePressure = _pfunc.getCapPercentage() / 00561 (1 - _pfunc.getCapPercentage()); 00562 for (u_int32_t j = 0 ; j < _sampleSize - 1 ; j += 2) 00563 { 00564 _current->randomize(); 00565 _current->compute(); 00566 _current->generateNeighbor(*_neighbor); 00567 if ( _current->getP() > 0 ) 00568 lambda0 = (_current->getF() / _current->getP()) * 00569 solutionScalePressure; 00570 if (lambda0 > cap) 00571 cap = lambda0; 00572 if (_neighbor->getP() > 0) 00573 lambda0 = (_neighbor->getF() / _neighbor->getP()) * 00574 solutionScalePressure; 00575 if (lambda0 > cap) 00576 cap = lambda0; 00577 double u = (_current->getF() + lambda1 * _current->getP()) - 00578 (_neighbor->getF() + lambda1 * _neighbor->getP()); 00579 sum += u > 0 ? u : (-1 * u); 00580 } 00581 00582 _pfunc.setPressureCap(cap); 00583 sum /= _sampleSize; 00584 _currentT = ((-1 * sum) / log(_acceptProb)); 00585 } 00586 00588 void tuneTemperature() 00589 { 00590 int acceptedWorse, uphill; 00591 u_int32_t iterations = 0; 00592 do 00593 { 00594 ++iterations; 00595 acceptedWorse = uphill = 0; 00596 for (u_int32_t count = 0 ; count < _maxIterations ; count++) 00597 { 00598 _current->generateNeighbor(*_neighbor); 00599 double delta, u; 00600 delta = (_neighbor->getF() + _lambda * _neighbor->getP()) - 00601 (_current->getF() + _lambda*_current->getP()); 00602 if (delta < 0) 00603 _current.swap(_neighbor); 00604 else 00605 { 00606 uphill++; 00607 u = 0 - delta / _currentT; 00608 double generated = Twister::generateDouble(); 00609 if (generated < exp(u)) 00610 { 00611 _current.swap(_neighbor); 00612 acceptedWorse++; 00613 } 00614 } 00615 if ( _current->getP() < _best->getP() ) 00616 *_best = *_current; 00617 else if ( (_current->getP() == _best->getP() ) && 00618 (_current->getF() < _best->getF()) ) 00619 *_best = *_current; 00620 } 00621 if ((static_cast<double>(acceptedWorse) / 00622 static_cast<double>(uphill)) < _acceptProb) 00623 _currentT = 1.5 * _currentT; 00624 } while ((static_cast<double>(acceptedWorse) / 00625 static_cast<double>(uphill)) < _acceptProb); 00626 } 00627 00629 void testNeigh() 00630 { 00631 double delta, u; 00632 delta = (_neighbor->getF() + _lambda * _neighbor->getP()) - 00633 (_current->getF() + _lambda * _current->getP()); 00634 if (delta < 0) 00635 _current.swap(_neighbor); 00636 else 00637 { 00638 u = 0 - delta / _currentT; 00639 if (Twister::generateDouble() < exp(u)) 00640 _current.swap(_neighbor); 00641 } 00642 } 00643 00645 void updateParam() 00646 { 00647 _currentT = _multiplierT * _currentT; 00648 _lambda = _pfunc(_iterations); 00649 } 00650 00651 SHARED_PTR<SolutionType> _best, _current, _neighbor; 00652 00653 u_int32_t _bestIter, _iterations, _maxIterations, _minIterations, 00654 _terminalBestIter; 00655 static const int _sampleSize = 10000; 00656 double _multiplierT, _acceptProb, _currentT; 00657 PenaltyFunc _pfunc; 00658 PenaltyType _lambda; 00659 }; 00660 00661 00667 template <class PenaltyFunc, class SolutionType> 00668 std::ostream& operator<<(std::ostream& os, 00669 const Annealer<PenaltyFunc, SolutionType>& engine) 00670 { 00671 return engine.dump(os); 00672 } 00673 00674 #endif 00675