Broken inplace integrand for HCubatureJL

[lxvm]

Hi,
I was writing a quadrature code with inplace integrands and I noticed that Integrals.jl has the same problem I noticed. Depending on how the quadrature library allocates arrays, it is possible that multiple integrand evaluations are made before they are summed together. When the result is overwritten to the same array this means some of the quadrature nodes are overwritten by the values of another. For example, integrating sin(x) on [-1,1] as below should give zero on the initial segment, but the inplace evaluator seems to integrate sin(abs(x)) on that segment.

julia> using Integrals

julia> prob = IntegralProblem((x,p) -> sin(p*x), -1, 1, 3.0)
IntegralProblem. In-place: false


julia> inplaceprob = IntegralProblem((y,x,p) -> y .= sin(p*x), -1.0, 1.0, 3.0)
IntegralProblem. In-place: true


julia> solve(prob, HCubatureJL(), maxiters=15)
u: 0.0

julia> solve(inplaceprob, HCubatureJL(), maxiters=15)
u: 1-element Vector{Float64}:
 1.3044979103237169

The problem is in the definition of the integrand wrapper

dx = zeros(eltype(lb), prob.nout)
f = x -> (prob.f(dx, x, prob.p); dx)

An issue with dx is that it may have the wrong type, since the user's function could return complex values at real nodes, for example. A fix is to allow the user to pass in the output array. The second issue in the example above can be fixed by returning a new array at each integrand evaluation, which can be done by returning one(eltype(lb))*dx. The downside to this last fix is that most of the array operations aren't inplace, but the only way to achieve this is if the library supports them (e.g. see quadgk! in QuadGK.jl for an implementation).

[ArnoStrouwen]

Yes, I think a breaking change is needed to fix this.
Making the user pass their own prototype container seems like the best solution to me also.
This would then also get rid of the nout keyword.

[ChrisRackauckas]

I think that's a good suggestion.

[lxvm]

Perhaps a combination of wrappers like InplaceIntegrand and BatchIntegrand could replace the nout and batch keywords and each algorithm can dispatch these integrands to the appropriate calls? For a batched and vectorized integrand we would want something like a BatchIntegrand(InplaceIntegrand()), although we could also make a separate wrapper

[lxvm]

Batch integrands could work similarly to how they will in QuadGK.jl after this PR gets merged JuliaMath/QuadGK.jl#80

[lxvm]

SciML PR SciML/SciMLBase.jl#497

[ChrisRackauckas]

Now that the SciMLBase PR has merged, Integrals.jl only supports SciMLBase v1 right now because the integrand interface is the v2 breaking change. This is the issue that has to be solved for the bump to v2 to occur.

[lxvm]

Yep, I'm working on a pr for it, but I have to change most of the library since scimlbase v2 breaks all the features. I also have to make a PR to scimlbase and I'll try to wrap them up today

[ChrisRackauckas]

Fixed by new interface 🎉 🎉 🎉 🎉 🎉