smdsearch.cc

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

         if( i == j )
            (*plex)[i][j] += edgeLengths[i];
      }
      
      if( edgeLengths[i] > delta ) delta = edgeLengths[i];
   }
   InitGeneralSimplex(plex);
   delete plex;
} // InitVariableLengthRightSimplex()

void SMDSearch::InitGeneralSimplex(const Matrix<double> *plex)
{
   functionCalls = 0;
   (*simplex) = (*plex);

   // zero out the valid bits
   for(int i = 0; i < dimensions; i++)
     simplexVBits[i] = 0;

   // NOTE: basePoint MUST be located in the last row of plex
   Vector<double> basePoint = (*plex).row(dimensions);

   // evaluate f(basePoint) and initialize it as the min
   int success;
   fcnCall(dimensions, (basePoint).begin(), simplexValues[dimensions], success);
   if(!success) cerr<<"Error with basePoint in initial simplex.\n";
   simplexVBits[dimensions] = 1;
   (*minPoint) = (basePoint);
   minValue = simplexValues[dimensions];
   currentIndex = minIndex = dimensions;

   // if we still haven't defined delta, go through the simplex and
   // define delta to be the length of the LONGEST simplex side
   double temp;
   if( delta < 0.0 ) {
     for( int j = 0; j < dimensions; j++ ) {
       for ( int k = j+1; k <= dimensions; k++ ) {
         temp = ( ((*simplex).row(j)) - ((*simplex).row(k)) ).l2norm();
         if( temp > delta ) delta = temp;
       } // inner for
     } // outer for
   } // outer if

   // if delta is still not defined, there is a definite problem
   if( delta < 0.0 )
     cout << "Error in simplex initialization: delta not set.\n";
   // cout << "\nminValue = " << minValue << " and minPoint = " << *minPoint << endl;
   //printSimplex();
} // InitGeneralSimplex()

void SMDSearch::ReadSimplexFile(istream& fp)
{
   if(fp == NULL) {
      cerr<<"No Input Stream in ReadSimplexFile()!\n";
      return; // There's no file handle!!
   }
   
   Vector<double> *basePoint = new Vector<double>(dimensions);
   Matrix<double> *plex = new Matrix<double>(dimensions+1,dimensions);
   for( int i = 0; i <= dimensions; i++ ) {
      for ( int j = 0; j < dimensions; j++ ) {
         fp >> (*plex)[i][j];
      } // inner for
   } // outer for
   (*basePoint) = (*plex).row(dimensions);
   InitGeneralSimplex(plex);
   delete basePoint;
   delete plex;
} // ReadSimplexFile()

// Query functions

int SMDSearch::GetFunctionCalls() const
{
   return functionCalls;
} // GetFunctionCalls()

void SMDSearch::GetMinPoint(Vector<double>* &minimum) const
{
   minimum = new Vector<double>((*minPoint));
} // GetMinPoint()

double SMDSearch::GetMinVal() const
{
   return minValue;
} // GetMinVal()

void SMDSearch::GetCurrentSimplex(Matrix<double>* &plex) const
{
   plex = new Matrix<double>((*simplex));
} // GetCurrentSimplex()

void SMDSearch::GetCurrentSimplexValues(double* &simValues) const
{
   simValues = new double[dimensions+1];
   for( int i = 0; i <= dimensions; i++ ) {
      simValues[i] = simplexValues[i];
   } // for
} // GetCurrentSimplexValues()

void SMDSearch::GetCurrentSimplexVBits(int* &simVBits) const
{
   simVBits = new int[dimensions+1];
   for( int i = 0; i <= dimensions; i++ ) {
      simVBits[i] = simplexVBits[i];
   } // for
} // GetCurrentSimplexValues()

int SMDSearch::GetVarNo() const
{
   return dimensions;
} // GetVarNo()

int SMDSearch::GetTolHit() const
{
   return toleranceHit;
} // GetTolHit()

// private functions

