libprima / prima

PRIMA is a package for solving general nonlinear optimization problems without using derivatives. It provides the reference implementation for Powell's derivative-free optimization methods, i.e., COBYLA, UOBYQA, NEWUOA, BOBYQA, and LINCOA. PRIMA means Reference Implementation for Powell's methods with Modernization and Amelioration, P for Powell.
http://libprima.net
BSD 3-Clause "New" or "Revised" License
291 stars 35 forks source link

Add callbacks #113

Closed nbelakovski closed 6 months ago

nbelakovski commented 7 months ago

I'm re-opening https://github.com/libprima/prima/pull/69, but since the original branch was closed I made a new branch on my own repo and copied the changes there. All credit of course goes to @jschueller as this is his work and not mine, I am merely shepherding it!

This is necessary for the upcoming work to integrate with SciPy, since the implementation of COBYLA that they have does provide a callback so obviously we don't want any regression in functionality upon switching to the PRIMA implementation of COBYLA.

nbelakovski commented 7 months ago

Also, @jschueller you're more than welcome to either make suggestions here or if you prefer to take over the work and make your own PR in favor of this one, that's fine too.

jschueller commented 7 months ago

great, I saw you rebased

while this PR was great to report progress,

my concern after that is that there is still no way to handle errors in the objective function callback itself (prima_obj/objcon) in order to handle a python exception for example, so I wondered if the objective function callback must return a "terminate" boolean too

what do you think ?

nbelakovski commented 7 months ago

Welcome back @jschueller !

It sounds like you're thinking of the following scenario (using Python as an example but it could apply to other languages as well):

def my_cost_function(x):
  try:
    ...some computation...
  catch err:
     ... what do we do here?...

And further it sounds like you're saying that if we add a terminate output to the objective/cost function, then in the catch block we could set it, and then once we're back in the Fortran code we can gracefully terminate the routine instead of crashing.

Seems like a reasonable addition, I wonder what @zaikunzhang thinks?

nbelakovski commented 7 months ago

Pushed updates now that the algo enum branch is merged.

nbelakovski commented 7 months ago

-"Fixed" the spelling issues -Addressed @zaikunzhang 's comments from the previous PR that the callback should report the "current best" values as opposed to simply the values at the current iteration (if the user is very interested in those they can access them via the objective function) -Added the callback at the end of the initialization -Added a no-op callback so that we could properly pass the callback to the initialization. This flows better with the result of the optional arguments in the [algorith]a.f90 files

nbelakovski commented 7 months ago

I spent most of today investing some of the build failures from CI. It's a long story, but the short version is that in my experiments with flang (the AMD optimizing C/C++ compiler, AOCC) I discovered that the only way to support this compiler is to make some rather ugly code changes. First of all, bobyqa.f90 and others cannot have the following code:

if (present(callback_fcn)) then
    callback_fcn_loc => callback_fcn
else
    callback_fcn_loc= => NO_OP_CALLBACK
end if

The problem is that while we can write call callback_fcn(...), we cannot copy the procedure pointer to another procedure pointer. Something gets lost in the copy and so call callback_fcn_loc(...) results in a segfault.

There is a workaround for this but it's not very pretty:

if (present(callback_fcn)) then
    call bobyqb(...., callback_fcn, ...)
else
    call bobyqb(..., NO_OP_CALLBACK, ....)
end if

The other problem is that we cannot call it from both bobyqa/initialize.f90 and then also from bobyqa/bobyqb.f90. I'm not sure why, something about too many levels. That one is solvable by simply not passing it to bobyqa/initialize.f90.

I don't know yet if these things will solve the issues with the other compilers as well (nvfortran and aflang are also having issues).

@zaikunzhang how would you feel about deciding that only gfortran and ifort will support interop with C? We can modify the CMake files so that they do not try to compile and run the C examples with flang/nvfortran/aflang, only with gfortran and ifort, and then we can keep the code a little bit cleaner. Does this make sense to you and does it sound sensible, or do you have thoughts on the topic you'd like to share?

nbelakovski commented 7 months ago

my concern after that is that there is still no way to handle errors in the objective function callback itself (prima_obj/objcon) in order to handle a python exception for example, so I wondered if the objective function callback must return a "terminate" boolean too

what do you think ?

Discussed this with @zaikunzhang offline and we think that it doesn't make much sense for the objective function to offer a terminate functionality. If it’s truly necessary, the terminate functionality provided by the callback can be used. The objective function can be called within the callback, and the value of terminate can be set appropriately.

I’ve pushed an update. Via discussion with @zaikunzhang we decided to go with the proposed code above that works with the other compilers. I've also modified the manner in which @jschueller implemented both the callback and the call to the C objective function.

The previous setup looked like this (I use bobyqa as an example but this applied equally to all algorithms):

  1. bobyqc_c.f90 captures the pointer to the callback function in a closure.
  2. The closure is passed to bobyqa.f90 and then to bobyqb.f90
  3. bobyqb call the closure, which then checks if the pointer is valid, and if so it calls evalcallback and passes the pointer and the other arguments
  4. evalcallback then casts the pointer to a callable fortran procedure, casts the arguments to their C equivalents, and calls the function.

The new setup is as follows:

  1. bobyqa_c.f90 check if the callback function pointer is valid, and if so casts it to a callable fortran procedure. This procedure is captured in a closure.
  2. The closure is passed to bobyqa.f90 and then to bobyqb.f90 (same as before)
  3. bobyqb call the closure, which casts the arguments to their C equivalents and calls the functions.

So two things changed:

  1. The casting of the function pointer is done during the initial setup, as opposed to every time the function is called.
  2. evalcallback and the closure were combined, so there's fewer functions to follow (this was essentially the main motivation, as the code was difficult to follow as it bounced around between many functions)

This was then also applied to evalcobj.

It looks like there are still some failing builds, but the errors appear to be happening on main as well (at least the ones that I tested), so they don't appear to be related to the code in this PR.

jschueller commented 7 months ago

so how would you handle a function that cannot be evaluated (an exception for example) on a new point proposed by the algorithm ?

nbelakovski commented 7 months ago

so how would you handle a function that cannot be evaluated (an exception for example) on a new point proposed by the algorithm ?

@zaikunzhang can speak in more detail on this, but from our conversations I understand that a common tactic in these situations is to return NaN. COBYLA and other algorithms are set up in such a way as to be able to recover from NaN by essentially going back and trying different points.

jschueller commented 7 months ago

ok, fine then

zaikunzhang commented 6 months ago

This is a milestone! Thank everyone!

nbelakovski commented 6 months ago

Woo hoo!