giordano / Cuba.jl

Library for multidimensional numerical integration with four independent algorithms: Vegas, Suave, Divonne, and Cuhre.
https://giordano.github.io/Cuba.jl/stable
MIT License
75 stars 9 forks source link

Improve user interface #3

Closed giordano closed 8 years ago

giordano commented 8 years ago

Currently, users have to write their own integrand functions using this awkward boilerplate:

function integrand(ndim::Cint, xx::Ptr{Cdouble}, ncomp::Cint, ff::Ptr{Cdouble},
                   userdata::Ptr{Void})
    # Take arrays from "xx" and "ff" pointers.
    x = pointer_to_array(xx, (ndim,))
    f = pointer_to_array(ff, (ncomp,))
    # Do calculations on "f" here
    #   ...
    # Store back the results to "ff"
    ff = pointer_from_objref(f)
return Cint(0)::Cint
end

It would be better to let users directly call integrator functions with something like this:

Vegas( (x)-> cos(x), ndim[; keywords])

Commit https://github.com/giordano/Cuba.jl/commit/e61c7f10410471f6fc74ec6df83bf4f322efe598 in branch ui enables this feature, but that implementation (using eval) is really _slow_ ever for simple integrand functions

stevengj commented 8 years ago

You shouldn't need to use eval for this. Just have a single integrand function, something like:

function generic_integrand!(ndim::Cint, xx::Ptr{Cdouble}, ncomp::Cint, ff::Ptr{Cdouble}, func_::Ptr{Void})
    x = pointer_to_array(xx, (ndim,))
    f = pointer_to_array(ff, (ncomp,))
    func! = unsafe_pointer_to_objref(func_)::Function
    func!(f, x) # call the user's function to compute f = func(x), in-place
    return Cint(0)
end
const c_generic_integrand! = cfunction(generic_integrand!, Cint, (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble})

and then call Cubu by passing c_generic_integrand! and the user's function as the Ptr{Void} pointer, as in the Passing closures via pass-through pointers on the Julia blog.

Note that Cint(0) already does a type-assertion in Julia 0.4, so Cint(0)::Cint should not be necessary.

Note also that ff = pointer_from_objref(f) does nothing. (You are assigning to the pointer, not the contents of the pointer ... this is just changing the value of a formal parameter inside the function and is not visible to the caller.) Nor is it necessary to do anything: by storing the result in f[i], you change the data that ff points to.

This requires's the user's function to act in-place on the output. It would be more natural to just let the user return an iterable result (e.g. a single result for ncomp == 1 or a tuple for ncomp > 1). You can do this by:

function generic_integrand(ndim::Cint, xx::Ptr{Cdouble}, ncomp::Cint, ff::Ptr{Cdouble}, func_::Ptr{Void})
    x = pointer_to_array(xx, (ndim,))
    func = unsafe_pointer_to_objref(func_)::Function
    y = func(f, x)
    length(y) != ncomp && return Cint(1) # error
    for i = 1:length(y)
        unsafe_store!(ff, y[i], i)
    end
    return Cint(0)
end
const c_generic_integrand = cfunction(generic_integrand, Cint, (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble})

(I haven't tested this, so I may have some typos and other errors, but this should give you the general idea.)

giordano commented 8 years ago

Thanks for the suggestions! Looking at Cubature.jl for vector-valued functions, I also thought about using a function that acts in-place on the output, but I didn't actually tried that yet.

I tested both solutions you suggested, and I found that while the second one is more natural, the first one is faster, and I prefer speed over ease of use, I think I'll go for the first way. In addition, that would be exactly the same kind of functions required by Cubature.jl, so probably many users are accustomed to them.