---

VariableStepsizeAnnealMinimizer.cpp

Go to the documentation of this file.

// VariableStepsizeAnnealMinimizer.cpp: implementation of the CVariableStepsizeAnnealMinimizer class.
//

#include "VariableStepsizeAnnealMinimizer.h"

// Construction/Destruction

CVariableStepsizeAnnealMinimizer::CVariableStepsizeAnnealMinimizer(int RNGSeed, CParameterFilter *pFilter, double T0, double lambda, double deltaX, double c)
:CSimulatedAnnealingStrategy(RNGSeed,pFilter)
{
        m_dTemperature = T0;
        m_dLambda = lambda;
        m_dDeltaX = deltaX;
        m_dC = c;
}

CVariableStepsizeAnnealMinimizer::~CVariableStepsizeAnnealMinimizer()
{

}

bool CVariableStepsizeAnnealMinimizer::AcceptMove(double Eold, double Enew)
{
        // probability has less simple form b/c of detailed balance correction
        // forced by shifting stepsize
        double normE = Enew/m_dTemperature;
        double Af = m_dDeltaX*((m_dC*normE + 1.0)/(normE + 1.0));
        normE = Eold/m_dTemperature;
        double Ai = m_dDeltaX*((m_dC*normE + 1.0)/(normE + 1.0));
        
        double logprob = -(Enew - Eold)/m_dTemperature;
        logprob += nParameters*log(Ai/Af);
        logprob -= 0.5*m_dSqrLogParDiff*(1.0/(Af*Af) - 1.0/(Ai*Ai));

        if(exp(logprob) > 1.0)
        {
                return true;
        }
        else if(m_pRNG->uniform() < exp(logprob))
        {
                return true;    
        }

        return false;
}

bool CVariableStepsizeAnnealMinimizer::Terminate()
{
        double Ebar = this->GetEBar();
        double ratio = (Ebar - m_dEBest)/Ebar;

        if(ratio < 1.0e-03)
        {
                return true;
        }
        
        return false;
}

bool CVariableStepsizeAnnealMinimizer::Equilibrated()
{
        if(m_iTrialCount > m_iNMCSteps)
        {
                return true;
        }

        return false;
}

void CVariableStepsizeAnnealMinimizer::Cool()
{
        
        double sigma = this->GetE2Bar() - (this->GetEBar())*(this->GetEBar());
        double expFac = -(m_dLambda*m_dTemperature)/sigma;
        m_dTemperature = m_dTemperature*exp(expFac);
        
        // increase the number of MC steps as we cool
        m_iNMCSteps = m_iNMCSteps*1.15;

        // write the best (temporary) parameters to a file and stdout
        fstream *pfbest = new fstream("satemp.par",ios::out);
        *pfbest << "##################################################################" << endl;
        *pfbest << "# BEST PARAMETERS (PARTIAL)" << endl;
        *pfbest << "##################################################################" << endl;
        for(int i = 0; i < nParameters; i++)
        {
                *pfbest << m_pdCurrentParameters[i] << endl;
        }
        delete pfbest;

}

void CVariableStepsizeAnnealMinimizer::InitializeTemperature(Minimizable *minimizable)
{
        // set number of equilibration steps
        m_iNMCSteps = 100;
        // save starting parameters
        double *savepar = new double[nParameters];
        m_pMO->ElementCopy(m_pdCurrentParameters,savepar,nParameters);
        double Eold = minimizable->ObjectiveFunction(m_pdCurrentParameters);
        double tempT = m_dTemperature;

        // run 50 steps and check the acceptance ratio, then increase
        // T if A < 0.8
        while(m_dAcceptanceRatio < 0.8)
        {
                m_dAcceptanceRatio = 0.0;
                m_dTemperature = tempT;
                for(int i = 0; i < 50; i++)
                {
                        this->GenerateMove();
                        double Enew = minimizable->ObjectiveFunction(m_pdTrialParameters);

                        if(this->AcceptMove(Eold,Enew))
                        {                       
                                Eold = Enew;
                                m_dAcceptanceRatio += 1.0;
                                m_pMO->ElementCopy(m_pdTrialParameters,m_pdCurrentParameters,nParameters);
                        }
                }
                m_dAcceptanceRatio = m_dAcceptanceRatio/50;
                tempT = tempT*1.5;
        }
        // restore original parameters
        m_pMO->ElementCopy(savepar,m_pdCurrentParameters,nParameters);
        
        cout << "T0 = " << m_dTemperature << endl;
        
        delete [] savepar;
        
        return;
}

void CVariableStepsizeAnnealMinimizer::GenerateMove()
{
        m_dSqrLogParDiff = 0.0;
        // compute the stepsize
        double normE = m_dECurrent/m_dTemperature;
        double Ai = m_dDeltaX*((m_dC*normE + 1.0)/(normE + 1.0));
        // now take a step w/ std dev. Ai
        m_pFilter->ForwardTransformation(m_pdCurrentParameters,nParameters);
        for(int i = 0; i < nParameters; i++)
        {
                double delta = m_pRNG->gaussian(Ai);
                m_pdTrialParameters[i] = m_pdCurrentParameters[i] + delta;
                m_dSqrLogParDiff += delta*delta;
        }
        m_pFilter->BackwardTransformation(m_pdTrialParameters,nParameters);
        m_pFilter->BackwardTransformation(m_pdCurrentParameters,nParameters);
}

---