NICE_Optimization/CombinatorialOptimizer.cpp

//////////////////////////////////////////////////////////////////////
//
//	CombinatorialOptimizer.cpp: Implementation of the 
//		CombinatorialOptimizer class.
//
//////////////////////////////////////////////////////////////////////

#include "optimization/CombinatorialOptimizer.h"
#include <stdlib.h>
#include <time.h>

#include "numbase.h" // ice random

#include <iostream>

using namespace OPTIMIZATION;

CombinatorialOptimizer::CombinatorialOptimizer(OptLogBase *loger): SuperClass(loger)
{
  m_mode = 1;
  m_T0 = 1;
  m_Tk = m_T0;
  m_alpha = 0.95;
  //m_stepLength = 1.0;
  m_fast = false;
}


CombinatorialOptimizer::CombinatorialOptimizer( const CombinatorialOptimizer &opt) : SuperClass(opt)
{
  m_mode			= opt.m_mode;
  m_T0			= opt.m_T0;
  m_Tk			= opt.m_Tk;
  m_alpha			= opt.m_alpha;
  //m_stepLength		= opt.m_stepLength;
  m_fast			= opt.m_fast;
  
}

/*
  operator=
*/
CombinatorialOptimizer & CombinatorialOptimizer::operator=(const CombinatorialOptimizer &opt)
{
    
  /*
      avoid self-copy
  */
  if(this != &opt)
  {
    
    /*
      =operator of SuperClass
    */
    SuperClass::operator=(opt);

      
    /*
      own values:
    */
    m_mode			= opt.m_mode;
    m_T0			= opt.m_T0;
    m_Tk			= opt.m_Tk;
    m_alpha			= opt.m_alpha;
    //m_stepLength	= opt.m_stepLength;
    m_fast			= opt.m_fast;
    
    
  }

    return *this;
}


CombinatorialOptimizer::~CombinatorialOptimizer()
{
}


void CombinatorialOptimizer::init()
{
  SuperClass::init();

  m_currentCostFunctionValue = evaluateCostFunction(m_parameters);

  /*
    seed rand
  */
  ice::Randomize();

  /*
    (re)set m_Tk
  */
  m_Tk = m_T0;



}

bool CombinatorialOptimizer::setMode(int mode)
{
  switch(mode)
  {
  case 1:
  case 2:
  case 3:
  case 4:
    m_mode = mode;
    return true;
  default:
    return false;
  }
}

bool CombinatorialOptimizer::setControlSeqParam(double T0,double alpha)
{
  if(T0 < 0)
    return false;

  m_T0 = T0;

  if(alpha < 0 || alpha > 1)
    return false;

  return true;

}

bool CombinatorialOptimizer::accepptPoint(double oldValue, double newValue)
{
  bool accept = false;
  //cout<<"old val: "<<oldValue<<endl;
    //    cout<<"new val: "<<newValue<<endl;
    //    cout<<endl;
  switch(m_mode)
  {
  case 1:
  {
    /*
      simulated annealing
    *****************************/
    
    if(newValue <= oldValue)
    {
      accept = true;
      break;
    }
        
    /*
      generate uniform distrubuted random number
    */
    double lambda = ice::RandomD();
    
    if(lambda <= exp(-(newValue - oldValue)/((m_Tk < 1.0e-20)? 1.0e-20 :m_Tk)))
    {
      accept =true;
    }
    break;
  }
  case 2:
    /*
      threshold acceptance
    *****************************/
    if((newValue - oldValue) < m_Tk)
    {
      accept = true;
    }
    break;
  case 3:
    /*
      great deluge algorithm
    *****************************/
    if(newValue <= m_Tk)
    {
      accept = true;
    }

    break;
  case 4:
    /*
      stochastic relaxation
    *****************************/
    if(newValue <= oldValue)
    {
      accept = true;
    }
  case 5:
    {
    /*
      annealing & theshold
    *****************************/
    
    /*
      accept better values
    */
    if(newValue <= oldValue)
    {
      accept = true;
      break;
    }
    
    /*
      dont accept values bigger than m_threshold
    */
    if(newValue > m_threshold )
    {
      accept = false;
      break;
    }

    /*
      generate uniform distrubuted random number
    */
    double lambda = ice::RandomD();
    
    if(lambda <= exp(-(newValue - oldValue)/((m_Tk < 1.0e-20)? 1.0e-20 :m_Tk)))
    {
      accept =true;
    }
    else
      accept = false;
    break;
    }
  default:
    accept = false;
  }

  return accept;
}

