sciml / integrals.jl Goto Github PK

A common interface for quadrature and numerical integration for the SciML scientific machine learning organization

Home Page: https://docs.sciml.ai/Integrals/stable/

License: MIT License

Julia 100.00%

quadrature integration automatic-differentiation algorithmic-differentiation differentiable-programming sciml scientific-machine-learning julia julia-language julialang numerical-integration

integrals.jl's Introduction

Integrals.jl

Integrals.jl is an instantiation of the SciML common IntegralProblem interface for the common numerical integration packages of Julia, including both those based upon quadrature as well as Monte-Carlo approaches. By using Integrals.jl, you get a single predictable interface where many of the arguments are standardized throughout the various integrator libraries. This can be useful for benchmarking or for library implementations, since libraries which internally use a quadrature can easily accept a integration method as an argument.

Tutorials and Documentation

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.

Examples

To perform one-dimensional quadrature, we can simply construct an IntegralProblem. The code below evaluates $\int_{-2}^5 \sin(xp)~\mathrm{d}x$ with $p = 1.7$. This argument $p$ is passed into the problem as the third argument of IntegralProblem.

using Integrals
f(x, p) = sin(x * p)
p = 1.7
domain = (-2, 5) # (lb, ub)
prob = IntegralProblem(f, domain, p)
sol = solve(prob, QuadGKJL())

For basic multidimensional quadrature we can construct and solve a IntegralProblem. Since we are using no arguments p in this example, we omit the third argument of IntegralProblem from above. The lower and upper bounds are now passed as vectors, with the ith elements of the bounds giving the interval of integration for x[i].

using Integrals
f(x, p) = sum(sin.(x))
domain = (ones(2), 3ones(2)) # (lb, ub)
prob = IntegralProblem(f, domain)
sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3)

If we would like to parallelize the computation, we can use the batch interface to compute multiple points at once. For example, here we do allocation-free multithreading with Cubature.jl:

using Integrals, Cubature, Base.Threads
function f(dx, x, p)
    Threads.@threads for i in 1:size(x, 2)
        dx[i] = sum(sin, @view(x[:, i]))
    end
end
domain = (ones(2), 3ones(2)) # (lb, ub)
prob = IntegralProblem(BatchIntegralFunction(f, zeros(0)), domain)
sol = solve(prob, CubatureJLh(), reltol = 1e-3, abstol = 1e-3)

If we would like to compare the results against Cuba.jl's Cuhre method, then the change is a one-argument change:

using Cuba
sol = solve(prob, CubaCuhre(), reltol = 1e-3, abstol = 1e-3)

integrals.jl's People

Contributors

Stargazers

Watchers

integrals.jl's Issues

Register

@JuliaRegistrator register()

Fix forward mode derivatives overloads

The non-Julia packages have forward mode overloads but it needs a fix still.

Finite transformations are not supported for non-finite bounds with inplace mass function.

Integrand definition for batch evaluations

Right now it looks like the integrand needs to be written specially to handle batch evaluation. Can we update the batch evaluation to work of off broadcast instead so the integrand doesn't need to be "batch aware"?

How does this impact GPU support?

Specifying maxiters and CubaSUAVE in solve causes julia to crash (segmentation fault)

Specifying maxiters with CubaSUAVE() causes julia to crash, e.g.

using Integrals, IntegralsCuba
f(x,p) = sum(sin.(x))
prob = IntegralProblem(f,ones(2),3ones(2))
sol = solve(prob,CubaSUAVE(),reltol=1e-3,abstol=1e-3, maxiters=10)

Also I'd like to check - is maxiters meant to be a hard limit for the number of function calls? I would have expected the below to have very different results, but they give the exact same value.

using Integrals
f(x,p) = sum(sin.(x))
prob = IntegralProblem(f,ones(2),3ones(2))
sol1 = solve(prob,HCubatureJL(),reltol=1e-8,abstol=1e-8, maxiters=1).u
sol2 = solve(prob,HCubatureJL(),reltol=1e-8,abstol=1e-8, maxiters=10).u
sol1 == sol2

Derivatives w.r.t. boundary for higher dimensional integrals

Just Stokes rules to do an integral on the boundary. Right now it's only written for 1 dimension, but both forward and reverse mode overloads should fix that.

Add adjoint definitions

Zygote.gradient failed with CubaCuhre and batch !=0

using Quadrature, ForwardDiff, FiniteDiff, Zygote, Cuba
f(x,p) = sum(sin.(x .* p))
lb = ones(3)
ub = 3ones(3)
p = [1.5,2.0,3.0]

prob = QuadratureProblem(f,lb,ub,p; batch=0)
sol = solve(prob,CubaCuhre(),reltol= 1e-4,abstol= 1e-4,maxiters=100)[1]

function testf(p)
    prob = QuadratureProblem(f,lb,ub,p, batch=0)
    solve(prob,CubaCuhre(),reltol= 1e-4,abstol= 1e-4,maxiters=100)[1]
end
dp1 = Zygote.gradient(testf,p)
dp2 = FiniteDiff.finite_difference_gradient(testf,p)
dp3 = ForwardDiff.gradient(testf,p)
# work fine

prob = QuadratureProblem(f,lb,ub,p; batch=10)
sol = solve(prob,CubaCuhre(),reltol= 1e-4,abstol= 1e-4,maxiters=100)[1]


function testf(p)
    prob = QuadratureProblem(f,lb,ub,p, batch=10)
    solve(prob,CubaCuhre(),reltol= 1e-4,abstol= 1e-4,maxiters=100)[1]
end
dp2 = FiniteDiff.finite_difference_gradient(testf,p)
dp3 = ForwardDiff.gradient(testf,p)
# work fine

dp1 = Zygote.gradient(testf,p)

ERROR: MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::ArrayInterface.StaticInt{N}) where {T<:Number, N} at /Users/kirill/.julia/packages/ArrayInterface/rw2kK/src/static.jl:18
  convert(::Type{R}, ::T) where {R<:Real, T<:ReverseDiff.TrackedReal} at /Users/kirill/.julia/packages/ReverseDiff/jFRo1/src/tracked.jl:255
  convert(::Type{T}, ::Unitful.Quantity) where T<:Real at /Users/kirill/.julia/packages/Unitful/1t88N/src/conversion.jl:145
  ...
