#ifndef EXPLOREPD1_H
#define EXPLOREPD1_H

#include <cassert>
#include <iostream>
using namespace std

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

template < typename EXP, typename LNM >
class explorepd1 : public EXP
{
public:

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

  /** The multiplier on hi[] for partial derivative calc. */
  T delta;

  /** x0 in line minimization. */
  T* linex0;
  /** Direction in line minimization. */
  T* linedi;
  /** The initial left t interval point. */
  T linet0;
  /** The initial right t interval point. */
  T linet1;

  /** The line minimizer. */
  LNM linemin;

  /** How many times to evaluate the line minimizer per iteration. */
  uint lineiterations;

  /** Construct a N dimensional direct search object on FN */
  explorepd1
  (
    uintc N_, 
    uintc lineiterations_,
    T const h0step=20.0, 
    uintc indexmax_=100
  );

  /** FN can be a reference. */
  explorepd1
  (
    FN fn_,
    uintc N_, 
    uintc lineiterations_,
    T const h0step=20.0, 
    uintc indexmax_=100
  );

  /** Cleanup. */
  ~explorepd1();

  void reset( T const * x0 );

  /** Iterate towards the minimum. */
  void operator ++ ();

};

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

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

template < typename EXP, typename LNM >
explorepd1<EXP,LNM>::explorepd1
(
  uintc N_, 
  uintc lineiterations_,
  T const h0step, 
  uintc indexmax_
)
  : EXP(N_,h0step,indexmax_),
  delta(0.01),
  linex0(new T[N_]), linedi(new T[N_]),
  linet0(-2.0), linet1(2.0), 
  linemin(N_,EXP::fn,EXP::xi,linex0,linedi),
  lineiterations(lineiterations_)
{ 
}

template < typename EXP, typename LNM >
explorepd1<EXP,LNM>::explorepd1
(
  FN fn_,
  uintc N_, 
  uintc lineiterations_,
  T const h0step, 
  uintc indexmax_
)
  : EXP(fn_,N_,h0step,indexmax_), 
  delta(0.01),
  linex0(new T[N_]), linedi(new T[N_]),
  linet0(-2.0), linet1(2.0), 
  linemin(N_,EXP::fn_,EXP::xi,linex0,linedi),
  lineiterations(lineiterations_)
{ 
}

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



template < typename EXP, typename LNM >
void explorepd1<EXP,LNM>::operator ++ ()
{
  //EXP::operator++();
  EXP::move();
//cout << printvecfunc(EXP::fn.xi,EXP::N) << endl;
//cout << SHOW(EXP::fmin) << endl;


  T x;
  T fnb;
  T fna;

  T h;

  for (uint i=0; i<EXP::N; ++i)
  {
    linedi[i] = 0;
    linex0[i] = EXP::fn.xi[i];

  }

  for (uint i=0; i<EXP::N; ++i)
  {
    h = EXP::hi[i]*delta;
    x = EXP::fn.xi[i];
    EXP::fn.xi[i] -= h;
    EXP::fn(fna);
    EXP::fn.xi[i] = x + h;
    EXP::fn(fnb);
    linex0[i] = x;
    linedi[i] = (fnb-fna)/h*0.5;

    //EXP::fn.xi = x;


    // Minimize on the line.
    linemin.reset(linet0,linet1);
    for (uint k=0; k<lineiterations; ++k)
    {
      ++linemin;
    }

    linedi[i] = 0;
    for (uint k=0; k<EXP::N; ++k)
      linex0[k] = EXP::fn.xi[k];

  }


}


#endif