User-defined callback routines

[ivan-pi]

Before we dive into the implementation of root-solvers, optimization routines, and quadrature routines (see #87 and #112 ) I think it is worth discussing the "standard" way of defining callback routines (or in other words user-defined functions) which are potentially parametrized by some variables.

Some of the possibilities are described at https://www.fortran90.org/src/best-practices.html#callbacks

The NAG library uses a very low level approach, where the user function carries along both integer and real user workspaces:

Function f(x, iuser, ruser)
Real (Kind=nag_wp)	:: 	f
Integer, Intent (Inout)	:: 	iuser(*)
Real (Kind=nag_wp), Intent (In)	:: 	x
Real (Kind=nag_wp), Intent (Inout)	:: 	ruser(*)

Do we want to support a high level API were the users extend a base class (see related issue #27 ) that carries an evaluate method, as suggested already in the comment by @rweed #87 (comment)?

In the ideal case, the n-D root finding routines, the optimization routines, and quadrature routines could share the same callback interface.

What should we do for routines that require both the function value and its derivative? Should we write multiple functions, as in

real function quad(x)
  real, intent(in) :: x
  quad = 2.0*x**2 + 3.0
end function
real function dquad(x)
  real, intent(in) :: x
  dquad = 4.0*x
end function

or would a subroutine interface `like

subroutine quadd(f,df,x)
  real, intent(in) :: x
  real, intent(out) :: f
  real, intent(out), optional :: df
  f = 2.0*x**2 + 3.0
  if (present(df)) df = 4.0*x
end subroutine

be preferable?

The downside of such a subroutine interface is that it cannot be nested in other expressions.I know that personally, from performance reasons, I sometimes use the second approach, as typically the function and the derivative share some common sub-expressions that can be reused, potentially saving a few CPU operations.

It is of course possible to support several interfaces using "adaptor" routines, but this introduces extra calls, de-referencing issues, and just makes things more complicated.

[ivan-pi]

To provide a more concrete example of low vs high level callbacks I will borrow an example from the NLopt library (written in C) that defines a Fortran 77 interface. In the original low level interface the user needs to provide a function of the form:

  subroutine square(val, n, x, grad, need_gradient, f_data)
  double precision val, x(n), grad(n)
  integer n, need_gradient
  if (need_gradient .ne. 0) then ! gradient is only needed by certain optimization methods
     grad(1) = 0.0
     grad(2) = 0.5 / dsqrt(x(2))
  endif
  val = dsqrt(x(2)) ! value of f(x) 
  end 

The f_data argument can be used to pass through a single variable containing any data and any type to your subroutine. (it is not even type-checked!!!)

In my refactored version of the NLopt Fortran interface (https://github.com/ivan-pi/nlopt) the user is required to extend an abstract base class nlopt_user_func carrying a deferred method eval:

module functor_mod
    use iso_c_binding, only: c_int, c_double
    use nlopt, only: nlopt_user_func

    implicit none

    type, extends(nlopt_user_func) :: square
    contains
        procedure, public :: eval => eval_square
    end type

contains

    real(c_double) function eval_square(this,n,x,grad)
        class(square), intent(in) :: this
        integer(c_int), intent(in) :: n
        real(c_double), intent(in) :: x(n)
        real(c_double), intent(out), optional :: grad(n)
        if (present(grad)) then
            grad(1) = 0.0_c_double
            grad(2) = 0.5_c_double/sqrt(x(2))
        end if
        eval_square = sqrt(x(2))
    end function
end module functor_mod

A benefit of the "functor" approach is that it allows the user to encapsulate any parameters of the function in the extended derived type, instead of resorting to non-type-safe methods or array workspaces.

[nshaffer]

In my own codes, Fortran's contains mechanism usually makes this a non-issue. As a sketch, consider a 1-D Newton-like rootfinder with the interface

function root(f, df, x0, ...)
    interface
        real function f(x)
            real, intent(in) :: x
        end function f
        real function df(x)
            real, intent(in) :: x
        end function df
    end interface
    real, intent(in) :: x0
    ...
end function root

If I need to call this inside a subroutine foo, then I define f and df in the contains section of foo. That way, f and df can refer (dynamically!) to any variables accessible within foo. This is usually all the flexibility that I need. The library defines what interfaces it needs, the user provides functions or subroutines to match, and the correctness of the user code can be determined at run time. The main thing you lose in this approach is that the contained procedures are internal to the subprogram they're defined in. I have not found this to actually matter in my own codes, but my experience may not be everyone's.

[ivan-pi]

Your example reminds me a bit of how one uses lambda functions to call the solve_ivp routine in Scipy which expects a function with the interface

def f(t,y):
    dy = ...
    return dy

If you have a function that requires some arguments, e.g. fun(t,y,args), you just wrap it in a lambda function:

sol = solve_ivp(fun=lambda t, y: fun(t, y, *args), ...)

Building upon your rootfinder example, you could use the functions externally from the calling subroutine by defining them first in a module

module funcs
implicit none
contains
real function f(x,a,b,c)
  real, intent(in) :: x, a, b, c
  f = a*x**2 + b*x + c
end function
real function df(x,a,b)
  real, intent(in) :: x, a, b
  df = 2*a*x + b
end function
end module

And then internally in the routine calling the root function only define wrapper functions to conform to the interface:

subroutine call_root(x0,y,x,a,b,c)
  use funcs
  real, intent(in) :: x0
  real, intent(out) :: x, y
  real, intent(in) :: a, b, c
  x = root(local_f, local_df, x0, ...)
  y = f(x,a,b,c)
contains
  real function local_f(x)
    real, intent(in) :: x
    local_f = f(x,a,b,c)
  end function f
  real function local_df(x)
    real, intent(in) :: x
    local_df = df(x,a,b)
  end function df  
end subroutine

This introduces some boilerplate code to make sure the interfaces conform. Since we are focusing on ease of use, I guess we can neglect the performance penalty incurred by having two function calls (I doubt we can make the callback routine interfaces pure).

I suppose the question I am asking is should we support the different callback options (e.g. simple scalar callbacks, callbacks using array workspaces, extended derived types, reverse communication interfaces...) from the start or do we stick to a single interface and force the user to worry about it?

[certik]

@nshaffer you described the "nested functions" approach: https://www.fortran90.org/src/best-practices.html#iv-nested-functions. It works well, and it is my method of choice also, but there is one caveat that we recently discovered and discussed at length with @pbrady, @nncarlson and @zjibben privately:

https://nullprogram.com/blog/2019/11/15/

There are two approaches how the nested functions (that can access the parent function scope) are implemented in a compiler. GFortran uses the so called "trampoline" as described in the article, with the security implications and a possible performance penalty. The NAG compiler avoids the trampoline, but introduces intermediate wrapper functions which compute the local environment. Again, there is a performance penalty. The question is how much of a performance penalty there is. We need to benchmark it.

@milancurcic this is another candidate for our benchmarking repository.

P.S. We also discussed that if the compiler has access to the source code, it could inline the integrator into user's code, which will eliminate the nested function (it gets inlined) and this whole performance issue disappears. However, I don't know if any compiler actually does that. That would be the optimal solution.

[ivan-pi]

P.S. We also discussed that if the compiler has access to the source code, it could inline the integrator into user's code, which will eliminate the nested function (it gets inlined) and this whole performance issue disappears.

This sounds like exactly what the Julia JIT compiler is doing. As stated here by Chris Rackauckas:

There's a lot of interesting things that come into play for small ODEs. Some things are JIT compiler tricks, like JIT compiling f calls into the solver (emphasis added) (C FFI is rather expensive in comparison to what a JIT can do https://github.com/dyu/ffi-overhead which is bring that to zero in this case)