anneal.h

/**************************************************************************
 * Copyright (c) 2004-2011 T. M. Murali                                   *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by   *
 * the Free Software Foundation, either version 3 of the License, or      *
 * (at your option) any later version.                                    *
 *                                                                        *
 * Biorithm is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 * GNU General Public License for more details.                           *
 *                                                                        *
 * You should have received a copy of the GNU General Public License      *
 * along with Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#ifndef _ANNEAL_H
#define _ANNEAL_H

#include <stdlib.h>
#include <iostream>
#include <vector>

using namespace std;

typedef double MyNT;

// first two constants taken from pathway perturbation paper.
const MyNT DEFAULT_START_TEMPERATURE = 100;
const MyNT DEFAULT_END_TEMPERATURE = 1E-5;
const unsigned int DEFAULT_NUM_ITERATIONS = 10000;
// #iterations = k => DET = DST / (DDT)^k => DDT = (DST/DET)^{1/k}.
const MyNT DEFAULT_DROP_TEMPERATURE = pow(DEFAULT_START_TEMPERATURE/DEFAULT_END_TEMPERATURE, 1.0/DEFAULT_NUM_ITERATIONS);


class MyAbstractAnnealState
{
private:
  MyNT _score;

public:
  MyAbstractAnnealState()
      : _score(0)
    {}
  
  virtual ~MyAbstractAnnealState()
    {}

  MyAbstractAnnealState &operator=(const MyAbstractAnnealState &rhs);
  
  virtual void computeRandomNeighbour(MyAbstractAnnealState &neighbour) = 0;


  virtual MyNT computeScore() const = 0;

  virtual void print(ostream &ostr) const = 0;

  virtual void undo() = 0;

};



template <typename AnnealState >//, typename AnnealMethod>
class MySimulatedAnnealer {
protected:                                              // member data

  // variables related to the annealing policy.
  MyNT _startTemperature;
  MyNT _endTemperature;
  // how much to decrease the temperature each time.
  MyNT _dropTemperature;
  unsigned int _numMaxIterations;
  unsigned int _numMaxIterationsWithoutImprovement;

  MyNT _bestScore;
  MyNT  _currentTemperature;
  unsigned int _numTotalIterations;
  unsigned int _numIterationsWithScoreDecrease;
  unsigned int _numIterationsWithoutImprovement;
  
  AnnealState _startState;
  AnnealState _currentState;
  // vector< AnnealState > _stateHistory;
  // vector< AnnealState > _neighbouringStates;
  
//   AnnealMethod &_method;

  
  
public:
  MySimulatedAnnealer() 
      :
      _startTemperature(DEFAULT_START_TEMPERATURE),
      _endTemperature(DEFAULT_END_TEMPERATURE),
      _dropTemperature(DEFAULT_DROP_TEMPERATURE),
      _numMaxIterations(DEFAULT_NUM_ITERATIONS),
      _numMaxIterationsWithoutImprovement(100),
      _bestScore(0), _currentTemperature(_startTemperature),
      _numTotalIterations(0), _numIterationsWithScoreDecrease(0),
      _numIterationsWithoutImprovement(0),
      _startState(), _currentState()// ,
        // _stateHistory(), _neighbouringStates()
//         _method(m)
    {
    }

  MySimulatedAnnealer(const AnnealState &s) 
      : _startTemperature(DEFAULT_START_TEMPERATURE),
      _endTemperature(DEFAULT_END_TEMPERATURE),
      _dropTemperature(DEFAULT_DROP_TEMPERATURE),
      _numMaxIterations(DEFAULT_NUM_ITERATIONS),
      _numMaxIterationsWithoutImprovement(100),
      _bestScore(0), _currentTemperature(_startTemperature),
        _numTotalIterations(0), _numIterationsWithScoreDecrease(0),
      _numIterationsWithoutImprovement(0),
      _startState(s), _currentState()//,
        // _stateHistory(), _neighbouringStates()
//         _method(m)
    {
    }

  virtual ~MySimulatedAnnealer()
    { 
    }

  bool acceptNeighbour(MyNT score)
    {
      if (score > _bestScore)
        return(true);
      // compute probability.
      MyNT prob = (score - _bestScore)/_currentTemperature;
      prob = exp(prob);
      MyNT sample = random()*1.0/RAND_MAX;
      if (sample <= prob)
        return(true);
      return(false);
    }
  
  MyNT computeBestScore(const vector< AnnealState > &neighbours,
                        AnnealState &bestNeighbour)
    {
      typename vector< AnnealState >::const_iterator aitr = neighbours.begin();
      MyNT score = aitr->computeScore();
      bestNeighbour = *aitr;
      MyNT currentScore;
      
      for (; aitr != neighbours.end(); aitr++)
        {
          currentScore = aitr->computeScore();
          if (currentScore > score)
            {
              score = currentScore;
              bestNeighbour = *aitr;
            }
        }
      return(score);
    }
  
  void decreaseTemperature()
    {
      _currentTemperature /= _dropTemperature;
    }
  
  void getFinalState(AnnealState &finalState) const
    {
      finalState = _currentState;// _stateHistory.back();
    }
  
  void initialise() 
    {
      _bestScore = 0;
      _currentTemperature = DEFAULT_START_TEMPERATURE;
      _numTotalIterations = 0;
      _numIterationsWithoutImprovement = 0;
      _currentState = _startState;
      // _stateHistory.push_back(_startState);
    }

  bool isDone()
    {
      if (_numMaxIterations == _numTotalIterations)
      // if (100 == _numIterationsWithoutImprovement)
        return(true);
      return(false);
    }

  void setNumTotalIterations(unsigned int num)
    {
      _numMaxIterations = num;
      // i should update _dropTemperature as well.
      _dropTemperature = pow(_startTemperature/_endTemperature, 1.0/_numMaxIterations);
    }
  
  void run() 
  {
    MyNT score;
    initialise();
    AnnealState neighbour;// = NULL;
    
    // compute its score.
    score = _startState.computeScore();
    
    while (true)
      {
        // compute the neighbours of the current state.
        // _neighbouringStates.clear();
        // _stateHistory.back().computeNeighbours(_neighbouringStates);
        // if (0 == _neighbouringStates.size())
        //   break;
        // score = computeBestScore(_neighbouringStates, bestNeighbour);

        // compute a random neighbour of the current state.
        _currentState.computeRandomNeighbour(neighbour);
        // _stateHistory.back().computeRandomNeighbour(neighbour);
        score = neighbour.computeScore();

        if (acceptNeighbour(score))
          {
            _currentState = neighbour;
            // _stateHistory.push_back(neighbour);
            if (_bestScore > score)
              // this is an iteration where the score decreased.
              _numIterationsWithScoreDecrease++;
            _bestScore = score;
            _numIterationsWithoutImprovement = 0;
          }
        else
          {
            // Reverse the changes made by the neighbour.
            neighbour.undo();
            _numIterationsWithoutImprovement++;
          }
        // I can delete neighbour in any case. If accepted, I have put
        // a copy of it at the back of _stateHistory.
        // delete neighbour;
        decreaseTemperature();
        _numTotalIterations++;
//#ifdef DEBUG
        if (0 == (_numTotalIterations%10000))
          {
            cout << "\t" << _numTotalIterations
                 << "\t" << _numIterationsWithScoreDecrease
                 << "\t" << _numIterationsWithScoreDecrease*1.0/_numTotalIterations;
            _currentState.print(cout);
          }
//#endif // DEBUG

        if (isDone())
          {
#ifdef DEBUG
            _currentState.print(cout);
            // _stateHistory.back().print(cout);
#endif // DEBUG
            break;
          }
      }
  }

  void setInitialState(AnnealState &s)
    {
      _startState = s;
    }
  
};


#endif // _ANNEAL_H