Annealer.h

#ifndef ANNEALER_H
#define ANNEALER_H

#include "../Twister/Twister.h"
#include "Penalties.h"
#include <iosfwd>
#include <string>
#include <sstream>
#include <iterator>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <sys/types.h>

#if defined(__GNUC__) && (__GNUC__ >= 4)
#include <tr1/memory>
#define SHARED_PTR std::tr1::shared_ptr
#else
#include <boost/shared_ptr.hpp>
#define SHARED_PTR boost::shared_ptr
#endif


template <class PenaltyFunc, class SolutionType>
class Annealer {
public:
    typedef typename PenaltyFunc::ReturnType PenaltyType;
    typedef unsigned long u_int32_t;

    Annealer(PenaltyFunc& pfunc, SolutionType& sol, double multT, 
        double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter):
        _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)),
        _bestIter(0), _iterations(0), _terminalBestIter(0),
        _pfunc(pfunc)
    {
        setParameters(multT, accept, tBI, minIter, maxIter);
    }

    Annealer(SolutionType& solution, double multT, 
        double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter) :
        _best(new SolutionType(solution)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)),
        _bestIter(0), _iterations(0), _terminalBestIter(0)
    {
        setParameters(multT, accept, tBI, minIter, maxIter);
    }
    
    Annealer(PenaltyFunc& pfunc, SolutionType& sol) :
        _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)),
        _bestIter(0), _iterations(0), _terminalBestIter(0),
        _pfunc(pfunc)
    {
    }

    virtual ~Annealer() {}
    
    PenaltyFunc& getPenaltyFunc() { return _pfunc; }

    const SolutionType& best() const { return *_best; }
    
    const SolutionType& current() const { return *_current; }
    
    void solve()
    {
        *_current = *_best;
        _best->setP(1000000);
        initializeParam();
        tuneTemperature();
        _iterations = 0;
        while ((_iterations <= _minIterations) || 
        (_bestIter < _terminalBestIter))
        {
            ++_iterations;
            for (u_int32_t count = 0; count < _maxIterations ; ++count)
            {
                _current->generateNeighbor(*_neighbor);
                testNeigh();
                if(_current->getP() < _best->getP())
                {
                        (*_best) = (*_current);
                        _bestIter = 1;
                }
                else if(_current->getP() == _best->getP() && _current->getF() < _best->getF())
                {
                        (*_best) = (*_current);
                        _bestIter = 1;
                }
            }
            ++_bestIter;
            updateParam();
        }
    }

    std::string solution() const
    {
        std::stringstream ss;
        ss << (*_best);
        std::string result = ss.str();
        return result;
    }
    
    void setParameters(double multT, double accept, u_int32_t tBI,
        u_int32_t minIterations, u_int32_t maxIterations)
    {
        _multiplierT = multT;
        _acceptProb = accept;
        _terminalBestIter = tBI;
        _minIterations = minIterations;
        _maxIterations = maxIterations;
    }
    
    void setSolutionParameters(const char* foo)
    {
        _best = SHARED_PTR<SolutionType>(new SolutionType(foo));
        _current = SHARED_PTR<SolutionType>(new SolutionType(*_best));
        _neighbor = SHARED_PTR<SolutionType>(new SolutionType(*_best));
    }
    
    std::ostream& dump(std::ostream& os) const
    {
        os << "Best solution:      " << (_best->getF()) << "\t" 
                   << (_best->getP()) << std::endl;
        os << "Best iter:          " << _bestIter << std::endl;
        os << "Iterations:         " << _iterations << std::endl;
        os << "Count limit:        " << _maxIterations << std::endl;
        os << "Minimum iterations: " << _minIterations << std::endl;
        os << "Sample size:        " << _sampleSize << std::endl;
        os << "Multiplier:         " << _multiplierT << std::endl;
        os << "Accept probability: " << _acceptProb << std::endl;
        os << "Terminal best iter: " << _terminalBestIter << std::endl;
        os << "Pressure:           " << _lambda << std::endl;
        return os;
    }
    
    double cost() const { return _best->getF(); }
    
    double penalty() const { return _best->getP(); }
    u_int32_t bestIter() const { return _bestIter; }

    u_int32_t iterations() const { return _iterations; }
    
    u_int32_t maxIterations() const { return _maxIterations; }
    
    u_int32_t minIterations() const { return _minIterations; }
    double multiplier() const { return _multiplierT; }
    
    double probability() const { return _acceptProb; }
    PenaltyType getLambda() const { return _lambda; }
    