void SMDSearch::CreateRefSimplex()
{ 
  // copy the known flip point over
  for( int i = 0; i < dimensions; i++ )
    (*refSimplex)[currentIndex][i] = (*simplex)[currentIndex][i];
  refSimplexValues[currentIndex] = simplexValues[currentIndex];
  refSimplexVBits[currentIndex] = simplexVBits[currentIndex];
  refCurrentIndex = currentIndex;

  // reflect the remaining points   
  for( int j = 0; j <= dimensions; j++ ) {
    if( j != currentIndex ) {
      refSimplexVBits[j] = 0;
      (*scratch) = ( (*simplex).row(currentIndex) * 2.0 ) - (*simplex).row(j);
      for( int k = 0; k < dimensions; k++ ) 
         (*refSimplex)[j][k] = (*scratch)[k];
    } // if
  } // outer for
} // CreateRefSimplex()

void SMDSearch::SwitchSimplices()
{
  // this allows us to remove the need to delete and
  // reallocate memory by simply swapping pointers
  // and using the same two "simplex memory slots"
  // for the entire search

  Matrix<double> *tmp1 = simplex;
  double         *tmp2 = simplexValues;
  int            *tmp3 = simplexVBits;
  int             tmp4 = currentIndex;

  simplex = refSimplex;
  simplexValues = refSimplexValues;
  simplexVBits = refSimplexVBits;
  currentIndex = refCurrentIndex;

  refSimplex = tmp1;
  refSimplexValues = tmp2;
  refSimplexVBits = tmp3;
  refCurrentIndex = tmp4;
} // SwitchSimplices()

void SMDSearch::ShrinkSimplex()
{
  if(DEBUG) cout << "Shrinking Simplex.\n\n";

   delta *= sigma;
   currentIndex = minIndex;
   Vector<double> *lowestPt = scratch;
   *lowestPt = (*simplex).row(minIndex);
   Vector<double> *tempPt = scratch2;
 
   for( int i = 0; i <= dimensions; i++ ) {
      if( i != minIndex ) {
         *tempPt = (*simplex).row(i);
         (*tempPt) = (*tempPt) + ( sigma * ( (*lowestPt)-(*tempPt) ) );
         for( int j = 0; j < dimensions; j++ ) {
            (*simplex)[i][j] = (*tempPt)[j];
         } // inner for

         simplexVBits[i] = 0;         

      } // if
   } // outer for
} // ShrinkSimplex()

int SMDSearch::GetAnotherIndex(int& index, int*& validBits)
{
  if ( !validBits[index] ) return 1;

  int initialIndex = index;
   
  do {
    index++;
    if( index > dimensions ) index = 0;
  } while ( ( index != initialIndex) && 
            ( validBits[index] ) );
    
  if( index == initialIndex )
    return 0;
  else
    return 1;
} // GetAnotherIndex()

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

void SMDSearch::printSimplex() const
{
  cout << "Primary Simplex:\n";

  for( int i = 0; i <= dimensions; i++ ) {
     cout << "Point: ";
     for ( int j = 0; j < dimensions; j++ ) {
       cout << (*simplex)[i][j] << " ";
     } // inner for
     cout << "   Value: " << simplexValues[i];
     
     if( simplexVBits[i] )
       cout << "   Valid\n";
     else
       cout << "   Invalid\n";
  } // outer for

  cout << "FCalls: " << functionCalls 
       << "   Delta: " << delta << "\n\n";
} // printSimplex()

void SMDSearch::printRefSimplex() const
{
  cout << "Reflection Simplex:\n";

  for( int i = 0; i <= dimensions; i++ ) {
     cout << "Point: ";
     for ( int j = 0; j < dimensions; j++ ) {
       cout << (*refSimplex)[i][j] << " ";
     } // inner for
     cout << "   Value: " << refSimplexValues[i];
     
     if( refSimplexVBits[i] )
       cout << "   Valid\n";
     else
       cout << "   Invalid\n";
  } // outer for

  cout << "FCalls: " << functionCalls 
       << "   Delta: " << delta << "\n\n";
} // printRefSimplex()