StochasticSensitivityAnalysis.cpp Source File

00001 // StochasticSensitivityAnalysis.cpp: implementation of the CStochasticSensitivityAnalysis class.
00002 //
00004 
00005 #include "StochasticSensitivityAnalysis.h"
00006 
00008 // Construction/Destruction
00010 
00011 CStochasticSensitivityAnalysis::CStochasticSensitivityAnalysis()
00012 {
00013 
00014 }
00015 
00016 CStochasticSensitivityAnalysis::CStochasticSensitivityAnalysis(bool readFlag, bool hessianFlag, int nEnsembleSize, int seed, CParameterFilter *pFilter, double funcAccuracy, double temperature)
00017 {
00018         m_bSVDFlag = false;
00019         m_bReadFlag = readFlag;
00020         m_bHessianFlag = hessianFlag;
00021         m_iRNGSeed = seed;
00022         m_iNEnsembleSize = nEnsembleSize;
00023         m_dTemperature = temperature;
00024         m_dRescale = 1.0/5.0;
00025         m_dAcceptanceRatio = 0.0;
00026         m_dRWValue = 0.0;
00027         m_dMNormSquared = 0.0;
00028         m_dFuncAccuracy = funcAccuracy;
00029         m_pFilter = pFilter;
00030         m_pRNG = new Rand(seed);
00031 }
00032 
00033 CStochasticSensitivityAnalysis::~CStochasticSensitivityAnalysis()
00034 {
00035         delete m_pRNG;
00036         delete [] m_pdSingularValues;
00037         for(int i = 0; i < m_vpEnsemble.size(); i++)
00038         {
00039                 delete m_vpEnsemble[i];
00040         }
00041 }
00042 
00043 void CStochasticSensitivityAnalysis::DoAnalysis(double *pFittedParameters, Minimizable *pMinimizable)
00044 {
00045         // cost difference between trial move and current move
00046         double deltaCost = 0.0;
00047         double quadratic = 0.0;
00048         double newCost = 0.0;
00049         // set the number of parameters
00050         m_iNParameters = pMinimizable->GetNParameters();
00051 
00052         // downcast to a sum-of-squares minimizable (this routine isn't designed to handle
00053         // anything else)
00054         NLLSMinimizable *pNLLS = dynamic_cast<NLLSMinimizable *>(pMinimizable);
00055         if(0 == pNLLS) {throw std::runtime_error("Minimizable is not a sum of squares!");}
00056 
00057         m_iNResiduals = pNLLS->GetNResiduals();
00058 
00059         // allocate memory
00060         m_pdHessian = new double *[m_iNParameters];
00061         m_pdSamplingMatrix = new double *[m_iNParameters];
00062         m_pdMatrixV = new double *[m_iNParameters];
00063         m_pdHessian[0] = new double[m_iNParameters*m_iNParameters];
00064         m_pdSamplingMatrix[0] = new double[m_iNParameters*m_iNParameters];
00065         m_pdMatrixV[0] = new double[m_iNParameters*m_iNParameters];
00066         for(int nPar = 1; nPar < m_iNParameters; nPar++)
00067         {
00068                 m_pdHessian[nPar] = &(m_pdHessian[0][nPar*m_iNParameters]);
00069                 m_pdSamplingMatrix[nPar] = &(m_pdSamplingMatrix[0][nPar*m_iNParameters]);
00070                 m_pdMatrixV[nPar] = &(m_pdMatrixV[0][nPar*m_iNParameters]);
00071         }
00072         m_pdJacobian = new double *[m_iNResiduals];
00073         m_pdJacobian[0] = new double[m_iNResiduals*m_iNParameters];
00074         for(int nR = 1; nR < m_iNResiduals; nR++)
00075         {
00076                 m_pdJacobian[nR] = &(m_pdJacobian[0][nR*m_iNParameters]);
00077         }
00078 
00079         // allocate space for current and trial parameters
00080         double *pParCurrent = new double[m_iNParameters];
00081         double *pParTrial = new double[m_iNParameters];
00082         double *pParDelta = new double[m_iNParameters];
00083         m_pdSingularValues = new double[m_iNParameters];
00084 
00085         // open a file to write the ratio of the actual change in cost to that 
00086         // predicted by the Hessian at the fitted values
00087         fstream *pfRatio = new fstream("ratio.par",ios::out);
00088         *pfRatio << "##################################################################" << endl;
00089         *pfRatio << "# delta(C)/(1/2) k.H0.k" << endl;
00090         *pfRatio << "# TEMP = " << m_dTemperature << endl;
00091         *pfRatio << "##################################################################" << endl;
00092 
00093         // set the current parameters equal to the best parameters
00094         for(int i = 0; i < m_iNParameters; i++)
00095         {
00096                 pParCurrent[i] = pFittedParameters[i];
00097         }
00098         // get the cost at the minimum; if you're using a synchronized integrator, this
00099         // has the additional benefit of setting up the list of stepsizes before the 
00100         // derivative is computed
00101         // OLD
00102         // double currentCost = pNLLS->ObjectiveFunction(pFittedParameters);
00103         // NEW
00104         double currentCost = pNLLS->F(pFittedParameters,m_dTemperature);
00105         // add the minimum to the ensemble
00106         AddParametersToEnsemble(pFittedParameters,pNLLS->GetResiduals(),0.0,currentCost);
00107         
00108         cout << "Cost at minimum is " << currentCost << endl;
00109 
00110         if(m_bReadFlag)
00111         {
00112                 ReadHessian();
00113         }
00114         else
00115         {
00116                 // compute the Hessian at the minimum
00117                 if(m_bHessianFlag)
00118                 {
00119                         ComputeTrueHessian(pFittedParameters,pNLLS);
00120                 }
00121                 else
00122                 {
00123                         ComputeApproximateHessian(pFittedParameters,pNLLS);
00124                 }
00125         }
00126 
00127         // now sample away
00128         double currentEnsembleSize = 1;
00129         int nTrials = 0;
00130         int nAccept = 0;
00131         while(currentEnsembleSize < m_iNEnsembleSize)
00132         {
00133                 // sample
00134                 quadratic = SampleGaussian(pParDelta);
00135                 //quadratic = SampleWithoutExplicitHessian(pParDelta);
00136 
00137                 // construct the trial move
00138                 ConstructTrialMove(pParTrial,pParCurrent,pParDelta);
00139                 
00140                 // OLD
00141                 // evaluate the cost at the new parameter values
00142                 //newCost = pNLLS->ObjectiveFunction(pParTrial);
00143                 //deltaCost = newCost - currentCost;
00144 
00145                 // NEW
00146                 newCost = pNLLS->F(pParTrial,m_dTemperature);
00147                 deltaCost = newCost - currentCost;
00148 
00149                 // write the cost ratio
00150                 *pfRatio << deltaCost/quadratic << endl;
00151 
00152                 // if we don't accept, just repeat
00153                 if(AcceptMove(quadratic,deltaCost))
00154                 {
00155                         nAccept++;
00156                         m_dAcceptanceRatio += 1.0;
00157                         cout << "Move accepted . . ." << endl;
00158                         // copy parameters
00159                         for(int i = 0; i < m_iNParameters; i++)
00160                         {
00161                                 pParCurrent[i] = pParTrial[i];
00162                                 currentCost = newCost;
00163                         }
00164                         currentEnsembleSize++;
00165                         // increment the RW value
00166                         m_dRWValue += m_dRescale*m_dMNormSquared;
00167                         // add another member to the ensemble
00168                         AddParametersToEnsemble(pParCurrent,pNLLS->GetResiduals(),quadratic,currentCost);
00169                         // notify observers
00170                         Notify();
00171                 }
00172         
00173                 nTrials++;
00174         }
00175 
00176         cout << "Acceptance Ratio : " << ((double) nAccept)/((double) nTrials) << endl;
00177 
00178         // temporary arrays
00179         double *minEigen = new double[m_iNParameters];
00180         double *eigenP = new double[m_iNParameters];
00181 
00182         // before returning, write costs/parameters/etc. to a file
00183         fstream *pfPars = new fstream("ensemble.par",ios::out);
00184         fstream *pfQuad = new fstream("quadratic.par",ios::out);
00185         fstream *pfCost = new fstream("cost.par",ios::out);
00186         fstream *pfResid = new fstream("residuals.par",ios::out);
00187         fstream *pfEigP = new fstream("eigenpars.par",ios::out);
00188         *pfPars << "##################################################################" << endl;
00189         *pfQuad << "##################################################################" << endl;
00190         *pfCost << "##################################################################" << endl;
00191         *pfResid << "##################################################################" << endl;
00192         *pfEigP << "##################################################################" << endl;
00193         *pfPars << "# PARAMETER ENSEMBLE - POUND IS SEPARATOR" << endl;
00194         *pfPars << "# TEMP = " << m_dTemperature << endl;
00195         *pfQuad << "# VALUE OF (1/2)k.H.k FOR ALL MEMBERS OF ENSEMBLE" << endl;
00196         *pfQuad << "# TEMP = " << m_dTemperature << endl;
00197         *pfCost << "# ACTUAL COST (NOT DELTA) FOR ALL MEMBERS OF ENSEMBLE" << endl;
00198         *pfCost << "# TEMP = " << m_dTemperature << endl;
00199         *pfResid << "# RESIDUALS FOR ALL MEMBERS OF ENSEMBLE - POUND IS SEPARATOR" << endl;
00200         *pfResid << "# TEMP = " << m_dTemperature << endl;
00201         *pfEigP << "# EIGENPARAMETER ENSEMBLE - DOUBLE NEWLINE IS SEPARATOR" << endl;
00202         *pfEigP << "# TEMP = " << m_dTemperature << endl;
00203         *pfResid << "##################################################################" << endl;
00204         *pfPars << "##################################################################" << endl;
00205         *pfQuad << "##################################################################" << endl;
00206         *pfCost << "##################################################################" << endl;
00207         *pfEigP << "##################################################################" << endl;
00208         for(int ensNum = 0; ensNum < m_vpEnsemble.size(); ensNum++)
00209         {
00210                 *pfQuad << setprecision(10) << m_vpEnsemble[ensNum]->GetQuadratic() << endl;
00211                 *pfCost << setprecision(10) << m_vpEnsemble[ensNum]->GetCost() << endl;
00212                 for(int j = 0; j < m_vpEnsemble[ensNum]->GetParameterSize(); j++)
00213                 {
00214                         *pfPars << m_vpEnsemble[ensNum]->GetParameter(j) << endl;
00215                 }
00216                 for(int j = 0; j < m_vpEnsemble[ensNum]->GetResidualsSize(); j++)
00217                 {
00218                         *pfResid << m_vpEnsemble[ensNum]->GetResidual(j) << endl;
00219                 }
00220                 // single pound symbol for the parameter file and residuals file
00221                 *pfPars << "#" << endl; 
00222                 *pfResid << "#" << endl;
00223         }
00224         
00225         // now do the eigenparameters
00226         // first the minimum (0th member of the ensemble)
00227         for(int m = 0; m < m_iNParameters; m++)
00228         {
00229                 minEigen[m] = 0.0;
00230         }
00231         // compute eigenparameters
00232         for(int j = 0; j < m_iNParameters; j++)
00233         {
00234                 double parameter = log(m_vpEnsemble[0]->GetParameter(j));
00235                 for(int k = 0; k < m_iNParameters; k++)
00236                 {
00237                         minEigen[k] += parameter*(m_pdMatrixV[j][k]);
00238                 }
00239         }
00240         // now print it (first line is just zeroes
00241         for(int j = 0; j < m_iNParameters; j++)
00242         {
00243                 *pfEigP << 0.0 << endl;
00244         }
00245         *pfEigP << endl << endl;
00246 
00247         
00248         for(int i = 1; i < m_vpEnsemble.size(); i++)
00249         {
00250                 // zero the eigenparameters
00251                 for(int m = 0; m < m_iNParameters; m++)
00252                 {
00253                         eigenP[m] = 0.0;
00254                 }
00255                 // compute eigenparameters
00256                 for(int j = 0; j < m_iNParameters; j++)
00257                 {
00258                         double parameter = log(m_vpEnsemble[i]->GetParameter(j));
00259                         for(int k = 0; k < m_iNParameters; k++)
00260                         {
00261                                 eigenP[k] += parameter*(m_pdMatrixV[j][k]);
00262                         }
00263                 }
00264                 // now print the eigenparameters
00265                 for(int j = 0; j < m_iNParameters; j++)
00266                 {
00267                         *pfEigP << eigenP[j] - minEigen[j] << endl;
00268                 }
00269                 *pfEigP << endl << endl;
00270         }       
00271         
00272         delete pfPars;
00273         delete pfQuad;
00274         delete pfCost;
00275         delete pfRatio;
00276         delete pfResid;
00277 
00278         delete [] minEigen;
00279         delete [] eigenP;
00280 
00281         // free memory
00282         delete [] pParCurrent;
00283         delete [] pParTrial;
00284         delete [] pParDelta;
00285         delete [] m_pdHessian[0];
00286         delete [] m_pdHessian;
00287         delete [] m_pdSamplingMatrix[0];
00288         delete [] m_pdSamplingMatrix;
00289         delete [] m_pdMatrixV[0];
00290         delete [] m_pdMatrixV;
00291         delete [] m_pdJacobian[0];
00292         delete [] m_pdJacobian;
00293 
00294         return;
00295 }
00296 
00297 void CStochasticSensitivityAnalysis::ComputeApproximateHessian(double *pParameters, NLLSMinimizable *pNLLS)
00298 {
00299         double transparameter = 0.0;
00300         double saveparameter = 0.0;
00301         double parplus = 0.0;
00302         double parminus = 0.0;
00303         double *fxPlusPlus = new double[m_iNResiduals];
00304         double *fxPlus = new double[m_iNResiduals];
00305         double *fx = new double[m_iNResiduals];
00306         double *fxMinus = new double[m_iNResiduals];
00307         double *fxMinusMinus = new double[m_iNResiduals];
00308 
00309         // compute residuals/cost at f(x)
00310         double chiSq = pNLLS->ObjectiveFunction(pParameters);
00311         for(int k = 0; k < m_iNResiduals; k++)
00312         {
00313                 fx[k] = (pNLLS->GetResiduals())[k];
00314         }
00315 
00316         // USE FOR CENTERED DIFFERENCE (5 pt as well?)
00317         double stepscale = pow(m_dFuncAccuracy,0.3333333);
00318         
00319         // 3 PT CENTERED DIFFERENCE
00320         // loop over parameters
00321         for(int i = 0; i < m_iNParameters; i++)
00322         {
00323                 // transform the parameter and save its old value
00324                 saveparameter = pParameters[i];
00325                 transparameter = m_pFilter->Operator(pParameters[i]);
00326 
00327                 double h = stepscale;
00328                 if(fabs(transparameter) > DBL_EPSILON*1000)
00329                 {
00330                         h = stepscale*fabs(transparameter);
00331                 }
00332 
00333                 parplus = transparameter + h;
00334                 parminus = transparameter - h;
00335                         
00336                 // forward step
00337                 pParameters[i] = m_pFilter->OperatorInverse(transparameter + h);
00338                 double chiSqPlus = pNLLS->ObjectiveFunction(pParameters);
00339                 for(int k = 0; k < m_iNResiduals; k++)
00340                 {
00341                         fxPlus[k] = (pNLLS->GetResiduals())[k];
00342                 }
00343 
00344                 // backward step
00345                 pParameters[i] = m_pFilter->OperatorInverse(transparameter - h);
00346                 double chiSqMinus = pNLLS->ObjectiveFunction(pParameters);
00347                 for(int k = 0; k < m_iNResiduals; k++)
00348                 {
00349                         fxMinus[k] = (pNLLS->GetResiduals())[k];
00350                 }
00351                 
00352                 // reset old parameter
00353                 pParameters[i] = saveparameter;
00354 
00355                 // fill in jacobian
00356                 for(int j = 0; j < m_iNResiduals; j++)
00357                 {
00358                         double delta = (fxPlus[j] - fxMinus[j])/(2*h);
00359                         m_pdJacobian[j][i] = delta;
00360                 }
00361         }
00362         
00363         // compute (1/2) of the scaled hessian (it is symmetric), and fill in 
00364         // the other half
00365         for(int i = 0; i < m_iNParameters; i++)
00366         {
00367                 for(int k = i; k < m_iNParameters; k++)
00368                 {
00369                         m_pdHessian[i][k] = 0.0;
00370                         for(int j = 0; j < m_iNResiduals; j++)
00371                         {
00372                                 m_pdHessian[i][k] += m_pdJacobian[j][i]*m_pdJacobian[j][k];
00373                         }
00374                         m_pdHessian[k][i] = m_pdHessian[i][k];
00375                 }
00376         }
00377 
00378         delete [] fxPlus;
00379         delete [] fx;
00380         delete [] fxMinus;
00381 
00382         cout << "Hessian calculated . . ." << endl;
00383         // write the approximate Hessian to a file for later use
00384         fstream *pfHessian = new fstream("approxH.par",ios::out);
00385         *pfHessian << "##################################################################" << endl;
00386         *pfHessian << "# APPROXIMATE LM HESSIAN AT MINIMUM" << endl;
00387         *pfHessian << "##################################################################" << endl;
00388         for(int row = 0; row < m_iNParameters; row++)
00389         {
00390                 for(int col = 0; col < m_iNParameters; col++)
00391                 {
00392                         *pfHessian << setprecision(10) << m_pdHessian[row][col] << endl;
00393                 }
00394         }
00395         delete pfHessian;
00396 
00397         cout << "Approximate Hessian written to file \"approxH.par\"" << endl;
00398 
00399         // also write the jacobian to a file for later use
00400         fstream *pfJacobian = new fstream("lmjacobian.par",ios::out); 
00401         *pfJacobian << "##################################################################" << endl;
00402         *pfJacobian << "# LEVENBERG-MARQUARDT JACOBIAN AT MINIMUM" << endl;
00403         *pfJacobian << "##################################################################" << endl;
00404         for(int row = 0; row < m_iNResiduals; row++)
00405         {
00406                 for(int col = 0; col < m_iNParameters; col++)
00407                 {
00408                         *pfJacobian << setprecision(10) << m_pdJacobian[row][col] << endl;
00409                 }
00410         }
00411         delete pfJacobian;
00412 
00413         // Hessian has been recalculated, so set the SVD flag to true
00414         m_bSVDFlag = true;
00415 
00416         return;
00417 }
00418 
00419 
00420 void CStochasticSensitivityAnalysis::ComputeTrueHessian(double *pParameters, NLLSMinimizable *pNLLS)
00421 {
00422         double transpi = 0.0;
00423         double transpj = 0.0;
00424         double savepi = 0.0;
00425         double savepj = 0.0;
00426         double chiSq,chiSqp,chiSqm;
00427         double chiSqpp,chiSqmm,chiSqpm,chiSqmp;
00428 
00429         // compute residuals/cost at f(x)
00430         chiSq = pNLLS->ObjectiveFunction(pParameters);
00431 
00432         // DON'T FOOL WITH THIS UNLESS YOU CHANGE DIFFERENCING FORMULA
00433         double stepscale = pow(m_dFuncAccuracy,0.3333333);
00434 
00435         // compute all (nParameters*(nParameters + 1))/2 unique hessian elements
00436         for(int i = 0; i < m_iNParameters; i++)
00437         {
00438                 for(int j = i; j < m_iNParameters; j++)
00439                 {
00440                         // transform parameters i and j and save their old values
00441                         savepi = pParameters[i];
00442                         savepj = pParameters[j];
00443                         transpi = m_pFilter->Operator(pParameters[i]);
00444                         transpj = m_pFilter->Operator(pParameters[j]);
00445 
00446                         // stepsizes - these are optimal for these formulas
00447                         double stepscalei = stepscale;
00448                         double stepscalej = stepscale;
00449                         double hi = stepscalei;
00450                         double hj = stepscalej;
00451                         if(fabs(transpi) > DBL_EPSILON)
00452                         {
00453                                 hi = stepscalei*fabs(transpi);
00454                         }
00455                         if(fabs(transpj) > DBL_EPSILON)
00456                         {
00457                                 hj = stepscalej*fabs(transpj);
00458                         }
00459 
00460                         // lowest order formulas for mixed second derivative
00461 
00462                         if(i == j)
00463                         {
00464                                 pParameters[i] = m_pFilter->OperatorInverse(transpi + hi);
00465                                 chiSqp = pNLLS->ObjectiveFunction(pParameters);
00466                                 pParameters[i] = savepi;
00467                                 pParameters[i] = m_pFilter->OperatorInverse(transpi - hi);
00468                                 chiSqm = pNLLS->ObjectiveFunction(pParameters);
00469                                 pParameters[i] = savepi;
00470 
00471                                 m_pdHessian[i][i] = chiSqp - 2*chiSq + chiSqm;
00472                                 m_pdHessian[i][i] = 0.5*m_pdHessian[i][i]/(hi*hi);
00473                         }
00474                         else
00475                         {                       
00476                                 // f(xi + hi, xj + h)
00477                                 pParameters[i] = m_pFilter->OperatorInverse(transpi + hi);
00478                                 pParameters[j] = m_pFilter->OperatorInverse(transpj + hj);
00479                                 chiSqpp = pNLLS->ObjectiveFunction(pParameters);
00480                                 pParameters[i] = savepi;
00481                                 pParameters[j] = savepj;
00482 
00483                                 // f(xi + hi, xj - hj)
00484                                 pParameters[i] = m_pFilter->OperatorInverse(transpi + hi);
00485                                 pParameters[j] = m_pFilter->OperatorInverse(transpj - hj);
00486                                 chiSqpm = pNLLS->ObjectiveFunction(pParameters);
00487                                 pParameters[i] = savepi;
00488                                 pParameters[j] = savepj;
00489 
00490                                 // f(xi - hi, xj + hj)
00491                                 pParameters[i] = m_pFilter->OperatorInverse(transpi - hi);
00492                                 pParameters[j] = m_pFilter->OperatorInverse(transpj + hj);
00493                                 chiSqmp = pNLLS->ObjectiveFunction(pParameters);
00494                                 pParameters[i] = savepi;
00495                                 pParameters[j] = savepj;
00496         
00497                                 // f(xi - hi, xj - hj)
00498                                 pParameters[i] = m_pFilter->OperatorInverse(transpi - hi);
00499                                 pParameters[j] = m_pFilter->OperatorInverse(transpj - hj);
00500                                 chiSqmm = pNLLS->ObjectiveFunction(pParameters);
00501                                 pParameters[i] = savepi;
00502                                 pParameters[j] = savepj;
00503 
00504 
00505                                 // compute hessian element and fill in its symmetric equivalent
00506                                 m_pdHessian[i][j] = chiSqpp - chiSqpm - chiSqmp + chiSqmm;
00507                                 m_pdHessian[i][j] = 0.5*m_pdHessian[i][j]/(4*hi*hj);
00508                                 m_pdHessian[j][i] = m_pdHessian[i][j];
00509                         }
00510 
00511                 } // end j loop
00512         } // end i loop
00513         
00514         cout << "True hessian calculated . . ." << endl;
00515         
00516         // write the approximate Hessian to a file for later use
00517         fstream *pfHessian = new fstream("trueH.par",ios::out);
00518         *pfHessian << "##################################################################" << endl;
00519         *pfHessian << "# TRUE HESSIAN AT MINIMUM" << endl;
00520         *pfHessian << "##################################################################" << endl;
00521         for(int row = 0; row < m_iNParameters; row++)
00522         {
00523                 for(int col = 0; col < m_iNParameters; col++)
00524                 {
00525                         *pfHessian << setprecision(10) << m_pdHessian[row][col] << endl;
00526                 }
00527         }
00528         delete pfHessian;
00529 
00530         cout << "True Hessian written to file \"trueH.par\"" << endl;
00531         // Hessian has been recalculated, so set the SVD flag to true
00532         m_bSVDFlag = true;
00533 
00534         return;
00535 }
00536 
00537 void CStochasticSensitivityAnalysis::ReadHessian()
00538 {
00539         char linebuf[1000];
00540         // open the file H.par for reading
00541         ifstream *pfHessian = new ifstream("H.par");
00542         // the file has three comment lines (begin with '#'), so read them in 
00543         // and throw them out
00544         for(int i = 0; i < 3; i++)
00545         {
00546                 pfHessian->getline(linebuf,sizeof(linebuf));
00547         }
00548         // now read in the NParameters*NParameters elements
00549         for(int row = 0; row < m_iNParameters; row++)
00550         {
00551                 for(int col = 0; col < m_iNParameters; col++)
00552                 {
00553                         pfHessian->getline(linebuf,sizeof(linebuf));
00554                         std::istringstream *dataStream = new std::istringstream(linebuf);
00555                         *dataStream >> m_pdHessian[row][col];
00556                 }
00557         }
00558         delete pfHessian;
00559         cout << "Hessian read in from file \"H.par\"" << endl;
00560         // set SVD flag to true so Hessian is decomposed in sampling routine
00561         m_bSVDFlag = true;
00562         return;
00563 }
00564 
00565 double CStochasticSensitivityAnalysis::SampleWithoutExplicitHessian(double *pParDelta)
00566 {
00567         double quadratic = 0.0;
00568         double *pRandVec = new double[m_iNParameters];
00569 
00570         // do the decomposition if the flag tells us to do so (basically once during the
00571         // entire analysis)
00572         if(m_bSVDFlag)
00573         {
00574                 cout << "Decomposing Jacobian to find matrix square root . . ." << endl;
00575                 
00576                 // decompose the hessian
00577                 char JOBU = 'A';
00578                 char JOBVT = 'A';
00579                 int M = m_iNResiduals;
00580                 int N = m_iNParameters;
00581                 int minMN = __min(M,N);
00582                 int maxMN = __max(M,N);
00583 
00584                 // matrix to decompose (J) in column-major order
00585                 double *A = new double[m_iNParameters*m_iNResiduals];
00586                 int counter = 0;
00587                 for(int j = 0; j < m_iNParameters; j++)
00588                 {
00589                         for(int i = 0; i < m_iNResiduals; i++)
00590                         {
00591                                 A[counter] = m_pdJacobian[i][j];
00592                                 counter++;
00593                         }
00594                 }
00595                 
00596                 double *D = new double[minMN];
00597                 double *E = new double[minMN];
00598                 double *TAUQ = new double[minMN];
00599                 double *TAUP = new double[minMN];
00600 
00601                 double *PT = new double[maxMN*maxMN];
00602                 double *Q = new double[maxMN*maxMN];
00603                 double *ATEMP = new double[M*N];
00604 
00605                 int LDA = M;
00606                 int LDU = M;
00607                 int LDVT = M;
00608                 int LWORK = 5*maxMN;
00609                 double *WORK = new double[LWORK];
00610                 int INFO = 0;
00611 
00612                 // the input/return arguments to the mkl SVD functions are different depending
00613                 // upon whether the problem is over- or underdetermined
00614 
00615                 if( M >= N )
00616                 {
00617                         cout << "Using decomposition procedure for M >= N . . ." << endl;
00618                         // put the matrix into band diagonal form
00619                         DGEBRD(&M,&N,A,&LDA,D,E,TAUQ,TAUP,WORK,&LWORK,&INFO);
00620                         // XXX debug
00621                         if(INFO != 0)
00622                         {
00623                                 cout << "ERROR : INFO = " << INFO << " in DGEBRD" << endl;
00624                         }
00625                         // now obtain the N leading columns of Q and copy them into our array
00626                         for(int l = 0; l < M*N; l++)
00627                         {
00628                                 ATEMP[l] = A[l];
00629                         }
00630                         char VECT = 'Q';
00631                         int K = N;
00632                         DORGBR(&VECT,&M,&N,&K,ATEMP,&LDA,TAUQ,WORK,&LWORK,&INFO);
00633                         // XXX debug
00634                         if(INFO != 0)
00635                         {
00636                                 cout << "ERROR : INFO = " << INFO << " in DORGBR (U)" << endl;
00637                         }
00638                         for(int l = 0; l < M*N; l++)
00639                         {
00640                                 Q[l] = ATEMP[l];
00641                         }
00642                 
00643                         // now obtain PT
00644                         for(int l = 0; l < M*N; l++)
00645                         {
00646                                 ATEMP[l] = A[l];
00647                         }
00648                         VECT = 'P';
00649                         K  = M;
00650                         DORGBR(&VECT,&N,&N,&K,ATEMP,&LDA,TAUP,WORK,&LWORK,&INFO);
00651                         if(INFO != 0)
00652                         {
00653                                 cout << "ERROR : INFO = " << INFO << " in DORGBR (VT)" << endl;
00654                         }
00655                         for(int l = 0; l < N*N; l++)
00656                         {
00657                                 PT[l] = ATEMP[l];
00658                         }
00659 
00660                         // obtain the SVD decomposition of A, using the matrices found above
00661                         char UPLO = 'U';
00662                         int NCVT = N;
00663                         int NRU = M;
00664                         int NCC = 0;
00665                         double *C = 0;
00666                         int LDC = 1;
00667                         LDVT = N;
00668                         LDU = M;
00669                         DBDSQR(&UPLO,&N,&NCVT,&NRU,&NCC,D,E,PT,&LDVT,Q,&LDU,C,&LDC,WORK,&INFO);
00670                         if(INFO != 0)
00671                         {
00672                                 cout << "ERROR : INFO = " << INFO << " in DBDSQR (VT)" << endl;
00673                         }
00674 
00675                         cout << "Matrix decomposed. Condition number is " << D[0]/D[N-1] << endl;
00676                                                 
00677                         // scroll through the singular values and find the ones whose inverses will be 
00678                         // huge and set them to zero; also, load up the array of singular values that 
00679                         // we store.  M >= N, so the number of singular values returned equals the
00680                         // number of parameters
00681                         double cutoff = D[0]*1.0e-06;
00682                         for(int i = 0; i < N; i++)
00683                         {
00684                                 m_pdSingularValues[i] = D[i];
00685                                 if(D[i] < cutoff) {D[i] = 0;}
00686                                 else {D[i] = 1.0/D[i];}
00687                         }
00688                         // D[i] now equals 1/D[i]
00689 
00690                         // now fill in the sampling matrix ("square root" of the Hessian)
00691                         for(int i = 0; i < m_iNParameters; i++)
00692                         {
00693                                 for(int j = 0; j < m_iNParameters; j++)
00694                                 {
00695                                         m_pdSamplingMatrix[i][j] = PT[i*m_iNParameters + j];
00696                                         m_pdMatrixV[i][j] = PT[i*m_iNParameters + j];
00697                                 }
00698                         }                       
00699                         
00700                         // open up a file to which we dump the eigenvectors (cols. of V)
00701                         fstream *pfEigs = new fstream("matrixV.par",ios::out);
00702                         *pfEigs << "##################################################################" << endl;
00703                         *pfEigs << "# EIGENVECTORS OF H - \"#\" IS SEPARATOR" << endl;
00704                         *pfEigs << "##################################################################" << endl;
00705                         for(int i = 0; i < m_iNParameters; i++)
00706                         {
00707                                 for(int j = 0; j < m_iNParameters; j++)
00708                                 {
00709                                         *pfEigs << setprecision(10) << m_pdSamplingMatrix[j][i] << endl;
00710                                 }
00711                                 *pfEigs << "#" << endl;
00712                         }
00713                         delete pfEigs;
00714 
00715                         for(int i = 0; i < m_iNParameters; i++)
00716                         {
00717                                 for(int j = 0; j < m_iNParameters; j++)
00718                                 {
00719                                         m_pdSamplingMatrix[i][j] = m_pdSamplingMatrix[i][j]*D[j];
00720                                 }
00721                         }
00722                 
00723                 }
00724                 else
00725                 {
00726                         cout << "Using decomposition procedure for M < N . . ." << endl;
00727                         // put the matrix into band diagonal form
00728                         DGEBRD(&M,&N,A,&LDA,D,E,TAUQ,TAUP,WORK,&LWORK,&INFO);
00729                         // XXX debug
00730                         if(INFO != 0)
00731                         {
00732                                 cout << "ERROR : INFO = " << INFO << " in DGEBRD" << endl;
00733                         }
00734                         // now obtain the whole M X M matrix Q and copy it into our array
00735                         for(int l = 0; l < M*N; l++)
00736                         {
00737                                 ATEMP[l] = A[l];
00738                         }
00739                         char VECT = 'Q';
00740                         int K = N;
00741                         DORGBR(&VECT,&M,&M,&K,ATEMP,&LDA,TAUQ,WORK,&LWORK,&INFO);
00742                         // XXX debug
00743                         if(INFO != 0)
00744                         {
00745                                 cout << "ERROR : INFO = " << INFO << " in DORGBR (U)" << endl;
00746                         }
00747                         for(int l = 0; l < M*M; l++)
00748                         {
00749                                 Q[l] = ATEMP[l];
00750                         }
00751                 
00752                         // now obtain the m leading rows of PT
00753                         for(int l = 0; l < M*N; l++)
00754                         {
00755                                 ATEMP[l] = A[l];
00756                         }
00757                         VECT = 'P';
00758                         K  = M;
00759                         DORGBR(&VECT,&M,&N,&K,ATEMP,&LDA,TAUP,WORK,&LWORK,&INFO);
00760                         if(INFO != 0)
00761                         {
00762                                 cout << "ERROR : INFO = " << INFO << " in DORGBR (VT)" << endl;
00763                         }
00764                         for(int l = 0; l < M*N; l++)
00765                         {
00766                                 PT[l] = ATEMP[l];
00767                         }
00768 
00769                         // obtain the SVD decomposition of A, using the matrices found above
00770                         char UPLO = 'L';
00771                         int NCVT = N;
00772                         int NRU = M;
00773                         int NCC = 0;
00774                         double *C = 0;
00775                         int LDC = 1;
00776                         LDU = M;
00777                         LDVT = M;
00778                         DBDSQR(&UPLO,&M,&NCVT,&NRU,&NCC,D,E,PT,&LDVT,Q,&LDU,C,&LDC,WORK,&INFO);
00779                         if(INFO != 0)
00780                         {
00781                                 cout << "ERROR : INFO = " << INFO << " in DBDSQR" << endl;
00782                         }
00783 
00784                         // in this case (M < N) there is are N - M guaranteed 0 singular values, so
00785                         // we have to remember that when computing matrix products and such
00786                         cout << "Matrix decomposed. Submatrix condition number is " << D[0]/D[M-1] << endl;
00787                         // scroll through the singular values and find the ones whose inverses will be 
00788                         // huge and set them to zero; also, load up the array of singular values that 
00789                         // we store.  the second loop fills in the SVs that are exactly zero                    
00790                         double cutoff = D[0]*1.0e-06;
00791                         for(int i = 0; i < M; i++)
00792                         {
00793                                 m_pdSingularValues[i] = D[i];
00794                                 if(D[i] < cutoff) {D[i] = 0;}
00795                                 else {D[i] = 1.0/D[i];}
00796                         }
00797                         for(int i = M; i < N; i++)
00798                         {
00799                                 m_pdSingularValues[i] = 0.0;
00800                         }
00801                         // D[i] now equals 1/D[i], except for the zeroes at the end
00802                         
00803                         // now fill in the sampling matrix ("square root" of the Hessian); the matrix
00804                         // will have N - M zero columns, corresponding to the N - M zero rows of 
00805                         // PT
00806                         for(int i = 0; i < m_iNParameters; i++)
00807                         {
00808                                 for(int j = 0; j < m_iNResiduals; j++)
00809                                 {
00810                                         m_pdSamplingMatrix[i][j] = PT[i*m_iNResiduals + j];
00811                                         m_pdMatrixV[i][j] = PT[i*m_iNResiduals + j];
00812                                 }
00813                                 for(int j = m_iNResiduals; j < m_iNParameters; j++)
00814                                 {
00815                                         m_pdSamplingMatrix[i][j] = 0.0;
00816                                         m_pdMatrixV[i][j] = 0.0;
00817                                 }
00818                         }
00819 
00820                         // open up a file to which we dump the eigenvectors (cols. of V)
00821                         fstream *pfEigs = new fstream("matrixV.par",ios::out);
00822                         *pfEigs << "##################################################################" << endl;
00823                         *pfEigs << "# EIGENVECTORS OF H - \"#\" IS SEPARATOR" << endl;
00824                         *pfEigs << "##################################################################" << endl;
00825                         for(int i = 0; i < m_iNParameters; i++)
00826                         {
00827                                 for(int j = 0; j < m_iNParameters; j++)
00828                                 {
00829                                         *pfEigs << setprecision(10) << m_pdSamplingMatrix[j][i] << endl;
00830                                 }
00831                                 *pfEigs << "#" << endl;
00832                         }
00833                         delete pfEigs;
00834 
00835                         for(int i = 0; i < m_iNParameters; i++)
00836                         {
00837                                 for(int j = 0; j < m_iNResiduals; j++)
00838                                 {
00839                                         m_pdSamplingMatrix[i][j] = m_pdSamplingMatrix[i][j]*D[j];
00840                                 }
00841                         }
00842                 }
00843                 
00844                 delete [] A;
00845                 delete [] ATEMP;
00846                 delete [] TAUQ;
00847                 delete [] TAUP;
00848                 delete [] Q;
00849                 delete [] PT;
00850                 delete [] D;
00851                 delete [] E;
00852                 delete [] WORK;
00853 
00854                 
00855                 // don't compute the matrix square root every time
00856                 m_bSVDFlag = false;
00857                 // dump the singular values for the matrix to a file
00858                 fstream *svd = new fstream("singval.par",ios::out);
00859                 *svd << "##################################################################" << endl;
00860                 *svd << "# SINGULAR VALUES OF J AT FITTED PARAMETERS" << endl;
00861                 *svd << "##################################################################" << endl;
00862                 for(int i = 0; i < m_iNParameters; i++)
00863                 {
00864                         *svd << setprecision(10) << m_pdSingularValues[i] << endl;
00865                 }
00866                 delete svd;
00867                 
00868         }
00869         // sample the gaussian in either case
00870 
00871         // draw random numbers r ~ N(0,1)
00872         for(int i = 0; i < m_iNParameters; i++)
00873         {
00874                 pRandVec[i] = m_pRNG->gaussian(1.0);
00875         }
00876         
00877         // now rescale by the sampling matrix, including the appropriate temperature
00878         // factor
00879         // RIGHT MULTIPLICATION
00880         for(int i = 0; i < m_iNParameters; i++)
00881         {
00882                 pParDelta[i] = 0.0;
00883                 for(int j = 0; j < m_iNParameters; j++)
00884                 {
00885                         pParDelta[i] += sqrt(m_dRescale)*m_pdSamplingMatrix[i][j]*pRandVec[j];
00886                 }
00887         }
00888                 
00889         // calculate and return the value of the quadratic form, simply (1/2)xT.JT.J.x
00890         // = (1/2) (Jx)T (Jx)
00891         double *tempVec = new double[m_iNResiduals];    
00892 
00893         for(int i = 0; i < m_iNResiduals; i++)
00894         {
00895                 tempVec[i] = 0.0;
00896                 for(int j = 0; j < m_iNParameters; j++)
00897                 {
00898                         tempVec[i] += m_pdJacobian[i][j]*pParDelta[j];
00899                 }
00900                 quadratic += 0.5*tempVec[i]*tempVec[i];
00901         }
00902         
00903         delete [] tempVec;
00904         delete [] pRandVec;
00905 
00906         return quadratic;
00907 }
00908 
00909 double CStochasticSensitivityAnalysis::SampleGaussian(double *pParDelta)
00910 {
00911         double quadratic = 0.0;
00912         double *pRandVec = new double[m_iNParameters];
00913 
00914         // do the decomposition if the flag tells us to do so (basically once during the
00915         // entire analysis)
00916         if(m_bSVDFlag)
00917         {
00918                 cout << "Decomposing Hessian to find matrix square root . . ." << endl;
00919                 
00920                 // decompose the hessian
00921                 char JOBU = 'A';
00922                 char JOBVT = 'A';
00923                 int M = m_iNParameters;
00924                 int N = m_iNParameters;
00925         
00926                 // matrix A in column-major order
00927                 double *A = new double[m_iNParameters*m_iNParameters];
00928                 int counter = 0;
00929                 for(int j = 0; j < m_iNParameters; j++)
00930                 {
00931                         pParDelta[j] = 0.0;
00932                         for(int i = 0; i < m_iNParameters; i++)
00933                         {
00934                                 A[counter] = m_pdHessian[i][j];
00935                                 counter++;
00936                         }
00937                 }
00938         
00939                 double *D = new double[m_iNParameters];
00940                 double *E = new double[m_iNParameters];
00941                 double *TAUQ = new double[m_iNParameters];
00942                 double *TAUP = new double[m_iNParameters];
00943 
00944                 double *PT = new double[m_iNParameters*m_iNParameters];
00945                 double *Q = new double[m_iNParameters*m_iNParameters];
00946                 double *ATEMP = new double[m_iNParameters*m_iNParameters];
00947 
00948                 int LDA = m_iNParameters;
00949                 int LDU = m_iNParameters;
00950                 int LDVT = m_iNParameters;
00951                 int LWORK = 5*M;
00952                 double *WORK = new double[LWORK];
00953                 int INFO = 0;
00954 
00955                 // put the matrix into band diagonal form
00956                 DGEBRD(&M,&N,A,&LDA,D,E,TAUQ,TAUP,WORK,&LWORK,&INFO);   
00957                 // now obtain Q and copy it into our array
00958                 for(int l = 0; l < m_iNParameters*m_iNParameters; l++)
00959                 {
00960                         ATEMP[l] = A[l];
00961                 }
00962                 char VECT = 'Q';
00963                 int K = M;
00964                 DORGBR(&VECT,&M,&N,&K,ATEMP,&LDA,TAUQ,WORK,&LWORK,&INFO);
00965                 for(int l = 0; l < m_iNParameters*m_iNParameters; l++)
00966                 {
00967                         Q[l] = ATEMP[l];
00968                 }
00969                 // now obtain PT
00970                 for(int l = 0; l < m_iNParameters*m_iNParameters; l++)
00971                 {
00972                         ATEMP[l] = A[l];
00973                 }
00974                 VECT = 'P';
00975                 DORGBR(&VECT,&M,&N,&K,ATEMP,&LDA,TAUP,WORK,&LWORK,&INFO);
00976                 for(int l = 0; l < m_iNParameters*m_iNParameters; l++)
00977                 {
00978                         PT[l] = ATEMP[l];
00979                 }
00980 
00981                 // obtain the SVD decomposition of A, using the matrices found above
00982                 char UPLO = 'U';
00983                 int NCVT = m_iNParameters;
00984                 int NRU = m_iNParameters;
00985                 int NCC = 0;
00986                 double *C = 0;
00987                 int LDC = 1;
00988                 LDVT = m_iNParameters;
00989                 LDU = m_iNParameters;
00990                 DBDSQR(&UPLO,&N,&NCVT,&NRU,&NCC,D,E,PT,&LDVT,Q,&LDU,C,&LDC,WORK,&INFO);
00991                 cout << "Matrix decomposed. Condition number is " << D[0]/D[m_iNParameters-1] << endl;
00992 
00993                 // scroll through the singular values and find the ones whose inverses will be 
00994                 // huge and set them to zero; also, load up the array of singular values that 
00995                 // we store
00996                 //double cutoff = D[0]*1.0e-06;
00997                 double cutoff = 1.0;
00998                 for(int i = 0; i < m_iNParameters; i++)
00999                 {
01000                         m_pdSingularValues[i] = D[i];
01001                         if(1.0/D[i] > cutoff) {D[i] = cutoff;}
01002                         else {D[i] = 1.0/D[i];}
01003                 }
01004                 // D[i] now equals 1/D[i]
01005 
01006                 // now fill in the sampling matrix ("square root" of the Hessian)
01007                 for(int i = 0; i < m_iNParameters; i++)
01008                 {
01009                         for(int j = 0; j < m_iNParameters; j++)
01010                         {
01011                                 m_pdSamplingMatrix[i][j] = PT[i*m_iNParameters + j];
01012                                 m_pdMatrixV[i][j] = PT[i*m_iNParameters + j];
01013                         }
01014                 }
01015 
01016                 // open up a file to which we dump the eigenvectors (cols. of V)
01017                 fstream *pfEigs = new fstream("matrixV.par",ios::out);
01018                 *pfEigs << "##################################################################" << endl;
01019                 *pfEigs << "# EIGENVECTORS OF H - \"#\" IS SEPARATOR" << endl;
01020                 *pfEigs << "##################################################################" << endl;
01021                 for(int i = 0; i < m_iNParameters; i++)
01022                 {
01023                         for(int j = 0; j < m_iNParameters; j++)
01024                         {
01025                                 *pfEigs << setprecision(10) << m_pdSamplingMatrix[j][i] << endl;
01026                         }
01027                         *pfEigs << "#" << endl;
01028                 }
01029                 delete pfEigs;
01030                 
01031                 for(int i = 0; i < m_iNParameters; i++)
01032                 {
01033                         for(int j = 0; j < m_iNParameters; j++)
01034                         {
01035                                 m_pdSamplingMatrix[i][j] = m_pdSamplingMatrix[i][j]*sqrt(D[j]);
01036                         }
01037                 }
01038 
01039                 // Square of the norm of the last row of the sampling matrix
01040                 for(int i = 0; i < m_iNParameters; i++)
01041                 {
01042                         m_dMNormSquared += m_pdSamplingMatrix[m_iNParameters-1][i]*m_pdSamplingMatrix[m_iNParameters-1][i];
01043                 }
01044                 
01045                 delete [] A;
01046                 delete [] ATEMP;
01047                 delete [] TAUQ;
01048                 delete [] TAUP;
01049                 delete [] Q;
01050                 delete [] PT;
01051                 delete [] D;
01052                 delete [] E;
01053                 delete [] WORK;
01054 
01055                 // don't compute the matrix square root every time
01056                 m_bSVDFlag = false;
01057                 // dump the singular values for the matrix to a file
01058                 fstream *svd = new fstream("singval.par",ios::out);
01059                 *svd << "##################################################################" << endl;
01060                 *svd << "# SINGULAR VALUES OF H AT FITTED PARAMETERS" << endl;
01061                 *svd << "##################################################################" << endl;
01062                 for(int i = 0; i < m_iNParameters; i++)
01063                 {
01064                         *svd << setprecision(10) << m_pdSingularValues[i] << endl;
01065                 }
01066                 delete svd;
01067         }
01068         // sample the gaussian in either case
01069 
01070         // draw random numbers r ~ N(0,1)
01071         for(int i = 0; i < m_iNParameters; i++)
01072         {
01073                 pRandVec[i] = m_pRNG->gaussian(1.0);
01074         }
01075         
01076         // now rescale by the sampling matrix, including the appropriate
01077         // reweighting by the control parameter.
01078         
01079         // RIGHT MULTIPLICATION
01080         for(int i = 0; i < m_iNParameters; i++)
01081         {
01082                 pParDelta[i] = 0.0;
01083                 for(int j = 0; j < m_iNParameters; j++)
01084                 {
01085                         pParDelta[i] += sqrt(m_dRescale)*m_pdSamplingMatrix[i][j]*pRandVec[j];
01086                 }
01087         }
01088         
01089         
01090         // calculate and return the value of the quadratic form
01091         double *tempPar = new double[m_iNParameters];
01092         for(int i = 0; i < m_iNParameters; i++)
01093         {
01094                 tempPar[i] = 0.0;
01095                 for(int j = 0; j < m_iNParameters; j++)
01096                 {
01097                         tempPar[i] += m_pdHessian[i][j]*pParDelta[j];
01098                 }
01099         }
01100         for(int i = 0; i < m_iNParameters; i++)
01101         {
01102                 quadratic += pParDelta[i]*tempPar[i];
01103         }
01104         delete [] tempPar;
01105         delete [] pRandVec;
01106 
01107         return quadratic;
01108 }
01109 
01110 
01111 bool CStochasticSensitivityAnalysis::AcceptMove(double quadratic, double deltaCost)
01112 {
01113         double arg = (-1.0*deltaCost)/m_dTemperature;
01114 
01115         if(arg >= 0.0)
01116         {
01117                 // accept the move
01118                 return true;
01119         }
01120         else
01121         {
01122                 double p = m_pRNG->uniform();
01123                 if(p < exp(arg))
01124                 {
01125                         return true;
01126                 }
01127                 else
01128                 {
01129                         return false;
01130                 }
01131         }
01132 }
01133 
01134 void CStochasticSensitivityAnalysis::AddParametersToEnsemble(double *pParSave, double *resid, double quadratic, double cost)
01135 {
01136         m_vpEnsemble.push_back(new CEnsembleMember(m_iNParameters,pParSave,m_iNResiduals,resid,quadratic,cost));
01137 }
01138 
01139 void CStochasticSensitivityAnalysis::ConstructTrialMove(double *pParTrial, double *pParCurrent, double *pParDelta)
01140 {
01141         m_pFilter->ForwardTransformation(pParCurrent,m_iNParameters);
01142         for(int j = 0; j < m_iNParameters; j++)
01143         {
01144                 pParTrial[j] = pParCurrent[j] + pParDelta[j];
01145         }
01146         m_pFilter->BackwardTransformation(pParCurrent,m_iNParameters);
01147         m_pFilter->BackwardTransformation(pParTrial,m_iNParameters);
01148 }
01149