matrix_type CombinatorialOptimizer::generatePoint()
{
    matrix_type newPoint(m_parameters);

  for(int i = 0; i < static_cast<int>(m_numberOfParameters); i++)
  {
    if(m_scales(i,0) > 0.0 )
    {
      if(m_fast == true)
      {
        //newPoint(i,0) = m_parameters(i,0) + m_stepLength * m_Tk * ice::GaussRandom(1.0) ;
        newPoint(i,0) = m_parameters(i,0) + m_scales(i,0) * m_Tk * ice::GaussRandom(1.0) ;
      }
      else
      {
        newPoint(i,0) = m_parameters(i,0) + m_scales(i,0) * ice::GaussRandom(1.0) ;
      }
    }
  }

        //cout<<"new Point: "<<newPoint<<endl;
        
  return newPoint;
}	

//bool CombinatorialOptimizer::setSteplength(double steplength)
//{
//	if(steplength <= 0)
//		return false;
//
//	m_stepLength = steplength;
//
//	return true;
//}

int CombinatorialOptimizer::optimize()
{

  // init
  init();


  /*
    start time criteria
  */
  m_startTime = clock();



  bool abort = false;

  matrix_type newPoint;
  double newValue;
  
  /* 
    do the optimization
  */
  while(abort == false)
  {
    /*
      increase iteration counter
    */
    m_numIter++;
    
    /*
      Check iteration limit
    */
    if(m_maxNumIterActive)
    {
      if(m_numIter >= m_maxNumIter)
      {
        /* set according return status and return */
        m_returnReason = SUCCESS_MAXITER;
        abort = true;
        continue;
      }
    }

    /*
      update control seq.
    */
    m_Tk *= m_alpha;


    /*
      generate new point
    */
    newPoint = generatePoint();
        

    /*
      evaluate cost function
    */
    newValue = evaluateCostFunction(newPoint);

    
    /*
      accept new point
    */
    if( accepptPoint(m_currentCostFunctionValue,newValue) )
    {
            m_parameters = newPoint;
      m_currentCostFunctionValue = newValue;
      
    }
                
                if(m_verbose)
                {
                    std::cout << "CombinatorialOptimizer :: parameter: ";
                    for(int r = 0; r < static_cast<int>(m_numberOfParameters); r++)
                    {
                        std::cout<< m_parameters(r,0) << " ";
                    }
                    std::cout << m_currentCostFunctionValue  << std::endl;
                }
    
    /*
      Check if it is in bounds, paramTol, funcTol, NumIter, gradienttol, maxSeconds
    */
    if(!checkParameters(m_parameters))
    {
      /* set according return status and the last parameters and return */
      m_returnReason = ERROR_XOUTOFBOUNDS;
      abort = true;
    }

    /*
      Check Optimization Timelimit, if active
    */
    if(m_maxSecondsActive)
    {
      m_currentTime = clock();

      /* time limit exceeded ? */
      if(((float)(m_currentTime - m_startTime )/CLOCKS_PER_SEC) >= m_maxSeconds )
      {
        /* set according return status and the last parameters and return */
        m_returnReason = SUCCESS_TIMELIMIT;
        abort = true;
        continue;
      }
    }

  }

  return m_returnReason;

}