derivativeMinimizer.hpp

Go to the documentation of this file.

#ifndef JULIAN_DERIVATIVEMINIMIZER_HPP
#define JULIAN_DERIVATIVEMINIMIZER_HPP

#include <gsl/gsl_multimin.h>
#include <mathematics/numericalAlgorithms/GslMultiminFunctionFdfAdapter.hpp>

namespace julian {

  enum class DerivativeMinimizer {FLETCHER_REEVES,         
                  POLAK_RIBIERE,           
                  BROYDEN_FLETCHER,        
                  STEEPEST_DESCENT         
  };

  template<typename F, typename dF, typename FdF>
  std::vector<double> derivativeMinimizer(F f,
                      dF df,
                      FdF fdf,
                      std::vector<double> x_init,
                      DerivativeMinimizer type,
                      double step_size = 0.1,
                      double tol = 1e-4,
                      double abs = 1e-6,
                      unsigned int max_iter = 200) {

    int n = x_init.size();
    GslMultiminFunctionFdfAdapter<decltype(f), decltype(df), decltype(fdf)> Func(f, df, fdf, n);
    
    size_t iter = 0;
    int status;

    const gsl_multimin_fdfminimizer_type *T;

    switch(type) {
    case DerivativeMinimizer::FLETCHER_REEVES:
      T = gsl_multimin_fdfminimizer_conjugate_fr;
      break;
    case DerivativeMinimizer::POLAK_RIBIERE:
      T = gsl_multimin_fdfminimizer_conjugate_pr;
      break;
    case DerivativeMinimizer::BROYDEN_FLETCHER:
      T = gsl_multimin_fdfminimizer_vector_bfgs2;
      break;
    case DerivativeMinimizer::STEEPEST_DESCENT:
      T = gsl_multimin_fdfminimizer_steepest_descent;
      break;
    }
    
    gsl_multimin_fdfminimizer *s;
    gsl_vector *x;
    x = gsl_vector_alloc (n);
    for (int i = 0; i < n; ++i) {
      gsl_vector_set(x, i, x_init.at(i));
    }
    s = gsl_multimin_fdfminimizer_alloc (T, n);
    gsl_multimin_fdfminimizer_set (s, &Func, x, step_size, tol);

    do {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate (s);
      if (status)
        break;
      status = gsl_multimin_test_gradient (s->gradient, abs);
    } while (status == GSL_CONTINUE && iter < max_iter);

    std::vector<double> ret;
    for (unsigned int i = 0; i < n; ++i) {
      ret.push_back(gsl_vector_get(s -> x, i));
    }

    gsl_multimin_fdfminimizer_free (s);
    gsl_vector_free (x);
    
    return ret;
  }

  template<typename F>
  std::vector<double> derivativeMinimizer(F f,
                      std::vector<double> x_init,
                      DerivativeMinimizer type,
                      double step_size = 0.1,
                      double tol = 1e-4,
                      double increment = 1e-5,
                      double abs = 1e-6,
                      unsigned int max_iter = 200) {
    int n = x_init.size();
    
    auto inverse_increment = 1.0 / increment;
    auto df = [f,n,increment,inverse_increment](std::vector<double>  x)->std::vector<double> {
      std::vector<double> ret;
      double val = f(x);
      for (int i = 0; i < n; ++i) {
    x[i] += increment;
    ret.push_back( (f(x) -val) * inverse_increment);
    x[i] -= increment;
      }      
      return ret;
    };

    auto fdf = [f,n,increment,inverse_increment](std::vector<double>  x)-> std::pair<double, std::vector<double> > {
      std::vector<double> ret;
      double val = f(x);
      auto inverse_increment = 1.0 / increment;
      for (int i = 0; i < n; ++i) {
    x[i] += increment;
    ret.push_back( (f(x) -val) * inverse_increment);
    x[i] -= increment;
      }      
      return std::make_pair(val,ret);
    };
    
    GslMultiminFunctionFdfAdapter<decltype(f), decltype(df), decltype(fdf)> Func(f, df, fdf, n);
    
    size_t iter = 0;
    int status ;

    const gsl_multimin_fdfminimizer_type *T;

    switch(type) {
    case DerivativeMinimizer::FLETCHER_REEVES:
      T = gsl_multimin_fdfminimizer_conjugate_fr;
      break;
    case DerivativeMinimizer::POLAK_RIBIERE:
      T = gsl_multimin_fdfminimizer_conjugate_pr;
      break;
    case DerivativeMinimizer::BROYDEN_FLETCHER:
      T = gsl_multimin_fdfminimizer_vector_bfgs2;
      break;
    case DerivativeMinimizer::STEEPEST_DESCENT:
      T = gsl_multimin_fdfminimizer_steepest_descent;
      break;
    }
    
    gsl_multimin_fdfminimizer *s;
    gsl_vector *x;
    
    x = gsl_vector_alloc (n);

    for (int i = 0; i < n; ++i) {
      gsl_vector_set(x, i, x_init.at(i));
    }

    s = gsl_multimin_fdfminimizer_alloc (T, n);
    gsl_multimin_fdfminimizer_set (s, &Func, x, step_size, tol);

    do {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate (s);
      if (status) {
        break;
      }
      status = gsl_multimin_test_gradient (s->gradient, abs);
    } while (status == GSL_CONTINUE && iter < max_iter);

    std::vector<double> ret;
    for (int i = 0; i < n; ++i) {
      ret.push_back(gsl_vector_get(s -> x, i));
    }
    
    gsl_multimin_fdfminimizer_free (s);
    gsl_vector_free (x);
    
    return ret;
  }
}  // namespace julian
#endif