#ifndef EXPLORELINE_H
#define EXPLORELINE_H

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

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


/*!
  \brief explore and the apply line minimization.

  This algorithm explores to a new position, then
  takes the previous position to construct a direction
  and minimizes on that direction.

  If the line minimization fails it backtracks and uses
  the explorer.

  The client controls iteration. Calling operator ++ does one iteration.
  
  The client can configures the number of line iterations in the
  constructor.  This is instead of a tollerance parameter and is
  a design choice.  

  The default interval parameters linet0 and linet1 can be changed
  by writing to them directly.

\par Example
\verbatim
  typedef exploreh<parab2,double*,double> exp2;
  typedef lineoptimizergold<parab2&,double*,double*,double> lnmin;
  exploreline< exp2, lnmin > w(3,lineiter);
  double x0[3] = { 0.0, 0.0, 0.0 };
  w.reset(x0);
  for (uint i=0; i<k; ++i)
    ++w;
  cout << w.xi[0] << " " << w.xi[1] << " " << w.xi[2] << endl;
\endverbatim
*/
;  //It has been years (10+) since I had to put such a hack into a compiler.
template < typename EXP, typename LNM >
class exploreline : public EXP
{
public:

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

  /** Previous position. */
  T* xi0;
  /** Temporary to save explore move. */
  T* yi0;
  /** Temporary to save the hi values. */
  //T* hi0;
  /** 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 */
  exploreline
  (
    uintc N_, 
    uintc lineiterations_,
    T const h0step=20.0, 
    uintc indexmax_=100
  );

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

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

  void reset( T const * x0 );

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

};


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

template < typename EXP, typename LNM >
void exploreline<EXP,LNM>::reset( T const * x0 )
{
//cout << "exploreline<EXP,LNM>::reset(T const *)" << endl;
  EXP::reset(x0);
  EXP::move();
  for (uint i=0; i<EXP::N; ++i)
    xi0[i] = EXP::xi[i];
  EXP::move();
}

template < typename EXP, typename LNM >
exploreline<EXP,LNM>::exploreline
(
  uintc N_, 
  uintc lineiterations_,
  T const h0step, 
  uintc indexmax_
)
  : EXP(N_,h0step,indexmax_),
  xi0(new T[N_]), yi0(new T[N_]), //hi0(new T[N_]), 
  linex0(new T[N_]), linedi(new T[N_]),
  linet0(-1.0), linet1(1.0), 
  linemin(N_,EXP::fn,EXP::xi,linex0,linedi),
  lineiterations(lineiterations_)
{ 
}

template < typename EXP, typename LNM >
exploreline<EXP,LNM>::exploreline
(
  FN fn_,
  uintc N_, 
  uintc lineiterations_,
  T const h0step, 
  uintc indexmax_
)
  : EXP(fn_,N_,h0step,indexmax_), 
  xi0(new T[N_]), yi0(new T[N_]), //hi0(new T[N_]), 
  linex0(new T[N_]), linedi(new T[N_]),
  linet0(-1.0), linet1(1.0), 
  linemin(N_,EXP::fn_,EXP::xi,linex0,linedi),
  lineiterations(lineiterations_)
{ 
}

template < typename EXP, typename LNM >
exploreline<EXP,LNM>::~exploreline()
{
  delete[] xi0;
  delete[] yi0;
  //delete[] hi0;
  delete[] linex0;
  delete[] linedi;
}


/*
Note yi0 can be eliminated. After setting the line,
 xi0[] = xi[] ...
Because a move is made at the end.  
 This algorithm has been very problamatic and perhaps
 I am going to have to wait for a better understanding.
*/

template < typename EXP, typename LNM >
void exploreline<EXP,LNM>::operator ++ ()
{
  // Preserve xi
  for (uint i=0; i<EXP::N; ++i)
    yi0[i] = EXP::xi[i];
  T fval0 = EXP::fmin;

//cout << "xi0: " << printvecfunc(xi0,EXP::N) << endl;
//cout << "xi:  " << printvecfunc(EXP::xi,EXP::N) << endl;
  

  // Set the new direction and starting point.
  for (uint i=0; i<EXP::N; ++i)
  {
    linex0[i] = EXP::xi[i];
    linedi[i] = EXP::xi[i] - xi0[i];
  }

//cout << "linex0: " << printvecfunc(linex0,EXP::N) << endl;
//cout << "linedi: " << printvecfunc(linedi,EXP::N) << endl;
//cout << SHOW(EXP::fmin) << endl;

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

//linemin.printstate();
//cout << SHOW(EXP::fmin) << endl;
  }

  T fval = linemin.minimumrealize();
//cout << SHOW(fval) << endl;
  if (fval < fval0)
  {
    EXP::fmin = fval;
//cout << "***" << endl;
  }
  else
  {
    cout << "error:  line minimization failed." << endl;
    cout << "  backtracking and defaulting to explore." << endl;
    for (uint i=0; i<EXP::N; ++i)
      EXP::xi[i] = yi0[i];
    EXP::fmin = fval0;

    //EXP::move();
  }
  EXP::move();

  for (uint i=0; i<EXP::N; ++i)
    xi0[i] = yi0[i];



//cout << "xi=" << printvecfunc(EXP::xi,EXP::N) << endl;

}

#endif