---

NelderMeadSimplexMinimizer.cpp

Go to the documentation of this file.

#include "NelderMeadSimplexMinimizer.h"

CNelderMeadSimplexMinimizer::CNelderMeadSimplexMinimizer()
{
}

CNelderMeadSimplexMinimizer::CNelderMeadSimplexMinimizer(CParameterFilter *pFilter, int nIterations, double lambda, double tol)
{
        m_pFilter = pFilter;
        m_dLambda = 1.0;
        m_dTol = 1.0e-06;
        m_pMO = new CMatrixOperations();
        m_iNIterations = nIterations;
}

CNelderMeadSimplexMinimizer::~CNelderMeadSimplexMinimizer(void)
{
}

void CNelderMeadSimplexMinimizer::InitializeSimplex(double *parameters)
{
        int i;
        m_pFilter->ForwardTransformation(parameters,nParameters);
        // fill in the initial point
        m_pMO->ElementCopy(parameters,m_vpdSimplex[0],nParameters);
        // fill in the rest
        for(i = 1; i < m_vpdSimplex.size(); i++)
        {
                m_pMO->ElementCopy(parameters,m_vpdSimplex[i],nParameters);
                m_vpdSimplex[i][i-1] += m_dLambda;
        }
        m_pFilter->BackwardTransformation(parameters,nParameters);
}

void CNelderMeadSimplexMinimizer::EvaluateSimplex()
{
        int i;
        for(i = 0; i < m_vpdSimplex.size(); i++)
        {
                m_pFilter->BackwardTransformation(m_vpdSimplex[i],nParameters);
                m_pdSimplexCosts[i] = m_pMinimizable->ObjectiveFunction(m_vpdSimplex[i]);
                m_pFilter->ForwardTransformation(m_vpdSimplex[i],nParameters);
        }
}

void CNelderMeadSimplexMinimizer::ComputePointSum()
{
        int i,j;

        for(j = 0; j < nParameters; j++)
        {
                double sum = 0.0;
                for(i = 0; i < m_vpdSimplex.size(); i++)
                {
                        sum += m_vpdSimplex[i][j];
                }
                m_pdPointSum[j] = sum;
        }
}

double CNelderMeadSimplexMinimizer::TryPoint(double exfac)
{
        int i,j;
        double *pTrial = new double[nParameters];
        double fac1 = (1.0 - exfac)/nParameters;
        double fac2 = fac1 - exfac;
        for(j = 0; j < nParameters; j++)
        {
                // extrapolate
                pTrial[j] = m_pdPointSum[j]*fac1 - m_vpdSimplex[m_iHi][j]*fac2;
        }
        // try the point
        m_pFilter->BackwardTransformation(pTrial,nParameters);
        double costTrial = m_pMinimizable->ObjectiveFunction(pTrial);
        m_pFilter->ForwardTransformation(pTrial,nParameters);
        // if it's less than the highest point, accept the point
        if(costTrial < m_pdSimplexCosts[m_iHi])
        {
                m_pdSimplexCosts[m_iHi] = costTrial;
                for(j = 0; j < nParameters; j++)
                {
                        m_pdPointSum[j] += pTrial[j] - m_vpdSimplex[m_iHi][j];
                        m_vpdSimplex[m_iHi][j] = pTrial[j];
                }
        }
        delete [] pTrial;

        return costTrial;
}

