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

Go to the documentation of this file.
```00001 // LeastSquaresADMinimizer.cpp: implementation of the LeastSquaresADMinimizer class.
00002 //
00004
00006
00007
00009 // Construction/Destruction
00011
00013 {
00014
00015 }
00016
00018 {
00019
00020 }
00021
00022
00024 // Compute approximate Hessian for use in Least Squares algs.
00025 // The hessian is symmetric so we only need compute 1/2 of it then
00026 // copy the rest.  This routine uses the ADOL-C Automatic Differentiation
00027 // package to compute the Jacobian.
00029
00030 void CLeastSquaresADMinimizer::ComputeDerivativeInformation(double *parameters, NLLSMinimizable *nlls)
00031 {
00032         int i,j,k,nR;
00033         double transparameter = 0.0;
00034         double saveparameter = 0.0;
00035         double *fxPlus = new double[m_iNResiduals];
00036         double *fx = new double[m_iNResiduals];
00037         double *fxMinus = new double[m_iNResiduals];
00038         double **jacobian = new double *[m_iNResiduals];
00039         jacobian[0] = new double[m_iNResiduals*nParameters];
00040         for(nR = 1; nR < m_iNResiduals; nR++)
00041         {
00042                 jacobian[nR] = &(jacobian[0][nR*nParameters]);
00043         }
00044
00045         double chiSq = nlls->ObjectiveFunction(parameters);
00046         for(k = 0; k < m_iNResiduals; k++)
00047         {
00048                 fx[k] = (nlls->GetResiduals())[k];
00049         }
00050
00051         // loop over parameters
00052         for(i = 0; i < nParameters; i++)
00053         {
00054                 // transform the parameter and save its old value
00055                 saveparameter = parameters[i];
00056                 transparameter = m_pFilter->Operator(parameters[i]);
00057
00058                 double h = pow(m_dFuncAccuracy,0.333333)*fabs(transparameter);
00059
00060                 // if statement warns of tiny stepsizes
00061                 if(h < DBL_EPSILON)
00062                 {
00063                         cout << "WARNING:Parmeter step small in ImprovedLevenbergMarquardtMinimizer::ComputeDerivativeInformation !" << endl;
00064                         if(transparameter < DBL_MIN)
00065                         {
00066                                 cout << "WARNING:Zero parameter detected, truncating . . ." << endl;
00067                                 // pretend like the parameter was really small
00068                                 h = pow(m_dFuncAccuracy,0.333333)*fabs(m_pFilter->Operator(parameters[i] + DBL_EPSILON));
00069                         }
00070                 }
00071
00072
00073                 // forward step
00074                 parameters[i] = m_pFilter->OperatorInverse(transparameter + h);
00075                 double chiSqPlus = nlls->ObjectiveFunction(parameters);
00076                 for(k = 0; k < m_iNResiduals; k++)
00077                 {
00078                         fxPlus[k] = (nlls->GetResiduals())[k];
00079                 }
00080                 // backward step
00081                 parameters[i] = m_pFilter->OperatorInverse(transparameter - h);
00082                 double chiSqMinus = nlls->ObjectiveFunction(parameters);
00083                 for(k = 0; k < m_iNResiduals; k++)
00084                 {
00085                         fxMinus[k] = (nlls->GetResiduals())[k];
00086                 }
00087
00088                 parameters[i] = saveparameter;
00089
00090                 // fill in jacobian
00091                 for(j = 0; j < m_iNResiduals; j++)
00092                 {
00093                         double delta = (fxPlus[j] - fxMinus[j])/(2*h);
00094                         jacobian[j][i] = delta;
00095                 }
00096         }
00097
00098         // compute -J^T f and store the ForceMatrix in the process
00099         for(i = 0; i < nParameters; i++)
00100         {
00102                 for(int j = 0; j < m_iNResiduals; j++)
00103                 {
00105                         m_pdForceMatrix[i][j] = -2.0*jacobian[j][i]*fx[j];
00106                 }
00107         }
00108
00109         // write the force matrix
00110         Notify();
00111
00112         // compute the scaled hessian (it is symmetric)
00113         for(i = 0; i < nParameters; i++)
00114         {
00115                 for(k = i; k < nParameters; k++)
00116                 {
00117                         m_pdAlpha[i][k] = 0.0;
00118                         for(j = 0; j < m_iNResiduals; j++)
00119                         {
00120                                 m_pdAlpha[i][k] += 2.0*jacobian[j][i]*jacobian[j][k];
00121                         }
00122                         m_pdAlpha[k][i] = m_pdAlpha[i][k];
00123                 }
00124         }
00125
00126         // Save the diagonal elements
00127         this->SaveDiagonal();
00128
00129         delete [] fxPlus;
00130         delete [] fx;
00131         delete [] fxMinus;
00132         delete [] jacobian[0];
00133         delete [] jacobian;
00134
00135         cout << "Hessian calculated . . ." << endl;
00136
00137         return;
00138 }
00139
00141 // Saves the current diagonal elements of the Hessian in the member
00142 // array m_pdDiagonal
00144