Stacktrace:
 [1] setindex! at ./array.jl:849 [inlined]
 [2] macro expansion at ./multidimensional.jl:802 [inlined]
 [3] macro expansion at ./cartesian.jl:64 [inlined]
 [4] macro expansion at ./multidimensional.jl:797 [inlined]
 [5] _unsafe_setindex!(::IndexLinear, ::Array{Float64,2}, ::Array{Array{Float64,1},1}, ::Base.Slice{Base.OneTo{Int64}}, ::Int64) at ./multidimensional.jl:789
 [6] _setindex! at ./multidimensional.jl:785 [inlined]
 [7] setindex!(::Array{Float64,2}, ::Array{Array{Float64,1},1}, ::Function, ::Int64) at ./abstractarray.jl:1153
 [8] (::Quadrature.var"#46#57"{QuadratureProblem{false,Array{Float64,1},typeof(f),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}})(::Array{Float64,2}, ::Array{Float64,1}) at /Users/kirill/.julia/packages/Quadrature/ZmWGy/src/Quadrature.jl:524
 [9] __solvebp_call(::QuadratureProblem{false,Array{Float64,1},Quadrature.var"#46#57"{QuadratureProblem{false,Array{Float64,1},typeof(f),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}},Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:reltol, :abstol, :maxiters),Tuple{Float64,Float64,Int64}}}}, ::CubaCuhre, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/kirill/.julia/packages/Quadrature/ZmWGy/src/Quadrature.jl:437
 [10] (::Quadrature.var"#quadrature_adjoint#52"{Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:reltol, :abstol, :maxiters),Tuple{Float64,Float64,Int64}}},QuadratureProblem{false,Array{Float64,1},typeof(f),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},CubaCuhre,Quadrature.ReCallVJP{Quadrature.ZygoteVJP},Array{Float64,1},Array{Float64,1},Array{Float64,1},Tuple{}})(::Array{Float64,1}) at /Users/kirill/.julia/packages/Quadrature/ZmWGy/src/Quadrature.jl:546
 [11] #65#back at /Users/kirill/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:55 [inlined]
 [12] #150 at /Users/kirill/.julia/packages/Zygote/c0awc/src/lib/lib.jl:191 [inlined]
 [13] (::Zygote.var"#1693#back#152"{Zygote.var"#150#151"{Quadrature.var"#65#back#64"{Quadrature.var"#quadrature_adjoint#52"{Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:reltol, :abstol, :maxiters),Tuple{Float64,Float64,Int64}}},QuadratureProblem{false,Array{Float64,1},typeof(f),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},CubaCuhre,Quadrature.ReCallVJP{Quadrature.ZygoteVJP},Array{Float64,1},Array{Float64,1},Array{Float64,1},Tuple{}}},Tuple{NTuple{8,Nothing},Tuple{}}}})(::Array{Float64,1}) at /Users/kirill/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49
 [14] #solve#10 at /Users/kirill/.julia/packages/Quadrature/ZmWGy/src/Quadrature.jl:149 [inlined]
 [15] (::typeof(∂(#solve#10)))(::Array{Float64,1}) at /Users/kirill/.julia/packages/Zygote/c0awc/src/compiler/interface2.jl:0
 [16] #150 at /Users/kirill/.julia/packages/Zygote/c0awc/src/lib/lib.jl:191 [inlined]
 [17] (::Zygote.var"#1693#back#152"{Zygote.var"#150#151"{typeof(∂(#solve#10)),Tuple{NTuple{5,Nothing},Tuple{}}}})(::Array{Float64,1}) at /Users/kirill/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49
 [18] (::typeof(∂(solve##kw)))(::Array{Float64,1}) at /Users/kirill/.julia/packages/Zygote/c0awc/src/compiler/interface2.jl:0
 [19] testf at ./none:3 [inlined]
 [20] (::typeof(∂(testf)))(::Float64) at /Users/kirill/.julia/packages/Zygote/c0awc/src/compiler/interface2.jl:0
 [21] (::Zygote.var"#41#42"{typeof(∂(testf))})(::Float64) at /Users/kirill/.julia/packages/Zygote/c0awc/src/compiler/interface.jl:45
 [22] gradient(::Function, ::Array{Float64,1}) at /Users/kirill/.julia/packages/Zygote/c0awc/src/compiler/interface.jl:54
 [23] top-level scope at none:1

Docs

Firstly, thanks for making this fantastic package! It's great to have the various quadrature / cubature packages accessible via a common interface.

The docs are currently a little harder to parse than they could be though. In particular:

The first example in the readme is kind of hard to parse because it's not explained what each of the arguments does until further down in the README. It would be helpful if there was at least quickly suggested that the first argument is the function you want to integrate, the second lower bound(s), and the third upper bound(s).
the second example is quite hard to parse, and I'm still not quite sure a) what format I can expect x to have and b) what format I can expect dx to have, as this doesn't appear to be documented. An explanation of precisely what I should assume about the structure of x, and dx (presumably as a function of x?) would be greatly appreciated.
assumptions about f when a parameter is not needed. It would be helpful to explicitly state that you still need to have an argument in f that is the parameters, assuming I've understood this correctly.
1, 2 and 3 would be greatly appreciated in docstring-form, in conjunction with the information in the QuadratureProblem section of the README.
solve doesn't appear to be properly documented anywhere, although I'll grant you that it's pretty obvious how to use it from the example. Regardless, at least the presence of a docstring would be great.
It's unclear what to do with the output type of a successful solve call. I think the u field is the thing I'm after, but I don't really know? It looks like this is defined in DiffEqBase, so maybe the docs already exist for it? If so, it would be great if they could be linked. Again, maybe a docstring for the type?

Thanks again for the package!

Inf Integrand handling

Many of the quadrature algorithms do not handle Inf in the integrand well.

QuadGK throws a DomainError, 2D or higher HCubature returns NaN. Cubature has similar issues but at different dimensions (see JuliaMath/HCubature.jl#29). Some Cuba methods seem to work, i.e. return Inf, while others segfault after a significant amount of computation time.

I'm not quite sure how to approach this issue. Seems like it should be handled upstream by the various algorithms, but this may be a big task to get cleaned up across the board.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Potential error in adjoint definition?

Solving IDEs with infinite bounds with NeuralPDE is not working. Upon further investigation (see SciML/NeuralPDE.jl#401), we realized this might be an error in the adjoint definition in Quadrature.jl.

The program outputs a LoadError: MethodError: no method matching (::Zygote.var"#pairs_namedtuple#351"{(), NamedTuple{(), Tuple{}}})(::Tuple{}) whenever the bounds are infinite. Quadrature.jl handles all of the logic of transforming the infinite bounds into finite bounds.

I will investigate this, but I am new to Quadrature.jl so if anybody has any guesses of what may be wrong here, LMK.

Setting limits of integration in solve to avoid new QuadratureProblem

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!

Incorrect integration results depending on integration bounds, using quadgk directly works

This comes from obtaining electron densities in a density of states given a chemical potential eta. Here's a MWE which shows the issue

using Integrals
using QuadGK
using Plots

function integrand(E::Real, eta::Real)
    return exp(-E^2/2) / (sqrt(2π) * (exp(E-eta) + 1))
end

function using_quadgk(eta::Real, lim::Real)
    iq(x) = integrand(x, eta)
    res, err = quadgk(iq, -lim, lim)
    return res
end

function using_integrals(eta::Real, lim::Real)
    ip = IntegralProblem(integrand, -lim, lim, eta)
    return solve(ip, QuadGKJL())[1]
end

er = LinRange(-30, 2, 200)
lim = 100
plot(er, using_quadgk.(er, lim), yscale=:log10, label="QuadGK")
plot!(er, using_integrals.(er, lim), yscale=:log10, label="Integrals")

lim = 1000
plot(er, using_quadgk.(er, lim), yscale=:log10, label="QuadGK")
plot!(er, using_integrals.(er, lim), yscale=:log10, label="Integrals")

The function should be smooth, so the calculation from Integrals.jl doesn't make sense. I've also tested other solver algorithms, but they give me also wrong results in varying degree.

Julia 1.6.6
Integrals v3.0.2

Zygote producing wrong result for ndim >1 & batch > 0

Only 1st element in the gradient is wrong and it is scaled exactly by the batchsize. NOTE: that the batchsize is determined by the internal algorithm and batch=10 is just a hint. In this following example the batchsize is 17. 136/17=8, where 8 is the correct answer

using Quadrature, Cuba, Cubature, Zygote, FiniteDiff, ForwardDiff, Test

f(x,p) = x[1,:]*p[1].+p[2]*p[3]

lb =[1.0,1.0]
ub = [3.0,3.0]
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(f,lb,ub,p)

function testf3(lb,ub,p; f=f)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=1)
    solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3)[1]
end

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3 # passes
@test_broken dp2 ≈ dp3 # Fail  [136,16,12] ≈ [9,16,12]

NNPDENS tests are failing

I am working on this PR SciML/NeuralPDE.jl#409. It is small, and doesn't touch NNPDENS at all, but about 2 days ago the CI started to fail in the NNPDENS tests with an error:

ERROR: Not implemented: convert TrackedArray{…,Vector{Float64}} to TrackedArray{…,Vector{Float32}}

You can find the whole stacktrace here https://github.com/SciML/NeuralPDE.jl/pull/409/checks?check_run_id=3817269601

Running the tests locally, I also see an error Display Error: ERROR: type TerminalPDEProblem has no field u0 coming from line 20:

https://github.com/SciML/NeuralPDE.jl/blob/23253b90eee3aeeb029782fb12cea631b90bb769/test/NNPDENS_tests.jl#L20

I am not sure why this is happening. I doubt it has anything to do with my PR, but if it does then feel free to close this issue.

EDIT: Fail, I opened this in Quadrature instead of NeuralPDE. Sorry!

Improve error when not specifying solver

Me, a common user, do not know the names of common solvers by heart AT ALL. I would therefore be tempted to try the following:

julia> using Integrals

julia> f(x, p) = sin(x)
f (generic function with 1 method)

julia> prob = IntegralProblem(f, 0, 1)
IntegralProblem. In-place: false


julia> solve(prob)
ERROR: MethodError: no method matching init(::IntegralProblem{false, SciMLBase.NullParameters, typeof(f), Int64, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
Stacktrace:
 [1] solve(args::IntegralProblem{false, SciMLBase.NullParameters, typeof(f), Int64, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ CommonSolve ~/.julia/packages/CommonSolve/TGRvG/src/CommonSolve.jl:23
 [2] solve(args::IntegralProblem{false, SciMLBase.NullParameters, typeof(f), Int64, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
   @ CommonSolve ~/.julia/packages/CommonSolve/TGRvG/src/CommonSolve.jl:23
 [3] top-level scope
   @ REPL[10]:1

That is not a very helpful error, and does not indicate the issue. I should have done the following:

julia> solve(prob, HCubatureJL())
u: 0.4596976941318603

So I would like 1 of the following 2 things:

A reasonable default integral algorithm
An improved error message that lists some common integral algorithms and demonstrates a proper function call of solve.

Batch Single dim derivative test is broken

This test is broken

https://github.com/SciML/Quadrature.jl/blob/f913f39d9766a5d15f3829c51155ab3b655269d3/test/derivative_tests.jl#L137

The problem occurs with dp2. It throws a LoadError: DimensionMismatch("variable with size(x) == (1, 15) cannot have a gradient with size(dx) == (15,)")

https://github.com/SciML/Quadrature.jl/blob/f913f39d9766a5d15f3829c51155ab3b655269d3/test/derivative_tests.jl#L133

Which causes the test to break.

For convenience, here is an MWE one can use to test.

using Quadrature, Cuba, Cubature, Zygote, FiniteDiff, ForwardDiff
using Test

### Batch Single dim
f(x,p) = x*p[1].+p[2]*p[3]

lb =1.0
ub = 3.0
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(f,lb,ub,p)

function testf3(lb,ub,p; f=f)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=1)
    solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3)[1]
end

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1] # LoadError: DimensionMismatch("variable with size(x) == (1, 15) cannot have a gradient with size(dx) == (15,)")
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3 #passes
@test dp2 ≈ dp3 # THIS IS BROKEN

Derivative with VEGAS: ERROR: AssertionError: prob.nout == 1

using Quadrature, ForwardDiff, FiniteDiff, Zygote, Cuba

f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]
prob = QuadratureProblem(f,lb,ub,p)
solve(prob,VEGAS(),reltol=1e-6,abstol=1e-6,maxiters =1000)[1]

function testf(p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,VEGAS(),reltol=1e-6,abstol=1e-6,maxiters =1000)[1])
end
dp1 = Zygote.gradient(testf,p)
#julia> ERROR: AssertionError: prob.nout == 1
dp2 = FiniteDiff.finite_difference_gradient(testf,p)
dp3 = ForwardDiff.gradient(testf,p)
#julia> ERROR: AssertionError: prob.nout == 1

Zygote tutorial broken

https://docs.sciml.ai/dev/modules/Integrals/tutorials/differentiating_integrals/

using Integrals, ForwardDiff, FiniteDiff, Zygote, IntegralsCuba
f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]

function testf(p)
    prob = IntegralProblem(f,lb,ub,p)
    sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end
dp1 = Zygote.gradient(testf,p)
dp2 = FiniteDiff.finite_difference_gradient(testf,p)
dp3 = ForwardDiff.gradient(testf,p)
dp1[1] ≈ dp2 ≈ dp3

julia> Zygote.gradient(testf,p)
ERROR: BoundsError: attempt to access 0-element Vector{Any} at index []
...

Derivative for Integration Limits Broken

Derivatives for the integration limits appear to be broken. For example, consider

julia> g(z) = solve(QuadratureProblem((x, p) -> exp(p * x), -Inf, z[1], z[2]), QuadGKJL(), reltol=1e-4).u
g (generic function with 1 method)

julia> g([1.5, 2.0])
10.042768453453952

julia> FiniteDifferences.grad(central_fdm(5, 1), g, [1.5, 2.0])[1]
2-element Array{Float64,1}:
 20.085536906907155
 10.042768652810555

Then

julia> Zygote.gradient(g, [1.5, 2.0])
ERROR: MethodError: no method matching getindex(::Float64, ::Tuple{})
Closest candidates are:
  getindex(::Number) at number.jl:75
  getindex(::Number, ::Integer) at number.jl:77
  getindex(::Number, ::Integer...) at number.jl:82
  ...
Stacktrace:
 [1] getindex(::DiffEqBase.QuadratureSolution{Float64,0,Float64,Float64,QuadratureProblem{false,Float64,var"#15#16",Float64,Float64,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},QuadGKJL,Nothing}) at /Users/wessel/.julia/packages/DiffEqBase/V7P18/src/solutions/solution_interface.jl:5
 [2] _broadcast_getindex at ./broadcast.jl:578 [inlined]

and

julia> ForwardDiff.gradient(g, [1.5, 2.0])
ERROR: StackOverflowError:
Stacktrace:
 [1] cachedrule(::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g),Float64},Float64,2}}, ::Int64) at /Users/wessel/.julia/packages/QuadGK/zzx9s/src/gausskronrod.jl:249 (repeats 79984 times)

Removing the dependence on the integration limit makes the error go away for ForwardDiff:

julia> g(z) = solve(QuadratureProblem((x, p) -> exp(p * x), -Inf, 1.5, z[2]), QuadGKJL(), reltol=1e-4).u
g (generic function with 1 method)

julia> FiniteDifferences.grad(central_fdm(5, 1), g, [1.5, 2.0])[1]
2-element Array{Float64,1}:
  0.0
 10.042768652810555

julia> ForwardDiff.gradient(g, [1.5, 2.0])
2-element Array{Float64,1}:
  0.0
 10.04276865281165

Zygote.gradient still breaks, unfortunately.

DomainIIntegrals.jl

I wasn't aware of this package previously. provides and interesting interface. https://github.com/JuliaApproximation/DomainIntegrals.jl

Quadrature with GPU

Hello,

I have a simple example which computes the mean of dirichlet distribution using quadrature. When I set ind_gpu = 1 to run on the gpu I get an error. When I set ind_gpu = 0, the code runs just fine and produces a reasonable answer for the mean of the distribution. Code and error message below:

Code:

Pkg.add(Pkg.PackageSpec(;name="Distributions", version="0.24.15"))
using Quadrature, Cuba, Cubature, Base.Threads
using Distributions, Random 
using DataFrames, CSV
using Flux, CUDA

## User Inputs
N = 3
tol = 2
ind_gpu = 1 # indicator whether to use gpu
alg = CubaDivonne() #CubatureJLh() # CubaDivonne() #works for CubaCuhre, CubaDivonne CubaSUAVE, fails for cubavega

# Setting up Variables
if ind_gpu == 1
    α0 = 10 .* (1 .- rand(N)) |> gpu
else
    α0 = 10 .* (1 .- rand(N))
end
reltol_val = 10.0^(-tol)
abstol_val = 10.0^(-tol)

# Setting up function
dist_dirichlet_pdf(x,p) = Distributions.pdf(Distributions.Dirichlet(p),x)
function f_dirichlet(dx,x,p)
    Threads.@threads for i in 1:N
        dx[i] = (dist_dirichlet_pdf([x;1.00-sum(x)],p) .* [x;1.00-sum(x)])[i]
    end
end

# Solving Integral
prob = QuadratureProblem(f_dirichlet,zeros(N-1),ones(N-1), α0, nout = N)
time_start = time()
mem_usage = @allocated sol_mean = Quadrature.solve(prob,alg,reltol=reltol_val,abstol=abstol_val)
total_mem = mem_usage/1000/2^20

# Checking Answer
mean_dirichlet(p) = p./sum(p)
display(mean_dirichlet(α0))
display(sol_mean)
'''

Error:

TaskFailedException

nested task error: InvalidIRError: compiling kernel partial_mapreduce_grid(typeof(identity), typeof(Base.add_sum), Float32, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, Val{true}, CuDeviceMatrix{Float32, 1}, Base.Broadcast.Broadcasted{CUDA.CuArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(SpecialFunctions.loggamma), Tuple{CuDeviceVector{Float32, 1}}}) resulted in invalid LLVM IR
Reason: unsupported call through a literal pointer (call to lgammaf_r)
Stacktrace:
 [1] logabsgamma
   @ ~/.julia/packages/SpecialFunctions/mFAQ4/src/gamma.jl:627
 [2] loggamma
   @ ~/.julia/packages/SpecialFunctions/mFAQ4/src/gamma.jl:670
 [3] _broadcast_getindex_evalf
   @ broadcast.jl:648
 [4] _broadcast_getindex
   @ broadcast.jl:621
 [5] getindex
   @ broadcast.jl:575
 [6] _map_getindex
   @ ~/.julia/packages/CUDA/M4jkK/src/mapreduce.jl:80
 [7] partial_mapreduce_grid
   @ ~/.julia/packages/CUDA/M4jkK/src/mapreduce.jl:117
Stacktrace:
  [1] check_ir(job::GPUCompiler.CompilerJob{GPUCompiler.PTXCompilerTarget, CUDA.CUDACompilerParams, GPUCompiler.FunctionSpec{typeof(CUDA.partial_mapreduce_grid), Tuple{typeof(identity), typeof(Base.add_sum), Float32, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, Val{true}, CuDeviceMatrix{Float32, 1}, Base.Broadcast.Broadcasted{CUDA.CuArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(SpecialFunctions.loggamma), Tuple{CuDeviceVector{Float32, 1}}}}}}, args::LLVM.Module)
    @ GPUCompiler ~/.julia/packages/GPUCompiler/XwWPj/src/validation.jl:123
  [2] macro expansion
    @ ~/.julia/packages/GPUCompiler/XwWPj/src/driver.jl:288 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/TimerOutputs/4QAIk/src/TimerOutput.jl:206 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/GPUCompiler/XwWPj/src/driver.jl:286 [inlined]
  [5] emit_asm(job::GPUCompiler.CompilerJob, ir::LLVM.Module, kernel::LLVM.Function; strip::Bool, validate::Bool, format::LLVM.API.LLVMCodeGenFileType)
    @ GPUCompiler ~/.julia/packages/GPUCompiler/XwWPj/src/utils.jl:62
  [6] cufunction_compile(job::GPUCompiler.CompilerJob)
    @ CUDA ~/.julia/packages/CUDA/M4jkK/src/compiler/execution.jl:306
  [7] check_cache
    @ ~/.julia/packages/GPUCompiler/XwWPj/src/cache.jl:44 [inlined]
  [8] cached_compilation
    @ ~/.julia/packages/CUDA/M4jkK/src/mapreduce.jl:87 [inlined]
  [9] cached_compilation(cache::Dict{UInt64, Any}, job::GPUCompiler.CompilerJob{GPUCompiler.PTXCompilerTarget, CUDA.CUDACompilerParams, GPUCompiler.FunctionSpec{typeof(CUDA.partial_mapreduce_grid), Tuple{typeof(identity), typeof(Base.add_sum), Float32, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, Val{true}, CuDeviceMatrix{Float32, 1}, Base.Broadcast.Broadcasted{CUDA.CuArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(SpecialFunctions.loggamma), Tuple{CuDeviceVector{Float32, 1}}}}}}, compiler::typeof(CUDA.cufunction_compile), linker::typeof(CUDA.cufunction_link))
    @ GPUCompiler ~/.julia/packages/GPUCompiler/XwWPj/src/cache.jl:0
 [10] cufunction(f::typeof(CUDA.partial_mapreduce_grid), tt::Type{Tuple{typeof(identity), typeof(Base.add_sum), Float32, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, Val{true}, CuDeviceMatrix{Float32, 1}, Base.Broadcast.Broadcasted{CUDA.CuArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(SpecialFunctions.loggamma), Tuple{CuDeviceVector{Float32, 1}}}}}; name::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ CUDA ~/.julia/packages/CUDA/M4jkK/src/compiler/execution.jl:294
 [11] cufunction
    @ ~/.julia/packages/CUDA/M4jkK/src/compiler/execution.jl:288 [inlined]
 [12] macro expansion
    @ ~/.julia/packages/CUDA/M4jkK/src/compiler/execution.jl:102 [inlined]
 [13] mapreducedim!(f::typeof(identity), op::typeof(Base.add_sum), R::CuArray{Float32, 1}, A::Base.Broadcast.Broadcasted{CUDA.CuArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(SpecialFunctions.loggamma), Tuple{CuArray{Float32, 1}}}; init::Float32)
    @ CUDA ~/.julia/packages/CUDA/M4jkK/src/mapreduce.jl:192
 [14] #_mapreduce#19
    @ ~/.julia/packages/GPUArrays/bjw3g/src/host/mapreduce.jl:62 [inlined]
 [15] #mapreduce#17
    @ ~/.julia/packages/GPUArrays/bjw3g/src/host/mapreduce.jl:28 [inlined]
 [16] mapreduce
    @ ~/.julia/packages/GPUArrays/bjw3g/src/host/mapreduce.jl:28 [inlined]
 [17] #_sum#682
    @ ./reducedim.jl:878 [inlined]
 [18] _sum
    @ ./reducedim.jl:878 [inlined]
 [19] #sum#680
    @ ./reducedim.jl:874 [inlined]
 [20] sum
    @ ./reducedim.jl:874 [inlined]
 [21] (Dirichlet{Float32, Ts, S} where {Ts<:AbstractVector{Float32}, S<:Real})(alpha::CuArray{Float32, 1}; check_args::Bool)
    @ Distributions ~/.julia/packages/Distributions/cNe2C/src/multivariate/dirichlet.jl:33
 [22] #Dirichlet#143
    @ ~/.julia/packages/Distributions/cNe2C/src/multivariate/dirichlet.jl:39 [inlined]
 [23] Dirichlet
    @ ~/.julia/packages/Distributions/cNe2C/src/multivariate/dirichlet.jl:39 [inlined]
 [24] dist_dirichlet_pdf(x::Vector{Float64}, p::CuArray{Float32, 1})
    @ Main ./In[26]:26
 [25] macro expansion
    @ ./In[26]:29 [inlined]
 [26] (::var"#483#threadsfor_fun#14"{Vector{Float64}, Vector{Float64}, CuArray{Float32, 1}, UnitRange{Int64}})(onethread::Bool)
    @ Main ./threadingconstructs.jl:81
 [27] (::var"#483#threadsfor_fun#14"{Vector{Float64}, Vector{Float64}, CuArray{Float32, 1}, UnitRange{Int64}})()
    @ Main ./threadingconstructs.jl:48

Stacktrace:
[1] wait
@ ./task.jl:317 [inlined]
[2] threading_run(func::Function)
@ Base.Threads ./threadingconstructs.jl:34
[3] macro expansion
@ ./threadingconstructs.jl:93 [inlined]
[4] f_dirichlet(dx::Vector{Float64}, x::Vector{Float64}, p::CuArray{Float32, 1})
@ Main ./In[26]:28
[5] (::Quadrature.var"#77#93"{QuadratureProblem{true, CuArray{Float32, 1}, typeof(f_dirichlet), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, CuArray{Float32, 1}})(x::Vector{Float64}, dx::Vector{Float64})
@ Quadrature ~/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:403
[6] generic_integrand!(ndim::Int32, x_::Ptr{Float64}, ncomp::Int32, f_::Ptr{Float64}, func!::Quadrature.var"#77#93"{QuadratureProblem{true, CuArray{Float32, 1}, typeof(f_dirichlet), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, CuArray{Float32, 1}})
@ Cuba ~/.julia/packages/Cuba/KIQTB/src/Cuba.jl:92
[7] dointegrate!
@ ~/.julia/packages/Cuba/KIQTB/src/divonne.jl:52 [inlined]
[8] dointegrate
@ ~/.julia/packages/Cuba/KIQTB/src/Cuba.jl:195 [inlined]
[9] divonne(integrand::Quadrature.var"#77#93"{QuadratureProblem{true, CuArray{Float32, 1}, typeof(f_dirichlet), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, CuArray{Float32, 1}}, ndim::Int64, ncomp::Int64; nvec::Int64, rtol::Float64, atol::Float64, flags::Int64, seed::Int64, minevals::Int64, maxevals::Int64, key1::Int64, key2::Int64, key3::Int64, maxpass::Int64, border::Float64, maxchisq::Float64, mindeviation::Float64, ngiven::Int64, ldxgiven::Int64, xgiven::Matrix{Float64}, nextra::Int64, peakfinder::Ptr{Nothing}, statefile::String, spin::Ptr{Nothing}, reltol::Missing, abstol::Missing)
@ Cuba ~/.julia/packages/Cuba/KIQTB/src/divonne.jl:145
[10] __solvebp_call(::QuadratureProblem{true, CuArray{Float32, 1}, typeof(f_dirichlet), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaDivonne, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Vector{Float64}, ::Vector{Float64}, ::CuArray{Float32, 1}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Quadrature ~/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:463
[11] #__solvebp#11
@ ~/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:153 [inlined]
[12] solve(::QuadratureProblem{true, CuArray{Float32, 1}, typeof(f_dirichlet), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaDivonne; sensealg::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
@ Quadrature ~/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:149
[13] top-level scope
@ ./timing.jl:321 [inlined]
[14] top-level scope
@ ./In[26]:0
[15] eval
@ ./boot.jl:360 [inlined]
[16] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1094
'''

Specify initial evaluation points [Feature Req.]

For DiffEqUncertainty expectation() applications, distributions that are narrow relative to the support can lead to incorrect results via Quadrature as the integrand is not sampled at points w/ non-zero (numerically) joint pdf values.

B/c we know the pdf, it would advantageous to "seed" or initialize any adaptive quadrature methods w/ the mean and/or random samples from the distribution.

E.g. for a simple linear system, u'=p*u, with uncertain IC,

u0_dist = [truncated(Normal(3.0,2.0),-1000,1000)]

will produce expectations of 0

while

u0_dist = [truncated(Normal(3.0,2.0),3-1000,3+1000)]

produces the correct result. The reason being that the midpoint of the integration domain is used as initial quadrature points in most algorithms supported.

Adaptive Quasi-Monte Carlo

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.127.8178&rep=rep1&type=pdf

I want it.

batch with scalar bounds

From the documentation:

batch: The preferred number of points to batch. This allows user-side parallelization of the integrand. If batch != 0, then each x[:,i] is a different point of the integral to calculate, and the output should be nout x batchsize. Note that batch is a suggestion for the number of points, and it is not necessarily true that batch is the same as batchsize in all algorithms.

That is not true for scalar bounds

using Quadrature
using Cubature
function f(x,p)
    x
end
function f_batched( dx_batched,x_batched,p)
    for i in 1:size(x_batched,1)
        dx_batched[i] = f(x_batched[i],p)
    end
    return nothing
end
prob = QuadratureProblem(f_batched,0.0,1.0,nout=1,batch=2)
sol = solve(prob,CubatureJLh(),reltol=1e-3,abstol=1e-3)

Infinite interval change of variables

Add automatic support for infinite and semi-infinite integral domains
https://github.com/stevengj/cubature/blob/master/README.md#infinite-intervals

Differentiation for a limit of an integral

I have a function with an integral, where the variable is the upper limit of the intgral. Using different integration methods and different automatic differentiation methods for differentiation for the upper limit leads to a StackOverflowError. If I understand the error messages right, the problem is in the different integration methods and the dual-type.

I tried to define with DiffRules.jl a rule for this case (which is simply the integrand at the upper limit), but I still got the error message. Because of the simple rule for this case, it would not be necessary to use the integration method itself and problems with the dual-type could be avoided.

eliminate need for nout parameter in most cases?

nout: The output size of the function f. Defaults to 1, i.e., a scalar integral output.

It seems like this nout parameter should not be needed for all algorithms.

In QuadGK and HCubature, you can integrate functions whose output is in any normed vector space, and the type is deteremined the first time the integrand is evaluated.

I get why you might need nout for some of the external C libraries, but it seems like you should be able to tell the user that it is not needed for certain algorithms.

FastGaussQuadrature.jl support?

I was looking though quadrature libraries and spotted https://github.com/JuliaApproximation/FastGaussQuadrature.jl

It might be nice to also hook that up to this interface.
(I don't actually have any need for this, was just speculating)

Differentiation handling of discontinuities

https://people.csail.mit.edu/sbangaru/projects/teg-2021/teg-2021.pdf

Add nout > 1 tests

QuadratureProblem is type unstable?

Julia 1.7.0 with Quadrature v1.12.0 and Cuba v2.2.0, calling @code_warntype on following function (taken from the documentation) shows that QuadratureProblem and the result of solve get typed as Any.

using Quadrature, Cuba

f(x,p) = sum(sin.(x .* p))

lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]

function testf(p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end

@code_warntype testf(p)

Variables
  #self#::Core.Const(testf)
  p::Vector{Float64}
  prob::Any

Body::Any
1 ─       (prob = Main.QuadratureProblem(Main.f, Main.lb, Main.ub, p))
│   %2  = Main.CubaCuhre()::Core.Const(CubaCuhre())
│   %3  = (:reltol, :abstol)::Core.Const((:reltol, :abstol))
│   %4  = Core.apply_type(Core.NamedTuple, %3)::Core.Const(NamedTuple{(:reltol, :abstol), T} where T<:Tuple)
│   %5  = Core.tuple(1.0e-6, 1.0e-6)::Core.Const((1.0e-6, 1.0e-6))
│   %6  = (%4)(%5)::Core.Const((reltol = 1.0e-6, abstol = 1.0e-6))
│   %7  = Core.kwfunc(Main.solve)::Core.Const(CommonSolve.var"#solve##kw"())
│   %8  = prob::Any
│   %9  = (%7)(%6, Main.solve, %8, %2)::Any
│   %10 = Base.getindex(%9, 1)::Any
│   %11 = Main.sin(%10)::Any
└──       return %11

Cuba.jl methods not dispatching properly

Using QuadGKJL() works fine, but all Cuba methods fail with.

using Quadrature, Cubature, Cuba

cost(x,p) = sin(x[1])
qprob = QuadratureProblem(cost,0,π; nout=1, batch=0)
sol = solve(qprob, QuadGKJL(), reltol=1e-3,abstol=1e-3) #OK

sol = solve(qprob, CubaVegas(), reltol=1e-3,abstol=1e-3)
sol = solve(qprob, CubaSUAVE(), reltol=1e-3,abstol=1e-3)
sol = solve(qprob, CubaDivonne(), reltol=1e-3,abstol=1e-3)
sol = solve(qprob, CubaCuhre(), reltol=1e-3,abstol=1e-3)

type QuadratureProblem has no field tspan
getproperty(::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Int64,Irrational{:π},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::Symbol) at Base.jl:33
get_concrete_tspan(::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Int64,Irrational{:π},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol,Symbol},NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}}) at solve.jl:135
get_concrete_problem(::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Int64,Irrational{:π},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol,Symbol},NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}}) at solve.jl:110
solve(::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Int64,Irrational{:π},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::CubaCuhre; kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol,Symbol},NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}}) at solve.jl:54
(::DiffEqBase.var"#solve##kw")(::NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}, ::typeof(solve), ::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Int64,Irrational{:π},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::CubaCuhre) at solve.jl:54
top-level scope at cuba_vegas.jl:11

