Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members  

PeriodicQuenchMinimizer.cpp

Go to the documentation of this file.
00001 #include "PeriodicQuenchMinimizer.h"
00002 
00003 CPeriodicQuenchMinimizer::CPeriodicQuenchMinimizer(int RNGSeed, CParameterFilter *pFilter, bool LMflag, int nQuench, int nSteps, double tol, double T, double deltaX, double c)
00004 {
00005         m_iRNGSeed = RNGSeed;
00006         m_pFilter = pFilter;
00007         m_pRNG = new Rand(m_iRNGSeed);
00008         m_dTemperature = T;
00009         m_dDeltaX = deltaX;
00010         m_dC = c;
00011         m_iNQuenches = nQuench;
00012         m_iNMCSteps = nSteps;
00013         // if true, the quencher is levenberg-marquardt
00014         if(LMflag = true)
00015         {
00016                 m_pQuencher = new CImprovedLevenbergMarquardtMinimizer(m_pFilter,false,1.0e-03,true,tol,true,tol,true,tol,500,1.0e-06);
00017         }
00018         else // do CG
00019         {
00020                 m_pQuencher = new CConjugateGradientMinimizer(false,m_pFilter,500,1.0e-06,3,tol,1.0e-06);
00021         }
00022         m_pMO = new CMatrixOperations();
00023 }
00024 
00025 CPeriodicQuenchMinimizer::~CPeriodicQuenchMinimizer(void)
00026 {
00027         delete m_pRNG;
00028         delete m_pQuencher;
00029         delete m_pMO;
00030 }
00031 
00032 double CPeriodicQuenchMinimizer::Minimize(double *parameters, Minimizable *minimizable)
00033 {
00034         int i,j = 0;
00035         nParameters = minimizable->GetNParameters();
00036         m_pdCurrentParameters = new double[nParameters];
00037         m_pdTrialParameters = new double[nParameters];
00038         double *pSavedParameters = new double[nParameters];
00039 
00040         // copy the current parameters into starting parameters
00041         m_pMO->ElementCopy(parameters,m_pdCurrentParameters,nParameters);
00042 
00043         // get the initial free energy
00044         double Eold = minimizable->F(parameters,m_dTemperature);
00045         cout << "Initial energy of " << Eold << endl;
00046         m_dECurrent = Eold;
00047         m_dAcceptanceRatio = 0;
00048         
00049         cout << "Equilibrating . . ." << endl;
00050         // now equilibrate
00051         for(i = 0; i < m_iNMCSteps; i++)
00052         {
00053                 this->GenerateMove();
00054                 double Enew = minimizable->F(m_pdTrialParameters,m_dTemperature);
00055 
00056                 if(this->AcceptMove(Eold,Enew))
00057                 {
00058                         Eold = Enew;
00059                         m_dAcceptanceRatio += 1.0;
00060                         m_pMO->ElementCopy(m_pdTrialParameters,m_pdCurrentParameters,nParameters);
00061                 }
00062         }
00063         cout << "Acceptance ratio during equilibration : " << m_dAcceptanceRatio/m_iNMCSteps << endl;
00064 
00065 
00066         for(j = 0; j < m_iNQuenches; j++)
00067         {
00068                 cout << "Starting iteration " << j << " of " << m_iNQuenches << endl;
00069                 m_dAcceptanceRatio = 0.0;
00070                 for(i = 0; i < m_iNMCSteps; i++)
00071                 {
00072                         this->GenerateMove();
00073                         double Enew = minimizable->ObjectiveFunction(m_pdTrialParameters);
00074 
00075                         if(this->AcceptMove(Eold,Enew))
00076                         {                       
00077                                 Eold = Enew;
00078                                 m_dAcceptanceRatio += 1.0;
00079                                 m_pMO->ElementCopy(m_pdTrialParameters,m_pdCurrentParameters,nParameters);
00080 
00081                         } // end AcceptMove if
00082                 }
00083                 cout << "Acceptance ratio of thermal step is " << m_dAcceptanceRatio/m_iNMCSteps << endl;
00084                 // now do a quench
00085                 // save the current parameters
00086                 m_pMO->ElementCopy(m_pdCurrentParameters,pSavedParameters,nParameters);
00087                 cout << "Quenching . . . " << endl;
00088                 double Equench = m_pQuencher->Minimize(m_pdCurrentParameters,minimizable);
00089                 // report the energy of the quenched parameters
00090                 cout << "Energy of quenched parameters is " << Equench << endl;
00091                 // now write them to a file
00092                 this->WriteQuenchedParameters(j,Equench);
00093                 // restore old parameters
00094                 m_pMO->ElementCopy(pSavedParameters,m_pdCurrentParameters,nParameters);
00095         }
00096 
00097         delete [] m_pdCurrentParameters;
00098         delete [] m_pdTrialParameters;
00099         delete [] pSavedParameters;
00100 
00101         return Eold;
00102 }
00103 
00104 void CPeriodicQuenchMinimizer::WriteQuenchedParameters(int iter, double E)
00105 {
00106         // dump whatever is in m_pdCurrentParameters to a file with the appropriate
00107         // number tag
00108         char fname[20];
00109         sprintf(fname,"quench%i.par",iter);
00110         // write the best (temporary) parameters to a file and stdout
00111         fstream *pfquench = new fstream(fname,ios::out);
00112         *pfquench << "##################################################################" << endl;
00113         *pfquench << "# QUENCHED PARAMETERS FROM ITERATION " << iter << endl;
00114         *pfquench << "# ENERGY : " << E << endl;
00115         *pfquench << "##################################################################" << endl;
00116         for(int i = 0; i < nParameters; i++)
00117         {
00118                 *pfquench << m_pdCurrentParameters[i] << endl;
00119         }
00120         delete pfquench;
00121         return;
00122 }
00123 
00124 void CPeriodicQuenchMinimizer::GenerateMove()
00125 {
00126         m_dSqrLogParDiff = 0.0;
00127         // compute the stepsize
00128         double normE = m_dECurrent/m_dTemperature;
00129         double Ai = m_dDeltaX*((m_dC*normE + 1.0)/(normE + 1.0));
00130         // now take a step w/ std dev. Ai
00131         m_pFilter->ForwardTransformation(m_pdCurrentParameters,nParameters);
00132         for(int i = 0; i < nParameters; i++)
00133         {
00134                 double delta = m_pRNG->gaussian(Ai);
00135                 m_pdTrialParameters[i] = m_pdCurrentParameters[i] + delta;
00136                 m_dSqrLogParDiff += delta*delta;
00137         }
00138         m_pFilter->BackwardTransformation(m_pdTrialParameters,nParameters);
00139         m_pFilter->BackwardTransformation(m_pdCurrentParameters,nParameters);
00140 }
00141 
00142 bool CPeriodicQuenchMinimizer::AcceptMove(double Eold, double Enew)
00143 {
00144         // probability has less simple form b/c of detailed balance correction
00145         // forced by shifting stepsize
00146         double normE = Enew/m_dTemperature;
00147         double Af = m_dDeltaX*((m_dC*normE + 1.0)/(normE + 1.0));
00148         normE = Eold/m_dTemperature;
00149         double Ai = m_dDeltaX*((m_dC*normE + 1.0)/(normE + 1.0));
00150         
00151         double logprob = -(Enew - Eold)/m_dTemperature;
00152         logprob += nParameters*log(Ai/Af);
00153         logprob -= 0.5*m_dSqrLogParDiff*(1.0/(Af*Af) - 1.0/(Ai*Ai));
00154 
00155         if(exp(logprob) > 1.0)
00156         {
00157                 return true;
00158         }
00159         else if(m_pRNG->uniform() < exp(logprob))
00160         {
00161                 return true;    
00162         }
00163 
00164         return false;
00165 }

Generated on Mon Nov 3 09:38:08 2003 by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002