Cache for QuadGKJL

[lxvm]

Hi, this pr makes the IntegralCache mechanism nontrivial and adds a working cache for the QuadGKJL algorithm. What makes this cache somewhat tricky is that it needs to be rebuilt eagerly, i.e. whenever the integrand or parameters are reset, we allocate a new cache since these changes could change the return type of the integrand. Consider the example below

using Integrals
f(x, p) = sin(x * p) # real-valued integrand
prob = IntegralProblem(f, 0, 1, 14.0)
cache = init(prob, QuadGKJL()) # cache for real values
solve!(cache) # u: 0.061661627270869046
g(x, p) = cis(x * p) # complex-valued integrand
cache = Integrals.set_f(cache, g) # new cache for complex values
solve!(cache) # u: 0.07075766826391927 + 0.061661627270869046im

If we could let the caller specify that the cache need not be changed, since they guarantee that the return type of the integral is unchanged, I think this would let us reuse the allocated cache across multiple calls. Otherwise I think we cannot safely reuse an immutable, concretely-typed cache without some guarantee of type stability when changing the integrand or input parameters.

In summary, the following cache-resetting methods were added

• set_f(cache, f, [nout]): rebuilds cache
• set_lb(cache, lb): does not rebuild cache
• set_ub(cache, ub): does not rebuild cache
• set_p(cache, p): rebuilds cache

I also made an effort to get the AD to work with these changes. ForwardDiff.jl should still be working, however I still run into #99 with Zygote.jl.

I would appreciate any feedback on the design and will gladly implement those changes!

[ChrisRackauckas]

In general it looks good. Any new interface like this needs to make sure it's documented.

We might want to consider making the cache mutable. This was recently changed in LinearSolve. It was done there mostly for correctness reasons SciML/LinearSolve.jl#303, but it means that Integrals is the only thing on an immutable cache now. There may still be reasons for it, but it is the odd one out now.

[lxvm]

Sounds good, I added documentation and tests. I can see why a mutable cache gives more flexibility, though for integral problems all of the information needed to build the cache is contingent on the types of the limits and the integrand. To show this, I changed the QuadGKJL cache to use Base.promote_op to infer the types it needs to allocate the cache. (In the case that the inferred type is not concrete, I don't allocate any cache since QuadGK.jl may have to deal with a type instability anyway.)

To be thorough, I can list what I think are pros and cons of mutable caches are for different algorithms that I know can currently use caches (there may be more)

1. FastGaussQuadrature.jl

• About: The algorithm is a quadrature rule, which we can store in the cache
• Pros: Mutability means we can update the cache (at solve time) after we update the algorithm (at init time)
• Cons: If solve!(cache) becomes type unstable because of changes to the cache type at solve time, then the caller will pay. In the case of nested integrals, i.e. an IntegralProblem that solves an inner IntegralProblem, this can lead to slow compilation.

2. QuadGK.jl

• About: The algorithm can cache a heap consisting of quadrature rule evaluations
• Pros: Allocating the heap at solve time would make it possible to let quadgk allocate the heap (though there is currently no mechanism to extract the heap unless it is the you pass one in)
• Cons: If we pass in a heap with the wrong types, then we may break quadgk (this is why I only allocate a cache at init time if Base.promote_op succeeds)

3. HCubature.jl

• About: Although hcubature_buffer is still waiting to be released, the points are identical to QuadGK.jl

In summary, I believe constructing the cache can always be done at init time since the cache requires a new allocation that the algorithm would otherwise make itself, and so building a cache won't modify the problem inputs. If there are algorithms for which this is not true, then I would be compelled to make the cache mutable

Additionally, I noticed that Setfield.@set! leads to inference failing on the return type of the cache for the set_X functions and I was wondering if there is a common way to fix it?

[ChrisRackauckas]

@ArnoStrouwen is there a reason this never merged?

[ArnoStrouwen]

@lxvm have a look at JuliaFormatter.jl, use it to get this PR into SciML style. Hopefully the tests will all pass and then we can merge.

[lxvm]

@ArnoStrouwen I updated the format , but it seems that codecov is failing on the new set_X functions for updating the cache. Do we still think an immutable cache is the best solution given that a concretely-typed cache is good for performance and changing the integrand in the cache would change the cache type?

I also wanted to request a code review on whether or not the user should be able to control allocation of a segment buffer for QuadGK.

[ChrisRackauckas]

I think an immutable cache style with an init function and solve!, matching the rest of SciML, would be the most desirable since it would then be uniform to all of the other solvers and not an odd ball, and we have seen a few cases where that improves performance and correctness more.

changing the integrand in the cache would change the cache type?

That's expected. I think for integrals the main use case for the cache would be changing parameters? If we setup function wrappers this will not be the case anyways, since then all functions are the same type.

[lxvm]

Ok, so to be consistent with the other SciML libraries I can make the cache mutable, specialize setproperty! for the cache type if necessary, and add all of this to the documentation.

[lxvm]

@ChrisRackauckas the cache is now mutable and I documented that its fields cannot change type when updated. I'm happy with the PR as-is, although I was wondering what you meant by setting up function wrappers? I assume you mean f = x -> prob.f(x, p), but if you meant using https://github.com/yuyichao/FunctionWrappers.jl then my question about using Base.promote_op to infer the integrand's type could be answered.

[ChrisRackauckas]

but if you meant using https://github.com/yuyichao/FunctionWrappers.jl then my question about using Base.promote_op to infer the integrand's type could be answered.

I meant this, and instead of promote_op it would need Core.return_type, or we change the interface to have that passed in (preferred) like is being discussed in #169

[ChrisRackauckas]

This is fantastic. Thanks so much for taking it to completion!