search

JuliaMath / HCubature.jl

pure-Julia multidimensional h-adaptive integration
Other
148 stars 23 forks source link

Inconsistent results for Inf integrands #29

Open agerlach opened 3 years ago

agerlach commented 3 years ago

I am looking at some SciML problems where, depending on the stability of the ODE, the integrand used w/ HCubature.jl may return Inf. I am noticing inconsistent behavior across hquadrature and hcubature for single and multi-dim domains. 

using HCubature

g(x) = x[1] > 0.0 ? Inf : x[1]

hquadrature(g, -1.0, 1.0)    #Inf
hcubature(g, -ones(1), ones(1)) #Inf
hcubature(g, -ones(2), ones(2)) #NaN

I see similar behavior in Cubature.jl as well, but Cubature.hcubature returns Inf up to dim = 4 and then NaN, while Cubature.pcubature returns Inf up dim = 17 (haven't test passed that).

Is this the expected behavior? I assume this is due to some Inf*0 that is propagating down.

Thanks.

stevengj commented 3 years ago

I'm not sure off the top of my head what might give a NaN, but probably only a small patch is required to return Inf.

agerlach commented 3 years ago

Sign changes in coefficients g.w leads to I=NaN due to Inf-Inf https://github.com/JuliaMath/HCubature.jl/blob/66074a8b75f71066624517ea9a89b4f28e963b0f/src/genz-malik.jl#L150

I am happy to put together a PR to address this corner case, but it is unclear to me what the appropriate behavior should be. My initial thought is to short-circuit the algorithm and return Inf if f ever returns Inf.

agerlach commented 3 years ago

On second thought short-circuiting would not be a good approach b/c integrand could be vector-valued mixing Inf and finite values.