protected:
    void initializeParam()
    {
        double lambda1 = 0, sum = 0;
        for (u_int32_t j = 0 ; j < _sampleSize - 1 ; j += 2)
        {
            _current->randomize();
            _current->compute();
            _current->generateNeighbor(*_neighbor);
            double u = (_current->getF() + lambda1 * _current->getP()) - 
                (_neighbor->getF() + lambda1 * _neighbor->getP());
               sum += u > 0 ? u : (-1 * u);
        }
        
        sum /= _sampleSize;
        _currentT = ((-1 * sum) / log(_acceptProb));
    }
    
    void tuneTemperature()
    {
        int acceptedWorse, uphill;
        do
        {
            acceptedWorse = uphill = 0;
            for (u_int32_t count = 0 ; count < _maxIterations ; count++)
            {
                _current->generateNeighbor(*_neighbor);
                double delta, u;
                delta = (_neighbor->getF() + _lambda * _neighbor->getP()) -
                        (_current->getF() + _lambda*_current->getP());
                if (delta < 0)
                    _current.swap(_neighbor);
                else
                {
                    uphill++;
                    u = 0 - delta / _currentT;
                    if (Twister::generateDouble() < exp(u))
                    {
                        _current.swap(_neighbor);
                        acceptedWorse++;
                    }
                }
                if ( _current->getP() < _best->getP() )
                    *_best = *_current;
                else if ( (_current->getP() == _best->getP() ) && 
                                                  (_current->getF() < _best->getF()) )
                    *_best = *_current;
            }
            if ((static_cast<double>(acceptedWorse) / 
                                 static_cast<double>(uphill)) < _acceptProb)
                _currentT = 1.5 * _currentT;
        } while ((static_cast<double>(acceptedWorse) / 
                                  static_cast<double>(uphill)) < _acceptProb);
    }

    void testNeigh()
    {
        double delta, u;
        delta = (_neighbor->getF() + _lambda * _neighbor->getP()) - 
                (_current->getF() + _lambda * _current->getP());
        if (delta < 0)
            _current.swap(_neighbor);
        else
        {
            u = 0 - delta / _currentT;
            if (Twister::generateDouble() < exp(u))
                _current.swap(_neighbor);
        }
    }
    
    void updateParam()
    {
        _currentT = _multiplierT * _currentT;
        _lambda = _pfunc(_iterations);
    }
        
    SHARED_PTR<SolutionType> _best, _current, _neighbor;
    
    u_int32_t _bestIter, _iterations, _maxIterations, _minIterations,
        _terminalBestIter;
    static const int _sampleSize = 10000;
    double _multiplierT, _acceptProb, _currentT;
    PenaltyFunc _pfunc;
    PenaltyType _lambda;
};




template <class SolutionType>
class Annealer<Compressed, SolutionType> {
public:

    typedef Compressed PenaltyFunc;
    typedef typename PenaltyFunc::ReturnType PenaltyType;
    typedef unsigned long u_int32_t;

    Annealer(PenaltyFunc& pfunc, SolutionType& sol, double multT, 
        double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter):
        _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)),
        _bestIter(0), _iterations(0), _terminalBestIter(0),
        _pfunc(pfunc)
    {
        setParameters(multT, accept, tBI, minIter, maxIter);
    }

    Annealer(SolutionType& solution, double multT, 
        double accept, u_int32_t tBI, u_int32_t minIter, u_int32_t maxIter) :
        _best(new SolutionType(solution)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)),
        _bestIter(0), _iterations(0), _terminalBestIter(0)
    {
        setParameters(multT, accept, tBI, minIter, maxIter);
    }
    
    Annealer(PenaltyFunc& pfunc, SolutionType& sol) :
        _best(new SolutionType(sol)), _current(new SolutionType(*_best)), _neighbor(new SolutionType(*_best)),
        _bestIter(0), _iterations(0), _terminalBestIter(0),
        _pfunc(pfunc)
    {
    }

    virtual ~Annealer() {}
    
    PenaltyFunc& getPenaltyFunc() { return _pfunc; }

    const SolutionType& best() const { return *_best; }
    
    const SolutionType& current() const { return *_current; }
    
    void solve()
    {
        *_current = *_best;
        _best->setP(1000000);
        initializeParam();
        tuneTemperature();
        _iterations = 0;
        while ((_iterations <= _minIterations) || 
        (_bestIter < _terminalBestIter))
        {
            ++_iterations;
            for (u_int32_t count = 0; count < _maxIterations ; ++count)
            {
                _current->generateNeighbor(*_neighbor);
                testNeigh();
                if(_current->getP() < _best->getP())
                {
                        (*_best) = (*_current);
                        _bestIter = 1;
                }
                else if(_current->getP() == _best->getP() && _current->getF() < _best->getF())
                {
                        (*_best) = (*_current);
                        _bestIter = 1;
                }
            }
            ++_bestIter;
            updateParam();
        }
    }

    std::string solution() const
    {
        std::stringstream ss;
        ss << (*_best);
        std::string result = ss.str();
        return result;
    }
    
    void setParameters(double multT, double accept, u_int32_t tBI,
        u_int32_t minIterations, u_int32_t maxIterations)
    {
        _multiplierT = multT;
        _acceptProb = accept;
        _terminalBestIter = tBI;
        _minIterations = minIterations;
        _maxIterations = maxIterations;
    }
    
    void setSolutionParameters(const char* foo)
    {
        _best = SHARED_PTR<SolutionType>(new SolutionType(foo));
        _current = SHARED_PTR<SolutionType>(new SolutionType(*_best));
        _neighbor = SHARED_PTR<SolutionType>(new SolutionType(*_best));
    }
    
    std::ostream& dump(std::ostream& os) const
    {
        os << "Best solution:      " << (_best->getF()) << "\t" 
                   << (_best->getP()) << std::endl;
        os << "Best iter:          " << _bestIter << std::endl;
        os << "Iterations:         " << _iterations << std::endl;
        os << "Count limit:        " << _maxIterations << std::endl;
        os << "Minimum iterations: " << _minIterations << std::endl;
        os << "Sample size:        " << _sampleSize << std::endl;
        os << "Multiplier:         " << _multiplierT << std::endl;
        os << "Accept probability: " << _acceptProb << std::endl;
        os << "Terminal best iter: " << _terminalBestIter << std::endl;
        os << "Pressure:           " << _lambda << std::endl;
        return os;
    }
    
    double cost() const { return _best->getF(); }
    
    double penalty() const { return _best->getP(); }
    u_int32_t bestIter() const { return _bestIter; }

    u_int32_t iterations() const { return _iterations; }
    
    u_int32_t maxIterations() const { return _maxIterations; }
    
    u_int32_t minIterations() const { return _minIterations; }
    double multiplier() const { return _multiplierT; }
    
    double probability() const { return _acceptProb; }
    PenaltyType getLambda() const { return _lambda; }
    
