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
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
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
00041 m_pMO->ElementCopy(parameters,m_pdCurrentParameters,nParameters);
00042
00043
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
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 }
00082 }
00083 cout << "Acceptance ratio of thermal step is " << m_dAcceptanceRatio/m_iNMCSteps << endl;
00084
00085
00086 m_pMO->ElementCopy(m_pdCurrentParameters,pSavedParameters,nParameters);
00087 cout << "Quenching . . . " << endl;
00088 double Equench = m_pQuencher->Minimize(m_pdCurrentParameters,minimizable);
00089
00090 cout << "Energy of quenched parameters is " << Equench << endl;
00091
00092 this->WriteQuenchedParameters(j,Equench);
00093
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
00107
00108 char fname[20];
00109 sprintf(fname,"quench%i.par",iter);
00110
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
00128 double normE = m_dECurrent/m_dTemperature;
00129 double Ai = m_dDeltaX*((m_dC*normE + 1.0)/(normE + 1.0));
00130
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
00145
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 }