smdsearch.cc

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

/*SMDSearch.cc
 *Declarations of Sequential version of Torczon's Multi-Directional Search
 *Adam Gurson College of William & Mary 2000
 *
 * slightly modified by Anne Shepherd (pls), 8/00
 */

#include "SMDSearch.h"
#include <iostream.h>
#include <iomanip.h>

// constructors & destructors

SMDSearch::SMDSearch(int dim)
{
   dimensions = dim;
   simplex = new Matrix<double>(dimensions+1,dimensions,0.0);
   simplexValues = new double[dimensions+1];
   simplexVBits = new int[dimensions+1];
   refSimplex = new Matrix<double>(dimensions+1,dimensions,0.0);
   refSimplexValues = new double[dimensions+1];
   refSimplexVBits = new int[dimensions+1];
   minPoint = new Vector<double>(dimensions,0.0);
   delta = -1.0;
   sigma = 0.5;
   functionCalls = 0;
   scratch = new Vector<double>(dimensions,0.0);
   scratch2 = new Vector<double>(dimensions,0.0);
} // SMDSearch() (default)

SMDSearch::SMDSearch(int dim, double Sigma)
{
   dimensions = dim;
   simplex = new Matrix<double>(dimensions+1,dimensions,0.0);
   simplexValues = new double[dimensions+1];
   simplexVBits = new int[dimensions+1];
   refSimplex = new Matrix<double>(dimensions+1,dimensions,0.0);
   refSimplexValues = new double[dimensions+1];
   refSimplexVBits = new int[dimensions+1];
   minPoint = new Vector<double>(dimensions,0.0);
   delta = -1.0;
   sigma = Sigma;
   functionCalls = 0;
   scratch = new Vector<double>(dimensions,0.0);
   scratch2 = new Vector<double>(dimensions,0.0);
} // SMDSearch() (special)

/*Slightly modified by Anne Shepherd, 8/00, to initialize some
  fields that were left out in the original
*/
SMDSearch::SMDSearch(const SMDSearch& Original)
{
   dimensions = Original.GetVarNo();
   Original.GetCurrentSimplex(simplex);
   Original.GetCurrentSimplexValues(simplexValues);
   Original.GetCurrentSimplexVBits(simplexVBits);
   refSimplex = new Matrix<double>(*(Original.refSimplex)); //added --pls
   refSimplexValues = new double[dimensions+1];             //added --pls
   refSimplexVBits = new int[dimensions+1];                 //added --pls
   for(int i = 0; i <= dimensions; i++) {                    //added --pls
     refSimplexValues[i] = Original.refSimplexValues[i];    
   } //for
   for(int j = 0; j <= dimensions; j++) {                    //added --pls
     refSimplexVBits[j] = Original.refSimplexVBits[j];
   } //for
   minPoint = new Vector<double>(*(Original.minPoint));
   minValue = Original.minValue;
   delta = Original.delta;
   sigma = Original.sigma;
   functionCalls = Original.functionCalls;
} // SMDSearch() (copy constructor)

SMDSearch::~SMDSearch()
{
   delete simplex;
   delete [] simplexValues;
   delete [] simplexVBits;
   delete refSimplex;
   delete [] refSimplexValues;
   delete [] refSimplexVBits;
   delete minPoint;
   delete scratch;
   delete scratch2;
   //NOTE: Matrix and Vector classes have their own destructors
} // ~SMDSearch

// algorithmic routines

void SMDSearch::ExploratoryMoves()
{
   int done;
   int lastMinIndex = minIndex;
   toleranceHit = 0;

   do {
     done = 0;
     CreateRefSimplex();

     if(DEBUG) {
       printSimplex();
       printRefSimplex();
     }

     // Go through the Reflection Simplex First
     refCurrentIndex = lastMinIndex;
     while( !done && GetAnotherIndex(refCurrentIndex, refSimplexVBits) ) {
       CalculateRefFunctionValue(refCurrentIndex);
       refSimplexVBits[refCurrentIndex] = 1;

       if(DEBUG) printRefSimplex();

       if( refSimplexValues[refCurrentIndex] < minValue ) {
         (*minPoint) = (*refSimplex).row(refCurrentIndex);
         minValue = refSimplexValues[refCurrentIndex];
         lastMinIndex = minIndex;
         minIndex = refCurrentIndex;
         SwitchSimplices();
         done = 1;
       } // if

       /* Changed to fix the problem of maxCalls == -1.
          --pls, 8/8/00
       */
       //if( functionCalls >= maxCalls ) return; 
       if( maxCalls > -1                            // pls
           && functionCalls >= maxCalls ) return;
     } // while (reflection search)

     // Go through the Primary Simplex Next
     while( !done && GetAnotherIndex(currentIndex, simplexVBits) ) {
       CalculateFunctionValue(currentIndex); 
       simplexVBits[currentIndex] = 1;
       // NOTE: currentIndex initialized in InitGeneralSimplex()

       if(DEBUG) printSimplex();

       if( simplexValues[currentIndex] < minValue ) {
         (*minPoint) = (*simplex).row(currentIndex);
         minValue = simplexValues[currentIndex];
         lastMinIndex = minIndex;
         minIndex = currentIndex;
         done = 1;
       } // if

       /* Changed to fix the problem of maxCalls == -1.
          --pls, 8/8/00
       */
       //if( functionCalls >= maxCalls ) return; 
       if( maxCalls > -1                            // pls
           && functionCalls >= maxCalls ) return;
     } // while (primary search)
     
     // Still there's no new min now, shrink
     if( !done )
       ShrinkSimplex();

   } while (!Stop());   // while stopping criteria is not satisfied
} // ExploratoryMoves()

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

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

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

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

   if(delta < stoppingStepLength) {
      toleranceHit = 1;
      return true;
   }
   else
      return false;
} // Stop()

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

// Simplex-altering functions

void SMDSearch::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++) {
         (*plex)[i][j] = (*plex)[0][j] + q;
      } // inner for 1
      j = i - 1; // is this line necessary (redundant)
      (*plex)[i][j] = (*plex)[0][j] + p;
      for(j = i; j < dimensions; j++) {
            (*plex)[i][j] = (*plex)[0][j] + q;
      } // inner for 2
   } // outer for

   delta = edgeLength;
   InitGeneralSimplex(plex);
   delete plex;
} // InitRegularTriangularSimplex()

void SMDSearch::InitFixedLengthRightSimplex(const Vector<double> *basePoint,
                                            const double edgeLength)
{
  // to take advantage of code reuse, this function simply turns
  // edgeLength into a vector of dimensions length, and then
  // calls InitVariableLengthRightSimplex()

  //  cerr << *basePoint;
   double* edgeLengths = new double[dimensions];
   for( int i = 0; i < dimensions; i++ ) {
      edgeLengths[i] = edgeLength;
   }
   InitVariableLengthRightSimplex(basePoint,edgeLengths);
   delete [] edgeLengths;
} // InitFixedLengthRightSimplex()

void SMDSearch::InitVariableLengthRightSimplex(const Vector<double> *basePoint,
                                               const double* edgeLengths)
{
   Matrix<double> *plex = new Matrix<double>(dimensions+1,dimensions,0.0);
   for( int i = 0; i < dimensions; i++ ) {
      // we're building the basePoint component-by-component into
      // the (n+1)st row
      (*plex)[dimensions][i] = (*basePoint)[i];

      // now fill in the ith row with the proper point
      for( int j = 0; j < dimensions; j++ ) {
         (*plex)[i][j] = (*basePoint)[j];