#ifndef MINPATTERSEARCH_H
#define MINPATTERSEARCH_H

#include <minexpdim.h>

/*!
\brief Pattern search with order 1 approximator.
*/
template< typename X >
class minpatternsearch : public minexpdim<X>
{
public:

  /** Consecutive X values are stored. */
  funchistory<X> fh;

  /** Allocate the fn with new. */
  minpatternsearch
  (
    funcstate<X>& fn_, 
    X const h0step=20.0,
    uintc indexmax_=500,
    boolc mem_=true
  );

  /** Set the current position as the new minimum. */
  virtual void reset();
  /** 1st order approximator to minimum. */
  virtual void operator ++ ();


};

//----------------------------------------------------
//  Implementation


template < typename X >  
minpatternsearch<X>::minpatternsearch
(
  funcstate<X>& fn_, 
  X const h0step,
  uintc indexmax_,
  boolc mem_
)
  : minexpdim<X>(fn_,h0step,indexmax_,mem_),
  fh(fn_)
{


}


template < typename X >
void minpatternsearch<X>::operator ++ ()
{
//cout << "minpatternsearch<X>::operator ++" << endl;
  if (!minexpdim<X>::valid)
    return;

  if (minexpdim<X>::index > minexpdim<X>::indexmax)
  {
    minexpdim<X>::valid=false;
    return;
  }

  assert( fh.xik.size() >= 2 );

  uintc dim = minexpdim<X>::fn->dim;

  // Algorithm on vectors
  // X* = X[k] + (X[k] - X[k-1])
  X* XK0 = fh[0];
  X* XK1 = fh[1];

//cout << "before pattern move" << endl;
//cout << "fmin=" << minexpdim<X>::fmin << " xi[dim]=" << minexpdim<X>::xi[dim] << endl;

  // TODO preserve hi, 
  //  have funcstate hi variable,
  //  funchistory have hipush, hipop, hidel_back.
  //  flag for optionally turning on preserving h.

//cout << "fmin=" << minexpdim<X>::fmin << " xi[dim]=" << minexpdim<X>::xi[dim] << endl;
  assert(minexpdim<X>::fmin==minexpdim<X>::xi[dim]);

  for (uint i=0; i<dim; ++i)
  {
    minexpdim<X>::xi[i] = XK0[i]*2.0-XK1[i]; 
  }
  // Was the pattern + exploratory move successful
  if (minexpdim<X>::moveOrder1())
  {
    minexpdim<X>::xi[dim] = minexpdim<X>::fmin;
    assert(minexpdim<X>::xi[dim]<XK0[dim]);
    fh.del_back();
    fh.push();
  }
  else
  {
    fh.restore();
//cout << "fmin=" << minexpdim<X>::fmin << " xi[dim]=" << minexpdim<X>::xi[dim] << endl;
    assert(minexpdim<X>::fmin==minexpdim<X>::xi[dim]);

    bool res=minexpdim<X>::moveOrder1(5); 
    assert(res); 
    if (res)
    {
      minexpdim<X>::xi[dim] = minexpdim<X>::fmin;
      fh.del_back();
      fh.push();
    }
  }
}

template < typename X >
void minpatternsearch<X>::reset()
{
  minexpdim<X>::reset();

  static uint maxloop=20;
  for (uint j=0; j<2; ++j)
  {
    minexpdim<X>::moveOrder1(maxloop);
    assert(minexpdim<X>::hasminimized==true);
    fh.push();
    minexpdim<X>::hasminimized=false;
  }
}


#endif