/*  HJSearch.cc
 *  implementation of the HJSearch class to find a minimal objective function
 *  solution.
 *  For a good description of the Hooke and Jeeves search algorithm
 *  I recommend Non-Linear Optimization Techniques by Box, Davies,
 *  and Swann, 1969. 
 *  Liz Dolan, The College of William and 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 "HJSearch.h"


HJSearch::HJSearch(int numberOfVariables): PatternSearch(numberOfVariables)
{
  step = initialStepLength;  //as provided by the user or default
  factor = 0.5;              //factor by which the step length is reduced
}

HJSearch::HJSearch(const HJSearch & Original): PatternSearch(Original)
{
  step = Original.step;
  factor = Original.factor;    //factor by which the step length is reduced
}

HJSearch::~HJSearch()
{
}

void HJSearch::CopySearch(const HJSearch & Original)
{
  PatternSearch::CopySearch(Original);
  step = Original.step;
  factor = Original.factor;
}

void HJSearch::CleanSlate(int dimensions)
//reinitialize all values
{
  PatternSearch::CleanSlate(dimensions);
  step = initialStepLength;
}

void HJSearch::ExploratoryMoves()
{
  int dimens = GetVarNo();
  Vector<double> currentPoint(dimens);
  Vector<double> lastImprovingPoint(dimens); //last base point
  Vector<double> storage(dimens); //for intermediate storage to reduce rounding error
  Vector<double> direction(dimens, 0.0); //direction of pattern extended step
  double value;                   //objective function value
  double positiveValue;           //obj.fun. value in the positive step
  double negativeValue;           // " for the negative step
  double lastImprovingValue;      //obj.fun.value of last base point
  int success = 0;
  GetMinPoint(currentPoint);      //initialize to the user initial point
  GetMinPoint(lastImprovingPoint);
  GetMinPoint(storage);
  GetMinVal(value);
  GetMinVal(lastImprovingValue);
  bool foundImprove = false;
  do
    {
      for(int iteration=0;iteration < dimens;iteration++)
	{
	  currentPoint[iteration] += step; 
	  fcnCall(dimens, currentPoint.begin(), positiveValue, success);
	  if(success!=1)
	    {
	      positiveValue = value + 1.0;
	      //if the call returned unsuccessfully, set positiveValue
              //to a value that will not be improving
	    }
	  if(positiveValue < value)
	    {
	      value = positiveValue; 
	      foundImprove = true;
	      //continue search in other dimensions from here
	    }//if positive is better
	    
	  if(!foundImprove)
	    {
	      currentPoint = storage;
	      currentPoint[iteration] -= (step);
	      fcnCall(dimens, currentPoint.begin(), negativeValue, success);
	      if(success!=1)
		{
		  negativeValue = value + 1.0;
		  //same kludge as in positive case
		}
	      if(negativeValue < value)
		{
		  value = negativeValue;
		  foundImprove = true;
		  //continue search in other dimensions from here
		}//if negative direction is better
	    }//if we need to check the negative
	  if(!foundImprove)
	    {//reset to original position
	      currentPoint = storage;
	    }//if neither direction gave improvement
	  else
	    {
	      storage = currentPoint;
	    }
	  foundImprove = false;      //reset for next iteration
	}//for
      direction =  currentPoint - lastImprovingPoint;
      //direction now holds the extended pattern step vector
      if(value < lastImprovingValue)
	{
	  //check whether the "new" point is within factor*step of the old
	  if(isnear(lastImprovingPoint, currentPoint, factor*step))
	    {
	      currentPoint = lastImprovingPoint;
	      value = lastImprovingValue;
	      storage = currentPoint;	   
	    }
	  else   //some step yielded improvement
	    {
	      lastImprovingValue = value;
	      ReplaceMinimum(currentPoint, value);
	      lastImprovingPoint = currentPoint;

	      //take the pattern extending step and find its value
	      currentPoint = direction + currentPoint;
	      storage = currentPoint;
	      fcnCall(dimens, currentPoint.begin(), value, success);
	    }	  
	}
      else
	{
	  if(isnear(currentPoint, lastImprovingPoint, factor*step))
	    {
	      step = step * factor;
	    }
	  else
	    {
	      //this case can only occur after an unsuccessful
	      //search about a pattern-step-located point
	      //move back to the point that was improving from the
	      //search about the last base point
	      GetMinPoint(currentPoint);
	      GetMinVal(value);
	      storage = currentPoint;
	    }
	}
    } while(!Stop(step));//while we haven't stopped()
}//ExploratoryMoves


bool HJSearch::Stop(double pace)
//Makes certain search has not exceeded maxCalls or 
//is stepping at less than stoppingStepLength
//If stop is true, resets step to last used length
{
  if(maxCalls > -1)
    if(GetFunctionCalls() > maxCalls)
      	return true;
  if(pace < stoppingStepLength)
      return true;
  return false;
}//Stop

double HJSearch::GetStepLength()
{
  return step;
}