Setting limits of integration in solve to avoid new QuadratureProblem about integrals.jl HOT 7 CLOSED

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?

from integrals.jl.

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

from integrals.jl.

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

from integrals.jl.

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

from integrals.jl.

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

from integrals.jl.

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

from integrals.jl.

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

from integrals.jl.