nmsearch.cc

在linux下面实现的单纯性算法的源代码

/*NMSearch.cc
 *class implementation of Nelder Mead Simplex Search
 *Adam Gurson College of William & Mary 1999
 *
 * modified slightly by Anne Shepherd (pls), 8/00
 */

#include "NMSearch.h"
#include "iostream.h"

// constructors & destructors

NMSearch::NMSearch(int dim)
{
   dimensions = dim;
   functionCalls = 0;
   simplex = NULL;
   simplexValues = NULL;
   centroid = new Vector<double>(dimensions,0.0);
   reflectionPt = new Vector<double>(dimensions,0.0);
   expansionPt = new Vector<double>(dimensions,0.0);
   contractionPt = new Vector<double>(dimensions,0.0);
   alpha = 1.0;
   beta = 0.5;
   gamma = 2.0;
   sigma = 0.5;
   scratch = new Vector<double>(dimensions,0.0);
   scratch2 = new Vector<double>(dimensions,0.0);
} // NMSearch() (default)

NMSearch::NMSearch(int dim, double Alpha, double Beta,
                            double Gamma, double Sigma)
{
   dimensions = dim;
   functionCalls = 0;
   simplex = NULL;
   simplexValues = NULL;
   centroid = new Vector<double>(dimensions,0.0);
   reflectionPt = new Vector<double>(dimensions,0.0);
   expansionPt = new Vector<double>(dimensions,0.0);
   contractionPt = new Vector<double>(dimensions,0.0);
   alpha = Alpha;
   beta = Beta;
   gamma = Gamma;
   sigma = Sigma;
   scratch = new Vector<double>(dimensions,0.0);
   scratch2 = new Vector<double>(dimensions,0.0);
} // NMSearch() (special)

NMSearch::NMSearch(const NMSearch& Original)
{
   dimensions = Original.GetVarNo();
   Original.GetCurrentSimplex(simplex);
   Original.GetCurrentSimplexValues(simplexValues);
   alpha = Original.alpha;
   beta = Original.beta;
   gamma = Original.gamma;
   sigma = Original.sigma;
   minIndex = Original.minIndex;
   maxIndex = Original.maxIndex;
   if(centroid != NULL) delete centroid;
   centroid = new Vector<double>(*(Original.centroid));
   if(reflectionPt != NULL) delete reflectionPt;
   reflectionPt = new Vector<double>(*(Original.reflectionPt));
   reflectionPtValue = Original.reflectionPtValue;
   if(expansionPt != NULL) delete expansionPt;
   expansionPt = new Vector<double>(*(Original.expansionPt));
   expansionPtValue = Original.expansionPtValue;
   if(contractionPt != NULL) delete contractionPt;
   contractionPt = new Vector<double>(*(Original.contractionPt));
   contractionPtValue = Original.contractionPtValue;
   functionCalls = Original.functionCalls;
} // NMSearch() (copy constructor)

NMSearch::~NMSearch()
{
   if(simplex != NULL) delete simplex;
   if(simplexValues != NULL) delete [] simplexValues;
   delete centroid;
   delete reflectionPt;
   delete expansionPt;
   delete contractionPt;
   delete scratch;
   delete scratch2;
   //NOTE: Matrix and Vector classes have their own destructors
} // ~NMSearch

// algorithmic routines

