auto wrapper take 2

[lxvm]

Alternative to @gsl_function. See the discussion in #130

TODO

• Implement a similar wrapper for gsl_function_fdf
• Implement a similar wrapper for gsl_function_vec
• Implement a similar wrapper for gsl_multiroot_function
• Implement a similar wrapper for gsl_multiroot_function_fdf
• Decide if/how to wrap the above in an api and add tests
• Make the @gsl_function... macros null-ops Fix and deprecate @gsl_function macros Keep macros for compatibility
• Update documentation

[yuyichao]

Ref https://github.com/JuliaMath/GSL.jl/tree/yyc/gsl_func

If you want to finish up the ones for other function structures that's fine. Otherwise I can finish it some time later.

[lxvm]

Thanks @yuyichao , I'd be happy to finish them up. I can add them as a pr to your branch just to keep the history more simple.

Also, should the changes to src/gen/direct_wrappers/ be made at the level of gen/makewrappers.jl so that they don't get over-written in the future?

[yuyichao]

So a few things,

The way the @gsl_function and @gsl_function_fdf works before is to do the conversion to the C types on the caller side, this is in general impossible to do without allocation and the proper mechanism to mark things as live during the ccall. This is what's fixed by using cconvert and unsafe_convert.

However, cconvert and unsafe_convert aren't special or magical on their own. They are only specially because ccall lowers to use them. It's ccall that keeps things properly alive and not cconvert, nor do any other local variable. As such, no local trick can extend the scope of caller side C types conversion beyond what @gsl_function and @gsl_function_fdf currently have to offer. Trying to use cconvert is therefore simply always wrong. The current implementation of the macro only works on the cases they do because that's the only case that can possibly work and it's a good thing that an error is thrown in the cases that can't work. It's of course possible to do some conversion and get some C types that might happen to work when the GC and the compiler doesn't feel like working too hard but in general, any attempt to make caller side conversion work in all cases are wrong.

So this leaves two possibility for what the macros can be used for:

1. To keep working to do caller side conversion in cases that can work
2. To simply do nothing and to pass this to the GSL functions to do the proper conversion (and then the macro can just be deprecated and removed later).

