Djinni  2.2
Annealer.h
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