@jacobwilliams and others, nlopt looks like a great optimization tool, but its fortran interface currently is Fortran 77. Could be updated to use iso_c_binding.
2 Likes
Yeah, some of that is Fortran 77 code that was converted to C and then a Fortran 77 wrapper written for it. If a Fortran version of something exists, the right way is to modernize it and use it. No need for C to be involved. Improvements can be made to the modernized Fortran code.
I’m also not a fan of writing Fortran 90 interfaces to a pile of Fortran 77 subroutines, as though they are our ancient legacy that must never be changed. 
But maybe I’ve mentioned that before. 
2 Likes
The native NLopt Fortran bindings are really only the bare minimum, i.e. they contain only the enumerator values, not even explicit interfaces, which makes them rather fragile to use.
As Emanuele @epagone pointed out, I already went through all the work to create interfaces to the actual C API of NLopt and wrap them nicely in a derived type to make the overall user experience more pleasant (announced in NLopt Fortran bindings). I also put some work into making the Fortran bindings compatible with all NLopt versions since 2.0.0, because not all package managers ship the latest version. Feedback and contributions are of course very welcome.
If you want to see an actual application of the Fortran bindings to NLopt, I recently pushed dftd4/dftd4-fit, which glues together a couple of fpm projects to a small command line optimization driver.
3 Likes
Not to undermine Sebastian’s large effort, however, I would like to mention that I also attempted to create a modern Fortran interface back in 2019. I even tried to submit it as pull request ([WIP] A Modern Fortran Interface by ivan-pi · Pull Request #233 · stevengj/nlopt · GitHub) to the original NLopt repository. Preceding my own attempt, GitHub member aitap described how to wrap NLopt in a blog post (Wrapping C opaque pointers in Fortran 2003 in type-strict way).
The factors blocking my PR at the time were CMake integration in the upstream project and also some struggles with memory leaks and interfaces in the object-oriented Fortran wrapper (Actual and dummy procedure argument mismatch - Error 7062 - Intel Communities).
Since some NLopt users already depend upon the “old” (standard-violating…) interface, I think it would be difficult to have a modern interface co-exist in the original NLopt development repository. I also had the feeling the maintainers (and creators) of NLopt aren’t particularly interested in modern Fortran as a contemporary programming language.
Hence, I ultimately closed my PR and switched to using @awvwgk’s wrappers instead of my own (incomplete) version. AFAICT the nlopt-f interface is completely type-safe and agnostic of the fact the underlying library is in C. It only assumes that you have nlopt installed somewhere.
My takeaway from this whole experience is that fragmentation of the community can really hinder progress. Another lesson I learned is the importance of “finishing” a project to a level where others can use it safely.
I would add a note for future Fortran optimization libraries (cc @zaikunzhang), if they could provide a nlopt-f callback compatible interface, it would make switching between optimization libraries much easier for consumers. The current interface expected in nlopt-f is:
function nlopt_func_interface(x, gradient, func_data) result(f)
import :: c_int, c_double, c_ptr
implicit none
real(c_double), intent(in) :: x(:)
real(c_double), intent(inout), optional :: gradient(:)
class(*), intent(in), optional :: func_data
real(c_double) :: f
end function
If you already have an optimization library with different callback conventions for the objective function and parameters, you can employ the adaptor pattern to offer consumers a type-safe interface. The same pattern is in fact used in nlopt-f to interface with the original library in C.
For even wider consumption potentially through NLopt itself (if the maintainers agree upon it), your optimization algorithm (as a Fortran subprogram) would need to export an interface similar to:
nlopt_result cobyla_minimize(unsigned n, nlopt_func f, void *f_data,
unsigned m, nlopt_constraint *fc,
unsigned p, nlopt_constraint *h,
const double *lb, const double *ub, /* bounds */
double *x, /* in: initial guess, out: minimizer */
double *minf,
nlopt_stopping *stop,
const double *dx);
This however means you agree to implement all the possible stopping criteria expected by the NLopt API (see nlopt_stopping or nlopt_result), as well as some further conditions upon the constraints (e.g. NLopt expects the objective function will never be evaluated outside of the bounds) among other things.
6 Likes
Good point, I opened an issue for discussing the ad hoc choice of callback in nlopt-f: Callback function design · Issue #16 · grimme-lab/nlopt-f · GitHub. If we can agree on a community standard for callbacks in optimization libraries, I would be happy to support it in nlopt-f as well.
1 Like
At least in the fpm ecosystem, a way of formalizing such community standards for callbacks would be to offer them in a module of abstract interfaces.
(As a side question, when did Fortran introduce abstract interfaces?)
You see this pattern carried to an extreme in the SciMLBase.jl sub-package from the Julia SciML project. The hierarchy starts with an abstract scientific machine learning problem:
abstract type SciMLProblem end
All types of problems, differential equations, PDE’s, linear systems, nonlinear systems, etc. are all sub-types of this. It is really impressive how they connected all the pieces with a huge hierarchy. For each abstract sub-problem, there will be a method looking something like
function solve(prob::AbstractProblem,args...,kwargs...)
# .
end
Using multiple-dispatch, the method will specialize based upon the problem, algorithms, parameters, and so forth.
Here is a concrete look of how it ends up on the consumer level:
using DifferentialEquations
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
For an optimization problem, you again have the same kind of interface:
using GalacticOptim, Optim
rosenbrock(x,p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
p = [1.0,100.0]
prob = OptimizationProblem(rosenbrock,x0,p)
sol = solve(prob,NelderMead())
The pattern works very well for them. It requires a shift in mindset from the classic Fortran interfaces with a subprogram (algorithm) that expects a callback function. In the Julia pattern you encapsulate both the problem and the algorithm. With some diligence and collaboration, a similar feat could be done in Fortran, building upon the same patterns exploited by Julia programmers.
3 Likes
Thank you @ivanpribec for the kind reminder and informative elaboration. I fully agree that it is important for libraries to provide compatible interfaces so that users can easily switch from one library to another. This is particularly relevant for optimization due to the inevitability of passing (scalar- or vector-valued) functions.
For the moment, my development is strictly focused on the modernization of Powell’s Derivative-free Optimization solvers, which were originally coded in Fortran 77 by M. J. D. Powell FRS (Wikipedia). These solvers are also included in NLopt, an example being the cobyla solver mentioned in your post. Currently, I concentrate more on the code (that does the real computation) instead of the interfaces. The current interfaces follow essentially the design by Powell with modernization.
For example, the interface I adopt for cobyla at this stage is as follows. I will investigate more carefully the interfaces of existing libraries (particularly NLopt-f and stdlib) when I am about to finish the modernization of the code.
Click to see the interface
subroutine cobyla(calcfc, m, x, f, &
& cstrv, constr, &
& f0, constr0, &
& nf, rhobeg, rhoend, ftarget, ctol, maxfun, iprint, &
& xhist, fhist, chist, conhist, maxhist, info)
!--------------------------------------------------------------------------------------------------!
! Among all the arguments, only CALCFC, M, X, and F are obligatory. The others are OPTIONAL and you can
! neglect them unless you are familiar with the algorithm. If you do not specify an optional input,
! it will be assigned the default value detailed below. For instance, we may write
!
! call cobyla(calcfc, m, x, f)
!
! or
!
! call cobyla(calcfc, m, x, f, rhobeg = 0.5D0, rhoend = 1.0D-3, maxfun = 100)
!
! See examples/cobyla_exmp.f90 for a concrete example.
!
! A detailed introduction to the arguments is as follows.
! N.B.: RP and IK are defined in the module CONSTS_MOD. See consts.F90 under the directory name
! "common". By default, RP = kind(0.0D0) and IK = kind(0), with REAL(RP) being the double-precision
! real, and INTEGER(IK) being the default integer. For ADVANCED USERS, RP and IK can be defined by
! setting __REAL_PRECISION__ and __INTEGER_KIND__ in common/ppf.h. Use the default if unsure.
!
! CALCFC
! Input, subroutine.
! CALCFC(X, F, CONSTR) should evaluate the objective function and constraints at the given
! REAL(RP) vector X; it should set the objective function value to the REAL(RP) scalar F and the
! constraint value to the REAL(RP) vector CONSTR. It must be provided by the user, and its
! definition must conform to the following interface:
! !-------------------------------------------------------------------------!
! !subroutine calcfc(x, f, constr)
! !real(RP), intent(in) :: x(:)
! !real(RP), intent(out) :: f
! !real(RP), intent(out) :: constr(:)
! !end subroutine calcfc
! !-------------------------------------------------------------------------!
! Besides, the size of CONSTR must be M, which is the second compulsory argument (see below).
!
! M
! Input, INTEGER(IK) scalar.
! M must be set to the number of constraints, namely the size (length) of CONSTR(X).
!
! X
! Input and output, REAL(RP) vector.
! As an input, X should be an N-dimensional vector that contains the starting point, N being the
! dimension of the problem. As an output, X will be set to an approximate minimizer.
!
! F
! Output, REAL(RP) scalar.
! F will be set to the objective function value of X at exit.
!
! CSTRV
! Output, REAL(RP) scalar.
! CSTRV will be set to the constraint violation of X at exit, i.e., MAXVAL([-CONSTR(X), 0]).
!
! CONSTR
! Output, ALLOCATABLE REAL(RP) vector.
! CONSTR will be set to the constraint value of X at exit.
!
! F0
! Input, REAL(RP) scalar.
! F0, if present, should be set to the objective function value of the starting X.
!
! CONSTR0
! Input, REAL(RP) vector.
! CONSTR0, if present, should be set to the constraint value of the starting X; in addition,
! SIZE(CONSTR0) must be M, or the solver will abort.
!
! NF
! Output, INTEGER(IK) scalar.
! NF will be set to the number of calls of CALCFC at exit.
!
! RHOBEG, RHOEND
! Inputs, REAL(RP) scalars, default: RHOBEG = 1, RHOEND = 10^-6. RHOBEG and RHOEND must be set to
! the initial and final values of a trust-region radius, both being positive and RHOEND <= RHOBEG.
! Typically RHOBEG should be about one tenth of the greatest expected change to a variable, and
! RHOEND should indicate the accuracy that is required in the final values of the variables.
!
! FTARGET
! Input, REAL(RP) scalar, default: -Inf.
! FTARGET is the target function value. The algorithm will terminate when a feasible point with a
! function value <= FTARGET is found.
!
! CTOL
! Input, REAL(RP) scalar, default: machine epsilon.
! CTOL is the tolerance of constraint violation. Any X with MAXVAL(-CONSTR(X)) <= CTOL is
! considered as feasible.
!
! MAXFUN
! Input, INTEGER(IK) scalar, default: MAXFUN_DIM_DFT*N with MAXFUN_DIM_DFT defined in the module
! CONSTS_MOD (see common/consts.F90). MAXFUN is the maximal number of calls of CALCFC.
!
! NPT
! Input, INTEGER(IK) scalar, default: 2N + 1.
! NPT is the number of interpolation conditions for each trust region model. Its value must be in
! the interval [N+2, (N+1)(N+2)/2].
!
! IPRINT
! Input, INTEGER(IK) scalar, default: 0.
! The value of IPRINT should be set to 0, 1, -1, 2, -2, 3, or -3, which controls how much
! information will be printed during the computation:
! 0: there will be no printing;
! 1: a message will be printed to the screen at the return, showing the best vector of variables
! found and its objective function value;
! 2: in addition to 1, each new value of RHO is printed to the screen, with the best vector of
! variables so far and its objective function value; each new value of CPEN is also printed;
! 3: in addition to 2, each function evaluation with its variables will be printed to the screen;
! -1, -2, -3: the same information as 1, 2, 3 will be printed, not to the screen but to a file
! named COBYLA_output.txt; the file will be created if it does not exist; the new output will
! be appended to the end of this file if it already exists. Note that IPRINT = -3 can be costly
! in terms of time and space.
!
! XHIST, FHIST, CHIST, CONHIST, MAXHIST
! XHIST: Output, ALLOCATABLE rank 2 REAL(RP) array;
! FHIST: Output, ALLOCATABLE rank 1 REAL(RP) array;
! CHIST: Output, ALLOCATABLE rank 1 REAL(RP) array;
! CONHIST: Output, ALLOCATABLE rank 2 REAL(RP) array;
! MAXHIST: Input, INTEGER(IK) scalar, default: MAXFUN
! XHIST, if present, will output the history of iterates; FHIST, if present, will output the
! history function values; CHIST, if present, will output the history of constraint violations;
! CONHIST, if present, will output the history of constraint values; MAXHIST should be a
! nonnegative integer, and XHIST/FHIST/CHIST/CONHIST will output only the history of the last
! MAXHIST iterations. Therefore, MAXHIST= 0 means XHIST/FHIST/CONHIST/CHIST will output nothing,
! while setting MAXHIST = MAXFUN requests XHIST/FHIST/CHIST/CONHIST to output all the history.
! If XHIST is present, its size at exit will be (N, min(NF, MAXHIST)); if FHIST is present, its
! size at exit will be min(NF, MAXHIST); if CHIST is present, its size at exit will be
! min(NF, MAXHIST); if CONHIST is present, its size at exit will be (M, min(NF, MAXHIST)).
!
! Important Notice:
! Setting MAXHIST to a large value can be costly in terms of memory for large problems.
! For instance, if N = 1000 and MAXHIST = 100, 000, XHIST will take up to 1 GB if we use double
! precision. MAXHIST will be reset to a smaller value if the memory needed exceeds MAXMEMORY
! defined in CONSTS_MOD (see consts.F90 under the directory named "common"; default: 2GB).
! Use *HIST with caution!!! (N.B.: the algorithm is NOT designed for large problems).
!
! MAXFILT
! Input, INTEGER(IK) scalar.
! MAXFILT is a nonnegative integer indicating the maximal length of the "filter" used for
! selecting the returned solution; default: 2000 (a value lower than 200 is not recommended)
!
! INFO
! Output, INTEGER(IK) scalar.
! INFO is the exit flag. It will be set to one of the following values defined in the module
! INFO_MOD (see common/info.F90):
! SMALL_TR_RADIUS: the lower bound for the trust region radius is reached;
! FTARGET_ACHIEVED: the target function value is reached;
! MAXFUN_REACHED: the objective function has been evaluated MAXFUN times;
! MAXTR_REACHED: the trust region iteration has been performed MAXTR times,
! the value of MAXTR being 4*MAXFUN, which is UNLIKELY to reach;
! NAN_INF_X: NaN or Inf occurs in x;
! NAN_MODEL: NaN occurs in the model;
! DAMAGING_ROUNDING: rounding errors are becoming damaging.
! !--------------------------------------------------------------------------!
! The following cases should NEVER occur unless there is a bug, because the
! code will try to continue in the corresponding scenarios.
! NAN_INF_F: the objective function returns NaN or +Inf
! TRSUBP_FAILED: a trust region step failed to reduce the quadratic model
! !--------------------------------------------------------------------------!
!--------------------------------------------------------------------------------------------------!
implicit none
! Compulsory arguments
procedure(OBJCON) :: calcfc
real(RP), intent(inout) :: x(:)
real(RP), intent(out) :: f
integer(IK), intent(in) :: m
! Optional inputs
integer(IK), intent(in), optional :: iprint
integer(IK), intent(in), optional :: maxfilt
integer(IK), intent(in), optional :: maxfun
integer(IK), intent(in), optional :: maxhist
real(RP), intent(in), optional :: constr0(:)
real(RP), intent(in), optional :: ctol
real(RP), intent(in), optional :: f0
real(RP), intent(in), optional :: ftarget
real(RP), intent(in), optional :: rhobeg
real(RP), intent(in), optional :: rhoend
! Optional outputs
integer(IK), intent(out), optional :: info
integer(IK), intent(out), optional :: nf
real(RP), intent(out), allocatable, optional :: chist(:)
real(RP), intent(out), allocatable, optional :: conhist(:, :)
real(RP), intent(out), allocatable, optional :: constr(:)
real(RP), intent(out), allocatable, optional :: fhist(:)
real(RP), intent(out), allocatable, optional :: xhist(:, :)
real(RP), intent(out), optional :: cstrv
end subroutine cobyla
Regarding this interface, I do have a comment to make, and it may be relevant for NLopt @awvwgk as well. In Powell’s design, the objective function and constraint function is evaluated by a single subroutine calcfc. This may seem unnatural or even inconvenient at the first glance. I have never asked Powell why he chose this design, but I believe he had a reason to do so. My speculation is as follows. Powell designed his solvers for problems that cannot provide derivatives, aka derivative-free optimization problems. In such problems, the objective and constraints are commonly evaluated by sophisticated and expensive computation/simulation, often black boxes. It is not abnormal that the objective and constraints are evaluated simultaneously (e.g., during the same simulation) and cannot be decoupled trivially. Therefore, it is reasonable to have a single subroutine calcfc that calculates both of them. Surely, we may decouple calcfc into two (e.g., calfun for the objective and calcon for the constraints) by an additional layer of wrapping, but we must avoid re-evaluation of calcfc at the same point because such evaluations are expensive.
In the future, I will implement an interface that allows passing the objective and constraints by not only calcfc , but also [calfun, calcon] or [calcfc, calcon].
Thanks.
I’ve been using the NLopt port of cobyla in my own work.
The design where both the objective function and constraints are evaluated together does make sense also for efficiency. Parts of the computations may be reused for both purposes. I’m not familiar enough with optimization to say, whether the constraints and objective function must always be computed simultaneously.
In a similar manner, the NLopt interface already combines the objective function and its gradient. If the algorithm only wants to update the objective function, but reuse the gradient from the previous iteration, it can do so easily, by passing a null pointer for the gradient callback. This does require the user, to check for the null pointer in his callback:
if (present(gradient)) then
gradient(1) = 3.0_wp * a * (a*x(1) + b)**2
gradient(2) = -1.0_wp
endif
f = (a*x(1) + b)**3 - x(2)
In this example a good compiler can reuse the temporary variable a*x(1) + b. This might not be the case, if the user has to provide two callback functions. For a derivative-free optimizer, you just skip the gradient part in the callback completely.
2 Likes
Thank you @ivanpribec for the explanation. It is very helpful to hear how practitioners use and regard these algorithms.
However, there are some significant misconceptions here. COBYLA does not use derivatives (hence the name “derivative-free optimization”). That’s the whole point of the algorithm.
If you have (approximate) (sub-)derivatives available, it is wrong to use derivative-free algorithms. Otherwise, it is like fighting a tiger by bare hands while you have machine guns with you. As a researcher in this particular area, this is always the first thing I tell to those who want to try derivative-free algorithms. For more information about this (quite special) category of algorithms, I would like to refer to the Introduction section of the PDFO package and a talk I gave some time ago.
Powell’s COBYLA algorithm (I mean, the algorithm in the mathematical sense) uses the objective and constraint values always simultaneously. So there is no efficiency problem to evaluate them always simultaneously.
What you described seems to be a much-adapted version of COBYLA, but again, it is hardly a good idea to use it when you have any kind of (approximate) (sub-)derivatives available. There are so many much more efficient algorithms and packages.
Thank you.
1 Like
I don’t think there is any misconception. My reply was only concerning the software interface. The example code with the objective function and gradient was just to highlight the mechanism by which gradient information is “requested” from the user. Putting the gradient- and function in the same callback is good because it means the compiler has a chance to optimize sub-expression reuse.
NLopt provides several algorithms for both constrained/unconstrained optimization, and derivative-free/derivative-available. In all cases, the low-level C library expects a function with the prototype:
double f(unsigned n, const double* x, double* grad, void* f_data);
Internally in NLopt, derivative-free algorithms will simply pass NULL for the grad array, and reuse the same interface. Users which do not know the gradient, will simply neglect the grad dummy variable entirely.
1 Like
I fully agree to decouple the evaluation of objective and constraints as much as possible, provided that we are talking about algorithms aiming at problems that can naturally evaluate the objective&constraints separately and can potentially provide derivatives. Things should not be coupled when they can be naturally decoupled.
Indeed, the current version of PDFO asks for separate evaluations of objective&constraints. We thought it was a convenient design. However, it turns out that Powell’s design is better for many problems that really need Powell’s algorithms, because objective&constraints are often impossible to be evaluated separately in those problems. That is why in future versions we plan to support more flexible ways of passing the objective&constraints.
1 Like
This weekend I’ve been playing around with a Fortran wrapper of CppNumericalSolvers, a C++ header-only library of numerical optimization methods with an interface striving to be close to the MATLAB’s fminsearch.
I haven’t published the code yet, but wanted to give an example of what the Fortran interface currently looks like. The following callbacks define the problem:
module rosenbrock_problem
use cppoptlib, only: wp => cppoptlib_wp
implicit none
public
contains
!> Objective F(x)
function rosenbrock(x,params) result(y)
real(wp), intent(in) :: x(:)
class(*), intent(in), optional :: params ! not used
real(wp) :: y, t1, t2
t1 = 1 - x(1)
t2 = x(2) - x(1)**2
y = t1**2 + 100*t2*t2
end function
!> Objective gradient, dF(x)/dx
subroutine rosenbrock_grad(x,grad,params)
real(wp), intent(in) :: x(:)
real(wp), intent(out) :: grad(:)
class(*), intent(in), optional :: params ! not used
grad(1) = -2*(1-x(1))+200*(x(2)-x(1)**2)*(-2*x(1))
grad(2) = 200*(x(2)-x(1)**2)
end subroutine
end module
The driver turned out quite pleasant:
program minimize
use cppoptlib, only: wp => cppoptlib_wp, problem_type, solver
use rosenbrock_problem
implicit none
type(problem_type) :: problem
real(wp) :: x(2), res
call problem%set(rosenbrock,rosenbrock_grad) ! or
! problem = problem_type(rosenbrock,rosenbrock_grad)
x = [-1._wp, 2._wp] ! initial guess
print *, "initial guess:", x
call problem%solve(res,x,solver=solver%bfgs)
print *, "objective: ", res
print *, "solution: ", x
end program
In contrast to NLopt which stores a pointer to an opaque C struct that stores the “state” of the numerical solver, in the present library the problem_type only stores pointers to the callback functions. Any C++ structures only live for the duration of the call to type-bound solve method.
It is also possible to use internal procedures instead of module procedures as the callback arguments:
program simple_problem
use cppoptlib, only: wp => cppoptlib_wp, problem_type
implicit none
real(wp) :: x(2), res
type(problem_type) :: problem
x = [-1, 2]
problem = problem_type(objective,gradient)
call problem%solve(res,x)
write(*,'(A,F0.6)') "f in argmin ", res
write(*,'(A)') "solver status "
write(*,'(A)') "solver criteria "
contains
real(wp) function objective(x,params)
real(wp), intent(in) :: x(:)
class(*), intent(in), optional :: params ! not used
objective = 5*x(1)**2 + 100*x(1)**2 + 5
end function
subroutine gradient(x,grad,params)
real(wp), intent(in) :: x(:)
real(wp), intent(out) :: grad(:)
class(*), intent(in), optional :: params ! not used
grad(1) = 2*5*x(1)
grad(2) = 2*100*x(2)
end subroutine
end program
I am currently investigating how to pass custom stopping criteria and printing callbacks. So far I quite like the interface. Having a “stateless” derived type, which only serves as a collection of callbacks seems to be a useful architectural pattern.
Say I wanted to solve the same parameterized problem multiple times in parallel (an embarrassingly parallel problem), presumably, I could achieve this easily with:
program parallel_minimization
use cppoptlib, only: wp => cppoptlib_wp, problem_type
implicit none
type :: problem_params
real(wp) :: a, b, c
end type
integer, parameter :: nproblems = 100
type(problem_type) :: problem
type(problem_params), allocatable :: params(:)
real(wp), allocatable :: obj(:), x(:,:)
integer :: i
problem = problem_type(objective,gradient)
allocate(params(nproblems),obj(nproblems),x(2,nproblems))
! [insert] ... initialize params from file or namelist
! we use a dynamic schedule to parallelize the loop because each solve might
! take a different number of function evaluations. we can fine-tune this
! with the chunk-size. to reduce over-head from scheduling, a guided schedule
! could also be used.
!$omp parallel do schedule(dynamic) default(private) shared(problem)
do i = 1, nproblems
call problem%solve(res(i),x(:,i),params(i))
end do
!$omp end parallel do
contains
! [insert] ... callbacks
end program
Does anyone know if sharing procedure pointers within a parallel section is a “legal” use? Are there any problems which could pop up? (Assume that the callbacks and type-bound solve routine reference no global data, i.e. they are thread-safe.)
Edit: a second embarassingly parallel usage case would be to start the solver from multiple initial guesses (maybe according to a Latin hypercube sampling). Picking the best solution (presumably the global minimum) is as simple as:
xbest = x(:,minloc(res,dim=1))
4 Likes
As long as the procedure is thread-safe accessing it in parallel shouldn’t be an issue.
How can you customize the problem to select for example a whether you have a minimization or maximization problem with this approach?
1 Like
- The C++ library only supports minimization. If you need maximization I presume you can multiply the objective with -1.
- Since the Fortran interface uses an “adaptor” between the Fortran callbacks and the C++ objective functor, you could replace
solvewith two methods,minimizeandmaximize, where the maximize adaptor would do the multiplication for you.
1 Like
I would avoid making the constructor type-bound, better create a new_problem procedure which creates a new instance of the type rather than the class, either as intent(out) argument or as return value from a function.
How do you decide which argument is part of the problem and which can be provided with solve. The initial guess is obviously associated with the solve procedure, the solver argument however seems misplaced here. Instead shouldn’t there be a BFGS solver which can consume a problem and return a result?
Great questions. Similar questions have been nagging me for a while now, should we encapsulate the problem, the solution, the solver, and how do we make them composable? There are similar dichotomies in writing object-oriented libraries for linear algebra, and ODE solvers. Do you make the matrix the derived type, or the factorization? Unfortunately, I haven’t yet reached a conclusion which is better.
Julia, which has arguably reached the highest degree of composability, follows a pattern where they have a solve method that gets passed a problem and a solver, something like:
# ODE Problem
using DifferentialEquations
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
# Optimization
using Optim
rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
result = optimize(rosenbrock, zeros(2), BFGS())
As you can see, for any abstract mathematical problem, the pattern they adopt is similar to:
result = action(problem, args, solver)
Since the language is dynamic this works great because the result type can carry all the output state. In Fortran I have the feeling it’s better to constrain output to intrisic types and arrays, perhaps merging solver statistics into a derived type to reduce the need for long lists of dummy variables.
If I understand correctly, you’re suggesting an interface like:
type(solver_type) :: solver
type(problem_type) :: problem
solver = bfgs() ! optionally pass solver settings like stop criteria
problem = problem_type(objective,gradient)
call solve(problem,initial_guess,solver) ! <-- here lies the dichotomy
! the solve method could be actually bound to the problem as in
call problem%solve(initial_guess,solver)
! or, consumed by the solver
call solver%solve(problem,initial_guess)
The original C++ library uses the second pattern, i.e.
cppoptlib::LbfgsSolver<FortranObjective<double>> solver;
solver.minimize(f, x);
we initialize a solver and use it to minimize the objective/problem (defined beforehand).
I guess you could call these two oppoising viewpoints solver-centric versus problem-centric. Julia libraries adopt the problem-centric viewpoint. The parameters enter this whole picture as a third design point. I consciously choose to avoid encapsulating the parameters (both problem- or solver-related) so that the library could be used easily in a parallel context.
Here is an overview of the current wrapper design.
The problem type only contains the callbacks:
!> A container for the callback objects
type :: problem_type
procedure(objective_callback), pointer, nopass :: objective => null()
procedure(gradient_callback), pointer, nopass :: gradient => null()
contains
procedure :: set => set_problem
procedure :: eval => eval_objective
procedure :: eval_grad => eval_gradient
procedure, pass(problem) :: solve => solve_problem
end type
The solvers are just integer “enumerators”
type :: solver_enum
integer(c_int) :: gd = 0 !> Gradient-descent
integer(c_int) :: newton = 1 !> Newton
integer(c_int) :: bfgs = 2 !> BFGS
integer(c_int) :: lbfgs = 3 !> L-BFGS
integer(c_int) :: lbfgsb = 4 !> L-BFGS-B
integer(c_int) :: cg = 5 !> Conjugate gradient
integer(c_int) :: cmaes = 6 !> CMAES
integer(c_int) :: neldermead = 7 !> Nelder-Mead
end type
type(solver_enum), parameter :: solver = solver_enum()
integer(c_int), parameter :: default_solver = solver%bfgs
The solve method prepares the adaptors and calls the C++ routine (hiding underneath a C interface)
subroutine solve_problem(problem,result,x,solver,params,istat)
class(problem_type), intent(in), target :: problem
real(wp), intent(out) :: result
real(wp), intent(inout), contiguous :: x(:)
class(*), intent(in), optional, target :: params
integer, intent(in), optional :: solver
integer, intent(out), optional :: istat
! TODO: investigate mimicking typing for solver argument
interface
subroutine cppoptlib_solve_d(solver,obj_fcn,grad_fcn, &
n,x,istat,udata) bind(c,name="cppoptlib_solve_d")
use, intrinsic :: iso_c_binding, only: c_double, &
c_int, c_ptr, c_funptr
integer(c_int), intent(in), value :: solver
type(c_funptr), value :: obj_fcn
type(c_funptr), value :: grad_fcn
integer(c_int), intent(in), value :: n
real(c_double), intent(inout) :: x(*)
integer(c_int), intent(out) :: istat
type(c_ptr), value :: udata
end subroutine
end interface
type(callback_adaptor), target :: cb
integer(c_int) :: istat_, solver_
solver_ = default_solver
if (present(solver)) then
solver_ = solver
end if
! The callback adaptor is essentially a container that stores
! pointers to the callback functions and parameters.
! It gets passed to the C routine as an void* pointer and later
! "cast" back in the actual callback routines on the Fortran side.
cb%obj => problem%objective
cb%grad => problem%gradient
if (present(params)) then
cb%params => params
end if
! TODO: bounded problems
! if (problem%bounded) ...
call cppoptlib_solve_d( &
solver=solver_, &
obj_fcn=c_funloc(cppoptlib_obj_callback_d), &
grad_fcn=c_funloc(cppoptlib_grad_callback_d), &
n=size(x), &
x=x, &
istat=istat_, &
udata=c_loc(cb))
! Handle status (error) flag
if (present(istat)) then
istat = istat_
else
if (istat_ < 0) then
print *, "cppoptlib: something went wrong"
return
end if
end if
! Evaluate objective for user convenience
result = problem%objective(x,params)
end subroutine
The C++ side first involves a class that wraps the Fortran callback adaptors:
template<typename T>
class FortranObjective : public cppoptlib::Problem<T> {
public:
using typename cppoptlib::Problem<T>::TVector;
using typename cppoptlib::Problem<T>::THessian;
T value(const TVector &x) {
// call Fortran objective callback
return obj_(x.size(),x.data(),params_);
}
void gradient(const TVector &x, TVector &grad) {
// call Fortran gradient callback
grad_(x.size(),x.data(),grad.data(),params_);
}
// Constructor
FortranObjective(obj_callback_d obj_fcn, grad_callback_d grad_fcn, void* params) :
obj_(obj_fcn), grad_(grad_fcn), params_(params) {}
private:
// TODO: handle callback type generically?
obj_callback_d obj_{nullptr};
grad_callback_d grad_{nullptr};
void* params_{nullptr};
};
The actual routine doing the work is effectively just a big switch to the various templated solver class instances, see below. I borrowed the pattern from the MATLAB wrapper. Compared with a native C++ driver program, the downside is that this leads to bloat of the binary executable as all solvers must be included. Perhaps link-time optimization could crop some of the dead branches, depending whether it can propagate the solver flag through the mixed-language interface?
void cppoptlib_solve_d(
int selected_solver,
obj_callback_d obj_fcn,
grad_callback_d grad_fcn,
int n,
double* x_,
int* ierr,
void* udata) {
// Create view of Fortran vector buffer
Eigen::Map<Eigen::VectorXd> x_view{x_, n};
// The minimize function cannot operate on a map, so we make a copy
auto x = x_view.eval();
FortranObjective<double> f{obj_fcn,grad_fcn,udata};
// Assume success
(*ierr) = 0;
enum solver_type {
GRADIENTDESCENT,
NEWTON,
BFGS,
LBFGS,
LBFGSB,
CONJUGATEDGRADIENTDESCENT,
CMAES,
NELDERMEAD};
// solve
// ----------------------------------------------------------
switch (selected_solver) {
case LBFGS: {
cppoptlib::LbfgsSolver<FortranObjective<double>> solver;
solver.minimize(f, x);
}
break;
case BFGS: {
cppoptlib::BfgsSolver<FortranObjective<double>> solver;
solver.minimize(f, x);
}
break;
case GRADIENTDESCENT: {
cppoptlib::GradientDescentSolver<FortranObjective<double>> solver;
solver.minimize(f, x);
}
break;
case CONJUGATEDGRADIENTDESCENT: {
cppoptlib::ConjugatedGradientDescentSolver<FortranObjective<double>> solver;
solver.minimize(f, x);
}
break;
case NEWTON: {
cppoptlib::NewtonDescentSolver<FortranObjective<double>> solver;
solver.minimize(f, x);
}
break;
case NELDERMEAD: {
cppoptlib::NelderMeadSolver<FortranObjective<double>> solver;
solver.minimize(f, x);
}
break;
default:
// TODO: Illegal input, print error
*ierr = -1;
break;
}
// over-write initial guess with result
x_view = x;
}