Quadrature.jl don't support GPU with batch != 0

While there is no required algorithms that would support vectorization (batch! = 0) and work with GPU

CUDA.allowscalar(false)
chain = FastChain(FastDense(2,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,1))
initθ = initial_params(chain) |>gpu
x_ = rand(2,10)|> gpu
chain(x_, initθ)


function g(x,p) 
  x = adapt(typeof(p),x)
  sum(abs2,chain(x, p), dims=1)
end

lb =[1.0,1.0]
ub = [3.0,3.0]
p = [2.0, 3.0, 4.0] |>gpu
prob = QuadratureProblem(g,lb,ub,p)

function testf3(p, p_)
    prob = QuadratureProblem(g,lb,ub,p, batch = 100, nout=1)
    solve(prob, CubatureJLp(); reltol=1e-3,abstol=1e-3)[1]
end

testf3(p, nothing)

scalar getindex is disallowed

Stacktrace:
 [1] #17 at /root/.julia/packages/Cubature/5zwuu/src/Cubature.jl:215 [inlined]
 [2] disable_sigint(::Cubature.var"#17#18"{Bool,Bool,Int64,Float64,Float64,Int64,Int32,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Cubature.IntegrandData{Quadrature.var"#79#91"{CuArray{Float32,1}}},Ptr{Nothing}}) at ./c.jl:446
 [3] cubature(::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::Quadrature.var"#79#91"{CuArray{Float32,1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Int32) at /root/.julia/packages/Cubature/5zwuu/src/Cubature.jl:169
 [4] #pcubature_v#24 at /root/.julia/packages/Cubature/5zwuu/src/Cubature.jl:230 [inlined]
 [5] __solvebp_call(::QuadratureProblem{true,CuArray{Float32,1},typeof(g),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::CubatureJLp, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Array{Float64,1}, ::Array{Float64,1}, ::CuArray{Float32,1}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /root/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:307
 [6] #__solvebp#11 at /root/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:153 [inlined]
 [7] solve(::QuadratureProblem{true,CuArray{Float32,1},typeof(g),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::CubatureJLp; sensealg::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol,Symbol},NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}}) at /root/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:149
 [8] testf3(::CuArray{Float32,1}, ::Nothing) at ./In[58]:24
 [9] top-level scope at In[58]:27
 [10] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

The same is for the automatic differentiation. While it isn't work with GPU.

dp1 = ForwardDiff.gradient(p->testf3(p, nothing),p)
MethodError: no method matching ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3}(::Float64, ::ForwardDiff.Partials{3,Float32})
Closest candidates are:
  ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3}(::Number) where {T, V, N} at /root/.julia/packages/ForwardDiff/sqhTO/src/dual.jl:73
  ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3}(::T) where T<:Number at boot.jl:716
  ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3}(::V, ::ForwardDiff.Partials{N,V}) where {T, V, N} at /root/.julia/packages/ForwardDiff/sqhTO/src/dual.jl:17
  ...

