wolfe_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_WOLFE_D_H
#define PRIMATOLOGIST_OPTIMIZATION_WOLFE_D_H

#include <primatologist/optimization/wolfe_decl.h>
#include <cmath>              // std::abs
#include <iostream>           // std::cout / std::endl
#include <iomanip>            // std::setw / std::setprecision
#include <string>             // std::string

namespace aims {

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

  template <typename O>
  WolfeLineSearch<O>::WolfeLineSearch( const Objective & f,
                                       const Vector & x,
                                       const Vector & p ):
    _f(f),
    _x(x),
    _p(p)
  {
    setC1();
    setC2();
    setStrongWolfe();
    setMaximize();
    setMaxIterations();
    setVerbose();
  }

  template <typename O>
  WolfeLineSearch<O>::WolfeLineSearch( const WolfeLineSearch<O> & other ):
    _f(other._f),
    _x(other._x),
    _p(other._p),
    _c1(other._c1),
    _c2(other._c2),
    _strong(other._strong),
    _maximize(other._maximize),
    _max_it(other._max_it),
    _verbose(other._verbose),
    _previous_value(other._previous_value),
    _previous_scalarprod(other._previous_scalarprod),
    _current_value(other._current_value)
  {}

  template <typename O>
  WolfeLineSearch<O>::~WolfeLineSearch()
  {}

  template <typename O>
  WolfeLineSearch<O> &
  WolfeLineSearch<O>::operator=( const WolfeLineSearch<O> & other )
  {
    if( this != &other )
    {
      _f = other._f;
      _x = other._x;
      _p = other._p;
      _c1 = other._c1;
      _c2 = other._c2;
      _strong = other._strong;
      _maximize = other._maximize;
      _max_it = other._max_it;
      _verbose = other._verbose;
      _previous_value = other._previous_value;
      _previous_scalarprod = other._previous_scalarprod;
      _current_value = other._current_value;
    }
    return *this;
  }

  //--------------------------------------------------------------------------
  //   SET PARAMETERS
  //--------------------------------------------------------------------------
  template <typename O>
  void WolfeLineSearch<O>::setC1( float c1 )
  {
    _c1 = c1;
  }

  template <typename O>
  void WolfeLineSearch<O>::setC2( float c2 )
  {
    _c2 = c2;
  }

  template <typename O>
  void WolfeLineSearch<O>::setPosition( const Vector & x )
  {
    _x = x;
  }

  template <typename O>
  void WolfeLineSearch<O>::setSearchDirection( const Vector & p )
  {
    _p = p;
  }

  template <typename O>
  void WolfeLineSearch<O>::setObjectiveFunction( const Objective & f )
  {
    _f = f;
  }

  template <typename O>
  void WolfeLineSearch<O>::setStrongWolfe( bool strong )
  {
    _strong = strong;
  }

  template <typename O>
  void WolfeLineSearch<O>::setVerbose( int v )
  {
    _verbose = v;
  }

  template <typename O>
  void WolfeLineSearch<O>::setMaximize( bool maximize )
  {
    _maximize = maximize;
  }

  template <typename O>
  void WolfeLineSearch<O>::setMaxIterations( int n )
  {
    _max_it = n;
  }

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

  template <typename O>
  float WolfeLineSearch<O>::c1() const
  {
    return _c1;
  }

  template <typename O>
  float WolfeLineSearch<O>::c2() const
  {
    return _c2;
  }

  template <typename O>
  const typename WolfeLineSearch<O>::Vector &
  WolfeLineSearch<O>::position() const
  {
    return _x;
  }

  template <typename O>
  const typename WolfeLineSearch<O>::Vector &
  WolfeLineSearch<O>::searchDirection() const
  {
    return _p;
  }

  template <typename O>
  const typename WolfeLineSearch<O>::Objective &
  WolfeLineSearch<O>::objectiveFunction() const
  {
    return _f;
  }

  template <typename O>
  bool WolfeLineSearch<O>::strongWolfe() const
  {
    return _strong;
  }

  template <typename O>
  int WolfeLineSearch<O>::verbose() const
  {
    return _verbose;
  }

  template <typename O>
  bool WolfeLineSearch<O>::maximize() const
  {
    return _maximize;
  }

  template <typename O>
  int WolfeLineSearch<O>::maxIterations() const
  {
    return _max_it;
  }

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

