#ifndef MINPATTERSEARCHORDER2_H
#define MINPATTERSEARCHORDER2_H

#include <minexpdim.h>


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

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

  /** Allocate the fn with new. */
  minpatternsearchorder2
  (
    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 >  
minpatternsearchorder2<X>::minpatternsearchorder2
(
  funcstate<X>& fn_, 
  X const h0step,
  uintc indexmax_,
  boolc mem_
)
  : minexpdim<X>(fn_,h0step,indexmax_,mem_),
  fh(fn_)
{
}

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

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

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

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

  // Algorithm on vectors
  // X* = 5/2*X[k] -2*X[k-1] + 1/2*X[k-2]
  X* XK0 = fh[0];
  X* XK1 = fh[1];
  X* XK2 = fh[2];

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

//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.5 + XK1[i]*(-2.0) + XK2[i]*0.5; 
  }
  // 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 minpatternsearchorder2<X>::reset()
{
  minexpdim<X>::reset();

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


#endif