---

LevenbergMarquardtMinimizer.cpp

Go to the documentation of this file.

// LevenbergMarquardtMinimizer.cpp: implementation of the CLevenbergMarquardtMinimizer class.
//

#include "LevenbergMarquardtMinimizer.h"

// Construction/Destruction

CLevenbergMarquardtMinimizer::CLevenbergMarquardtMinimizer()
{

}

CLevenbergMarquardtMinimizer::~CLevenbergMarquardtMinimizer()
{

}

// Adds to each diagonal element of the Hessian the scaleParameter, i.e. 
// after calling RescaleDiagonal(alpha) you get
//                              H[i][i] = H[i][i] + alpha; 

void CLevenbergMarquardtMinimizer::RescaleDiagonal(double scaleParameter)
{
        for(int i = 0; i < nParameters; i++)
        {
                m_pdAlpha[i][i] = m_pdAlpha[i][i] + scaleParameter;
        }
}

// Restores the diagonal elements of the Hessian to the contents of
// the array m_pdDiagonal

void CLevenbergMarquardtMinimizer::RestoreDiagonal()
{
        int i;
        for(i = 0; i < nParameters; i++)
        {
                m_pdAlpha[i][i] = m_pdDiagonal[i];
        }
}

// If parameters are required to be positive semidefinite, this function
// ensures they remain so before the objective function is evaluated

void CLevenbergMarquardtMinimizer::CheckBounds(double *parameters)
{
        for(int i = 0; i < nParameters; i++)
        {
                if(parameters[i] < 0.0) 
                {
                        parameters[i] = 0.0;
                        cout << "Parameter " << i << " hit a bound, truncating . . ." << endl;
                }
        }
}

void CLevenbergMarquardtMinimizer::ObtainTrialParameters(double *currentP, double *deltaP, double *trialP)
{
        int j;
        // forward transform the current parameters, then add and back transform
        m_pFilter->ForwardTransformation(currentP,nParameters);
        for(j = 0; j < nParameters; j++)
        {
                trialP[j] = currentP[j] + deltaP[j];
        }
        m_pFilter->BackwardTransformation(trialP,nParameters);
        m_pFilter->BackwardTransformation(currentP,nParameters);
        if(m_bPositiveSemiDef) {this->CheckBounds(trialP);}
/*      // XXX debug
        for(j = 0; j < nParameters; j++)
        {
                cout << trialP[j] << endl;
        }
*/
}

void CLevenbergMarquardtMinimizer::AcceptParameters(double *parameters, double *trialP, double newCost)
{
        int i;
        for(i = 0; i < nParameters; i++)
        {
                parameters[i] = trialP[i];
        }

        // write the accepted parameters to a file, because numerically 
        // computing the hessian is extremely expensive
        fstream *pftemp = new fstream("lstsqrtemp.par",ios::out);
        *pftemp << "##################################################################" << endl;
        *pftemp << "# BEST PARAMETERS" << endl;
        *pftemp << "# COST : " << newCost << endl;
        *pftemp << "##################################################################" << endl;
        for(i = 0; i < nParameters; i++)
        {
                *pftemp << parameters[i] << endl;
        }
        delete pftemp;
        // END WRITE
}

bool CLevenbergMarquardtMinimizer::ComputeChiSqTol(double oldCost, double newCost)
{
        if(m_bCheckChi == false) {return false;}
        double fracDeltaC = (oldCost - newCost)/oldCost;
        if(fabs(fracDeltaC) < m_dChiSqTol)
        {
                return true;
        }
        return false;
}

bool CLevenbergMarquardtMinimizer::ComputeGradTol()
{
        if(m_bCheckGrad == false) {return false;}
        double gradNorm = 0.0;
        for(int i = 0; i < nParameters; i++)
        {
                gradNorm += m_pdGrad[i]*m_pdGrad[i];
        }
        gradNorm = sqrt(gradNorm);
        if(gradNorm < m_dGradTol)
        {
                return true;
        }
        return false;
}

bool CLevenbergMarquardtMinimizer::ComputeParTol(double *oldP, double *deltaP)
{
        if(m_bCheckPar == false) {return false;}
        double parTol = 0.0;
        for(int i = 0; i < nParameters; i++)
        {
                double tempPar = m_pFilter->Operator(oldP[i]);
                parTol = fabs(deltaP[i])/(m_dTau + fabs(tempPar));
                if(parTol > m_dParTol)
                {
                        return false;
                }
        }
        return true;
}

---