  template <typename O>
  typename WolfeLineSearch<O>::Vector
  WolfeLineSearch<O>::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) << "a"
                << std::endl;
      std::cout << std::string(80, '-') << std::endl;
    }

    if( verbose() > 0 )
    {
      std::cout << std::setw(7) << "W"
                << std::setw(4) << 0
                << std::setw(15) << std::setprecision(5) << _previous_value
                << std::setw(15) << std::setprecision(5) << 0.;
    }

    if( !checkDirection() )
    {
      if( verbose() > 0 )
      {
        std::cout << " Not "
                  << std::string( maximize() ? "an ascent" : "a descent" )
                  << std::string( " direction." ) << std::endl;
      }
      return _x;
    }
    else if( verbose() > 0 )
      std::cout << std::endl;

    float a = 1.;
    bool ok = false;
    int i = 0;
    float a_min = 0., a_max;
    bool a_max_set = false;

    do {
      ++i;
      _current_value = objectiveFunction().value( _x + a * _p );
      if( ( maximize() &&_current_value <= _previous_value ) ||
          _current_value >= _previous_value )
      {
        // Because we know the direction is an ascent direction, the objective
        // function should not decrease. If this is the case, it is because
        // we took a step too big and found ourseleves farther than the
        // local maxima.
        // In this case, we thus decrease the step length.
        if( verbose() > 0 )
        {
          std::cout << std::setw(7) << "W"
                    << std::setw(4) << i
                    << std::setw(15) << std::setprecision(5) << _current_value
                    << std::setw(15) << std::setprecision(5) << a
                    << " B"
                    << std::endl;
        }
        ok = false;
        a_max = a;
        a = ( a_max - a_min ) / 2. + a_min;
      }
      else if( not wolfe1(a) or not wolfe2(a) )
      {
        if( verbose() > 0 )
        {
          std::cout << std::setw(7) << "W"
                    << std::setw(4) << i
                    << std::setw(15) << std::setprecision(5) << _current_value
                    << std::setw(15) << std::setprecision(5) << a
                    << " W"
                    << std::endl;
        }
        ok = false;
        a_min = a;
        if( a_max_set )
          a = ( a_max - a_min ) / 2. + a_min;
        else
          a *= 2.;
      }
      else
      {
        if( verbose() > 0 )
        {
          std::cout << std::setw(7) << "W"
                    << std::setw(4) << i
                    << std::setw(15) << std::setprecision(5) << _current_value
                    << std::setw(15) << std::setprecision(5) << a
                    << " OK"
                    << std::endl;
        }
        ok = true;
      }

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

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

    return _x + a * _p;
  }

  //--------------------------------------------------------------------------
  //   LINE SEARCH
  //--------------------------------------------------------------------------

  template <typename O>
  bool WolfeLineSearch<O>::checkDirection()
  {
    _previous_value = objectiveFunction().value( _x );
    _previous_scalarprod = (math::transpose( _p ) * objectiveFunction().derivative( _x )).at(0);
    if( maximize() )
      return _previous_scalarprod > 0;
    else
      return _previous_scalarprod < 0;
  }

  template <typename O>
  bool WolfeLineSearch<O>::wolfe1( float a ) const
  {
    float righthand = _previous_value + _c1 * a * _previous_scalarprod;
    if( maximize() )
      return _current_value >= righthand;
    else
      return _current_value <= righthand;
  }

  template <typename O>
  bool WolfeLineSearch<O>::wolfe2w( float a ) const
  {
    float current_scalarprod = (math::transpose(_p) * objectiveFunction().derivative( _x + a * _p )).at(0);
    float righthand = _c2 * _previous_scalarprod;
    if( maximize() )
      return current_scalarprod <= righthand;
    else
      return current_scalarprod >= righthand;
  }

  template <typename O>
  bool WolfeLineSearch<O>::wolfe2s( float a ) const
  {
    float current_scalarprod = (math::transpose(_p) * objectiveFunction().derivative( _x + a * _p )).at(0);
    float righthand = _c2 * _previous_scalarprod;
    return std::abs(current_scalarprod) <= std::abs(righthand);
  }

  template <typename O>
  bool WolfeLineSearch<O>::wolfe2( float a ) const
  {
    if( _strong )
      return wolfe2s(a);
    else
      return wolfe2w(a);
  }

} // namespace aims

#endif // PRIMATOLOGIST_OPTIMIZATION_WOLFE_D_H