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 )