/*  NLessSearch.cc 
 *  implementations of the derived class NLessSearch functions.
 *  The NLessSearch searches about a regular simplex (i.e. 
 *  minimal positive basis) until finding improvement in the
 *  objective function value.  Then the search relocates to the
 *  improving point and begins again.
 *  Liz Dolan, The College of William & Mary, 1999
 *
 *                                                                 
 *  The author of this software is Elizabeth  D. Dolan.                     
 *  Permission to use, copy, modify, and distribute this software  
 *  for any purpose without fee is hereby granted, provided that   
 *  this entire notice is included in all copies of any software   
 *  which is or includes a copy or modification of this software   
 *  and in all copies of the supporting documentation for such     
 *  software.  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT    
 *  ANY EXPRESS OR IMPLIED WARRANTY.  IN PARTICULAR, THE AUTHOR    
 *  OFFERS NO REPRESENTATION OR WARRANTY OF ANY KIND                   
 *  CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS    
 *  FITNESS FOR ANY PARTICULAR PURPOSE.                            
 */                                                                 

#include "NLessSearch.h" 

NLessSearch::NLessSearch(int numberOfVariables): PatternSearch(numberOfVariables)
{
}

NLessSearch::~NLessSearch()
{
}

void NLessSearch::ExploratoryMoves()
{
  double pace = 0.0;
  pace = GetStepLength();  //as given by the user or default
  CreatePattern(pace);     //creates a regular simplex
  //note that the actual use of pace has been deprecated.
  Vector<double> currentPoint(GetVarNo());
  Vector<double> nextPoint(GetVarNo());
  double value;
  double nextValue;
  int length;
  int success = 0;
  GetMinPoint(currentPoint);  //initialize min point
  GetMinVal(value);           //and obj. func. value
  GetPatternLength(length);
  //search the pattern in each direction until improvement is found,
  //then stop and begin a new iteration at the better point

  do
    {
      for( int i = 0; i < length; i++)
	{
	  NextPoint(i, currentPoint, nextPoint);
	  fcnCall(GetVarNo(), nextPoint.begin(), nextValue, success);
	  if(success==1)
	    {
	      if(nextValue < value)
		{
		  ReplaceMinimum(nextPoint, nextValue);
		  value = nextValue;
		  currentPoint = (nextPoint);
		  i = -1; //start the search over at the new point
		}//if there is improvement
	    }//if able to get a function value
	}//for now we know that there aren't better points around this one
	UpdatePattern();
      
    }while(!Stop());//while we haven't stopped()
}//ExploratoryMoves

void NLessSearch::CreatePattern(double pace)
/*
  For information on how to build a regular
  simplex, see pages 79-81 S.L.S. Jacoby, J.S. Kowalik 
  and J.T. Pizzo.
  Iterative Methods for Nonlinear Optimization Problems. 
  Prentice Hall, Inc., Englewood Cliffs, NJ. 1972.
*/
{
  int vars = GetVarNo();
  //refer to book for p and q
  double p = (sqrt(vars+1) -1 + vars)/(vars*sqrt(2));
  double q = (sqrt(vars+1) -1)/(vars*sqrt(2));
  Matrix<double> nlessPattern(vars, vars+1);  //vars + 1 = # of vectors in pattern
  InitializePattern(vars + 1, &nlessPattern); //initialize pattern to default
  Vector<double> basis(vars);  
  //basis is the first point used to create a simplex according to the algorithm
  if(vars > 0)
    {
      for(int j = 0; j < vars; j++) 
	{
	  basis[j] = (- (p+((vars-1)*q))/(vars+1));
	  nlessPattern[j][0] = basis[j];
	}

      for(int i = 1; i < vars + 1; i++)
	{
	  for(int k = 0; k < vars; k++)
	    {
	      nlessPattern[k][i] = basis[k] + q;
	    }
	  nlessPattern[i-1][i] = basis[i-1] + p;
	}//outer for
      //make sure that the vectors of the pattern not only point
      //in the right direction, but are also of desired length
      SizePattern(vars, nlessPattern, 1.0);
      InitializePattern(vars + 1, &nlessPattern);
    }//if there's anything to allocate
}//CreatePattern
  
void NLessSearch::UpdatePattern()
{
  ScalePattern(0.5);
}//UpdatePattern

void NLessSearch::SizePattern(int dimens, Matrix<double> &pat, double size)
{
  int length;
  GetPatternLength(length);  //length of the pattern
  Vector<double> compare(dimens);
  Vector<double> compareTo(dimens);
  double compDist; //distance from the centroid
  compareTo = pat.col(0); //initialize for testing
  for(int i = 0;i < length;i++)
    {
      compare = pat.col(i);
      compDist = 0;

      //find the distance between the vectors being compared
      compDist = compare.l2norm();

      /*
	calculate the angles and assure that they are equal
	between compareTo and compare for each compare n.e.
	compareTo
	see page 106 of Gilbert Strang, Linear Algebra and its
	Applications
      
	if(i!=0)
	{
          double angle = 0;
          for(int calc = 0; calc < dimens; calc++)
	    {
	      angle += (compare[calc] * compareTo[calc]);
	    }
          angle = angle/(compDist*compDist);
          angle = acos(angle);
        }
      */

      //resize to desired length, based on knowledge of
      //relationships in the ratios of similar triangles
      for(int j = 0; j < dimens; j++)
	{
	  pat[j][i] = size*compare[j]/compDist;
	}
      compDist = 0;
      
    }//for

}