lineopgold.h

Go to the documentation of this file.

#ifndef LINEOPGOLD_H
#define LINEOPGOLD_H

#include <vec.h>
#include <typeop.h>

/*
  \brief  Golden ratio line optimization.

  Find the minimum in an interval.
*/
template< typename EXP >
class lineopgold
{

public:


//  typedef typename typeop<EXP>::Tconst::Tconst Tconst
  typedef typename typeop<EXP>::Tbare::Ttype T;
  typedef typename typeop<EXP>::Tref EXP2;

  typedef T const Tconst;

  EXP2 exp;

  uintc N;

  T* di;

  lineopgold(EXP2 _exp); 
  ~lineopgold();

  T const hlarge() const;

  bool const minimize( T const hstep, uintc kiterations );

protected:

  T* X0;
  T X0fmin; 

  T ti[4];

  T fxk[4];

  static T goldratio;

#ifndef NDEBUG
  bool const verifyti( uintc ai0, uintc ai1, uintc ai2, uintc ai3) const;
#endif

};


template< typename EXP >
typename lineopgold<EXP>::T lineopgold<EXP>::goldratio = 0.61803398875;


template< typename EXP >
bool const lineopgold<EXP>::minimize( T const hstep, uintc kiterations )
{
  assert( kiterations>6); // Make minimization effective, min 4+k function evaluations.
  // Save the current state.
  uint k;
  for ( k=0; k<N; ++k)
    X0[k] = exp.xi[k];
  X0fmin = exp.fn.fmin;

  // Indexes to vectors.
  uint ai0=0;
  uint ai1=1;
  uint ai2=2;
  uint ai3=3;
  uint rej;

  // these are the 4 points about xi which are being sampled.
  T len = goldratio*hstep*2.0;
  ti[0] = -hstep;
  ti[3] = hstep;
  ti[1] = ti[3]-len;
  ti[2] = ti[0]+len; 

  // Evaluate the function at the 4 points.
  for ( k=0; k<N; ++k)
    exp.xi[k] = X0[k] + di[k]*ti[ai0]; 
  exp.fn(fxk[ai0]);

  for ( k=0; k<N; ++k)
    exp.xi[k] = X0[k] + di[k]*ti[ai3];
  exp.fn(fxk[ai3]);

  for ( k=0; k<N; ++k)
    exp.xi[k] = X0[k] + di[k]*ti[ai1];
  exp.fn(fxk[ai1]);

  for ( k=0; k<N; ++k)
    exp.xi[k] = X0[k] + di[k]*ti[ai2];
  exp.fn(fxk[ai2]);

  // Minimization
  for (uint i=0; i<kiterations; ++i)
  {
    len *= goldratio;

    // Test - reject x0 if successful.
    if (fxk[ai2]<fxk[ai1])
    {
      // Reject x0 
      rej = ai0;
      ai0 = ai1;
      ai1 = ai2;
      ai2 = rej;
      ti[ai2] = ti[ai3]-len;
      assert(verifyti(ai0,ai1,ai2,ai3));
      for ( k=0; k<N; ++k)
        exp.xi[k] = X0[k] + di[k]*ti[ai2];
      exp.fn(fxk[ai2]);
    }
    else
    {
      // Reject x3 
      rej = ai3;
      ai3 = ai2;
      ai2 = ai1;
      ai1 = rej;
      ti[ai1] = ti[ai0]+len;
      assert(verifyti(ai0,ai1,ai2,ai3));
      for ( k=0; k<N; ++k)
        exp.xi[k] = X0[k] + di[k]*ti[ai1];
      exp.fn(fxk[ai1]);
    }
  }

  // Find the new minimum.
  if (fxk[ai1]<fxk[ai2])
  {
    if (fxk[ai1]<X0fmin)
    {
      for ( k=0; k<N; ++k)
        exp.xi[k] = X0[k] + di[k]*ti[ai1];
      exp.fmin = fxk[ai1]; 

      return true;
    }
  }
  else
  {
    if (fxk[ai2]<X0fmin)
    {
      for ( k=0; k<N; ++k)
        exp.xi[k] = X0[k] + di[k]*ti[ai2];
      exp.fmin = fxk[ai2]; 

      return true;
    }
  }

  // Restore the original position, no change.
  for ( k=0; k<N; ++k)
    exp.xi[k] = X0[k];
  exp.fmin = X0fmin;

  return false;
}


template< typename EXP >
lineopgold<EXP>::lineopgold(EXP2 _exp) 
  : exp(_exp), N(exp.N), di(new T[N]),  
    X0(new T[N])
{
}


template< typename EXP >
lineopgold<EXP>::~lineopgold()
{
  delete[] di; 
  delete[] X0;
}


template< typename EXP >
typename lineopgold<EXP>::T const lineopgold<EXP>::hlarge() const
{
  T h = exp.hi[0];
  if (h<0)
    h *= -1.0;
  T y;
  for (uint i=1; i<N; ++i)
  {
    y = exp.hi[0];
    if (y<0)
      y += -1.0;

    if (y>h)
      h=y;
  }

  return h;
}


template< typename EXP >
bool const lineopgold<EXP>::verifyti
(
  uintc ai0, 
  uintc ai1, 
  uintc ai2, 
  uintc ai3
) const
{
  bool valid=true;

  if ( ti[ai1] < ti[ai0] )
    valid=false;
  if ( ti[ai2] < ti[ai0] )
    valid=false;
  if ( ti[ai3] < ti[ai0] )
    valid=false;

  if ( ti[ai2] < ti[ai1] )
    valid=false;
  if ( ti[ai3] < ti[ai1] )
    valid=false;

  if ( ti[ai3] < ti[ai2] )
    valid=false;

  if (valid==false)
  {
    cout << "error: ai indexes incorrect" << endl;
    cout << SHOW(ai0) << " " << SHOW(ai1) << " ";
    cout << SHOW(ai2) << " " << SHOW(ai3) << endl;

    cout << SHOW(ti[0]) << " " << SHOW(ti[1]) << " ";
    cout << SHOW(ti[2]) << " " << SHOW(ti[3]) << endl;
  }
     
  return valid;
}




#endif