void NMSearch::ExploratoryMoves()
{
   double secondHighestPtValue; // used for contraction/reflection decision
   toleranceHit = 0;

   FindMinMaxIndices();
   do {
      if(DEBUG) printSimplex();
      FindCentroid();
      secondHighestPtValue = simplexValues[SecondHighestPtIndex()];
    // reflection step
      FindReflectionPt();

      // stop if at maximum function calls and update the simplex
      /*changed  8/8/00 to fix the problem of maxCalls == -1  --pls
        formerly read if(functionCalls <= maxCalls)
      */
      if (maxCalls > -1
	  && functionCalls >= maxCalls) {
	FindMinMaxIndices();
	ReplaceSimplexPoint(maxIndex, *reflectionPt);
	simplexValues[maxIndex] = reflectionPtValue;
	FindMinMaxIndices(); 
	return;
      } // if using call budget
      
      // possibility 1
      if(simplexValues[minIndex] > reflectionPtValue) {
	FindExpansionPt(); // expansion step
	
	if (reflectionPtValue > expansionPtValue) {
	  ReplaceSimplexPoint(maxIndex, *expansionPt);
	  simplexValues[maxIndex] = expansionPtValue;
	} // inner if
	else {
	  ReplaceSimplexPoint(maxIndex, *reflectionPt);
	  simplexValues[maxIndex] = reflectionPtValue;
	} // else         
      } // if for possibility 1
      
    // possibility 2

      else if( (secondHighestPtValue > reflectionPtValue        ) &&
               (   reflectionPtValue >= simplexValues[minIndex]) ) {
         ReplaceSimplexPoint(maxIndex, *reflectionPt);
         simplexValues[maxIndex] = reflectionPtValue;
      } // else if for possibility 2

    // possibility 3
      else if( reflectionPtValue >= secondHighestPtValue ) {
         FindContractionPt(); // contraction step
         if(maxPrimePtId == 0) {
           if( contractionPtValue > maxPrimePtValue ) {
             ShrinkSimplex();
           } // inner if
           else {
             ReplaceSimplexPoint(maxIndex, *contractionPt);
             simplexValues[maxIndex] = contractionPtValue;
           } // inner else
         } // maxPrimePtId == 0
         else if(maxPrimePtId == 1) {
           if( contractionPtValue >= maxPrimePtValue ) {
             ShrinkSimplex();
           } // inner if
           else {
             ReplaceSimplexPoint(maxIndex, *contractionPt);
             simplexValues[maxIndex] = contractionPtValue;
           } // inner else
         } // maxPrimePtId == 1
      } // else if for possibility 3

    // if we haven't taken care of the current simplex, something's wrong
      else {
         cerr << "Error in ExploratoryMoves() - "
              << "Unaccounted for case.\nTerminating.\n";
         return;
      }
   FindMinMaxIndices();
   } while (!Stop());   // while stopping criteria is not satisfied
} // ExploratoryMoves()

void NMSearch::ReplaceSimplexPoint(int index, const Vector<double>& newPoint)
{
   for( int i = 0; i < dimensions; i++ ) {
      (*simplex)[index][i] = newPoint[i];
   } // for
} // ReplaceSimplexPoint()

void NMSearch::CalculateFunctionValue(int index)
{
   *scratch = (*simplex).row(index);
   int success;
   fcnCall(dimensions, (*scratch).begin(), simplexValues[index], success);
   if(!success) cerr<<"Error calculating point in CalculateFunctionValue().\n";
} // CalculateFunctionValue()

void NMSearch::SetAlpha(double newAlpha)
{
   alpha = newAlpha;
} // SetAlpha()

void NMSearch::SetBeta(double newBeta)
{
   beta = newBeta;
} // SetBeta()

void NMSearch::SetGamma(double newGamma)
{
   gamma = newGamma;
} // SetGamma()

void NMSearch::SetSigma(double newSigma)
{
   sigma = newSigma;
} // SetGamma()

bool NMSearch::Stop()
{
   if(maxCalls > -1) {
      if(functionCalls >= maxCalls)
         return true;
   }

   double mean = 0.0;

   for( int i = 0; i <= dimensions; i++) {
      if( i != minIndex ) {
         mean += simplexValues[i];
      } // if
   } //for 

   mean /= (double)dimensions;

   // Test for the suggested Nelder-Mead stopping criteria
   double total = 0.0;
   for( int i = 0; i <= dimensions; i++ ) {
      total += pow( simplexValues[i] - mean ,2);
   } //for
   total /= ((double)dimensions + 1.0);
   total = sqrt(total);
   
   
   // printSimplex();
   if(total < stoppingStepLength) {
      toleranceHit = 1;
      return true;
   }
   else
      return false;
} // Stop()

void NMSearch::fcnCall(int n, double *x, double& f, int& flag)
{
   fcn(n, x, f, flag);
   functionCalls++;
} // fcnCall()

// Simplex-altering functions

void NMSearch::InitRegularTriangularSimplex(const Vector<double> *basePoint,
                                            const double edgeLength)
{
  //  This routine constructs a regular simplex (i.e., one in which all of 
  //  the edges are of equal length) following an algorithm given by Jacoby,
  //  Kowalik, and Pizzo in "Iterative Methods for Nonlinear Optimization 
  //  Problems," Prentice-Hall (1972).  This algorithm also appears in 
  //  Spendley, Hext, and Himsworth, "Sequential Application of Simplex 
  //  Designs in Optimisation and Evolutionary Operation," Technometrics, 
  //  Vol. 4, No. 4, November 1962, pages 441--461.

   int i,j;
   double p, q, temp;
   Matrix<double> *plex = new Matrix<double>(dimensions+1,dimensions,0.0);
   for( int col = 0; col < dimensions; col++ ) {
      (*plex)[0][col] = (*basePoint)[col];
   }

   temp = dimensions + 1.0;
   q = ((sqrt(temp) - 1.0) / (dimensions * sqrt(2.0))) * edgeLength;
   p = q + ((1.0 / sqrt(2.0)) * edgeLength);

   for(i = 1; i <= dimensions; i++) { 
      for(j = 0; j <= i-2; j++) {