Simplex Optimization

Listing 1: The basic simplex template.

#if !defined ( SIMPLEX_H_INCLUDED )
#define SIMPLEX_H_INCLUDED
#include <cmath>
#include <cfloat>
#include <valarray>
#include <utility>
using namespace std;
typedef valarray<double> vd;
template < typename T > class simplex
{
public:
//---------------------------------------------------------------
   simplex ( const T t ) : m_nVar( t( ).size( ) ) // constructor
   {
      m_FunctionObject       = t;  // copy the function object in
                                   // for partial safety
      m_vInitialVector.resize( m_nVar );
      m_vInitialVector       = t( );
      m_vDeltaVector.resize  ( m_nVar );
      m_vFinalVector.resize  ( m_nVar );
      m_bDeltaVectorIsValid  = false;
      m_dFinalFOM            = DBL_MAX; // make huge starting FOM
      m_lMaximumIterations   = 5000;
      m_dSmallestAllowedShift = 1.0E-7;
       //if the change falls below this value, minimization stops
   }; 
//---------------------------------------------------------------
   std::pair<double, vd> fnMinimize( void )  
   {  //  this does all the work for simplex
      if ( !m_bDeltaVectorIsValid ) m_fnComputeDeltaVector( );
      m_fnSetupSimplex( );
      valarray<double> vShift(m_nVar);
      valarray<double> vPreviousResponseVector( m_nVar+1 );
      vPreviousResponseVector = DBL_MAX;
      m_lIterationsPerformed = 0;
      for ( size_t iterations=0; 
                 iterations<m_lMaximumIterations; iterations++ )
      {  double dResponse1, dResponse2;
         valarray<double> vCentroid( m_nVar );
         valarray<double> vTest1( m_nVar );
         vCentroid = 0.0;
         vTest1 = 0;
         fnFindMinMaxResponse( );
         vCentroid = fnCalculateSimplexCentroid( );
         vShift    = vCentroid - m_vSimplexVertex[ m_lWorstFitVertex ];
         vTest1    = vCentroid + vShift;
         dResponse1 = m_FunctionObject( vTest1 );
         //  if better than the best then try a bigger step
         if (dResponse1 < m_vVectorOfFOM[ m_lBestFitVertex ] )
         {
            valarray<double> vTest2( m_nVar );
            vTest2 = vTest1 + vShift;
            dResponse2 = m_FunctionObject ( vTest2 );
            //  if better than the best, use this one
            if ( dResponse2 < m_vVectorOfFOM[m_lBestFitVertex] )
            {  // doubled reflection shift
               m_vSimplexVertex[ m_lWorstFitVertex ] = vTest2;
               m_vVectorOfFOM[ m_lWorstFitVertex ] = dResponse2;
            }
            else
            {  // reflected the worst vertex
               m_vSimplexVertex[ m_lWorstFitVertex ] = vTest1;
               m_vVectorOfFOM[ m_lWorstFitVertex ] = dResponse1;
            }
         }
         // if better than the worst, then just replace the worst
         else if ( dResponse1 < m_vVectorOfFOM[ m_lWorstFitVertex ] )
         {  // simple replacement
            m_vSimplexVertex[ m_lWorstFitVertex ] = vTest1;
            m_vVectorOfFOM[ m_lWorstFitVertex ]     = dResponse1;
         }
         //  otherwise, it was worse than the worst we know
         else
         {
            //try the point midway between worst and the centroid
            vTest1 = vCentroid - 0.5 * vShift;
            dResponse1 = m_FunctionObject ( vTest1 );
            // if better than the worst we know, then just use it
            if (dResponse1 < m_vVectorOfFOM[m_lWorstFitVertex] )
            {  // replace with halfway between worst and centroid
               m_vSimplexVertex[ m_lWorstFitVertex ] = vTest1;
               m_vVectorOfFOM[ m_lWorstFitVertex ]  = dResponse1;
            }
            else
            //  nothing worked; contract toward the best point 
            {  //  just replace each point (except best) with the
              //  point halfway between best and itself
              size_t i;
              for ( i=0; i<m_nVar; i++ )
              {
                 if ( i != m_lBestFitVertex )
                 {
                    m_vSimplexVertex[ i ] =
                    ( m_vSimplexVertex[ i ] + 
                       m_vSimplexVertex[m_lBestFitVertex] ) /2.0;
                    m_vVectorOfFOM[ i ] = 
                       m_FunctionObject( m_vSimplexVertex[ i ] );
                 }
              }
            }
         }
         if ( iterations > 1 )
         {
            valarray<double> vTest = vPreviousResponseVector - 
              m_vVectorOfFOM;
            if ( sqrt( vTest.sum( ) ) <= m_dSmallestAllowedShift ) break;
         }
         vPreviousResponseVector = m_vVectorOfFOM;
         fnFindMinMaxResponse(  );
      } // iterations
      m_lIterationsPerformed = iterations;
      fnFindMinMaxResponse(  );
      m_vFinalVector = m_vSimplexVertex[ m_lBestFitVertex ];
      m_dFinalFOM    = m_vVectorOfFOM[ m_lBestFitVertex ];
      return ( make_pair( // return best FOM and best vertex
           m_dFinalFOM, m_vSimplexVertex[ m_lBestFitVertex ] ) );
   };
//---------------------------------------------------------------
private:
    const size_t     m_nVar;
    size_t           m_lMaximumIterations;
    double           m_dSmallestAllowedShift;
    valarray<double> m_vInitialVector;
    valarray<double> m_vFinalVector;
    valarray<double> m_vDeltaVector;
    valarray<double> m_vVectorOfFOM;
    valarray<vd>     m_vSimplexVertex; // vd is just a typedef
                                       // for the Microsoft compiler 
    bool             m_bDeltaVectorIsValid;
    double           m_dFinalFOM;
    size_t           m_lBestFitVertex;
    size_t           m_lWorstFitVertex;
    double           m_dBestFitValue;
    double           m_dWorstFitValue;
    size_t           m_lIterationsPerformed;
    // copy the function object in for safety    
    T                m_FunctionObject; 
//---------------------------------------------------------------
   void m_fnComputeDeltaVector ( void )
   {
      size_t i;
      for ( i=0; i<m_vDeltaVector.size( ); ++i )
      {
         if ( fabs( m_vInitialVector[i] ) > 100.0 * DBL_MIN )
         {
            m_vDeltaVector[i] = 0.01 * m_vInitialVector[i]; 
         }
         else
         {
            m_vDeltaVector = 0.1;
         }
      }
      m_bDeltaVectorIsValid = true;
   };
//---------------------------------------------------------------
   void m_fnSetupSimplex ( void )
   {
      size_t i;
      m_vVectorOfFOM.resize( m_nVar + 1 );
      m_vSimplexVertex.resize( m_nVar + 1 );
      valarray<double> tempVertex( m_nVar );
      for ( i=0; i<m_nVar+1; ++i )
      {
         tempVertex = m_vInitialVector;
         if ( i < m_nVar ) tempVertex[i] += m_vDeltaVector[i];
         m_vVectorOfFOM[i] = m_FunctionObject( tempVertex );
         m_vSimplexVertex[ i ] = tempVertex;
      }
      fnFindMinMaxResponse ( );
   };
//---------------------------------------------------------------
   void fnFindMinMaxResponse( void )
   {
      m_dBestFitValue  = m_vVectorOfFOM.min( );
      m_dWorstFitValue = m_vVectorOfFOM.max( );
      size_t i;
      m_lBestFitVertex = -1;
      for ( i=0; i<m_nVar+1; ++i )
      {
         if ( m_vVectorOfFOM[i] == m_dBestFitValue )
         {
            m_lBestFitVertex = i;
            break;
         }
      }
      for ( i=0; i<m_nVar+1; ++i )
      {
         if ( m_vVectorOfFOM[i] == m_dWorstFitValue )
         {
            m_lWorstFitVertex = i;
            break;
         }
      }
   };
//---------------------------------------------------------------
   valarray<double> fnCalculateSimplexCentroid( void ) const
   {
      return ( ( m_vSimplexVertex.sum( ) - 
          m_vSimplexVertex[ m_lWorstFitVertex ] ) / 
          double( m_nVar ) );
   }
};  //  class simplex
#endif  // if !defined ( SIMPLEX_H_INCLUDED )