#ifndef MINEXPDIM_H
#define MINEXPDIM_H

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

#include <funcstate.h>

/*!
  \brief Exploratory search minimization algorithm.

 Does not use derivatives.  Simplist and likely
 expensive minimization algorithm.
*/
template< typename X >
class minexpdim
{
public:

  /** Function being minimized. */
  funcstate<X>* fn;

  /** The step sizes. */
  X* hi;
  /** N dimensional array of initial step sizes. */
  X* h0; 
  /** From fn. */
  X* xi;
  /** The minimum. */
  X fmin;

// TODO dim from fn to avoid fn->dim references.

// TODO fn to reference rather than pointer.


  /** Is the memory managed? */
  bool mem;
  /** Local index. */
  uint index;
  /** Associated with ++ operator and reset() */
  bool valid;
  /** Set to true when there is a change in state with ++ opertion. */ 
  bool hasminimized;

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

  /** Use this object as an iterator. */
  boolc operator !() const
    { return valid; }
  /** Set the current position as the new minimum. */
  virtual void reset();
  /** 0th order approximator to minimum. */
  virtual void operator ++ ();

  /** Allocate the fn with new. */
  minexpdim
  (
    funcstate<X>& fn_, 
    X const h0step=20.0,
    uintc indexmax_=500,
    boolc mem_=true
  );

  /** Cleanup. */
  virtual ~minexpdim();

  /** Update the new state if successful. */
  boolc minmove();
  /** Update the new state if successful. */
  boolc exploredimension(uintc i);

  /** First order approxmation with N moves. */
  boolc moveOrder1();

  /** Move until a successful minimization. */
  boolc moveOrder1(uintc maxloop);

  /** To increase algorithm stability at cost 
      increase hi[i] on successive search. */
  bool magnitude_increase;

  /** Maximum hi[i]^2. */
  X himax();
  /** Maximum tollerance abs(hi[i]) */
  X himaxtol();

  void h0set(X h);
  void hiset(X h);
};

template< typename X >
class minexpdimN : public minexpdim<X>
{
public:

  /** Allocate the fn with new. */
  minexpdimN
  (
    funcstate<X>& fn_, 
    X const h0step=20.0,
    uintc indexmax_=1000,
    boolc mem_=true
  );

  /** 1st order approximator to minimum. */
  virtual void operator ++ ();

};


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


template < typename X >
boolc minexpdim<X>::moveOrder1()
{
  bool res=false;
  for (uint i=0; i<fn->dim; ++i)
  {
    minexpdim<X>::operator ++ ();
    if (hasminimized)
      res=true;
  }

  hasminimized=res;
  return res;
}


template < typename X >
boolc minexpdim<X>::moveOrder1(uintc maxloop)
{
  for (uint k=0; k<maxloop; ++k)
  {
    if(moveOrder1())
      return true;
  }

  return false;
}

template < typename X >
void minexpdim<X>::operator ++ ()
{
//cout << "minexpdim<X>::operator ++" << endl;
  if (!valid)
    return;

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

  exploredimension(index % fn->dim);

  ++index;
}

template < typename X >
boolc minexpdim<X>::minmove()
{
  (*fn)();
  if (xi[fn->dim]<fmin)
  {
    fmin=xi[fn->dim];
//cout << SHOW(fmin) << endl;
    hasminimized=true;
    return true;
  }
  assert(hasminimized==false);
  return false;
}

template < typename X >
boolc minexpdim<X>::exploredimension(uintc i)
{
  hasminimized=false;

  // Remembering last hi's direction
  xi[i] += hi[i];
  if (minmove())
  {
    if (magnitude_increase)
      hi[i] *= 1.61803398875;
    return true;
  }

  hi[i] *= -1;
  xi[i] += hi[i];
  xi[i] += hi[i];
  if (minmove())
  {
    if (magnitude_increase)
      hi[i] *= 1.61803398875;
    return true;
  }

  xi[i] -= hi[i];
  xi[fn->dim] = fmin;

  // Restore last good direction.
  hi[i] *= -1;

  // Unsuccessful search - decrease h
  // Golden Section Search.
  hi[i] *= 0.61803398875; 

  return false;
}

template < typename X >  
minexpdim<X>::minexpdim
(
  funcstate<X>& fn_, 
  X const h0step,
  uintc indexmax_,
  boolc mem_
)
  : fn(&fn_), xi(fn->xi), mem(mem_),
  index(0), valid(false), indexmax(indexmax_), magnitude_increase(false) 
{
  assert(fn!=0);
  assert(fn->dim>0);
  hi = new X[fn->dim];
  h0 = new X[fn->dim];

  h0set(h0step);
}

template < typename X >
minexpdim<X>::~minexpdim()
{
  if (mem)
  {
    delete[] hi;
    delete[] h0;
    delete fn;
  }
  hi=0;
  h0=0;
  xi=0;
}

template < typename X >
void minexpdim<X>::h0set(X h)
{
  for (uint i=0; i<fn->dim; ++i)
  {
    assert(h!=0);
    h0[i] = h;
  }
}

template < typename X >
void minexpdim<X>::hiset(X h)
{
  for (uint i=0; i<fn->dim; ++i)
  {
    assert(h!=0);
    hi[i] = h;
  }
}

template < typename X >
void minexpdim<X>::reset()
{
  index=0;
  for (uint i=0; i<fn->dim; ++i)
  {
    assert(h0[i]!=0.0);
    hi[i] = h0[i];
  }
  valid=true;

  fmin = (*fn)();

  hasminimized=false;
}

template < typename X >
X minexpdim<X>::himax()
{
  assert(fn->dim!=0);
  X h=hi[0]*hi[0];
  for (uint i=1; i<fn->dim; ++i)
  {
    X h2(hi[i]*hi[i]);
    if (h2>h)
      h=h2;
  }

  return h;
}

template < typename X >
X minexpdim<X>::himaxtol()
{
  assert(fn->dim!=0);
  X h=abs(hi[0]);
  for (uint i=1; i<fn->dim; ++i)
  {
    X h2(abs(hi[i]));
    if (h2>h)
      h=h2;
  }

  return h;
}

template < typename X >
void minexpdimN<X>::operator ++ ()
{
  bool res = minexpdim<X>::moveOrder1();
  minexpdim<X>::hasminimized=res;
}

template < typename X >  
minexpdimN<X>::minexpdimN
(
  funcstate<X>& fn_, 
  X const h0step,
  uintc indexmax_,
  boolc mem_
)
  : minexpdim<X>::minexpdim(fn_,h0step,indexmax_,mem_)
{
}


#endif