#ifndef EXPLOREH_H
#define EXPLOREH_H


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

typedef unsigned int uint;
typedef unsigned int const uintc;

#define SHOW(x) #x << '=' << (x)

/*!
  \brief N dimensional Pattern search minimization algorithm.

  This is a 0 order aproximator.  When all N variables are
  iterated over it becomes a 1st order approximator.

  F is the function with the following interface.
  void FN::operator( T & val ) 
  T * FN::xi
  The state xi is owned by FN.
  xi supports vector operations [].

  be owned by either the minimizer (the default) or the
  function F.

*/
template < typename FN, typename XI, typename T >
class exploreh
{
protected:
  /** Frees resources. */
  void clean();
  /** Allocate memory and state. Excludes x and h! */
  void init();
public:

  /** XI type information exported. */
  typedef FN FNtype;
  /** XI type information exported. */
  typedef XI XItype;
  /** T type information exported. */
  typedef T Ttype;

  /** The function to be minimized. */
  FN fn;

  /** Local index. */
  uint index;
  /** Associated with ++ operator and reset() */
  bool valid;
  /** Set to true when there is a change in state with ++ opertion. */ 
  bool hasStateChanged;
  /** The minimum associated with the current state x. */
  T fmin;
  /** The step sizes. */
  T * hi;
  /** Set all the h values to the new step length. */ 
  void h0Set( T const hnew);


  /** Update the new state if successful. */
  bool const evaluate()
  {
    T f1;
    fn(f1);
    if (f1<fmin)
    {
      fmin=f1;
//cout << SHOW(fmin) << endl;
      hasStateChanged=true;
      return true;
    }
    return false;
  }

  /** Move to local minimum from current position, ignores the 
      current global minimum. */
  void movelocal();

  /** Perterb N dimensions till the state is changed. */
  void move();

  /** The number of dimensions. */
  uintc N;

  /** N dimensional array of initial step sizes. */
  T * h0; 

  /** N dimensional array of the current state. */
  XI xi;  

  /** Set the pointers manually. eg if fn does not manage xi. */
  void xiSet(XI xi_)
    { xi = fn.xi = xi_; }

  /** Each index has a counter, if during iteration this is 
      exceeded valid is set to false. */
  uint indexmax;


  // Set with initial high h values to 
  // avoid problem of never converging with small h values.

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

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

  /** Clean up memory allocations. */
  ~exploreh();

  /** Use this object as an iterator. */
  bool const operator !() const
    { return valid; }

  /** Set the current position as the new minimum. */
  void reset();


  /** The current state is evaluated for a new minimum.
     The step sizes h are reset to h0.
     The index counters are set to 0.
     xi is set to x0.
  */
  void reset( T const * x0 );

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

  /** Display the minimum, hi, and xi. */
  ostream & print(ostream & os) const;

};


/*
 *---------------------------------------------------------------
 * Implementation
 */


/* silly
template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::resetposition()
{
  T f1;
  fn(f1);
  fmin=f1;
}
*/



template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::movelocal()
{

#ifndef NDEBUG
if (valid==false)
  cout << "error" << endl;
  assert(valid==true);
#endif

  // preserve fmin;
  T fmin0 = fmin;
  fn(fmin);
  //resetposition();
/*
  T f1;
  fn(f1);
  fmin=f1;
*/

  hasStateChanged=false;

  for ( ; (hasStateChanged==false)&&valid; )
  {
    operator ++ ();
  }

  // Restore fmin.
  fmin = fmin0;
}


template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::move()
{
#ifndef NDEBUG
if (valid==false)
  cout << "error:  exploreh.valid is false" << endl;
//  assert(valid==true);
#endif
  hasStateChanged=false;

  for ( ; (hasStateChanged==false)&&valid; )
  {
    operator ++ ();
  }

}


template < typename FN, typename XI, typename T >
ostream & exploreh<FN,XI,T>::print(ostream & os) const
{
  os << "fmin=" << fmin;

  os << " h: " << hi[0];
  for (uint i=1; i<N; ++i)
    os << " " << hi[i];
  os << " x: " << xi[0];
  for (uint i=1; i<N; ++i)
    os << " " << xi[i];

  os << endl;

  return os;
}


template < typename FN, typename XI, typename T >
exploreh<FN,XI,T>::exploreh(uintc N_, T const h0step, uintc indexmax_)
  : valid(false), N(N_), indexmax(indexmax_)
{
  init();
  h0Set(h0step);
}

template < typename FN, typename XI, typename T >
exploreh<FN,XI,T>::exploreh
(
  FN fn_,
  uintc N_, 
  T const h0step,
  uintc indexmax_
)
  : fn(fn_), valid(false), N(N_), indexmax(indexmax_)
{
  init();
  h0Set(h0step);
}

template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::clean()
{
  delete[] hi;
  delete[] h0;
  hi=0;
  h0=0;
}

template < typename FN, typename XI, typename T >
exploreh<FN,XI,T>::~exploreh()
{
  clean();
}

template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::init()
{
  assert(N!=0);

  index=0;
  hi = new T[N];
  h0 = new T[N];

  xi = fn.xi;
}


template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::h0Set( T const hnew)
{
  assert(N!=0);

  for (uint i=0; i<N; ++i)
    hi[i] = h0[i] = hnew;
}

template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::reset()
{
  index=0;
  for (uint i=0; i<N; ++i)
  {
    assert(h0[i]!=0.0);
    hi[i] = h0[i];
  }
  valid=true;

  fn(fmin);
}

template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::reset( T const * x0)
{
  index=0;
  for (uint i=0; i<N; ++i)
  {
    assert(h0[i]!=0.0);
    hi[i] = h0[i];
    xi[i] = x0[i];
  }
  valid=true;

  fn(fmin);
}



template < typename FN, typename XI, typename T >
void exploreh<FN,XI,T>::operator ++ ()
{
//cout << "exploreh<FN,XI,T>::operator ++" << endl;

  hasStateChanged=false;

  if (!valid)
    return;

  if (index>indexmax)
  {
    valid=false;
#ifndef NDEBUG
cout << "error: indexmax reached  " << SHOW(index) << endl;
#endif
    return;
  }

  ++index;

  for (uint i=0; i<N; ++i)
  {
//cout << "***" << i << endl;
    xi[i] += hi[i];
    if (evaluate())
      continue;

    hi[i] *= -1;
    xi[i] += hi[i];
    xi[i] += hi[i];
    if (evaluate())
      continue;

    xi[i] -= hi[i];

    // Golden Section Search.
    hi[i] *= 0.61803398875; 
  }
//cout << printvecfunc(xi,N) << endl;
}




#endif