powell02.h

Go to the documentation of this file.

#ifndef POWELL_H  
#define POWELL_H   
#include <cassert>
#include <iostream>
using namespace std

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

#include <cirbuffarr.h>

#include <linepathd1.h>

;
template< typename EXP, typename LNM >
class powell02
{
public:

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

  uintc N;

  EXP exp;

  cirbuffarr< T > xi0;

  cirbuffarr< T > di0;

  T* pi;

  T* linex0;
  T* linedi;
  T linet0;
  T linet1;

 
  linepathd1<FN,XI,XI,T> path;
  LNM linemin;

  uint lineiterations;

  powell02
  (
    uintc N_, 
    uintc lineiterations_,
    T const h0step=20.0, 
    uintc indexmax_=20
  );

  powell02
  (
    FN fn_,
    uintc N_, 
    uintc lineiterations_,
    T const h0step=20.0, 
    uintc indexmax_=20
  );

  ~powell02();

  void reset( T const * x0 );

  void operator ++ ();

  void pushbackadirection();

};


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

template < typename EXP, typename LNM >
void powell02<EXP,LNM>::pushbackadirection()
{
  static T * p0;
  static T * pN;
  static T * d0;

  ++di0;

  p0 = xi0[0];
  pN = xi0[N];
//cout << "p0: " << printvecfunc(p0,N) << endl;
//cout << "pN: " << printvecfunc(pN,N) << endl;
  d0 = di0[0];
  for (uint i=0; i<N; ++i)
    d0[i] = p0[i] - pN[i];
}


template < typename EXP, typename LNM >
powell02<EXP,LNM>::powell02
(
  uintc N_, 
  uintc lineiterations_,
  T const h0step, 
  uintc indexmax_
)
  : N(N_), exp(N_,h0step,indexmax_), xi0(N_*N_*N_,N_),
  di0(N_,N_),
  linex0(new T[N_]), linedi(new T[N_]),
  linet0(-1.0), linet1(1.0), 
  path(N_,exp.fn,exp.xi,linex0,linedi),
  linemin(path),
  lineiterations(lineiterations_)
{ 
}

template < typename EXP, typename LNM >
powell02<EXP,LNM>::powell02
(
  FN fn_,
  uintc N_, 
  uintc lineiterations_,
  T const h0step, 
  uintc indexmax_
)
  : N(N_), exp(fn_,N_,h0step,indexmax_), xi0(N_*N_*N_,N_), 
  di0(N_,N_),
  linex0(new T[N_]), linedi(new T[N_]),
  linet0(-1.0), linet1(1.0), 
  path(N_,EXP::fn,EXP::xi,linex0,linedi),
  linemin(path),
  lineiterations(lineiterations_)
{ 
}


template < typename EXP, typename LNM >
powell02<EXP,LNM>::~powell02()
{
  delete[] linex0;
  delete[] linedi;
}



template < typename EXP, typename LNM >
void powell02<EXP,LNM>::reset( T const * x0 )
{
//cout << "powell02<EXP,LNM>::reset" << endl;

  exp.reset(x0);
//cout << printvecfunc(exp.xi,N) << endl;
  xi0.push(exp.xi);
  for (uint i=1; i<N; ++i)
  {
    exp.move();
//cout << SHOW(exp.fn.counter) << endl;
    xi0.push(exp.xi);
//cout << printvecfunc(exp.xi,N) << endl;
  }

  for (uint i=0; i<N; ++i)
  {
    exp.move();
    xi0.push(exp.xi);
    pushbackadirection();
  }
//cout << "xi0" << endl;
//cout << xi0 << endl;
//cout << "di0" << endl;
//cout << di0 << endl;

}

template < typename EXP, typename LNM >
void powell02<EXP,LNM>::operator ++ ()
{
  T* p0 = & exp.xi[0];
  T* d0;

  T fval0=exp.fmin;

//cout << "powell ++ operator" << endl;

  for (uint i=0; i<N; ++i)
  {
    d0 = di0[N-i-1];

    for (uint k=0; k<N; ++k)
    {
      linex0[i] = p0[i];
      linedi[i] = d0[i];

      linemin.reset(linet0,linet1);
      for (uint j=0; j<lineiterations; ++j)
      {
        ++linemin;
      }
//cout << "***" <<  printvecfunc(exp.xi,N) << endl;
    }

  }
  //exp.move();

  if (fval0<exp.fmin)
  {
cout << "fmin failed" << endl;
    xi0.copyto(exp.xi);
    exp.move();
  }

  xi0.push(exp.xi);
  pushbackadirection();
}


#endif