The second option allows it to be used in more cases though the first one is more backward compatible. It's conceivable that someone could be doing that and manually construct other structures containing gsl_function (or pointers to it, and there are a few candidate types in GSL AFAICT). The first option will keep that working and the second will break that. I usually prefer to err on the compatibility side especially when it's easy to do so. (edit: OTOH, I do agree and anyone that was doing this probably know what they are doing and can fix their code easily so I don't really mind turning the macro into a no-op either, as long as no-one come to claim that they have an actual use of this kind.)

And then there's the question of gsl_function_fdf. Just like what's mentioned above, using cconvert in the macro is just wrong, always. Related, cconvert has a fixed semantics in ccall and passing in more than two arguments to it is also just wrong as well. It should take a type and an object to convert and no more.

The first step to fix gsl_function_fdf is actually to determine the user API. In my experience there are two possible (and non-conflicting) ways to do this. One is to just pass in multiple functions as f and df (and maybe fdf), whereas the other one is to pass in an object with user customizable methods defined on them, e.g. GSL.function_f, GSL.function_df, etc... (the fdf function can have a fallback method that calls f and df, i.e. GSL.function_fdf(obj, x) = GSL.function_f(obj, x), GSL.function_df(obj, x)).

The latter may look/feel nicer though it can be annoy to use with closures so I would recommend still use the multiple callback one as the API, the other one can be added later if needed. And I also think it's safe enough to just allow the user to pass the functions in as a tuple. Both (f, df) and (f, df, fdf) can be supported and the corresponding Base.cconvert(::Type{Ref{gsl_function_fdf}}, ::NTuple{2}) and Base.cconvert(::Type{Ref{gsl_function_fdf}}, ::NTuple{3}) would need to be implemented to do the actual conversion in ccall.

[lxvm]

Yes, just adding a new cconvert method to convert tuples (f, df) or (f, df, fdf) into gsl_function_fdf structs seems like the easiest, non-breaking change to make. I'll restore the macro definitions and focus on writing and testing this new api for the remaining function types

[lxvm]

Actually, for consistency I think the new api for gsl_function could be to pass a (f,) to the routine instead of f, since this has the advantage we don't have to modify any of the generated routines to specialize on f::Function.

[yuyichao]

Actually, for consistency I think the new api for gsl_function could be to pass a (f,) to the routine instead of f, since this has the advantage we don't have to modify any of the generated routines to specialize on f::Function.

That's just a very strange API.

Also, it seems that there are other API's that are inherently unsafe to even use the ccall conversion. root_*solver_set methods uses the function pointer passed in beyond the calling scope and requires a complete different way of handling it.

[lxvm]

Fair enough, I was thinking a tuple-based api could be a building block for something higher-level. In light of the unsafe apis you mentioned, it could make sense to define a higher-level "safe" api for which a solver can consistently use a function pointer. Although you mentioned it is wrong to use cconvert outside of ccall, if the references returned by cconvert are stored in a mutable container does that protect them from the GC?

I'm still trying to understand if this is similar to the 'PyCall' solution in this discussion, although what I had in mind was this

Possible api

mutable struct GSLSolverWrapper{T<:Tuple,S}
    isset::Bool       # whether the solver has been given a gsl_function
    ref               # store the output of the gsl_function cconvert here
    solver::Ptr{S}
    fs::T
end

# default behaves like a gsl_function in a ccall, using the tuple cconvert methods
GSLSolverWrapper(fs...) = GSLSolverWrapper(false, nothing, nothing, C_NULL, fs)
Base.cconvert(T::Type{Ref{gsl_function}}, s::GSLSolverWrapper) = Base.cconvert(T, s.fs)

# for calling a function with a specific solver, pass the desired solver type
function GSLSolverWrapper(solvertype::Ref{gsl_root_fsolver_type}, f)
    solver = root_fsolver_alloc(solvertype)
    gsl = GSLSolverWrapper(false, (gsl_function(C_NULL, C_NULL), Ref((f,))), solver, (f,))
    finalizer(f -> root_fsolver_free(f.solver), gsl)
    return gsl
end

# e.g. GSLSolverWrapper(gsl_root_fsolver_brent, sin) then implements
function GSL.root_fsolver_set(solver::GSLSolverWrapper{F,gsl_root_fsolver}, a, b) where F
    if !solver.isset
        solver.ref = Base.cconvert(Ref{gsl_function}, solver.fs) # we allocate/store our own gsl_function to bypass ccall mechanism
    end
    root_fsolver_set(solver.solver, solver.ref[1], a, b) # creates a pointer to solver.ref in solver.solver
    solver.isset = true
end
function GSL.root_fsolver_set(solver::GSLSolverWrapper{Tuple{F},gsl_root_fsolver}, f::F, a, b) where F
    # implicitly free references
    solver.ref = (gsl_function(C_NULL, C_NULL), Ref((f,)))
    solver.fs = (f,)
    solver.isset = false
    root_fsolver_set(solver, a, b)
end
function GSL.root_fsolver_iterate(solver::GSLSolverWrapper{F,gsl_root_fsolver}) where F
    !solver.isset && error()
    root_fsolver_iterate(solver.solver)
end
function GSL.root_fsolver_root(solver::GSLSolverWrapper{F,gsl_root_fsolver}) where F
    !solver.isset && error()
    root_fsolver_root(solver.solver)
end

which behaves like

example

julia> struct Linear
       a::Float64
       end

julia> (p::Linear)(x) = x-p.a

julia> f = GSLSolverWrapper(gsl_root_fsolver_brent, Linear(0.73));

julia> root_fsolver_set(f, -1.0,1.0)
true

julia> root_fsolver_root(f)
0.0

julia> root_fsolver_iterate(f)
0

julia> root_fsolver_root(f)
0.73

julia> root_fsolver_set(f, Linear(0.2), -1.0,1.0)
true

julia> root_fsolver_root(f)
0.0

julia> root_fsolver_iterate(f)
0

julia> root_fsolver_root(f)
0.19999999999999996

This would be a lot of work to implement for each solver type. I'm more interested if this could be a viable new api that could later be extended to a higher level interface. I'd also have to check if this can express other unsafe apis in the library.

[lxvm]

I've gone ahead and added cconvert, unsafe_convert pairs for the remaining macros. As you've mentioned, they can't be used safely with the root_*solver types given ccall's semantics so they are unusable without a new api. I agree this would have to be done carefully

Just a cranky idea for a workaround: could we use llvm tokens to keep an object alive such as a gsl_function after a call to root_*solver_set by registering a finalizer on the root_*solver in that same call?

[yuyichao]

if the references returned by cconvert are stored in a mutable container does that protect them from the GC?

Not automatically. That mutable container has to be protected from the GC.

although what I had in mind was this

The details are not all correct but yes that's the idea.

As you've mentioned, they can't be used safely with the root_*solver types given ccall's semantics so they are unusable without a new api.

I have not looked at the use of other function types but at least for gsl_function there's also cheb_init, min_fminimizer_set* that are unsafe.

These, and other similar cases should be restricted on the function signature to make sure they do not accept the new input type.

[yuyichao]

Please don't. This tuple is totally unnecessary and it also create more opportunity for allocation.

[yuyichao]

This is wrong. Also I suggest using different helper for fdf directly which avoid the creation of a closure unnecessarily.

See 8a5a1d1

Though this is currently a moot point since there's no user of this.

[yuyichao]

I don't think this is correct either.

[yuyichao]

And see https://github.com/JuliaMath/GSL.jl/compare/yyc/gsl_func, which I believe is all that can be done without introducing new API. For the two FDF variance, I don't believe there's anything that can use them so there's basically no public API for them right now. The wrapper type is added for future API. Note that in this case, not even the macro version currently used is safe since they don't garantee that the object will be live between calls to GSL.

Edit: Actually while trying to make the safe wrapper I think the wrapper type is unnecessary and I'll revert those.

[lxvm]

Closing in favor of #132

@yuyichao I really appreciate your guidance and assistance with this work and I've learned a lot. Thank you for going the extra mile and writing a new interface for the improved functionality! Is there anything I can help with, such as tests?

[yuyichao]

Yeah, tests would be nice. I'm pretty sure all of the four function types have tests already.

The cconvert/unsafe_convert implementation for gsl_multiroot_function has one use in multiroot_fdjacobian that doesn't have a test.

Out of the safe wrapper, only the solvers have tests but not the GSLCheb and GSLMinFMinimizer.