Stacktrace:
 [1] __solvebp(::QuadratureProblem{true,CuArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3},1},typeof(g),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::CubatureJLp, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Array{Float64,1}, ::Array{Float64,1}, ::CuArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3},1}; kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol,Symbol},NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}}) at /root/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:629
 [2] solve(::QuadratureProblem{true,CuArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3},1},typeof(g),Array{Float64,1},Array{Float64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::CubatureJLp; sensealg::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol,Symbol},NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}}) at /root/.julia/packages/Quadrature/NPUfc/src/Quadrature.jl:149
 [3] testf3(::CuArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3},1}, ::Nothing) at ./In[52]:26
 [4] #37 at ./In[54]:1 [inlined]
 [5] vector_mode_dual_eval at /root/.julia/packages/ForwardDiff/sqhTO/src/apiutils.jl:37 [inlined]
 [6] vector_mode_gradient(::var"#37#38", ::CuArray{Float32,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#37#38",Float32},Float32,3,CuArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3},1}}) at /root/.julia/packages/ForwardDiff/sqhTO/src/gradient.jl:106
 [7] gradient(::Function, ::CuArray{Float32,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#37#38",Float32},Float32,3,CuArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3},1}}, ::Val{true}) at /root/.julia/packages/ForwardDiff/sqhTO/src/gradient.jl:19
 [8] gradient(::Function, ::CuArray{Float32,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#37#38",Float32},Float32,3,CuArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#37#38",Float32},Float32,3},1}}) at /root/.julia/packages/ForwardDiff/sqhTO/src/gradient.jl:17 (repeats 2 times)
 [9] top-level scope at In[54]:1
 [10] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

Single dimensional batch tests are broken

[ Info: Dimension = 1, Alg = VEGAS(100, 1000), Integrand = 1
Converged in 1000 evaluations
nevals = 1000
[ Info: Dimension = 1, Alg = CubatureJLh(), Integrand = 1
[ Info: Dimension = 1, Alg = CubatureJLp(), Integrand = 1
[ Info: Dimension = 1, Alg = CubaVegas(), Integrand = 1
[ Info: Dimension = 1, Alg = CubaSUAVE(), Integrand = 1
[ Info: Dimension = 1, Alg = VEGAS(100, 1000), Integrand = 2
Converged in 177000 evaluations
nevals = 177000
[ Info: VEGAS(100, 1000) failed
(sol.u, _sol.u) = (4.000095651237495, 60.0)
[ Info: Dimension = 1, Alg = CubatureJLh(), Integrand = 2
[ Info: Dimension = 1, Alg = CubatureJLp(), Integrand = 2
[ Info: CubatureJLp() failed
(sol.u, _sol.u) = (12.0, 60.0)
[ Info: Dimension = 1, Alg = CubaVegas(), Integrand = 2
[ Info: CubaVegas() failed
(sol.u, _sol.u) = (4.001139501949519, 60.0)
[ Info: Dimension = 1, Alg = CubaSUAVE(), Integrand = 2
[ Info: CubaSUAVE() failed
(sol.u, _sol.u) = (3.999560368540787, 60.0)
[ Info: Dimension = 1, Alg = VEGAS(100, 1000), Integrand = 3
Converged in 101000 evaluations
nevals = 101000
[ Info: VEGAS(100, 1000) failed
(sol.u, _sol.u) = (1.5295789271704923, 20.900856510613945)
[ Info: Dimension = 1, Alg = CubatureJLh(), Integrand = 3
[ Info: Dimension = 1, Alg = CubatureJLp(), Integrand = 3
[ Info: CubatureJLp() failed
(sol.u, _sol.u) = (3.7837768393868902, 20.900856510613945)
[ Info: Dimension = 1, Alg = CubaVegas(), Integrand = 3
[ Info: CubaVegas() failed
(sol.u, _sol.u) = (1.5296718426011646, 20.900856510613945)
[ Info: Dimension = 1, Alg = CubaSUAVE(), Integrand = 3
[ Info: CubaSUAVE() failed
(sol.u, _sol.u) = (1.5302221213067564, 20.900856510613945)
Test Summary:                       | Pass  Broken  Total
Batched Single Dimension Integrands |    7       8     15

It looks like MonteCarloIntegration.jl and Cuba.jl are interpreting it the same way, while CubatureJLh and CubatureJLp are doing two entirely different things. I haven't looked into who's correct here.

CubaVegas, CubaSUAVE, CubaDivonne: Julia has exited

For small values of abstol, reltol, maxiters, we get the exit error in solve for CubaVegas, CubaSUAVE, CubaDivonne.

using Quadrature, Cuba

f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]

prob = QuadratureProblem(f,lb,ub,p)
algs = [CubaVegas(), CubaSUAVE(), CubaDivonne()]

for alg in algs
    solve(prob,alg,reltol=1e-4,abstol=1e-4, maxiters =200)[1]
end
#Julia has exited.

Additional keyword arguments are not passed to `Hcubature`

Additional keyword arguments seem to be not passed to Hcubature.

In the example below I try to pass maxevals=typemax(Int) and initdiv=50, however, if used via the Quadrature.jl interface they have no effect.

# Aim: \int _infty^infty exp(rosenbrock(x)) dx
rosenbrock(x) = -400*(8*x[1]^2 - x[2] - 0.5)^2 - (4*x[1] - 1)^2

# change of variable
f(t) = exp(rosenbrock(t./(1 .- t.^2))) * prod((1 .+ t.^2)/(1 .- t.^2).^2)

# Quadrature.jl
prob = QuadratureProblem((t,p) -> f(t), -ones(2), ones(2))
int1 = solve(prob, HCubatureJL(),
             reltol=1e-9, abstol=0,
             maxevals=typemax(Int), initdiv=50)

# HCubature.jl
int2 = HCubature.hcubature(f, -ones(2), ones(2),
                           rtol=1e-9, atol=0,
                           maxevals=typemax(Int), initdiv=50)

int1.u ≈ int2[1]  # substantially different!

Without the additional keyword arguments, both interface give the same result.

Zygote w/ nout > 1 and Batch errors

Currently not working. See various commented out tests, e.g.
https://github.com/SciML/Quadrature.jl/blob/2fa4a7d10e213a1206f759310264b411b304e731/test/derivative_tests.jl#L153

test fails on Julia 1.8

using Integrals, Zygote

julia> f(x, p) = sum(sin.(p[1] * x))
f (generic function with 1 method)

julia> lb = 1.0
1.0

julia> ub = 3.0
3.0

julia> p = 2.0
2.0

julia> function testf(lb, ub, p)
           prob = IntegralProblem(f, lb, ub, p)
           sin(solve(prob, QuadGKJL(), reltol = 1e-3, abstol = 1e-3)[1])
       end
testf (generic function with 1 method)

julia> dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p)
ERROR: BoundsError: attempt to access 0-element Vector{Any} at index []
Stacktrace:
  [1] throw_boundserror(A::Vector{Any}, I::Tuple{})
    @ Base ./abstractarray.jl:703

1D example

Hello,

great repo! I am trying to add one parameter to the 1D example, see self-explanatory code below

`using Quadrature, Cuba, Cubature, Zygote, FiniteDiff, ForwardDiff
using Test

One Dimensional

f(y,p) = sum(sin.(p[1].*y) + p[2])
lb = 1.0
ub = 3.0
p = [2.0,1.0]
prob = QuadratureProblem(f,lb,ub,p)
sol = solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3)

function testf(lb,ub,p)
prob = QuadratureProblem(f,lb,ub,p)
solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3)
end

dlb1,dub1,dp1 = Zygote.gradient(testf,lb,ub,p)`

However, I find the following error at the Zygote gradient call:

ERROR: Output should be scalar; gradients are not defined for output u: 1.31184143840124581

Also, is there any way to handle additional input parameters?

Thank you in advance.

maxiters kwarg not functioning

using Quadrature
function f(x,p,c) 
    c[1]+=1
    sum(sin.(x))
end

algs = [HCubatureJL(), CubatureJLh(), CubatureJLp(), CubaSUAVE(), CubaDivonne(), CubaCuhre()]

count_array = zeros(length(algs))
for (i,alg) ∈ enumerate(algs)
    try
        count = [0]
        prob = QuadratureProblem((x,p)->f(x,p,count),ones(2),3ones(2))
        sol = solve(prob,alg,reltol=1e-3,abstol=1e-3, maxiter=3)
        count_array[i] = count[1]
    catch
        count_array[i] = NaN
    end
end
count_array #[17, 18, 82, NaN, NaN, NaN]. NaN indicates invalid kwarg

Hessian throws error

I am trying to run the example from the Readme.md file but computing hessians instead of gradients, and getting an error for both Zygote and ForwardDiff.

using Quadrature, ForwardDiff, FiniteDiff, Zygote, Cuba
f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]

function testf(p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end
    
hp1 = Zygote.hessian(testf,p) # throws error
hp2 = FiniteDiff.finite_difference_hessian(testf,p) # ok
hp3 = ForwardDiff.hessian(testf,p) # throws error

Status `~/.julia/environments/v1.6/Project.toml`
[67601950] Quadrature v1.9.0
[f6369f11] ForwardDiff v0.10.18
[e88e6eb3] Zygote v0.6.14
[8a292aeb] Cuba v2.2.0
[6a86dc24] FiniteDiff v2.8.0

HCubature() not working for vector-valued integrands

using Quadrature, Cubature, Cuba

cost(x,p) = sin(x[1])

qprob = QuadratureProblem(cost,0,π; nout=2, batch=0)
sol = solve(qprob, HCubatureJL(), reltol=1e-3,abstol=1e-3)

AssertionError: prob.nout == 1
#solve#6 at Quadrature.jl:70 [inlined]
(::DiffEqBase.var"#solve##kw")(::NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}, ::typeof(solve), ::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Int64,Irrational{:π},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::HCubatureJL) at Quadrature.jl:62
top-level scope at cuba_vegas.jl:6

Is this from a unneeded @assert or is other logic in Quadrature.jl needed to support nout > 1? Also, Quadrature.jl docs does not list this as a limitation.

From HCubature.jl Docs:

Because hcubature is written purely in Julia, the integrand f(x) can return any vector-like object (technically, any type supporting +, -, * real, and norm: a Banach space). You can integrate real, complex, and matrix-valued integrands, for example.

`transformation_if_inf` is not type stable

Updatingtransformation_if_inf is sufficient to make HCubatureJL() type stable. Other methods are still type unstable.

CubaSUAVE returning wrong result, NaN

suave() from Cuba.jl is not handling large maxiters properly. See issue.

Currently Quadtrature.jl sets by default, maxiters = typemax(Int). I recommend changing the default to the default used by Cuba.jl, maxiters = 1000000

Single dim domain bounds as arrays

Single dim algorithms error when integration bounds are single element arrays.

using Quadrature, Cubature, Cuba

cost(x,p) = sin.(x)

qprob = QuadratureProblem(cost,0,π; nout=1, batch=0)
sol = solve(qprob, QuadGKJL(), reltol=1e-3,abstol=1e-3)

qprob2 = QuadratureProblem(cost,[0],[π]; nout=1, batch=0)
sol2 = solve(qprob2, QuadGKJL(), reltol=1e-3,abstol=1e-3)

QuadGKJL only accepts one-dimensional quadrature problems.
error(::String) at error.jl:33
solve(::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Array{Int64,1},Array{Irrational{:π},1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::QuadGKJL; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at Quadrature.jl:46
(::DiffEqBase.var"#solve##kw")(::NamedTuple{(:reltol, :abstol),Tuple{Float64,Float64}}, ::typeof(solve), ::QuadratureProblem{false,DiffEqBase.NullParameters,typeof(cost),Array{Int64,1},Array{Irrational{:π},1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::QuadGKJL) at Quadrature.jl:45
top-level scope at singledim_arraybounds.jl:9

Passing data to function

Hello,

I have a function that I would like to integrate repeatedly but which depends on data whose values might change. Is there a way to specify a function in Quadrature.jl such that I don't need to create a new QuadratureProblem each time the values of the data change?

For example, is there a way to pass f(x, p, d) = sum(sin.(x * p + d)) so that I could change the value of d in an outer loop and just call solve() instead of creating a new QuadratureProblem each time?

Apologies in advance if this is obvious and there is something trivial I am overlooking.

Thanks,
Jon

Quadrature Time to Integrate as function of N and tolerance

Hi, I'm trying to find the mean of a dirichlet distribution using Quadrature. However, I'm finding that the time to complete increases drastically as a function of N. For example for N =5 it takes less than 20 seconds, but N = 10 takes over an hour. I wonder if I'm doing anything wrong with this code or if these results should be expected.

I'm curious about your intuition with respect to N, the tolerance level, and the algorithm - CubaDivonne vs CubatureJLh().

Thank you, code is below:

Pkg.add(Pkg.PackageSpec(;name="Distributions", version="0.24.15"))
using Quadrature, Cuba, Cubature, Base.Threads
using Distributions, Random
using DataFrames, CSV
using Flux, CUDA

## User Inputs
N = 7
tol = 5
ind_gpu = 0 # indicator whether to use gpu
alg = CubatureJLh() # CubaDivonne() #works for CubaCuhre, CubaDivonne CubaSUAVE, fails for cubavega

# Setting up Variables
if ind_gpu == 1
    α0 = 10 .* (1 .- rand(N)) |> gpu
else
    α0 = 10 .* (1 .- rand(N))
end
reltol_val = 10.0^(-tol)
abstol_val = 10.0^(-tol)

# Setting up function
dist_dirichlet_pdf(x,p) = Distributions.pdf(Distributions.Dirichlet(p),x)
function f_dirichlet(dx,x,p)
    Threads.@threads for i in 1:N
        dx[i] = (dist_dirichlet_pdf([x;1.00-sum(x)],p) .* [x;1.00-sum(x)])[i]
    end
end

# Solving Integral
prob = QuadratureProblem(f_dirichlet,zeros(N-1),ones(N-1), α0, nout = N)
time_start = time()
mem_usage = @allocated sol_mean = Quadrature.solve(prob,alg,reltol=reltol_val,abstol=abstol_val)
total_mem = mem_usage/1000/2^20

# Checking Answer
mean_dirichlet(p) = p./sum(p)
display(mean_dirichlet(α0))
display(sol_mean)

test_passed = sum((abs.(sol_mean .- mean_dirichlet(α0)) .< 1e-4))

## Saving Meta Data
time_end = time()
total_time = time_end-time_start
'''

Register

@JuliaRegistrator register()

sciml / integrals.jl Goto Github PK

integrals.jl's Introduction

Integrals.jl

Tutorials and Documentation

Examples

integrals.jl's People

Contributors

Stargazers

Watchers

Forkers

integrals.jl's Issues

One Dimensional

Recommend Projects

Recommend Topics

Recommend Org