#ifndef PATTERNSEARCHD2_H
#define PATTERNSEARCHD2_H

#include <vec.h>
#include <cirbuffarr.h>
#include <typeop.h>
#include <print.h>

/*!
  \brief Minor adaption of the Pattern search algorithm by an K-dimensional taylor series pattern move.

  
  The traditional pattern move has dim=1 and h=1.

  I modified the algorithm for a variable step length.
  The pattern move was changed from a first order approximation to a 
  dim order approximation.

  Pattern move : 
  X[n+1] = X[n] + h*D(X[n-1]) + h^2/2*D^2(X[n-1]) + ...
              + h^dim/(dim!)D^dim(X[n-dim])
*/
template< typename EXP >
class patternsearchD2
{
public:

  typedef typename typeop<EXP>::Tbare::Ttype T;
  typedef typename typeop<EXP>::Tbare::XItype XI;

  /** The exploratory mover. */
  EXP exp;
  
  /** Pattern move step length. */
  T hstep;

  /** Dimension of the pattern move approximator. */
  uintc dim;

  /** Pattern move is X[n+1] = X[n]*coef[0] + X[n-1]*coef[1] +... */
  T* coef;

  /** The pattern move is configured with an exploratory mover and variable step length. */
  patternsearchD2( EXP exp_, uintc dim_, T hstep_=1.0);

  /** Memory cleanup. */
  ~patternsearchD2();

  /** The previous state. */
  cirbuffarr< T > xi0;

  /** The current state as a vec for arithmetic operations. */
  vec< XI, T > xivec;

  /** Reset from the current position and keep the state. */
  void reset();

  /** New position and new state. ie a full reset. */
  void reset( T const * x0 )
    { exp.reset(x0); reset(); }

  /** One exploratory move then pattern moves. */
  void operator ++ ();

  bool const operator ! () const
    { return exp.valid; }
  
};


/*
 *---------------------------------------------------------------
 * Implementation
 */

template< typename EXP >
patternsearchD2<EXP>::~patternsearchD2()
{
  delete[] coef;
};

template< typename EXP >
void patternsearchD2<EXP>::reset()
{ 
  exp.reset(); 
  xi0.push(exp.xi); 
  for (uint i=1; i<dim; ++i)
    exp.move(); 
    xi0.push(exp.xi); 
}


template< typename EXP >
void patternsearchD2<EXP>::operator ++ ()
{
  uintc N(exp.N);
  uint k;
  uint j;
  for (bool loop=true; loop;)
  {
    xi0.push(exp.xi);

    for (k=0; k<N; ++k)
    {
      exp.xi[k] = coef[0]*xi0[0][k];
      for (j=1; j<=dim; ++j)
        exp.xi[k] += coef[j]*xi0[j][k];
    }

//      xivec[k] = coef[0]*xi0[0][k] + coef[1]*xi0[1][k] 
//                 + coef[2]*xi0[2][k];

    exp.movelocal();

    loop = exp.evaluate();
  }
  xi0.copyto(exp.xi);
}


template< typename EXP >
patternsearchD2<EXP>::patternsearchD2( EXP exp_, uintc dim_, T hstep_ )
  : exp(exp_), hstep(hstep_), dim(dim_), coef(new T[dim_+1]),
    xi0(dim+3,exp.N), xivec(exp.xi,exp.N)
{
  assert(dim>0);

  uintc w(dim+1);

  T t0[w];
  vec< T*, T > t(t0,w);
  
  T xop[w];
  vec< T*, T > x(xop,w);
  x.identity(0);

  vec< T*, T > coefv(coef,w);

  coefv = x;
  T hfactor(1);
  for (uint k=1; k<w; ++k)
  {
    hfactor *= (hstep/(T)k);
    x.difference();
    t = x;
    t *= hfactor;
    coefv += t;
  }

  //cout << printvecfunc(coef,dim+1) << endl;

/*

  coef[0] = hstep*hstep*0.5 + hstep + 1.0;
  coef[1] = -hstep-hstep*hstep;
  coef[2] = hstep*hstep*0.5;
*/


}


#endif