27 #include <gsl/gsl_multimin.h>
28 #include <gsl/gsl_blas.h>
29 #include <queso/GslVector.h>
30 #include <queso/VectorSpace.h>
31 #include <queso/ScalarFunction.h>
32 #include <queso/GslOptimizer.h>
33 #include <queso/OptimizerMonitor.h>
41 double c_evaluate(
const gsl_vector * x,
void * context) {
49 for (
unsigned int i = 0; i < state.sizeLocal(); i++) {
50 state[i] = gsl_vector_get(x, i);
68 gsl_vector * derivative) {
77 for (
unsigned int i = 0; i < state.sizeLocal(); i++) {
78 state[i] = gsl_vector_get(x, i);
88 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
89 gsl_vector_set(derivative, i, GSL_NAN);
104 bool userComputedDerivative =
true;
105 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
108 if (gsl_isnan(deriv[i])) {
109 userComputedDerivative =
false;
114 if (userComputedDerivative) {
115 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
116 gsl_vector_set(derivative, i, -deriv[i]);
124 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
125 double tempState = state[i];
134 state[i] = tempState;
137 if (!gsl_isnan(fx) && !gsl_isnan(fxph)) {
138 gsl_vector_set(derivative, i, (fxph - fx) / h);
141 gsl_vector_set(derivative, i, GSL_NAN);
150 double * f, gsl_vector * derivative) {
160 m_objectiveFunction(objectiveFunction),
161 m_initialPoint(new
GslVector(objectiveFunction.domainSet().
162 vectorSpace().zeroVector())),
163 m_minimizer(new
GslVector(this->m_objectiveFunction.domainSet().
164 vectorSpace().zeroVector())),
165 m_solver_type(BFGS2),
166 m_fstep_size(this->m_objectiveFunction.domainSet().vectorSpace().zeroVector()),
189 std::cerr <<
"Minimization was given initial point outside of domain"
196 zeroVector().sizeLocal();
219 for (
unsigned int i = 0; i < initialPoint.
sizeLocal(); i++) {
237 bool gradient_needed =
false;
247 gradient_needed =
true;
263 return gradient_needed;
269 gsl_vector * x = gsl_vector_alloc(dim);
270 for (
unsigned int i = 0; i < dim; i++) {
275 const gsl_multimin_fdfminimizer_type* type = NULL;
280 type = gsl_multimin_fdfminimizer_conjugate_fr;
283 type = gsl_multimin_fdfminimizer_conjugate_pr;
286 type = gsl_multimin_fdfminimizer_vector_bfgs;
289 type = gsl_multimin_fdfminimizer_vector_bfgs2;
292 type = gsl_multimin_fdfminimizer_steepest_descent;
303 gsl_multimin_fdfminimizer * solver =
304 gsl_multimin_fdfminimizer_alloc(type, dim);
307 gsl_multimin_function_fdf minusLogPosterior;
308 minusLogPosterior.n = dim;
312 minusLogPosterior.params = (
void *)(
this);
321 status = gsl_multimin_fdfminimizer_iterate(solver);
326 std::cerr <<
"Error while GSL does optimisation. "
327 <<
"See below for GSL error type." << std::endl;
328 std::cerr <<
"Gsl error: " << gsl_strerror(status) << std::endl;
333 status = gsl_multimin_test_gradient(solver->gradient, this->getTolerance());
337 gsl_vector* x = gsl_multimin_fdfminimizer_x(solver);
338 std::vector<double> x_min(dim);
339 for(
unsigned int i = 0; i < dim; i++)
340 x_min[i] = gsl_vector_get(x,i);
342 double f = gsl_multimin_fdfminimizer_minimum(solver);
344 gsl_vector* grad = gsl_multimin_fdfminimizer_gradient(solver);
345 double grad_norm = gsl_blas_dnrm2(grad);
347 monitor->
append( x_min, f, grad_norm );
352 for (
unsigned int i = 0; i < dim; i++) {
353 (*m_minimizer)[i] = gsl_vector_get(solver->x, i);
357 gsl_multimin_fdfminimizer_free(solver);
366 gsl_vector* x = gsl_vector_alloc(dim);
367 for (
unsigned int i = 0; i < dim; i++) {
372 const gsl_multimin_fminimizer_type* type = NULL;
377 type = gsl_multimin_fminimizer_nmsimplex;
380 type = gsl_multimin_fminimizer_nmsimplex2;
383 type = gsl_multimin_fminimizer_nmsimplex2rand;
396 gsl_multimin_fminimizer* solver =
397 gsl_multimin_fminimizer_alloc(type, dim);
400 gsl_multimin_function minusLogPosterior;
401 minusLogPosterior.n = dim;
403 minusLogPosterior.params = (
void *)(
this);
406 gsl_vector* step_size = gsl_vector_alloc(dim);
408 for(
unsigned int i = 0; i < dim; i++) {
412 gsl_multimin_fminimizer_set(solver, &minusLogPosterior, x, step_size);
421 status = gsl_multimin_fminimizer_iterate(solver);
426 std::cerr <<
"Error while GSL does optimisation. "
427 <<
"See below for GSL error type." << std::endl;
428 std::cerr <<
"Gsl error: " << gsl_strerror(status) << std::endl;
433 size = gsl_multimin_fminimizer_size(solver);
435 status = gsl_multimin_test_size (size, this->
getTolerance());
439 gsl_vector* x = gsl_multimin_fminimizer_x(solver);
440 std::vector<double> x_min(dim);
441 for(
unsigned int i = 0; i < dim; i++)
442 x_min[i] = gsl_vector_get(x,i);
444 double f = gsl_multimin_fminimizer_minimum(solver);
446 monitor->
append( x_min, f, size );
453 for (
unsigned int i = 0; i < dim; i++) {
454 (*m_minimizer)[i] = gsl_vector_get(solver->x, i);
458 gsl_vector_free(step_size);
459 gsl_multimin_fminimizer_free(solver);
479 if( solver == std::string(
"fletcher_reeves_cg") )
481 else if( solver == std::string(
"polak_ribiere_cg") )
483 else if( solver == std::string(
"bfgs") )
485 else if( solver == std::string(
"bfgs2") )
487 else if( solver == std::string(
"steepest_decent") )
489 else if( solver == std::string(
"nelder_mead") )
491 else if( solver == std::string(
"nelder_mead2") )
493 else if( solver == std::string(
"nelder_mead2_rand") )
499 std::cerr <<
"Error: Invalid GslOptimizer solver name: " << solver << std::endl
500 <<
" Valids choices are: fletcher_reeves_cg" << std::endl
501 <<
" polak_ribiere_cg" << std::endl
502 <<
" bfgs" << std::endl
503 <<
" bfgs2" << std::endl
504 <<
" steepest_decent" << std::endl
505 <<
" nelder_mead" << std::endl
506 <<
" nelder_mead2" << std::endl
507 <<
" nelder_mead2_rand" << std::endl;
const BaseScalarFunction< GslVector, GslMatrix > & objectiveFunction() const
Returns the objective function.
void append(std::vector< double > &x_min, double objective, double norm)
Add current value of minimizer, objective, and error norm.
Object to monitor convergence of optimizers.
A base class for handling optimisation of scalar functions.
void cwSet(double value)
Component-wise sets all values to this with value.
const BaseScalarFunction< GslVector, GslMatrix > & m_objectiveFunction
virtual ~GslOptimizer()
Destructor.
unsigned int sizeLocal() const
Returns the length of this vector.
GslOptimizer(const BaseScalarFunction< GslVector, GslMatrix > &objectiveFunction)
Constructs an object that will maximize a scalar function.
A templated (base) class for handling scalar functions.
void c_evaluate_with_derivative(const gsl_vector *x, void *context, double *f, gsl_vector *derivative)
GslVector * m_initialPoint
GslVector m_fstep_size
For use in gradient-free algorithms.
Class for vector operations using GSL library.
virtual void minimize(OptimizerMonitor *monitor=NULL)
Minimize the objective function, starting at m_initialPoint.
void minimize_no_gradient(unsigned int dim, OptimizerMonitor *monitor)
double getFiniteDifferenceStepSize() const
Returns the step size used in the finite difference formula.
double getTolerance() const
Returns the tolerance used to test for an extremum in the optimizer.
void setInitialPoint(const GslVector &intialPoint)
Set the point at which the optimization starts.
double m_line_tol
Line minimization tolerance in gradient-based algorithms.
double c_evaluate(const gsl_vector *x, void *context)
void set_step_size(const GslVector &step_size)
Sets step size used in gradient-free solvers.
void set_solver_type(SolverType solver)
void c_evaluate_derivative(const gsl_vector *x, void *context, gsl_vector *derivative)
SolverType string_to_enum(std::string &solver)
A base class for handling optimisation of scalar functions.
bool solver_needs_gradient(SolverType solver)
Helper function.
const GslVector & minimizer() const
Return the point that minimizes the objective function.
double m_fdfstep_size
For use in gradient-based algorithms.
void minimize_with_gradient(unsigned int dim, OptimizerMonitor *monitor)
unsigned int getMaxIterations() const
Returns the maximum number of iterations the optimizer will do.