`nvec` parameter not implemented?

[afniedermayer]

I was trying to use the nvec parameter. (I think nvec could be used to implement parallelism, this is the way the Mathematica interface to cubalib implements parallelism.) However, I couldn't figure out how to do it. Looking at the code, it looks like the nvec parameter isn't passed to the integrand function: https://github.com/giordano/Cuba.jl/blob/master/src/Cuba.jl#L86 and https://github.com/giordano/Cuba.jl/blob/master/src/Cuba.jl#L96

[giordano]

I think you're right! Feel free to provide a patch, I'm not sure I'll be able to have a look very soon.

[giordano]

I quickly tried this:

function generic_integrand!(ndim::Cint, x_::Ptr{Cdouble}, ncomp::Cint,
                            f_::Ptr{Cdouble}, func!, nvec::Cint)
    # Get arrays from "x_" and "f_" pointers.
    x = unsafe_wrap(Array, x_, (ndim, nvec))
    f = unsafe_wrap(Array, f_, (ncomp, nvec))
    func!(x, f)
    return Cint(0)
end

# Return pointer for "integrand", to be passed as "integrand" argument to Cuba functions.
integrand_ptr{T}(integrand::T) = cfunction(generic_integrand!, Cint,
                                           (Ref{Cint}, # ndim
                                            Ptr{Cdouble}, # x
                                            Ref{Cint}, # ncomp
                                            Ptr{Cdouble}, # f
                                            Ref{T}, # userdata
                                            Ref{Cint})) # nvec

But I'm not sure this suffices.

[afniedermayer]

I've just tried the very same thing :-) I'm also not sure whether that's enough, I'll try to test it once I have more time...

Thanks for looking into this and more generally thanks for providing this package!

[giordano]

Uhm, with this change the the result of integration isn't even correct and with nvec > 1 the function is also slower:

julia> cuhre((x, f) -> f[] = x[], nvec = 1)
Component:
 1: 0.5 ± 4.291960664289122e-15 (prob.: 0.0)
Integrand evaluations: 195
Fail:                  0
Number of subregions:  2

julia> cuhre((x, f) -> f[] = x[], nvec = 2)
Component:
 1: 0.37544225534137515 ± 0.0030736506895798653 (prob.: 0.0)
Integrand evaluations: 1000025
Fail:                  1
Number of subregions:  7693

julia> cuhre((x, f) -> f[] = x[], nvec = 4)
Component:
 1: 0.12829062941575486 ± 0.004796612531625397 (prob.: 0.0)
Integrand evaluations: 1000025
Fail:                  1
Number of subregions:  7693

[giordano]

Looking better to how the x and f arrays are treated when nvec > 1, probably I was using the wrong function calls:

julia> cuhre((x, f) -> f[1,:] = x[1,:], nvec = 1)
Component:
 1: 0.5 ± 4.291960664289122e-15 (prob.: 0.0)
Integrand evaluations: 195
Fail:                  0
Number of subregions:  2

julia> cuhre((x, f) -> f[1,:] = x[1,:], nvec = 2)
Component:
 1: 0.5 ± 4.291960664289122e-15 (prob.: 0.0)
Integrand evaluations: 195
Fail:                  0
Number of subregions:  2

julia> cuhre((x, f) -> f[1,:] = x[1,:], nvec = 4)
Component:
 1: 0.5 ± 4.291960664289122e-15 (prob.: 0.0)
Integrand evaluations: 195
Fail:                  0
Number of subregions:  2

which now are indeed faster than nvec = 1:

julia> @benchmark cuhre((x, f) -> f[1,:] = x[1,:], nvec = 1)
BenchmarkTools.Trial: 
  memory estimate:  67.52 KiB
  allocs estimate:  1374
  --------------
  minimum time:     33.144 μs (0.00% GC)
  median time:      35.663 μs (0.00% GC)
  mean time:        48.058 μs (22.65% GC)
  maximum time:     3.921 ms (98.52% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark cuhre((x, f) -> f[1,:] = x[1,:], nvec = 2)
BenchmarkTools.Trial: 
  memory estimate:  34.52 KiB
  allocs estimate:  702
  --------------
  minimum time:     20.266 μs (0.00% GC)
  median time:      21.898 μs (0.00% GC)
  mean time:        28.492 μs (19.39% GC)
  maximum time:     4.049 ms (98.82% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark cuhre((x, f) -> f[1,:] = x[1,:], nvec = 4)
BenchmarkTools.Trial: 
  memory estimate:  18.77 KiB
  allocs estimate:  366
  --------------
  minimum time:     14.691 μs (0.00% GC)
  median time:      15.964 μs (0.00% GC)
  mean time:        19.161 μs (11.79% GC)
  maximum time:     2.996 ms (98.78% GC)
  --------------
  samples:          10000
  evals/sample:     1