---

PeriodicQuenchMinimizer.cpp

Go to the documentation of this file.

#include "PeriodicQuenchMinimizer.h"

CPeriodicQuenchMinimizer::CPeriodicQuenchMinimizer(int RNGSeed, CParameterFilter *pFilter, bool LMflag, int nQuench, int nSteps, double tol, double T, double deltaX, double c)
{
        m_iRNGSeed = RNGSeed;
        m_pFilter = pFilter;
        m_pRNG = new Rand(m_iRNGSeed);
        m_dTemperature = T;
        m_dDeltaX = deltaX;
        m_dC = c;
        m_iNQuenches = nQuench;
        m_iNMCSteps = nSteps;
        // if true, the quencher is levenberg-marquardt
        if(LMflag = true)
        {
                m_pQuencher = new CImprovedLevenbergMarquardtMinimizer(m_pFilter,false,1.0e-03,true,tol,true,tol,true,tol,500,1.0e-06);
        }
        else // do CG
        {
                m_pQuencher = new CConjugateGradientMinimizer(false,m_pFilter,500,1.0e-06,3,tol,1.0e-06);
        }
        m_pMO = new CMatrixOperations();
}

CPeriodicQuenchMinimizer::~CPeriodicQuenchMinimizer(void)
{
        delete m_pRNG;
        delete m_pQuencher;
        delete m_pMO;
}

double CPeriodicQuenchMinimizer::Minimize(double *parameters, Minimizable *minimizable)
{
        int i,j = 0;
        nParameters = minimizable->GetNParameters();
        m_pdCurrentParameters = new double[nParameters];
        m_pdTrialParameters = new double[nParameters];
        double *pSavedParameters = new double[nParameters];

        // copy the current parameters into starting parameters
        m_pMO->ElementCopy(parameters,m_pdCurrentParameters,nParameters);

        // get the initial free energy
        double Eold = minimizable->F(parameters,m_dTemperature);
        cout << "Initial energy of " << Eold << endl;
        m_dECurrent = Eold;
        m_dAcceptanceRatio = 0;
        
        cout << "Equilibrating . . ." << endl;
        // now equilibrate
        for(i = 0; i < m_iNMCSteps; i++)
        {
                this->GenerateMove();
                double Enew = minimizable->F(m_pdTrialParameters,m_dTemperature);

                if(this->AcceptMove(Eold,Enew))
                {
                        Eold = Enew;
                        m_dAcceptanceRatio += 1.0;
                        m_pMO->ElementCopy(m_pdTrialParameters,m_pdCurrentParameters,nParameters);
                }
        }
        cout << "Acceptance ratio during equilibration : " << m_dAcceptanceRatio/m_iNMCSteps << endl;


        for(j = 0; j < m_iNQuenches; j++)
        {
                cout << "Starting iteration " << j << " of " << m_iNQuenches << endl;
                m_dAcceptanceRatio = 0.0;
                for(i = 0; i < m_iNMCSteps; 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);

                        } // end AcceptMove if
                }
                cout << "Acceptance ratio of thermal step is " << m_dAcceptanceRatio/m_iNMCSteps << endl;
                // now do a quench
                // save the current parameters
                m_pMO->ElementCopy(m_pdCurrentParameters,pSavedParameters,nParameters);
                cout << "Quenching . . . " << endl;
                double Equench = m_pQuencher->Minimize(m_pdCurrentParameters,minimizable);
                // report the energy of the quenched parameters
                cout << "Energy of quenched parameters is " << Equench << endl;
                // now write them to a file
                this->WriteQuenchedParameters(j,Equench);
                // restore old parameters
                m_pMO->ElementCopy(pSavedParameters,m_pdCurrentParameters,nParameters);
        }

        delete [] m_pdCurrentParameters;
        delete [] m_pdTrialParameters;
        delete [] pSavedParameters;

        return Eold;
}

void CPeriodicQuenchMinimizer::WriteQuenchedParameters(int iter, double E)
{
        // dump whatever is in m_pdCurrentParameters to a file with the appropriate
        // number tag
        char fname[20];
        sprintf(fname,"quench%i.par",iter);
        // write the best (temporary) parameters to a file and stdout
        fstream *pfquench = new fstream(fname,ios::out);
        *pfquench << "##################################################################" << endl;
        *pfquench << "# QUENCHED PARAMETERS FROM ITERATION " << iter << endl;
        *pfquench << "# ENERGY : " << E << endl;
        *pfquench << "##################################################################" << endl;
        for(int i = 0; i < nParameters; i++)
        {
                *pfquench << m_pdCurrentParameters[i] << endl;
        }
        delete pfquench;
        return;
}

void CPeriodicQuenchMinimizer::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);
}

bool CPeriodicQuenchMinimizer::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;
}

---