nonDerivativeMinimizer.hpp

Go to the documentation of this file.

#ifndef JULIAN_NONDERIVATIVEMINIMIZER_HPP
#define JULIAN_NONDERIVATIVEMINIMIZER_HPP

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

namespace julian {

  enum class NonDerivativeMinimizer {NELDER_MEAD_SIMPLEX,        
                     RANDOM_NELDER_MEAD_SIMPLEX  
  };

  template<typename F>
  std::vector<double> nonDerivativeMinimizer(F f,std::vector<double> x_init, NonDerivativeMinimizer t, double abs = 1e-3, int number_of_iterations = 100) {
    int n = x_init.size();

    GslMultiminFunctionAdapter<decltype(f)> Func(f, n);

    gsl_multimin_function *gsl_f = static_cast<gsl_multimin_function*>(&Func); 

    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
    
    switch(t) {
    case NonDerivativeMinimizer::NELDER_MEAD_SIMPLEX:
      T = gsl_multimin_fminimizer_nmsimplex2;
      break;
    case NonDerivativeMinimizer::RANDOM_NELDER_MEAD_SIMPLEX:
      T = gsl_multimin_fminimizer_nmsimplex2rand;
      break;
    }
    
    gsl_multimin_fminimizer* s = gsl_multimin_fminimizer_alloc (T, n);
    
    gsl_vector *ss, *x;
    x = gsl_vector_alloc (n);
    ss = gsl_vector_alloc (n);
    gsl_vector_set_all (ss, 1.0);
    
    for (int i = 0; i < n; i++) {
      gsl_vector_set (x, i, x_init.at(i));
    }
    
    gsl_multimin_fminimizer_set (s, gsl_f, x, ss);
    int status = 0;
    int iter = 0;
    do {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);

      if (status)
        break;
      
      double size = gsl_multimin_fminimizer_size (s);
      status = gsl_multimin_test_size (size, abs);
      
    } while (status == GSL_CONTINUE && iter < number_of_iterations);

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