cuhre fails at high precisions (?)

[antoine-levitt]

The question mark in the title is because it seems to fail, but really I have no idea what the correct result is.

The attached code compares Cuba's cuhre with Cubature's hcubature. The results differ by about 1e-3, despite an abstol of 1e-6 and a reltol of 0, and no limit in maxevals. Both integrators report success. I'm filing this here (rather than at Cubature) because Cubature's results are more in line with what I expect, but I can't exclude an error in hcubature. Both integrators stick to their guns when increasing the precision. The integrand is an indicator function of a polygon

using Interpolations
using Cubature
using Cuba
using PyPlot

const xmax = 1.0
const reltol = 0.0
const abstol = 1e-6

const P1 = BSpline(Linear())

ε(x,y) = 3cos(2pi*x)*cos(2pi*y)+sin(4pi*x)*cos(4pi*y)
function test(p,N)
    x = 0:1/N:xmax-xmax/N
    y = x

    z = ε.(x,y')

    itp = interpolate(z, p, OnGrid())
    εp = scale(itp,x,y)

    function f(x,y,εF)
        εqxy = εp[x,y]
        εqxy < εF ? 1.0 : 0.0
    end

    function error(εF)
        res = cuhre(((x,out)->(out[1] = f(x[1],x[2],εF))), 2,1, reltol=reltol, abstol=abstol,maxevals=1e10)
        println()
        println(res[1][1])
        # println(res)
        res = hcubature(x->f(x[1],x[2],εF), (0,0), (xmax,xmax), reltol=reltol, abstol=abstol, maxevals=0)
        println(res[1])
    end

    error(1.95)
end

test(P1,10)

[giordano]

Honestly, I don't think there is anything I can do about this. If it's an issue it's most probably an upstream issue. However, you should really work out the integral to be sure what the expected result is, in order to claim there is a bug somewhere, but I acknowledge this is a non-trivial function.

Just to add two other data-points I computed the integral also with HCubature.jl and NIntegration.jl:

using Interpolations
using Cuba, Cubature, HCubature, NIntegration

const xmax = 1.0
const reltol = 0.0
const abstol = 1e-6

const P1 = BSpline(Linear())

ε(x,y) = 3cos(2pi*x)*cos(2pi*y)+sin(4pi*x)*cos(4pi*y)
function test(p,N)
    x = 0:1/N:xmax-xmax/N
    y = x

    z = ε.(x,y')

    itp = interpolate(z, p, OnGrid())
    εp = scale(itp,x,y)

    function f(x,y,εF)
        εqxy = εp[x,y]
        εqxy < εF ? 1.0 : 0.0
    end

    function error(εF)
        res_cuba = cuhre(((x,out)->(out[1] = f(x[1],x[2],εF))), reltol=reltol, abstol=abstol,maxevals=1e10)
        println(res_cuba[1][1], " ± ", res_cuba[2][1], " (Cuba.jl)")
        res_cubature = hcubature(x->f(x[1],x[2],εF), (0,0), (xmax,xmax), reltol=reltol, abstol=abstol, maxevals=0)
        println(res_cubature[1], " ± ", res_cubature[2], " (Cubature.jl)")
        res_hcubature = HCubature.hcubature(x->f(x[1],x[2],εF), (0,0), (xmax,xmax), rtol=reltol, atol=abstol, maxevals=0)
        println(res_hcubature[1], " ± ", res_hcubature[2], " (HCubature.jl)")
        res_nintegration = nintegrate((x,y,z)->f(x,y,εF), (0,0,0), (xmax,xmax,xmax), reltol=reltol, abstol=abstol, maxevals=1000000)
        println(res_hcubature[1], " ± ", res_hcubature[2], " (NIntegration.jl)")
    end

    error(1.95)
end

test(P1,10)

This is the output:

0.8840867311005485 ± 9.99974855348049e-7 (Cuba.jl)
0.8848940796576487 ± 9.999981972972212e-7 (Cubature.jl)
1.1532286744906772 ± 1.041558705481888 (HCubature.jl)
1.1532286744906772 ± 1.041558705481888 (NIntegration.jl)

If I didn't mess up anything with those package, they seem to provide the same highly inaccurate result, with the same large error.

[giordano]

The function is discontinuous, so the result most probably depends on sampling. I think that this should give an approximate evaluation of the integral, right?

julia> mean(f.(0:0.0001:1,(0:0.0001:1)',1.95))
0.8840487514092307

If this is correct, Cuba.jl gives the most accurate result among the packages tested above

[giordano]

This is the smallest sampling which doesn't eat all my memory up:

julia> mean(f.(0:0.00005:1,(0:0.00005:1)',1.95))
0.8840676910207287

The result seems to be converging to Cuba.jl answer.

[antoine-levitt]

Indeed, thanks for the detective work! Using a generator rather than f. to save memory, I get 0.8840827781 with 0:0.00001:1. So that would be a success for Cuba and a failure for Cubature. Maybe the takeaway is: don't integrate discontinuous functions with Cubature. CC @stevengj.

[giordano]

Well, I think that any integrator may fail with discontinuous function. Perhaps Cuba in this case has been lucky to get the right sampling.

Side note: if you know that a function has peaks and you also know their positions, the divonne integrator in Cuba.jl is particularly useful, because you can provide it with the positions of the peaks, in order to increase sampling around them.

[antoine-levitt]

Alright, thanks and sorry for the noise (in any case this has nothing to do with the julia bindings).

[giordano]

No problem, it has been interesting to look at this :-) I guess you may want to report the issue to the other packages, though.