27 #include <queso/Defines.h>
28 #include <queso/GslVector.h>
29 #include <queso/VectorSpace.h>
30 #include <queso/ScalarFunction.h>
31 #include <queso/GslOptimizer.h>
32 #include <queso/OptimizerMonitor.h>
34 #include <gsl/gsl_multimin.h>
35 #include <gsl/gsl_blas.h>
43 double c_evaluate(
const gsl_vector * x,
void * context) {
51 for (
unsigned int i = 0; i < state.sizeLocal(); i++) {
52 state[i] = gsl_vector_get(x, i);
70 gsl_vector * derivative) {
79 for (
unsigned int i = 0; i < state.sizeLocal(); i++) {
80 state[i] = gsl_vector_get(x, i);
90 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
91 gsl_vector_set(derivative, i, GSL_NAN);
106 bool userComputedDerivative =
true;
107 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
110 if (gsl_isnan(deriv[i])) {
111 userComputedDerivative =
false;
116 if (userComputedDerivative) {
117 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
118 gsl_vector_set(derivative, i, -deriv[i]);
126 for (
unsigned int i = 0; i < deriv.sizeLocal(); i++) {
127 double tempState = state[i];
136 state[i] = tempState;
139 if (!gsl_isnan(fx) && !gsl_isnan(fxph)) {
140 gsl_vector_set(derivative, i, (fxph - fx) / h);
143 gsl_vector_set(derivative, i, GSL_NAN);
152 double * f, gsl_vector * derivative) {
162 m_objectiveFunction(objectiveFunction),
163 m_initialPoint(new
GslVector(objectiveFunction.domainSet().
164 vectorSpace().zeroVector())),
165 m_minimizer(new
GslVector(this->m_objectiveFunction.domainSet().
166 vectorSpace().zeroVector())),
167 m_solver_type(BFGS2),
168 m_fstep_size(this->m_objectiveFunction.domainSet().vectorSpace().zeroVector()),
191 std::cerr <<
"Minimization was given initial point outside of domain"
198 zeroVector().sizeLocal();
221 for (
unsigned int i = 0; i < initialPoint.
sizeLocal(); i++) {
239 bool gradient_needed =
false;
249 gradient_needed =
true;
265 return gradient_needed;
271 gsl_vector * x = gsl_vector_alloc(dim);
272 for (
unsigned int i = 0; i <
dim; i++) {
277 const gsl_multimin_fdfminimizer_type* type = NULL;
282 type = gsl_multimin_fdfminimizer_conjugate_fr;
285 type = gsl_multimin_fdfminimizer_conjugate_pr;
288 type = gsl_multimin_fdfminimizer_vector_bfgs;
291 type = gsl_multimin_fdfminimizer_vector_bfgs2;
294 type = gsl_multimin_fdfminimizer_steepest_descent;
305 gsl_multimin_fdfminimizer * solver =
306 gsl_multimin_fdfminimizer_alloc(type, dim);
309 gsl_multimin_function_fdf minusLogPosterior;
310 minusLogPosterior.n =
dim;
314 minusLogPosterior.params = (
void *)(
this);
323 status = gsl_multimin_fdfminimizer_iterate(solver);
328 std::cerr <<
"Error while GSL does optimisation. "
329 <<
"See below for GSL error type." << std::endl;
330 std::cerr <<
"Gsl error: " << gsl_strerror(status) << std::endl;
335 status = gsl_multimin_test_gradient(solver->gradient, this->getTolerance());
339 gsl_vector* x = gsl_multimin_fdfminimizer_x(solver);
340 std::vector<double> x_min(dim);
341 for(
unsigned int i = 0; i <
dim; i++)
342 x_min[i] = gsl_vector_get(x,i);
344 double f = gsl_multimin_fdfminimizer_minimum(solver);
346 gsl_vector* grad = gsl_multimin_fdfminimizer_gradient(solver);
347 double grad_norm = gsl_blas_dnrm2(grad);
349 monitor->
append( x_min, f, grad_norm );
354 for (
unsigned int i = 0; i <
dim; i++) {
355 (*m_minimizer)[i] = gsl_vector_get(solver->x, i);
359 gsl_multimin_fdfminimizer_free(solver);
368 gsl_vector* x = gsl_vector_alloc(dim);
369 for (
unsigned int i = 0; i <
dim; i++) {
374 const gsl_multimin_fminimizer_type* type = NULL;
379 type = gsl_multimin_fminimizer_nmsimplex;
382 type = gsl_multimin_fminimizer_nmsimplex2;
385 type = gsl_multimin_fminimizer_nmsimplex2rand;
398 gsl_multimin_fminimizer* solver =
399 gsl_multimin_fminimizer_alloc(type, dim);
402 gsl_multimin_function minusLogPosterior;
403 minusLogPosterior.n =
dim;
405 minusLogPosterior.params = (
void *)(
this);
408 gsl_vector* step_size = gsl_vector_alloc(dim);
410 for(
unsigned int i = 0; i <
dim; i++) {
414 gsl_multimin_fminimizer_set(solver, &minusLogPosterior, x, step_size);
423 status = gsl_multimin_fminimizer_iterate(solver);
428 std::cerr <<
"Error while GSL does optimisation. "
429 <<
"See below for GSL error type." << std::endl;
430 std::cerr <<
"Gsl error: " << gsl_strerror(status) << std::endl;
435 size = gsl_multimin_fminimizer_size(solver);
437 status = gsl_multimin_test_size (size, this->
getTolerance());
441 gsl_vector* x = gsl_multimin_fminimizer_x(solver);
442 std::vector<double> x_min(dim);
443 for(
unsigned int i = 0; i <
dim; i++)
444 x_min[i] = gsl_vector_get(x,i);
446 double f = gsl_multimin_fminimizer_minimum(solver);
448 monitor->
append( x_min, f, size );
455 for (
unsigned int i = 0; i <
dim; i++) {
456 (*m_minimizer)[i] = gsl_vector_get(solver->x, i);
460 gsl_vector_free(step_size);
461 gsl_multimin_fminimizer_free(solver);
481 if( solver == std::string(
"fletcher_reeves_cg") )
483 else if( solver == std::string(
"polak_ribiere_cg") )
485 else if( solver == std::string(
"bfgs") )
487 else if( solver == std::string(
"bfgs2") )
489 else if( solver == std::string(
"steepest_decent") )
491 else if( solver == std::string(
"nelder_mead") )
493 else if( solver == std::string(
"nelder_mead2") )
495 else if( solver == std::string(
"nelder_mead2_rand") )
501 std::cerr <<
"Error: Invalid GslOptimizer solver name: " << solver << std::endl
502 <<
" Valids choices are: fletcher_reeves_cg" << std::endl
503 <<
" polak_ribiere_cg" << std::endl
504 <<
" bfgs" << std::endl
505 <<
" bfgs2" << std::endl
506 <<
" steepest_decent" << std::endl
507 <<
" nelder_mead" << std::endl
508 <<
" nelder_mead2" << std::endl
509 <<
" nelder_mead2_rand" << std::endl;
GslOptimizer(const BaseScalarFunction< GslVector, GslMatrix > &objectiveFunction)
Constructs an object that will maximize a scalar function.
void append(std::vector< double > &x_min, double objective, double norm)
Add current value of minimizer, objective, and error norm.
double m_line_tol
Line minimization tolerance in gradient-based algorithms.
void set_step_size(const GslVector &step_size)
Sets step size used in gradient-free solvers.
virtual ~GslOptimizer()
Destructor.
A base class for handling optimisation of scalar functions.
void setInitialPoint(const GslVector &intialPoint)
Set the point at which the optimization starts.
const BaseScalarFunction< GslVector, GslMatrix > & m_objectiveFunction
const BaseScalarFunction< GslVector, GslMatrix > & objectiveFunction() const
Returns the objective function.
A templated (base) class for handling scalar functions.
virtual void minimize(OptimizerMonitor *monitor=NULL)
Minimize the objective function, starting at m_initialPoint.
bool solver_needs_gradient(SolverType solver)
Helper function.
void minimize_with_gradient(unsigned int dim, OptimizerMonitor *monitor)
void cwSet(double value)
Component-wise sets all values to this with value.
A base class for handling optimisation of scalar functions.
double m_fdfstep_size
For use in gradient-based algorithms.
GslVector * m_initialPoint
unsigned int getMaxIterations() const
Returns the maximum number of iterations the optimizer will do.
double getTolerance() const
Returns the tolerance used to test for an extremum in the optimizer.
void set_solver_type(SolverType solver)
unsigned int sizeLocal() const
Returns the length of this vector.
void minimize_no_gradient(unsigned int dim, OptimizerMonitor *monitor)
const GslVector & minimizer() const
Return the point that minimizes the objective function.
void c_evaluate_derivative(const gsl_vector *x, void *context, gsl_vector *derivative)
void c_evaluate_with_derivative(const gsl_vector *x, void *context, double *f, gsl_vector *derivative)
double getFiniteDifferenceStepSize() const
Returns the step size used in the finite difference formula.
Class for vector operations using GSL library.
GslVector m_fstep_size
For use in gradient-free algorithms.
Object to monitor convergence of optimizers.
double c_evaluate(const gsl_vector *x, void *context)
SolverType string_to_enum(std::string &solver)