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

Go to the documentation of this file.

#ifndef OPENTISSUE_CORE_MATH_OPTIMIZATION_PROJECTED_STEEPEST_DESCENT_H
#define OPENTISSUE_CORE_MATH_OPTIMIZATION_PROJECTED_STEEPEST_DESCENT_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/big/big_types.h>
#include <OpenTissue/core/math/optimization/optimization_constants.h>
#include <OpenTissue/core/math/optimization/optimization_armijo_projected_backtracking.h>
#include <OpenTissue/core/math/optimization/optimization_stationary_point.h>

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

#include <cmath>
#include <stdexcept>
#include <cassert>

// Convexity parameter
extern double convex_lambda;

namespace OpenTissue
{
  namespace math
  {
    namespace optimization
    {

      // Hauberg: I have hacked this function a bit to support convex IK.
      template < 
        typename T
        , typename function_functor
        , typename gradient_functor
        , typename projection_operator
      >
      inline void projected_steepest_descent(
      function_functor  & f
      , gradient_functor & nabla_f
      , boost::numeric::ublas::vector<T> & x
      , projection_operator const & P
      , size_t const & max_iterations
      , T      const & absolute_tolerance
      , T      const & relative_tolerance
      , T      const & stagnation_tolerance
      , size_t       & status
      , size_t       & iteration
      , T            & error
      , T      const & alpha
      , T      const & beta
      , ublas::vector<T> * profiling = 0
      )  
      {
        using std::fabs;
        using std::min;
        using std::max;

        typedef          ublas::compressed_matrix<T>       matrix_type;
        typedef          ublas::vector<T>                  vector_type;
        typedef          T                                 real_type;
        typedef          OpenTissue::math::ValueTraits<T>  value_traits;

        if(max_iterations <= 0)
          throw std::invalid_argument("max_iterations must be larger than zero");
        if(absolute_tolerance < value_traits::zero() )
          throw std::invalid_argument("absolute_tolerance must be non-negative");
        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");
        if(profiling == &x)
          throw std::logic_error("profiling must not point to x-vector");

        error             = value_traits::infinity();
        iteration = 0;

        status = OK;

        size_t const m = x.size();
        if(m==0)
          return;

        status = ITERATING; // Indicate that we are iterating and have not converged

        if(profiling)
        {
          (*profiling).resize( max_iterations );
          (*profiling).clear();
        }

        // Declare temporary storage
        vector_type dx;
        vector_type x_old;
        vector_type x0;
        vector_type nabla_f_k;

        // Allocate space for temporaries
        dx.resize(m);
        x_old.resize(m);
        nabla_f_k.resize(m);
        x0.resize (m);
        x0.assign (x);
        
        // Initialize 

        x = P(x); // Make sure that the initial x-value is a feasible iterate!
        real_type             f_0   = f(x);
        ublas::noalias( nabla_f_k ) = nabla_f(x) - convex_lambda * (x - x0);

        // Iterate until convergence
        for (; iteration < max_iterations; ++iteration)
        {
          if(profiling)
            (*profiling)(iteration) = f_0;

          // Check for absolute convergence
          if(stationary_point( nabla_f_k, absolute_tolerance, error ) )
          {
            status = ABSOLUTE_CONVERGENCE;
            return;
          }

          // Obtain descent direction
          dx.assign( -nabla_f_k );

          // Perform a projected line-search along the descent direction
          x_old.assign( x );
          real_type f_tau = f_0;
          armijo_projected_backtracking(
            f
            , nabla_f_k
            , x_old
            , x
            , dx
            , relative_tolerance
            , stagnation_tolerance
            , alpha
            , beta
            , f_tau
            , status
            , P
            );

          if(status != OK)
            return;

          // Update values for next iteration
          f_0 = f_tau;
          ublas::noalias( nabla_f_k ) = nabla_f(x) - convex_lambda * (x - x0);

        }//end for loop
      }



      template < 
        typename T
        , typename function_functor
        , typename gradient_functor
        , typename projection_operator
      >
      inline void projected_steepest_descent(
      function_functor  & f
      , gradient_functor & nabla_f
      , boost::numeric::ublas::vector<T> & x
      , projection_operator const & P
      , size_t const & max_iterations
      , T      const & absolute_tolerance
      , T      const & relative_tolerance
      , T      const & stagnation_tolerance
      , size_t       & status
      , size_t       & iteration
      , T            & error
      , T      const & tau
      , ublas::vector<T> * profiling = 0
      )  
      {
        using std::fabs;
        using std::min;
        using std::max;

        typedef          ublas::compressed_matrix<T>       matrix_type;
        typedef          ublas::vector<T>                  vector_type;
        typedef          T                                 real_type;
        typedef          OpenTissue::math::ValueTraits<T>  value_traits;

        if(max_iterations <= 0)
          throw std::invalid_argument("max_iterations must be larger than zero");
        if(absolute_tolerance < value_traits::zero() )
          throw std::invalid_argument("absolute_tolerance must be non-negative");
        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(profiling == &x)
          throw std::logic_error("profiling must not point to x-vector");
        if(tau < value_traits::zero() )
          throw std::invalid_argument("step size must be a positive number");

        error             = value_traits::infinity();
        iteration = 0;

        status = OK;

        size_t const m = x.size();
        if(m==0)
          return;

        status = ITERATING; // Indicate that we are iterating and have not converged

        if(profiling)
        {
          (*profiling).resize( max_iterations );
          (*profiling).clear();
        }

        // Declare temporary storage
        vector_type dx;
        vector_type x_new;
        vector_type nabla_f_k;

        // Allocate space for temporaries
        dx.resize(m);
        x_new.resize(m);
        nabla_f_k.resize(m);

        // Initialize 

        x = P(x); // Make sure that the initial x-value is a feasible iterate!
        real_type             f_0   = f(x);
        ublas::noalias( nabla_f_k ) = nabla_f(x);

        // Iterate until convergence
        for (; iteration < max_iterations; ++iteration)
        {
          if(profiling)
            (*profiling)(iteration) = f_0;

          // Check for absolute convergence
          if(stationary_point( nabla_f_k, absolute_tolerance, error ) )
          {
            status = ABSOLUTE_CONVERGENCE;
            return;
          }

          dx.assign( -nabla_f_k );

          x_new = P(x + tau*dx);
          T f_1 = f(x_new);

          assert( is_number( f_1 ) || !"projected_steepest_descent(): internal error, NAN is encountered?");

          T gamma = ublas::inner_prod( nabla_f_k, x_new - x );
          assert( is_number( gamma ) || !"projected_steepest_descent(): internal error, NAN is encountered?");

          // Test if we have non descent direction
          if( gamma >= value_traits::zero() )
          {
            status = NON_DESCENT_DIRECTION;
            f_1 = f_0;
            x_new.assign( x );
          }

          if(stagnation( x, x_new, stagnation_tolerance ) )
          {
            status = STAGNATION;
          }

          if(relative_convergence(f_0,f_1, relative_tolerance) )
          {
            status = RELATIVE_CONVERGENCE;
          }

          x.assign( x_new );
          f_0 = f_1;
          nabla_f_k.assign( nabla_f(x) );

          if(status != ITERATING)
            return;

        }//end for loop
      }


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

// OPENTISSUE_CORE_MATH_OPTIMIZATION_PROJECTED_STEEPEST_DESCENT_H
#endif