#ifndef POWELL03_H
#define POWELL03_H

#include <print.h>

#include <linepathd1.h>


template< typename LNM, typename T >
class powell03
{
public:

  LNM opt;
  /** The number of iterations for one minimization. */
  uint optiterations;

  uintc dim;

  T * memoryblock;

  /** Positions */
  T * pos;
  /** Direction vectors */
  T * dir;
  /** function state */
  T * xi;

  /** The start of the direction stack. */
  uintc dir_current;

  T ** linepathx0;
  T ** linepathdi;

  /** O(fn) multiplier */
  uint optiterations;

  powell03(LNM opt_)
    : opt(opt_), optiterations(10), dim(opt.linepath.dim)
  {
    memoryblock = new T[dim*dim*2 + dim*2];
    uint memi=0;
    pos = & memoryblock[memi];
    memi += dim*dim;
    dir = & memoryblock[memi];
    memi += dim*dim;
    xi = & memoryblock[memi];
    memi += dim;
    u0 = & memoryblock[memi];
    memi += dim;
  
    opt.linepath.xi = xi;
    linepathx0 = & opt.linepath.xo;
    linepathdi = & opt.linepath.di;
  }

  ~powell03()
    { delete[] memoryblock; }

  /** Reset the direction vectors to be perpendicular unit vectors. */
  void dirUnitVectors();

protected:

  /** The k'th line is minimized with the new state in xi. */
  void minimize(uintc k)
  {
    *linepathx0 = pos+dim*k; 
    uintc k2 = (dir_current + k ) % dim;
    *linepathd1 = dir+dim*k2;

    opt.reset(0.0,5.0);
    for (uint i=0; i<optiterations; ++i)
      ++opt;
  }

};


template< typename LNM, typename T >
void powell03<LNM,T>::dirUnitVectors()
{
  for (uint i=0; i<dim*dim; ++i)
    dir[i]=0;
  for (uint i=0; i<dim; ++i)
    dir[i*N+i]=(T)1;
}
    





#endif