patternsearchD2.h

Go to the documentation of this file.

#ifndef PATTERNSEARCHD2_H
#define PATTERNSEARCHD2_H

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

template< typename EXP >
class patternsearchD2
{
public:

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

  EXP exp;
  
  T hstep;

  uintc dim;

  T* coef;

  patternsearchD2( EXP exp_, uintc dim_, T hstep_=1.0);

  ~patternsearchD2();

  cirbuffarr< T > xi0;

  vec< XI, T > xivec;

  void reset();

  void reset( T const * x0 )
    { exp.reset(x0); reset(); }

  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