double CNelderMeadSimplexMinimizer::Minimize(double *parameters, Minimizable *minimizable)
{
        int i,j,iteration;
        double cTry;
        double cBest;
        bool restartFlag = false;
        m_pMinimizable = minimizable;
        nParameters = m_pMinimizable->GetNParameters();

        
        // allocate memory
        for(i = 0; i < nParameters + 1; i++)
        {
                m_vpdSimplex.push_back(new double[nParameters]);
        }
        m_pdSimplexCosts = new double[nParameters + 1];
        m_pdPointSum = new double[nParameters];
        
        // initialize simplex
        this->InitializeSimplex(parameters);
        // evaluate costs
        this->EvaluateSimplex();
        // compute initial pointsum
        this->ComputePointSum();
        
        iteration = 0;
        // loop while not converged (outer loop) and less than max iterations
        while(iteration < m_iNIterations)
        {
                // find the highest, next highest, and lowest points
                m_iLo = 0;
                if(m_pdSimplexCosts[1] > m_pdSimplexCosts[0]) 
                {
                        m_iHi = 1;
                        m_iNextHi = 0;
                }
                else
                {
                        m_iHi = 0;
                        m_iNextHi = 1;
                }
                for(i = 0; i < nParameters + 1; i++)
                {
                        if(m_pdSimplexCosts[i] <= m_pdSimplexCosts[m_iLo]) {m_iLo = i;}
                        if(m_pdSimplexCosts[i] > m_pdSimplexCosts[m_iHi])
                        {
                                m_iNextHi = m_iHi;
                                m_iHi = i;
                        }
                        else if((m_pdSimplexCosts[i] > m_pdSimplexCosts[m_iNextHi]) && (i != m_iHi))
                        {
                                m_iHi = i;
                        }
                } // end hi, lo loop

                cout << "Best cost after iteration " << iteration << " : " << m_pdSimplexCosts[m_iLo] << endl; 
                
                // check the function tolerance to see if we jump out
                // XXX remember a restart!
                if(CheckFuncTol()) 
                {
                        if(restartFlag == true)
                        {
                                cout << "Relative Function Tolerance Achieved, Terminating . . ." << endl;
                                break;
                        }
                        else
                        {
                                cout << "Minimum suspected, restarting . . ." << endl;
                                m_pFilter->BackwardTransformation(m_vpdSimplex[m_iLo],nParameters);
                                this->InitializeSimplex(m_vpdSimplex[m_iLo]);
                                this->EvaluateSimplex();
                                this->ComputePointSum();
                                restartFlag = true;
                                continue;
                        }
                }
                // reflect from the high point - this starts a new iteration
                cTry = TryPoint(-1.0);
                if(cTry <= m_pdSimplexCosts[m_iLo])
                {
                        // XXX debug
                        cout << "\tTrying additional extrapolation . . . " << endl;
                        // result is better than the best point, so try an additional extrapolation
                        cTry = TryPoint(2.0);
                }
                else if(cTry >= m_pdSimplexCosts[m_iNextHi])
                {
                        // point is worse than the second highest, so try to find an intermediate
                        double cSave = m_pdSimplexCosts[m_iHi];
                        cout << "\tTrying intermediate extrapolation . . ." << endl;
                        cTry = TryPoint(0.5);
                        if(cTry >= cSave)
                        {
                                cout << "\tContracting around lowest point . . ." << endl;
                                // can't get rid of bad point, contract around the lowest point
                                for(i = 0; i < m_vpdSimplex.size(); i++)
                                {
                                        if(i != m_iLo)
                                        {
                                                for(j = 0; j < nParameters; j++)
                                                {
                                                        double Q = 0.5*(m_vpdSimplex[i][j] + m_vpdSimplex[m_iLo][j]);
                                                        m_vpdSimplex[i][j] = Q;
                                                        m_pdPointSum[j] = Q;
                                                }
                                                m_pFilter->BackwardTransformation(m_pdPointSum,nParameters);
                                                cTry = m_pMinimizable->ObjectiveFunction(m_pdPointSum);
                                                m_pFilter->ForwardTransformation(m_pdPointSum,nParameters);
                                        }
                                }

                                // recompute the PointSum
                                this->ComputePointSum();
                        }
                }
                
                iteration++;
                // check to see if we should write the best parameters we've found
                if((iteration % 20) == 0)
                {
                        // XXX debug
                        cout << "Writing parameters to file . . ." << endl;
                        // write best parameters
                        fstream *pftemp = new fstream("nmtemp.par",ios::out);
                        *pftemp << "##################################################################" << endl;
                        *pftemp << "# BEST PARAMETERS" << endl;
                        *pftemp << "# COST : " << m_pdSimplexCosts[m_iLo] << endl;
                        *pftemp << "##################################################################" << endl;
                        int pcount;
                        for(pcount = 0; pcount < nParameters; pcount++)
                        {
                                *pftemp << setprecision(10) << m_pFilter->OperatorInverse(m_vpdSimplex[m_iLo][pcount]) << endl;
                        }
                        delete pftemp;
                }
        }       // end while loop
        
        // put best parameters, unfiltered, into parameters vector
        for(j = 0; j < nParameters; j++)
        {
                parameters[j] = m_pFilter->OperatorInverse(m_vpdSimplex[m_iLo][j]);
        }
        cBest = m_pdSimplexCosts[m_iLo];
        // now clear memory
        delete [] m_pdPointSum;
        delete [] m_pdSimplexCosts;
        for(i = 0; i < m_vpdSimplex.size(); i++)
        {
                delete [] m_vpdSimplex[i];
        }

        return cBest;
}

bool CNelderMeadSimplexMinimizer::CheckFuncTol()
{
        double num = 2.0*(fabs(m_pdSimplexCosts[m_iHi] - m_pdSimplexCosts[m_iLo]));
        double den = fabs(m_pdSimplexCosts[m_iHi]) + fabs(m_pdSimplexCosts[m_iLo]) + 1.0e-10;
        double fracTol = num/den;
        
        if(fracTol < m_dTol)
        {
                return true;
        }
        return false;
}

---