newton_d.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */

#ifndef PRIMATOLOGIST_OPTIMIZATION_NEWTON_D_H
#define PRIMATOLOGIST_OPTIMIZATION_NEWTON_D_H

#include <primatologist/optimization/newton_decl.h>

namespace aims {

  //--------------------------------------------------------------------------
  //   CONSTRUCTOR / DESTRUCTOR / COPY
  //--------------------------------------------------------------------------

  template <typename O, template <typename> class L>
  Newton<O,L>::Newton( const Objective & f, const Vector & x ):
    _f(f), _x(x)
  {
    setGamma();
    setStopCriterion();
    setStopValue();
    setLineSearch();
    setMaximize();
    setMaxIterations();
    setVerbose();
  }

  template <typename O, template <typename> class L>
  Newton<O,L>::Newton( const Newton<O,L> & other ):
    _f(other._f),
    _x(other._x),
    _gamma(other._gamma),
    _ls(other._ls),
    _verbose(other._verbose),
    _maximize(other._maximize),
    _max_it(other._max_it)
  {}

  template <typename O, template <typename> class L>
  Newton<O,L>::~Newton()
  {}

  template <typename O, template <typename> class L>
  Newton<O,L> &
  Newton<O,L>::operator=( const Newton<O,L> & other )
  {
    if( this != &other )
    {
      _f = other._f;
      _x = other._x;
      _gamma = other._gamma;
      _ls = other._ls;
      _verbose = other._verbose;
      _maximize = other._maximize;
      _max_it = other._max_it;
    }
    return *this;
  }

  //--------------------------------------------------------------------------
  //   SET PARAMETERS
  //--------------------------------------------------------------------------

  template <typename O, template <typename> class L>
  void Newton<O,L>::setGamma( float gamma )
  {
    _gamma = gamma;
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setStopCriterion( const std::string & crit )
  {
    _crit = crit;
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setStopValue( float epsilon )
  {
    _epsilon = epsilon;
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setLineSearch( const LineSearch & ls )
  {
    _ls = ls;
    lineSearch().setObjectiveFunction( objectiveFunction() );
    lineSearch().setMaximize( maximize() );
    lineSearch().setVerbose( verbose() );
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setMaximize( bool maximize )
  {
    _maximize = maximize;
    lineSearch().setMaximize(maximize);
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setMaxIterations( int max )
  {
    _max_it = max;
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setVerbose( int v )
  {
    _verbose = v;
    lineSearch().setVerbose(v);
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setObjectiveFunction( const Objective & f )
  {
    _f = f;
    lineSearch().setObjectiveFunction(f);
  }

  template <typename O, template <typename> class L>
  void Newton<O,L>::setPosition( const Vector & x )
  {
    _x = x;
  }

  //--------------------------------------------------------------------------
  //   GET PARAMETERS
  //--------------------------------------------------------------------------

  template <typename O, template <typename> class L>
  float Newton<O,L>::gamma() const
  {
    return _gamma;
  }

  template <typename O, template <typename> class L>
  const std::string & Newton<O,L>::stopCriterion() const
  {
    return _crit;
  }

  template <typename O, template <typename> class L>
  float Newton<O,L>::stopValue() const
  {
    return _epsilon;
  }

  template <typename O, template <typename> class L>
  typename Newton<O,L>::LineSearch &
  Newton<O,L>::lineSearch()
  {
    return _ls;
  }

  template <typename O, template <typename> class L>
  bool Newton<O,L>::maximize() const
  {
    return _maximize;
  }

  template <typename O, template <typename> class L>
  int Newton<O,L>::maxIterations() const
  {
    return _max_it;
  }

  template <typename O, template <typename> class L>
  int Newton<O,L>::verbose() const
  {
    return _verbose;
  }

  template <typename O, template <typename> class L>
  const typename Newton<O,L>::Vector &
  Newton<O,L>::position() const
  {
    return _x;
  }

  template <typename O, template <typename> class L>
  const typename Newton<O,L>::Objective &
  Newton<O,L>::objectiveFunction() const
  {
    return _f;
  }

  //--------------------------------------------------------------------------
  //   EXECUTE
  //--------------------------------------------------------------------------

  template <typename O, template <typename> class L>
  typename Newton<O,L>::Vector Newton<O,L>::execute()
  {
    if( verbose() > 1 )
    {
      std::cout << std::string(80, '-') << std::endl;
      std::cout << std::setw(7) << "Search "
                << std::setw(4) << "i "
                << std::setw(15) << "f"
                << std::setw(15) << "crit"
                << std::endl;
      std::cout << std::string(80, '-') << std::endl;
    }

    int i =0;
    bool ok;
    Vector current_x = position();
    Vector previous_x;
    float current_value = objectiveFunction().value( position() );
    float previous_value;
    float crit_value;
    Vector dir;

    if( verbose() > 0 )
    {
      std::cout << std::setw(7) << "N"
                << std::setw(4) << i
                << std::setw(15) << std::setprecision(5) << current_value
                << std::setw(15) << std::setprecision(5) << ""
                << std::endl;
    }

    do
    {
      ++i;
      previous_x = current_x;
      previous_value = current_value;
      dir = gamma() * math::invert( objectiveFunction().hessian( previous_x ) )
                    * objectiveFunction().derivative( previous_x );
      if( not maximize() )
        dir *= -1.;
      lineSearch().setPosition( previous_x );
      lineSearch().setSearchDirection( dir );
      current_x = lineSearch().execute();
      current_value = objectiveFunction().value( current_x );

      crit_value = std::abs( current_value - previous_value );
      if( stopCriterion() == "gain" )
        crit_value /= std::abs( previous_value );
      ok = crit_value < stopValue();

      if( verbose() > 0 )
      {
        std::cout << std::setw(7) << "N"
                  << std::setw(4) << i
                  << std::setw(15) << std::setprecision(5) << current_value
                  << std::setw(15) << std::setprecision(5) << crit_value
                  << std::string( ok ? " OK" : " C" )
                  << std::endl;
      }

    } while( not ok and ( maxIterations() < 0 or i < maxIterations() ) );

    if( maxIterations() >= 0 and i >= maxIterations() )
      std::cout << "Maximum number of iterations reached." << std::endl;

    return current_x;
  }



} // namespace aims

#endif // PRIMATOLOGIST_OPTIMIZATION_NEWTON_D_H