Setting limits of integration in solve to avoid new QuadratureProblem

[fintzij]

Hello,

Is there a way to specify the limits of integration in the solve() call instead of when I create a new QuadratureProblem? The problem I am solving involves repeatedly integrating the same function, but over different intervals that are not necessarily known at compile time. I would like to avoid having to create a new instance of QuadratureProblem every time I get a new set of intervals. The context for this is a package for Bayesian inference for time to event models where sample paths are sampled as part of the MCMC algorithm (the limits of integration change at each step of the Markov chain) and the likelihood depends on the integrated hazard over each inter-event interval.

Looking at the source code, it seems like perhaps something like this is possible since one of the signatures for solve is: solve(prob::QuadratureProblem, ::Nothing, sensealg, lb, ub, p, args...; reltol, abstol, kwargs...) in Quadrature at /Users/fintzijr/.julia/packages/Quadrature/UCzIP/src/Quadrature.jl:139 However, I'm not entirely sure how to initialize a QuadratureProblem and solve it using this signature.

For example, how should I modify the following?

using Quadrature

# this works 
f(x,p) = p
qp = QuadratureProblem(f, 1.0, 3.2, 10.0)
solve(qp, QuadGKJL())

# I'd like to be able to do something like this
solve(qp, QuadGKJL(); lb = 0.0, ub = 5.0)

I've also tried extending Quadrature to define a QuadratureProblem that is a mutable struct so that I could modify qp.lb and qp.ub directly. I'm sure that there are very good reasons why instances of QuadratureProblem are immutable, so this is probably a very silly thing to do for other reasons and maybe it's good that it didn't work.

Thank you!

[ChrisRackauckas]

The Julia compiler will optimize out immutable constructions in local instances, so remake on QuadratureProblem to change the bounds is pretty lightweight. So I wouldn't worry about the performance impact here. Did you see that in your case it actually heap allocates and doesn't just perform the operation on the stack?

[fintzij]

Ah thanks!! I didn't know about the remake function, that's super helpful! For this example, it's about a 10x speedup vs. creating a new QuadratureProblem each time (timings below).

Regarding the second point about heap vs. stack, we're only solving one-dimensional integrals (we just happen to have a lot of them) of fairly well-behaved functions. Based on some small trial timings, it seemed like using GuadGKJL() would probably be the most performant algorithm, but I couldn't figure out how to do the integration in-place. Interestingly, the non-in-place version seems to outperform the in-place version when I use the HCubatureJL() backend. Maybe I'm doing something wrong?

Code for timings:

exph(x,p) = p
function dexph(dx,x,p) 
    dx .= p
end

qp = QuadratureProblem(exph, 1.0, 2.2, 10.0)
qp2 = QuadratureProblem{true}(dexph, 1.0, 2.2, 10.0)
solve(qp, HCubatureJL())
solve(qp2, HCubatureJL())

@time begin
    for k in 1:1000000
        # remake has minimal overhead
        # solve(qp, QuadGKJL())
        solve(remake(qp, lb = 1.0, ub = 2.2), QuadGKJL())

        # Making a new QuadratureProblem each time is slowwww
        # qp = QuadratureProblem(exph, 1.0, 2.2, 10.0)
        # solve(qp, QuadGKJL())

        # HCubatureJL also slow, even in-place
        # solve(remake(qp2, lb = 1.0, ub = 2.2), HCubatureJL())

        # HCubatureJL non-in-place is faster than in-place
        # solve(remake(qp, lb = 1.0, ub = 2.2), HCubatureJL())
    end
end

[ChrisRackauckas]

@YingboMa check this when you look at the stability stuff.

[YingboMa]

See #65 (comment) we cannot really make the current API fast because it's inherently expensive.

[fintzij]

Hi @ChrisRackauckas,

I am reopening this to ask whether it is possible to use remake to alter data that gets passed to the integrand or if there is another strategy that would make sense. I would like to pass an argument to the integrand, like the value of a covariate, that is not the variable of integration and is not a parameter that I want exposed to autodiff, but I also still want to avoid creating a new instance of (now) IntegralProblem each time the data changes since there is some overhead to this. I'll explain my confusion:

• The signature for the integrand is strict: the first argument is the variable(s) of integration, the second consists of parameters, which I'm taking to mean any argument to the integrand that is fixed in a particular instance when we solve the IntegralProblem.
• If we are using quadrature within an optimization problem, say to fit a model to data, the parameters in the integrand are varying as the optimization routine explores the parameter space. Clearly, we want variables that are truly parameters to be exposed to autodiff, but it wouldn't make sense to treat data as such.
• I thought the following might work, but it doesn't:

using Integrals

# define integrand
# p is a true parameter, d is data. I want to change the value of d from call to call
g(x,p; d = 1.0) = x .* p .* d
ip = IntegralProblem(g, 1.0, 3.2, 10.0)
solve(ip, QuadGKJL()) # returns 46.2

# can't just pass a new value for d via kwargs
# so try remaking ip with a new integrand
g2(x,p;d=2.0) = x .* p .* d
remake(ip, f=g2)
solve(ip, QuadGKJL()) # returns 46.2

• I thought maybe I could solve the problem with clever use of closures, but actually I'm not sure exactly how that would work without creating a new IntegralProblem each time the data changes and incurring the overhead from this. Here's as far as I got with that:

function g2(lb, ub, p0, d) 
    IntegralProblem(
        (x,p) -> x .* p .* d,
        lb,
        ub,
        p0)
end
ip = g2(1.0, 3.2, 10.0, 2.1) # g2 creates a new IntegralProblem, so nothing solved
solve(ip, QuadGKJL())

Thanks

[ChrisRackauckas]

Remake won't handle that well. We're working on an enhanced parameter interface to solve this kind of thing, SciML/DifferentialEquations.jl#881

[fintzij]

Thanks! I'll stay tuned, closing this issue now.