/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/math/optimization/optimization_armijo_projected_backtracking.h

Go to the documentation of this file.

#ifndef OPENTISSUE_CORE_MATH_OPTIMIZATION_ARMIJO_PROJECTED_BACKTRACKING_H
#define OPENTISSUE_CORE_MATH_OPTIMIZATION_ARMIJO_PROJECTED_BACKTRACKING_H
//
// OpenTissue Template Library
// - A generic toolbox for physics-based modeling and simulation.
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
//
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
//
#include <OpenTissue/configuration.h>

#include <OpenTissue/core/math/math_value_traits.h>
#include <OpenTissue/core/math/math_is_number.h>

#include <OpenTissue/core/math/big/big_types.h>
#include <OpenTissue/core/math/optimization/optimization_constants.h>
#include <OpenTissue/core/math/optimization/optimization_stagnation.h>
#include <OpenTissue/core/math/optimization/optimization_relative_convergence.h>

#include <stdexcept>
#include <cassert>

namespace OpenTissue
{
  namespace math
  {
    namespace optimization
    {

      template < typename T, typename function_type, typename projection_operator >
      inline T armijo_projected_backtracking(
        function_type const & f
        , ublas::vector<T> const & nabla_f
        , ublas::vector<T> const & x
        , ublas::vector<T> & x_tau
        , ublas::vector<T> const & dx
        , T      const & relative_tolerance
        , T      const & stagnation_tolerance
        , T      const & alpha
        , T      const & beta
        , T            & f_tau
        , size_t       & status
        , projection_operator const & P
        )
      {
        using std::fabs;
        using std::max;

        typedef OpenTissue::math::ValueTraits<T>  value_traits;

        T      const TOO_TINY              = ::boost::numeric_cast<T>(0.00001);

        status = OK;

        if(relative_tolerance < value_traits::zero() )
          throw std::invalid_argument("relative_tolerance must be non-negative");

        if(stagnation_tolerance < value_traits::zero() )
          throw std::invalid_argument("stagnation_tolerance must be non-negative");

        if (beta >= value_traits::one() )
          throw std::invalid_argument("Illegal beta value");

        if (alpha <= value_traits::zero() )
          throw std::invalid_argument("Illegal alpha value");

        if(beta<=alpha)
          throw std::invalid_argument("beta must be larger than alpha");

        size_t const m = x.size();

        // Test if we have a problem to work on
        if(m==0)
          return value_traits::zero();

        assert(x.size() == x_tau.size() || !"armijo_projected_backtracking(): x and x_tau have incompatible dimensions");






        T gamma = alpha*ublas::inner_prod( nabla_f, dx );
        assert( is_number( gamma ) || !"armijo_projected_backtracking(): internal error, NAN is encountered?");

        // Test if initial search direction is a descent direction
        if( gamma >= value_traits::zero() )
        {
          status = NON_DESCENT_DIRECTION;
          x_tau.assign( x );
          return value_traits::zero();
        }

        // Do line-search by performing Armijo back-tracking
        T const f_0 = f_tau;
        T       tau = value_traits::one();

        x_tau = P(x + dx);
        f_tau = f(x_tau);
        assert( is_number( f_tau ) || !"armijo_projected_backtracking(): internal error, NAN is encountered?");


        // Test for whether x_tau = x. If this is the
        // case then we can not possible make a line-search along
        // the direction dx. 
        if(stagnation( x, x_tau, value_traits::zero() ) )
        {
          status = DESCEND_DIRECTION_IN_NORMAL_CONE;
          f_tau = f_0;
          x_tau.assign( x );
          return value_traits::zero();
        }

        gamma = alpha*ublas::inner_prod( nabla_f, x_tau - x );
        assert( is_number( gamma ) || !"armijo_projected_backtracking(): internal error, NAN is encountered?");

        // Perform pojected back-tracking line-search
        while( ( f_tau > (f_0 + gamma*tau ) ) && tau > TOO_TINY )
        {
          tau *= beta;
          assert( is_number( tau ) || !"armijo_projected_backtracking(): internal error, NAN is encountered?");
          x_tau = P( x + dx*tau );
          gamma = alpha*ublas::inner_prod( nabla_f, x_tau - x );
          assert( is_number( gamma ) || !"armijo_projected_backtracking(): internal error, NAN is encountered?");
          f_tau = f(x_tau);
          assert( is_number( f_tau ) || !"armijo_projected_backtracking(): internal error, NAN is encountered?");
        }

        // Test if a new step length was computed
        if( (tau < TOO_TINY) && (f_tau > f_0))
          status = BACKTRACKING_FAILED;

        if(stagnation( x, x_tau, stagnation_tolerance ) )
          status = STAGNATION;

        if(relative_convergence(f_0,f_tau, relative_tolerance) )
          status = RELATIVE_CONVERGENCE;

        // Return the new step length to caller
        return tau;
      }

    } // namespace optimization
  } // namespace math
} // namespace OpenTissue

// OPENTISSUE_CORE_MATH_OPTIMIZATION_ARMIJO_PROJECTED_BACKTRACKING_H
#endif