exploreh.h

Go to the documentation of this file.

#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)

template < typename FN, typename XI, typename T >
class exploreh
{
protected:
  void clean();
  void init();
public:

  typedef FN FNtype;
  typedef XI XItype;
  typedef T Ttype;

  FN fn;

  uint index;
  bool valid;
  bool hasStateChanged;
  T fmin;
  T * hi;
  void h0Set( T const hnew);


  bool const evaluate()
  {
    T f1;
    fn(f1);
    if (f1<fmin)
    {
      fmin=f1;
//cout << SHOW(fmin) << endl;
      hasStateChanged=true;
      return true;
    }
    return false;
  }

  void movelocal();

  void move();

  uintc N;

  T * h0; 

  XI xi;  

  void xiSet(XI xi_)
    { xi = fn.xi = xi_; }

  uint indexmax;


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

  exploreh
  (
    uintc N_, 
    T const h0step=20.0, 
    uintc indexmax_=500
  );

  exploreh
  (
    FN fn_,
    uintc N_, 
    T const h0step=20.0, 
    uintc indexmax_=500
  );

  ~exploreh();

  bool const operator !() const
    { return valid; }

  void reset();


  void reset( T const * x0 );

  void operator ++ ();

  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