protected:
    void initializeParam()
    {
        double lambda1 = 0, sum = 0;
        double lambda0 = 0, cap = 0;
        double solutionScalePressure = _pfunc.getCapPercentage() /
             (1 - _pfunc.getCapPercentage());
        for (u_int32_t j = 0 ; j < _sampleSize - 1 ; j += 2)
        {
            _current->randomize();
            _current->compute();
            _current->generateNeighbor(*_neighbor);
            if ( _current->getP() > 0 )
                lambda0 = (_current->getF() / _current->getP()) * 
                    solutionScalePressure;
            if (lambda0 > cap)
                cap = lambda0;
            if (_neighbor->getP() > 0)
                lambda0 = (_neighbor->getF() / _neighbor->getP()) * 
                    solutionScalePressure;
            if (lambda0 > cap)
                cap = lambda0;
            double u = (_current->getF() + lambda1 * _current->getP()) - 
                (_neighbor->getF() + lambda1 * _neighbor->getP());
            sum += u > 0 ? u : (-1 * u);
        }
        
        _pfunc.setPressureCap(cap);
        sum /= _sampleSize;
        _currentT = ((-1 * sum) / log(_acceptProb));
    }
    
    void tuneTemperature()
    {
        int acceptedWorse, uphill;
        u_int32_t iterations = 0;
        do
        {
            ++iterations;
            acceptedWorse = uphill = 0;
            for (u_int32_t count = 0 ; count < _maxIterations ; count++)
            {
                _current->generateNeighbor(*_neighbor);
                double delta, u;
                delta = (_neighbor->getF() + _lambda * _neighbor->getP()) -
                        (_current->getF() + _lambda*_current->getP());
                if (delta < 0)
                    _current.swap(_neighbor);
                else
                {
                    uphill++;
                    u = 0 - delta / _currentT;
                    double generated = Twister::generateDouble();
                    if (generated < exp(u))
                    {
                        _current.swap(_neighbor);
                        acceptedWorse++;
                    }
                }
                if ( _current->getP() < _best->getP() )
                    *_best = *_current;
                else if ( (_current->getP() == _best->getP() ) && 
                                                  (_current->getF() < _best->getF()) )
                    *_best = *_current;
            }
            if ((static_cast<double>(acceptedWorse) / 
                                 static_cast<double>(uphill)) < _acceptProb)
                _currentT = 1.5 * _currentT;
        } while ((static_cast<double>(acceptedWorse) / 
                                  static_cast<double>(uphill)) < _acceptProb);
    }

    void testNeigh()
    {
        double delta, u;
        delta = (_neighbor->getF() + _lambda * _neighbor->getP()) - 
                (_current->getF() + _lambda * _current->getP());
        if (delta < 0)
            _current.swap(_neighbor);
        else
        {
            u = 0 - delta / _currentT;
            if (Twister::generateDouble() < exp(u))
                _current.swap(_neighbor);
        }
    }
    
    void updateParam()
    {
        _currentT = _multiplierT * _currentT;
        _lambda = _pfunc(_iterations);
    }
        
    SHARED_PTR<SolutionType> _best, _current, _neighbor;
    
    u_int32_t _bestIter, _iterations, _maxIterations, _minIterations,
        _terminalBestIter;
    static const int _sampleSize = 10000;
    double _multiplierT, _acceptProb, _currentT;
    PenaltyFunc _pfunc;
    PenaltyType _lambda;
};


template <class PenaltyFunc, class SolutionType>
std::ostream& operator<<(std::ostream& os, 
        const Annealer<PenaltyFunc, SolutionType>& engine)
{
    return engine.dump(os);
}

#endif