From f2e5b5c9f483ea0f2cd70e36d5767ddf58fa3b27 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Thu, 22 Feb 2024 05:27:14 +0100 Subject: [PATCH] reapply formatting --- .JuliaFormatter.toml | 3 +- docs/pages.jl | 2 +- docs/src/basics/FAQ.md | 7 +- docs/src/basics/IntegralFunction.md | 2 +- docs/src/basics/SampledIntegralProblem.md | 6 +- .../tutorials/differentiating_integrals.md | 2 +- ext/IntegralsArblibExt.jl | 24 +-- ext/IntegralsCubaExt.jl | 7 +- ext/IntegralsFastGaussQuadratureExt.jl | 2 +- ext/IntegralsForwardDiffExt.jl | 4 +- ext/IntegralsMCIntegrationExt.jl | 15 +- ext/IntegralsZygoteExt.jl | 7 +- src/Integrals.jl | 20 +- src/algorithms.jl | 3 +- src/algorithms_extension.jl | 37 ++-- src/algorithms_sampled.jl | 2 - src/infinity_handling.jl | 32 ++-- src/simpsons.jl | 10 +- test/derivative_tests.jl | 179 +++++++++++------- test/inf_integral_tests.jl | 158 ++++++++-------- test/interface_tests.jl | 26 +-- test/nested_ad_tests.jl | 2 +- test/runtests.jl | 1 - test/sampled_tests.jl | 2 +- 24 files changed, 314 insertions(+), 239 deletions(-) diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 9c793591..959ad88a 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1,2 +1,3 @@ style = "sciml" -format_markdown = true \ No newline at end of file +format_markdown = true +format_docstrings = true \ No newline at end of file diff --git a/docs/pages.jl b/docs/pages.jl index d8c8a7d6..594e50ef 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -6,5 +6,5 @@ pages = ["index.md", "basics/SampledIntegralProblem.md", "basics/solve.md", "basics/FAQ.md"], - "Solvers" => Any["solvers/IntegralSolvers.md"], + "Solvers" => Any["solvers/IntegralSolvers.md"] ] diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index baaff810..4cd33aaa 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -40,8 +40,9 @@ For badly-behaved integrands, such as (nearly) singular and highly oscillatory functions, most algorithms will fail to converge and either throw an error or silently return the incorrect result. In some cases Integrals.jl can provide an error code when things go wrong, but otherwise you can always check if the error -estimate for the integral is less than the requested tolerance, e.g. `sol.resid -< max(abstol, reltol*norm(sol.u))`. Sometimes using a larger tolerance or higher +estimate for the integral is less than the requested tolerance, e.g. +`sol.resid < max(abstol, reltol*norm(sol.u))`. +Sometimes using a larger tolerance or higher precision arithmetic may help. ## How can I integrate arbitrarily-spaced data? @@ -60,4 +61,4 @@ Otherwise, feel free to open an issue or pull request. ## Can I take derivatives with respect to the limits of integration? -Currently this is not implemented. \ No newline at end of file +Currently this is not implemented. diff --git a/docs/src/basics/IntegralFunction.md b/docs/src/basics/IntegralFunction.md index eee17094..58a67dba 100644 --- a/docs/src/basics/IntegralFunction.md +++ b/docs/src/basics/IntegralFunction.md @@ -3,4 +3,4 @@ ```@docs SciMLBase.IntegralFunction SciMLBase.BatchIntegralFunction -``` \ No newline at end of file +``` diff --git a/docs/src/basics/SampledIntegralProblem.md b/docs/src/basics/SampledIntegralProblem.md index fd2aa84f..f856c68a 100644 --- a/docs/src/basics/SampledIntegralProblem.md +++ b/docs/src/basics/SampledIntegralProblem.md @@ -12,7 +12,7 @@ Say, by some means we have generated a dataset `x` and `y`: ```@example 1 using Integrals # hide f = x -> x^2 -x = range(0, 1, length=20) +x = range(0, 1, length = 20) y = f.(x) ``` @@ -63,9 +63,9 @@ using Integrals # hide f1 = x -> x^2 f2 = x -> x^3 f3 = x -> x^4 -x = range(0, 1, length=20) +x = range(0, 1, length = 20) y = [f1.(x) f2.(x) f3.(x)] -problem = SampledIntegralProblem(y, x; dim=1) +problem = SampledIntegralProblem(y, x; dim = 1) method = SimpsonsRule() solve(problem, method) ``` diff --git a/docs/src/tutorials/differentiating_integrals.md b/docs/src/tutorials/differentiating_integrals.md index 763f8a09..e2c4b5ed 100644 --- a/docs/src/tutorials/differentiating_integrals.md +++ b/docs/src/tutorials/differentiating_integrals.md @@ -16,7 +16,7 @@ function testf(p) sin(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1]) end testf(p) -dp1 = Zygote.gradient(testf,p) +dp1 = Zygote.gradient(testf, p) dp2 = FiniteDiff.finite_difference_gradient(testf, p) dp3 = ForwardDiff.gradient(testf, p) dp1[1] ≈ dp2 ≈ dp3 diff --git a/ext/IntegralsArblibExt.jl b/ext/IntegralsArblibExt.jl index 37811bd7..e0c408a8 100644 --- a/ext/IntegralsArblibExt.jl +++ b/ext/IntegralsArblibExt.jl @@ -3,9 +3,9 @@ module IntegralsArblibExt using Arblib using Integrals -function Integrals.__solvebp_call(prob::IntegralProblem, alg::ArblibJL, sensealg, domain, p; - reltol = 1e-8, abstol = 1e-8, maxiters = nothing) - +function Integrals.__solvebp_call( + prob::IntegralProblem, alg::ArblibJL, sensealg, domain, p; + reltol = 1e-8, abstol = 1e-8, maxiters = nothing) lb_, ub_ = domain lb, ub = map(first, domain) if !isone(length(lb_)) || !isone(length(ub_)) @@ -17,16 +17,18 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::ArblibJL, sensealg res = Acb(0) y_ = similar(prob.f.integrand_prototype, typeof(res)) f_ = (y, x; kws...) -> (prob.f(y_, x, p; kws...); Arblib.set!(y, only(y_))) - val = Arblib.integrate!(f_, res, lb, ub, atol=abstol, rtol=reltol, - check_analytic=alg.check_analytic, take_prec=alg.take_prec, - warn_on_no_convergence=alg.warn_on_no_convergence, opts=alg.opts) - SciMLBase.build_solution(prob, alg, val, get_radius(val), retcode = ReturnCode.Success) + val = Arblib.integrate!(f_, res, lb, ub, atol = abstol, rtol = reltol, + check_analytic = alg.check_analytic, take_prec = alg.take_prec, + warn_on_no_convergence = alg.warn_on_no_convergence, opts = alg.opts) + SciMLBase.build_solution( + prob, alg, val, get_radius(val), retcode = ReturnCode.Success) else f_ = (x; kws...) -> only(prob.f(x, p; kws...)) - val = Arblib.integrate(f_, lb, ub, atol=abstol, rtol=reltol, - check_analytic=alg.check_analytic, take_prec=alg.take_prec, - warn_on_no_convergence=alg.warn_on_no_convergence, opts=alg.opts) - SciMLBase.build_solution(prob, alg, val, get_radius(val), retcode = ReturnCode.Success) + val = Arblib.integrate(f_, lb, ub, atol = abstol, rtol = reltol, + check_analytic = alg.check_analytic, take_prec = alg.take_prec, + warn_on_no_convergence = alg.warn_on_no_convergence, opts = alg.opts) + SciMLBase.build_solution( + prob, alg, val, get_radius(val), retcode = ReturnCode.Success) end end diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl index e391c650..f5c16977 100644 --- a/ext/IntegralsCubaExt.jl +++ b/ext/IntegralsCubaExt.jl @@ -2,8 +2,8 @@ module IntegralsCubaExt using Integrals, Cuba import Integrals: transformation_if_inf, - scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm, - CubaSUAVE, CubaDivonne, CubaCuhre + scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm, + CubaSUAVE, CubaDivonne, CubaCuhre function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm, sensealg, @@ -22,7 +22,8 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori nvec = min(maxiters, prob.f.max_batch) # nvec == 1 in Cuba will change vectors to matrices, so we won't support it when # batching - nvec > 1 || throw(ArgumentError("BatchIntegralFunction must take multiple batch points")) + nvec > 1 || + throw(ArgumentError("BatchIntegralFunction must take multiple batch points")) if mid isa Real _x = zeros(typeof(mid), nvec) diff --git a/ext/IntegralsFastGaussQuadratureExt.jl b/ext/IntegralsFastGaussQuadratureExt.jl index f3868bf0..1dc5c398 100644 --- a/ext/IntegralsFastGaussQuadratureExt.jl +++ b/ext/IntegralsFastGaussQuadratureExt.jl @@ -33,7 +33,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::Integrals.GaussLeg sensealg, domain, p; reltol = nothing, abstol = nothing, maxiters = nothing) where {C} - if !all(isone∘length, domain) + if !all(isone ∘ length, domain) error("GaussLegendre only accepts one-dimensional quadrature problems.") end @assert prob.f isa IntegralFunction diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 300b6063..49abe0b5 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -22,8 +22,8 @@ end # Manually split for the pushforward function Integrals.__solvebp(cache, alg, sensealg, domain, - p::Union{D,AbstractArray{<:D}}; - kwargs...) where {T, V, P, D<:ForwardDiff.Dual{T, V, P}} + p::Union{D, AbstractArray{<:D}}; + kwargs...) where {T, V, P, D <: ForwardDiff.Dual{T, V, P}} # we need the output type to avoid perturbation confusion while unwrapping nested duals # We compute a vector-valued integral of the primal and dual simultaneously diff --git a/ext/IntegralsMCIntegrationExt.jl b/ext/IntegralsMCIntegrationExt.jl index fb9a68cf..e4d67792 100644 --- a/ext/IntegralsMCIntegrationExt.jl +++ b/ext/IntegralsMCIntegrationExt.jl @@ -3,10 +3,10 @@ module IntegralsMCIntegrationExt using MCIntegration, Integrals function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, domain, p; - reltol = nothing, abstol = nothing, maxiters = 10) + reltol = nothing, abstol = nothing, maxiters = 10) lb, ub = domain mid = vec(collect((lb + ub) / 2)) - var = Continuous(vec([tuple(a,b) for (a,b) in zip(lb, ub)])) + var = Continuous(vec([tuple(a, b) for (a, b) in zip(lb, ub)])) if prob.f isa BatchIntegralFunction error("VEGASMC doesn't support batching. See https://github.com/numericalEFT/MCIntegration.jl/issues/29") @@ -16,7 +16,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, f_ = (x, f, c) -> begin n = 0 for v in x - mid[n+=1] = first(v) + mid[n += 1] = first(v) end prob.f(f0, mid, p) f .= vec(f0) @@ -26,21 +26,22 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, f_ = (x, c) -> begin n = 0 for v in x - mid[n+=1] = first(v) + mid[n += 1] = first(v) end fx = prob.f(mid, p) fx isa AbstractArray ? vec(fx) : fx end end dof = ones(Int, length(f0)) # each composite Continuous var gets 1 dof - res = integrate(f_; var, dof, inplace=isinplace(prob), type=eltype(f0), - solver=:vegasmc, niter=maxiters, verbose=-2, print=-2, alg.kws...) + res = integrate(f_; var, dof, inplace = isinplace(prob), type = eltype(f0), + solver = :vegasmc, niter = maxiters, verbose = -2, print = -2, alg.kws...) out, err, chi = if f0 isa Number map(only, (res.mean, res.stdev, res.chi2)) else map(a -> reshape(a, size(f0)), (res.mean, res.stdev, res.chi2)) end - SciMLBase.build_solution(prob, alg, out, err, chi=chi, retcode = ReturnCode.Success) + SciMLBase.build_solution( + prob, alg, out, err, chi = chi, retcode = ReturnCode.Success) end end diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index 90d6f882..1b389565 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -28,7 +28,8 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal # zygote doesn't support mutation, so we build an oop pullback if sensealg.vjp isa Integrals.ZygoteVJP if cache.f isa BatchIntegralFunction - dx = similar(cache.f.integrand_prototype, size(cache.f.integrand_prototype)[begin:end-1]..., 1) + dx = similar(cache.f.integrand_prototype, + size(cache.f.integrand_prototype)[begin:(end - 1)]..., 1) _f = x -> (cache.f(dx, x, p); dx) # TODO: let the user pass a batched jacobian so we can return a BatchIntegralFunction dfdp_ = function (x, p) @@ -104,7 +105,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal # TODO replace evaluation at endpoint (which anyone can do without Integrals.jl) # with integration of dfdx uing the same quadrature dlb = cache.f isa BatchIntegralFunction ? -batch_unwrap(_f([lb])) : -_f(lb) - dub = cache.f isa BatchIntegralFunction ? batch_unwrap(_f([ub])) : _f(ub) + dub = cache.f isa BatchIntegralFunction ? batch_unwrap(_f([ub])) : _f(ub) return (NoTangent(), NoTangent(), NoTangent(), @@ -128,7 +129,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal out, quadrature_adjoint end -batch_unwrap(x::AbstractArray) = dropdims(x; dims=ndims(x)) +batch_unwrap(x::AbstractArray) = dropdims(x; dims = ndims(x)) Zygote.@adjoint function Zygote.literal_getproperty(sol::SciMLBase.IntegralSolution, ::Val{:u}) diff --git a/src/Integrals.jl b/src/Integrals.jl index c5038278..91d47720 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -66,7 +66,6 @@ end # Give a layer to intercept with AD __solvebp(args...; kwargs...) = __solvebp_call(args...; kwargs...) - function quadgk_prob_types(f, lb::T, ub::T, p, nrm) where {T} DT = float(T) # we need to be careful to infer the same result as `evalrule` RT = Base.promote_op(*, DT, Base.promote_op(f, DT, typeof(p))) # kernel @@ -106,10 +105,10 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p BatchIntegrand((y, x) -> prob.f(y, x, p), similar(bu)) else fsize = size(bu)[begin:(end - 1)] - BatchIntegrand{Array{eltype(bu),ndims(bu)-1}}() do y, x + BatchIntegrand{Array{eltype(bu), ndims(bu) - 1}}() do y, x y_ = similar(bu, fsize..., length(y)) prob.f(y_, x, p) - map!(collect, y, eachslice(y_; dims=ndims(bu))) + map!(collect, y, eachslice(y_; dims = ndims(bu))) return nothing end end @@ -121,8 +120,8 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p f = if u isa AbstractVector BatchIntegrand((y, x) -> y .= prob.f(x, p), u) else - BatchIntegrand{Array{eltype(u),ndims(u)-1}}() do y, x - map!(collect, y, eachslice(prob.f(x, p); dims=ndims(u))) + BatchIntegrand{Array{eltype(u), ndims(u) - 1}}() do y, x + map!(collect, y, eachslice(prob.f(x, p); dims = ndims(u))) return nothing end end @@ -134,7 +133,8 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p if isinplace(prob) result = prob.f.integrand_prototype * mid # result may have different units than prototype f = (y, x) -> prob.f(y, x, p) - val, err = quadgk!(f, result, lb, ub, segbuf = cache.cacheval, maxevals = maxiters, + val, err = quadgk!( + f, result, lb, ub, segbuf = cache.cacheval, maxevals = maxiters, rtol = reltol, atol = abstol, order = alg.order, norm = alg.norm) SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) else @@ -189,8 +189,9 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p; # This is an ugly hack that is compatible with both f = x -> (prob.f(y, eltype(x) <: SubArray ? parent(first(x)) : x', p); vec(y)) else - y = prob.f(mid isa Number ? typeof(mid)[] : - Matrix{eltype(mid)}(undef, length(mid), 0), + y = prob.f( + mid isa Number ? typeof(mid)[] : + Matrix{eltype(mid)}(undef, length(mid), 0), p) f = x -> prob.f(eltype(x) <: SubArray ? parent(first(x)) : x', p) end @@ -221,7 +222,8 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p; SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success) end -export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule, TrapezoidalRule, SimpsonsRule +export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule, + TrapezoidalRule, SimpsonsRule export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre export CubatureJLh, CubatureJLp export ArblibJL diff --git a/src/algorithms.jl b/src/algorithms.jl index 9f1bef72..4222b994 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -127,9 +127,8 @@ function GaussLegendre(; n = 250, subintervals = 1, nodes = nothing, weights = n return GaussLegendre(nodes, weights, subintervals) end - """ - QuadratureRule(q; n=250) +QuadratureRule(q; n=250) Algorithm to construct and evaluate a quadrature rule `q` of `n` points computed from the inputs as `x, w = q(n)`. It assumes the nodes and weights are for the standard interval diff --git a/src/algorithms_extension.jl b/src/algorithms_extension.jl index b79eaa8d..995f712f 100644 --- a/src/algorithms_extension.jl +++ b/src/algorithms_extension.jl @@ -123,13 +123,15 @@ end function CubaVegas(; flags = 0, seed = 0, minevals = 0, nstart = 1000, nincrease = 500, gridno = 0) - isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaVegas requires `using Cuba`") + isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && + error("CubaVegas requires `using Cuba`") return CubaVegas(flags, seed, minevals, nstart, nincrease, gridno) end function CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000, nmin = 2, flatness = 25.0) - isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaSUAVE requires `using Cuba`") + isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && + error("CubaSUAVE requires `using Cuba`") return CubaSUAVE(flags, seed, minevals, nnew, nmin, flatness) end @@ -138,13 +140,15 @@ function CubaDivonne(; flags = 0, seed = 0, minevals = 0, maxchisq = 10.0, mindeviation = 0.25, xgiven = zeros(Cdouble, 0, 0), nextra = 0, peakfinder = C_NULL) - isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaDivonne requires `using Cuba`") + isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && + error("CubaDivonne requires `using Cuba`") return CubaDivonne(flags, seed, minevals, key1, key2, key3, maxpass, border, maxchisq, mindeviation, xgiven, nextra, peakfinder) end function CubaCuhre(; flags = 0, minevals = 0, key = 0) - isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaCuhre requires `using Cuba`") + isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && + error("CubaCuhre requires `using Cuba`") return CubaCuhre(flags, minevals, key) end @@ -174,8 +178,9 @@ publisher={Elsevier} struct CubatureJLh <: AbstractCubatureJLAlgorithm error_norm::Int32 end -function CubatureJLh(; error_norm=0) - isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && error("CubatureJLh requires `using Cubature`") +function CubatureJLh(; error_norm = 0) + isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && + error("CubatureJLh requires `using Cubature`") return CubatureJLh(error_norm) end @@ -193,13 +198,12 @@ Defaults to `Cubature.INDIVIDUAL`, other options are struct CubatureJLp <: AbstractCubatureJLAlgorithm error_norm::Int32 end -function CubatureJLp(; error_norm=0) - isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && error("CubatureJLp requires `using Cubature`") +function CubatureJLp(; error_norm = 0) + isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && + error("CubatureJLp requires `using Cubature`") return CubatureJLp(error_norm) end - - """ ArblibJL(; check_analytic=false, take_prec=false, warn_on_no_convergence=false, opts=C_NULL) @@ -221,8 +225,10 @@ struct ArblibJL{O} <: AbstractIntegralExtensionAlgorithm warn_on_no_convergence::Bool opts::O end -function ArblibJL(; check_analytic=false, take_prec=false, warn_on_no_convergence=false, opts=C_NULL) - isnothing(Base.get_extension(@__MODULE__, :IntegralsArblibExt)) && error("ArblibJL requires `using Arblib`") +function ArblibJL(; check_analytic = false, take_prec = false, + warn_on_no_convergence = false, opts = C_NULL) + isnothing(Base.get_extension(@__MODULE__, :IntegralsArblibExt)) && + error("ArblibJL requires `using Arblib`") return ArblibJL(check_analytic, take_prec, warn_on_no_convergence, opts) end @@ -232,14 +238,15 @@ end Markov-chain based Vegas algorithm from MCIntegration.jl Refer to -[`MCIntegration.integrate`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/#MCIntegration.integrate-Tuple{Function}) +[`MCIntegration.integrate`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/#MCIntegration.integrate-Tuple%7BFunction%7D) for documentation on the keywords, which are passed directly to the solver with a set of defaults that works for conforming integrands. """ -struct VEGASMC{K<:NamedTuple} <: AbstractIntegralExtensionAlgorithm +struct VEGASMC{K <: NamedTuple} <: AbstractIntegralExtensionAlgorithm kws::K end function VEGASMC(; kws...) - isnothing(Base.get_extension(@__MODULE__, :IntegralsMCIntegrationExt)) && error("VEGASMC requires `using MCIntegration`") + isnothing(Base.get_extension(@__MODULE__, :IntegralsMCIntegrationExt)) && + error("VEGASMC requires `using MCIntegration`") return VEGASMC(NamedTuple(kws)) end diff --git a/src/algorithms_sampled.jl b/src/algorithms_sampled.jl index 6cd5dd57..9e71237e 100644 --- a/src/algorithms_sampled.jl +++ b/src/algorithms_sampled.jl @@ -5,7 +5,6 @@ abstract type AbstractSampledIntegralAlgorithm <: SciMLBase.AbstractIntegralAlgo Struct for evaluating an integral via the trapezoidal rule. - Example with sampled data: ``` @@ -27,7 +26,6 @@ Struct for evaluating an integral via the Simpson's composite 1/3-3/8 rule over `AbstractRange`s (evenly spaced points) and Simpson's composite 1/3 rule for non-equidistant grids. - Example with equidistant data: ``` diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl index ef03f69b..0cf5c54b 100644 --- a/src/infinity_handling.jl +++ b/src/infinity_handling.jl @@ -37,7 +37,7 @@ function substitute_t(t::Number, lb::Number, ub::Number) end end function substitute_t(t::AbstractVector, lb::AbstractVector, ub::AbstractVector) - x = similar(t, typeof(one(eltype(t))*(first(lb)+first(ub)))) + x = similar(t, typeof(one(eltype(t)) * (first(lb) + first(ub)))) jac = one(eltype(t)) for i in eachindex(lb) x[i], dj = substitute_t(t[i], lb[i], ub[i]) @@ -58,7 +58,7 @@ function substitute_f(f::F, dt, t, p, lb, ub) where {F} end function substitute_t(t::AbstractVector, lb::Number, ub::Number) - x = similar(t, typeof(one(eltype(t))*(lb+ub))) + x = similar(t, typeof(one(eltype(t)) * (lb + ub))) jac = similar(x) for (i, ti) in enumerate(t) x[i], jac[i] = substitute_t(ti, lb, ub) @@ -66,10 +66,11 @@ function substitute_t(t::AbstractVector, lb::Number, ub::Number) return x, jac end function substitute_t(t::AbstractArray, lb::AbstractVector, ub::AbstractVector) - x = similar(t, typeof(one(eltype(t))*(first(lb)+first(ub)))) + x = similar(t, typeof(one(eltype(t)) * (first(lb) + first(ub)))) jac = similar(x, size(t, ndims(t))) for (i, it) in enumerate(axes(t)[end]) - x[axes(x)[begin:end-1]..., i], jac[i] = substitute_t(t[axes(t)[begin:end-1]..., it], lb, ub) + x[axes(x)[begin:(end - 1)]..., i], jac[i] = substitute_t( + t[axes(t)[begin:(end - 1)]..., it], lb, ub) end return x, jac end @@ -77,13 +78,13 @@ end function substitute_batchf(f::F, t, p, lb, ub) where {F} x, jac = substitute_t(t, lb, ub) r = f(x, p) - return r .* reshape(jac, ntuple(d -> d==ndims(r) ? length(jac) : 1, ndims(r))) + return r .* reshape(jac, ntuple(d -> d == ndims(r) ? length(jac) : 1, ndims(r))) end function substitute_batchf(f::F, dt, t, p, lb, ub) where {F} x, jac = substitute_t(t, lb, ub) f(dt, x, p) for (i, j) in zip(axes(dt)[end], jac) - for idt in CartesianIndices(axes(dt)[begin:end-1]) + for idt in CartesianIndices(axes(dt)[begin:(end - 1)]) dt[idt, i] *= j end end @@ -95,22 +96,31 @@ function transformation_if_inf(prob, ::Val{true}) f = prob.f bounds = map(substitute_bounds, lb, ub) lb_sub = lb isa Number ? first(bounds) : map(first, bounds) - ub_sub = ub isa Number ? last(bounds) : map(last, bounds) + ub_sub = ub isa Number ? last(bounds) : map(last, bounds) f_sub = if isinplace(prob) if f isa BatchIntegralFunction - BatchIntegralFunction{true}(let f=f.f; (dt, t, p) -> substitute_batchf(f, dt, t, p, lb, ub); end, + BatchIntegralFunction{true}( + let f = f.f + (dt, t, p) -> substitute_batchf(f, dt, t, p, lb, ub) + end, f.integrand_prototype, max_batch = f.max_batch) else - IntegralFunction{true}(let f=f.f; (dt, t, p) -> substitute_f(f, dt, t, p, lb, ub); end, + IntegralFunction{true}(let f = f.f + (dt, t, p) -> substitute_f(f, dt, t, p, lb, ub) + end, f.integrand_prototype) end else if f isa BatchIntegralFunction - BatchIntegralFunction{false}(let f=f.f; (t, p) -> substitute_batchf(f, t, p, lb, ub); end, + BatchIntegralFunction{false}(let f = f.f + (t, p) -> substitute_batchf(f, t, p, lb, ub) + end, f.integrand_prototype) else - IntegralFunction{false}(let f=f.f; (t, p) -> substitute_f(f, t, p, lb, ub); end, + IntegralFunction{false}(let f = f.f + (t, p) -> substitute_f(f, t, p, lb, ub) + end, f.integrand_prototype) end end diff --git a/src/simpsons.jl b/src/simpsons.jl index 7c45b0f0..06b5fb7f 100644 --- a/src/simpsons.jl +++ b/src/simpsons.jl @@ -49,11 +49,13 @@ end isodd(j) && return (x[begin + j + 1] - x[begin + j - 1])^3 / (x[begin + j] - x[begin + j - 1]) / (x[begin + j + 1] - x[begin + j]) / 6 - iseven(j) && + iseven(j) && return (x[begin + j] - x[begin + j - 2]) / 6 * - (2 - (x[begin + j - 1] - x[begin + j - 2]) / (x[begin + j] - x[begin + j - 1])) + - (x[begin + j + 2] - x[begin + j]) / 6 * - (2 - (x[begin + j + 2] - x[begin + j + 1]) / (x[begin + j + 1] - x[begin + j])) + (2 - + (x[begin + j - 1] - x[begin + j - 2]) / (x[begin + j] - x[begin + j - 1])) + + (x[begin + j + 2] - x[begin + j]) / 6 * + (2 - + (x[begin + j + 2] - x[begin + j + 1]) / (x[begin + j + 1] - x[begin + j])) end function find_weights(x::AbstractVector, ::SimpsonsRule) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index ffd9e81f..ba22e8f1 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -10,39 +10,39 @@ reltol = 1e-3 abstol = 1e-3 alg_req = Dict( - QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, + QuadratureRule(gausslegendre, n = 50) => ( + nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, allows_iip = false), - GaussLegendre(n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, + GaussLegendre(n = 50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, allows_iip = false), QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1, allows_iip = true), HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1, - max_dim = Inf, allows_iip = true), - # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, - # allows_iip = true), - CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, + max_dim = Inf, allows_iip = true), CubatureJLh() => ( + nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, allows_iip = true), CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1, - max_dim = Inf, allows_iip = true), - # CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, - # allows_iip = true), - # CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, - # allows_iip = true), - # CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2, - # max_dim = Inf, allows_iip = true), - # CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, - # allows_iip = true), -) + max_dim = Inf, allows_iip = true)) +# VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, +# allows_iip = true), +# CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, +# allows_iip = true), +# CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, +# allows_iip = true), +# CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2, +# max_dim = Inf, allows_iip = true), +# CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, +# allows_iip = true), # integrands should have same shape as parameters, independent of dimensionality -integrands =( - (x, p) -> map(q -> prod(y -> sin(y*q), x), p), +integrands = ( + (x, p) -> map(q -> prod(y -> sin(y * q), x), p), ) # function to turn the output into a scalar / test different tangent types scalarize_solution = ( sol -> sin(sum(sol)), - sol -> sin(sol[1]), + sol -> sin(sol[1]) ) # we will be able to use broadcasting for this after https://github.com/FluxML/Zygote.jl/pull/1488 @@ -58,7 +58,7 @@ function f_helper!(f, y, x, p) end # the Zygote implementation is inconsistent about 0-d so we hijack it -struct Scalar{T<:Real} <: Real +struct Scalar{T <: Real} <: Real x::T end Base.iterate(a::Scalar) = (a.x, nothing) @@ -67,13 +67,13 @@ Base.IteratorSize(::Type{Scalar{T}}) where {T} = Base.HasShape{0}() Base.eltype(::Type{Scalar{T}}) where {T} = T Base.length(a::Scalar) = 1 Base.size(::Scalar) = () -Base.:+(a::Scalar, b::Scalar) = Scalar(a.x+b.x) -Base.:*(a::Number, b::Scalar) = a*b.x -Base.:*(a::Scalar, b::Number) = a.x*b -Base.:*(a::Scalar, b::Scalar) = Scalar(a.x*b.x) +Base.:+(a::Scalar, b::Scalar) = Scalar(a.x + b.x) +Base.:*(a::Number, b::Scalar) = a * b.x +Base.:*(a::Scalar, b::Number) = a.x * b +Base.:*(a::Scalar, b::Scalar) = Scalar(a.x * b.x) Base.zero(a::Scalar) = Scalar(zero(a.x)) Base.map(f, a::Scalar) = map(f, a.x) -(::Type{T})(a::Scalar) where {T<:Real} = T(a.x) +(::Type{T})(a::Scalar) where {T <: Real} = T(a.x) struct ScalarAxes end # the implementation doesn't preserve singleton axes Base.axes(::Scalar) = ScalarAxes() Base.iterate(::ScalarAxes) = nothing @@ -85,7 +85,7 @@ Base.reshape(A::AbstractArray, ::ScalarAxes) = Scalar(only(A)) # p can't be either without both prs function batch_helper(f, x, p) t = f(zero(eltype(x)), zero(eltype(eltype(p)))) - typeof(t).([f(y, q) for q in p, y in eachslice(x; dims=ndims(x))]) + typeof(t).([f(y, q) for q in p, y in eachslice(x; dims = ndims(x))]) end function batch_helper!(f, y, x, p) @@ -101,7 +101,8 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) end testf(lb, ub, p) - dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p isa Number && f isa BatchIntegralFunction ? Scalar(p) : p) + dlb1, dub1, dp1 = Zygote.gradient( + testf, lb, ub, p isa Number && f isa BatchIntegralFunction ? Scalar(p) : p) f_lb = lb -> testf(lb, ub, p) f_ub = ub -> testf(lb, ub, p) @@ -113,16 +114,16 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) dub2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dub))(f_ub, ub) if lb isa Number - @test dlb1 ≈ dlb2 atol=abstol rtol=reltol - @test dub1 ≈ dub2 atol=abstol rtol=reltol + @test dlb1≈dlb2 atol=abstol rtol=reltol + @test dub1≈dub2 atol=abstol rtol=reltol else # TODO: implement multivariate limit derivatives in ZygoteExt - @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol - @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol + @test_broken dlb1≈dlb2 atol=abstol rtol=reltol + @test_broken dub1≈dub2 atol=abstol rtol=reltol end # TODO: implement limit derivatives in ForwardDiffExt - @test_broken dlb2 ≈ getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol - @test_broken dub2 ≈ getproperty(ForwardDiff, dub)(dfdub, ub) atol=abstol rtol=reltol + @test_broken dlb2≈getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol + @test_broken dub2≈getproperty(ForwardDiff, dub)(dfdub, ub) atol=abstol rtol=reltol f_p = p -> testf(lb, ub, p) @@ -131,15 +132,16 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) dp2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dp))(f_p, p) dp3 = getproperty(ForwardDiff, dp)(f_p, p) - @test dp1 ≈ dp2 atol=abstol rtol=reltol - @test dp2 ≈ dp3 atol=abstol rtol=reltol + @test dp1≈dp2 atol=abstol rtol=reltol + @test dp2≈dp3 atol=abstol rtol=reltol return end - ### One Dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution) + req.nout > 1 || continue req.min_dim <= 1 || continue @@ -148,16 +150,21 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize end ## One-dimensional nout -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.nout > 1 || continue req.min_dim <= 1 || continue @info "One-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout - do_tests(; f, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; + f, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### N-dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @@ -166,104 +173,134 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize end ### N-dimensional nout -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, + nout in 1:max_nout_test + req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @info "Multi-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout - do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim), + p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### One Dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution) + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= 1 || continue @info "One-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1)) - do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + do_tests(; f = fiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## One-dimensional nout -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= 1 || continue @info "One-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout)) - do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f = fiip, scalarize, lb = 1.0, ub = 3.0, + p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### N-dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @info "Multi-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1)) - do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + do_tests(; + f = fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### N-dimensional nout iip -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, + nout in 1:max_nout_test + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @info "Multi-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout)) - do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f = fiip, scalarize, lb = ones(dim), ub = 3ones(dim), + p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, One Dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution) + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= 1 || continue @info "Batched, one-dimensional, scalar, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) - do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + do_tests(; f = bf, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## Batch, One-dimensional nout -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= 1 || continue @info "Batched, one-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) - do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f = bf, scalarize, lb = 1.0, ub = 3.0, + p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, N-dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @info "Batched, multi-dimensional, scalar, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) - do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + do_tests(; + f = bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### Batch, N-dimensional nout -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, + nout in 1:max_nout_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @info "Batch, multi-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) - do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f = bf, scalarize, lb = ones(dim), ub = 3ones(dim), + p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, one-dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution) + req.allows_batch || continue req.allows_iip || continue req.nout > 1 || continue @@ -271,11 +308,13 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize @info "Batched, one-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0)) - do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + do_tests(; f = bfiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## Batch, one-dimensional nout -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.allows_batch || continue req.allows_iip || continue req.nout > 1 || continue @@ -283,11 +322,14 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize @info "Batched, one-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0)) - do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f = bfiip, scalarize, lb = 1.0, ub = 3.0, + p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, N-dimensional -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.allows_batch || continue req.allows_iip || continue req.nout > 1 || continue @@ -295,11 +337,15 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize @info "Batched, multi-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0)) - do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + do_tests(; f = bfiip, scalarize, lb = ones(dim), + ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### Batch, N-dimensional nout iip -for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), + (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, + nout in 1:max_nout_test + req.allows_batch || continue req.allows_iip || continue req.nout > 1 || continue @@ -307,5 +353,6 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize @info "Batched, multi-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0)) - do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f = bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), + p = [2.0i for i in 1:nout], alg, abstol, reltol) end diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl index c93ea51f..9164fd4d 100644 --- a/test/inf_integral_tests.jl +++ b/test/inf_integral_tests.jl @@ -5,116 +5,116 @@ abstol = 1e-3 # not all quadratures are compatible with infinities if they evaluate the endpoints alg_req = Dict( - QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, - allows_iip = false, allows_inf=true), - # GaussLegendre(n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, - # allows_iip = false, allows_inf=true), - QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1, - allows_iip = true, allows_inf=true), + QuadratureRule(gausslegendre, n = 50) => ( + nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, + allows_iip = false, allows_inf = true), QuadGKJL() => ( + nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1, + allows_iip = true, allows_inf = true), HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1, - max_dim = Inf, allows_iip = true, allows_inf=true), - # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, - # allows_iip = true), - CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, - max_dim = Inf, allows_iip = true, allows_inf=true), - # CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1, - # max_dim = Inf, allows_iip = true, allows_inf=false), - # CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, - # allows_iip = true), - # CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, - # allows_iip = true), - # CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2, - # max_dim = Inf, allows_iip = true), - # CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, - # allows_iip = true), + max_dim = Inf, allows_iip = true, allows_inf = true), CubatureJLh() => ( + nout = Inf, allows_batch = true, min_dim = 1, + max_dim = Inf, allows_iip = true, allows_inf = true) ) - +# GaussLegendre(n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, +# allows_iip = false, allows_inf=true), +# VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, +# allows_iip = true), +# CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1, +# max_dim = Inf, allows_iip = true, allows_inf=false), +# CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, +# allows_iip = true), +# CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, +# allows_iip = true), +# CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2, +# max_dim = Inf, allows_iip = true), +# CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, +# allows_iip = true), problems = ( (; # 1. multi-variate infinite limits: Gaussian - f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), - domain=(@SVector[-Inf, -Inf], @SVector[Inf, Inf]), - solution=1.00, + f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), + domain = (@SVector[-Inf, -Inf], @SVector[Inf, Inf]), + solution = 1.00 ), (; # 2. multi-variate flipped infinite limits: Gaussian - f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), - domain=(@SVector[Inf, Inf], @SVector[-Inf, -Inf]), - solution=1.00, + f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), + domain = (@SVector[Inf, Inf], @SVector[-Inf, -Inf]), + solution = 1.00 ), (; # 3. multi-variate mixed infinite/semi-infinite upper limit: Gaussian - f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), - domain=(@SVector[-Inf, 0], @SVector[Inf, Inf]), - solution=0.5, + f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), + domain = (@SVector[-Inf, 0], @SVector[Inf, Inf]), + solution = 0.5 ), (; # 4. multi-variate mixed infinite/semi-infinite lower limit: Gaussian - f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), - domain=(@SVector[-Inf, -Inf], @SVector[Inf, 0]), - solution=0.5, + f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x), + domain = (@SVector[-Inf, -Inf], @SVector[Inf, 0]), + solution = 0.5 ), (; # 5. single-variable infinite limit: Gaussian - f=(x, p) -> pdf(Normal(0.00, 1.00), x), - domain=(-Inf, Inf), - solution=1.0, + f = (x, p) -> pdf(Normal(0.00, 1.00), x), + domain = (-Inf, Inf), + solution = 1.0 ), (; # 6. single-variable flipped infinite limit: Gaussian - f=(x, p) -> pdf(Normal(0.00, 1.00), x), - domain=(Inf, -Inf), - solution=-1.0, + f = (x, p) -> pdf(Normal(0.00, 1.00), x), + domain = (Inf, -Inf), + solution = -1.0 ), (; # 7. single-variable semi-infinite upper limit: Gaussian - f=(x, p) -> pdf(Normal(0.00, 1.00), x), - domain=(0.00, Inf), - solution=0.5, + f = (x, p) -> pdf(Normal(0.00, 1.00), x), + domain = (0.00, Inf), + solution = 0.5 ), (; # 8. single-variable flipped, semi-infinite upper limit: Gaussian - f=(x, p) -> pdf(Normal(0.00, 1.00), x), - domain=(0.00, -Inf), - solution=-0.5, + f = (x, p) -> pdf(Normal(0.00, 1.00), x), + domain = (0.00, -Inf), + solution = -0.5 ), (; # 9. single-variable semi-infinite lower limit: Gaussian - f=(x, p) -> pdf(Normal(0.00, 1.00), x), - domain=(-Inf, 0.00), - solution=0.5, + f = (x, p) -> pdf(Normal(0.00, 1.00), x), + domain = (-Inf, 0.00), + solution = 0.5 ), (; # 10. single-variable flipped, semi-infinite lower limit: Gaussian - f=(x, p) -> pdf(Normal(0.00, 1.00), x), - domain=(Inf, 0.00), - solution=-0.5, + f = (x, p) -> pdf(Normal(0.00, 1.00), x), + domain = (Inf, 0.00), + solution = -0.5 ), (; # 11. single-variable infinite limit: Lorentzian - f=(x, p) -> 1 / (x^2 + 1), - domain=(-Inf, Inf), - solution=pi/1, + f = (x, p) -> 1 / (x^2 + 1), + domain = (-Inf, Inf), + solution = pi / 1 ), (; # 12. single-variable shifted, semi-infinite lower limit: Lorentzian - f=(x, p) -> 1 / ((x-2)^2 + 1), - domain=(-Inf, 2), - solution=pi/2, + f = (x, p) -> 1 / ((x - 2)^2 + 1), + domain = (-Inf, 2), + solution = pi / 2 ), (; # 13. single-variable shifted, semi-infinite upper limit: Lorentzian - f=(x, p) -> 1 / ((x-2)^2 + 1), - domain=(2, Inf), - solution=pi/2, + f = (x, p) -> 1 / ((x - 2)^2 + 1), + domain = (2, Inf), + solution = pi / 2 ), (; # 14. single-variable flipped, shifted, semi-infinite lower limit: Lorentzian - f=(x, p) -> 1 / ((x-2)^2 + 1), - domain=(Inf, 2), - solution=-pi/2, + f = (x, p) -> 1 / ((x - 2)^2 + 1), + domain = (Inf, 2), + solution = -pi / 2 ), (; # 15. single-variable flipped, shifted, semi-infinite upper limit: Lorentzian - f=(x, p) -> 1 / ((x-2)^2 + 1), - domain=(2, -Inf), - solution=-pi/2, + f = (x, p) -> 1 / ((x - 2)^2 + 1), + domain = (2, -Inf), + solution = -pi / 2 ), (; # 16. single-variable finite limits: constant - f=(x, p) -> 1.0, - domain=(1, 3), - solution=2, + f = (x, p) -> 1.0, + domain = (1, 3), + solution = 2 ), (; # 17. single-variable flipped, finite limits: constant - f=(x, p) -> 1.0, - domain=(3, 1), - solution=-2, - ), + f = (x, p) -> 1.0, + domain = (3, 1), + solution = -2 + ) ) function f_helper!(f, y, x, p) @@ -123,7 +123,7 @@ function f_helper!(f, y, x, p) end function batch_helper(f, x, p) - map(i -> f(x[axes(x)[begin:end-1]..., i], p), axes(x)[end]) + map(i -> f(x[axes(x)[begin:(end - 1)]..., i], p), axes(x)[end]) end function batch_helper!(f, y, x, p) @@ -134,7 +134,7 @@ end do_tests = function (; f, domain, alg, abstol, reltol, solution) prob = IntegralProblem(f, domain) sol = solve(prob, alg; reltol, abstol, do_inf_transformation = Val(true)) - @test abs(only(sol) - solution) < max(abstol, reltol*abs(solution)) + @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) cache = @test_nowarn @inferred init(prob, alg; do_inf_transformation = Val(true)) @test_nowarn @inferred solve!(cache) @test_nowarn @inferred solve(prob, alg; do_inf_transformation = Val(true)) @@ -159,7 +159,7 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob @info "iip infinity test" alg=nameof(typeof(alg)) problem=j fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(size(solution))) - do_tests(; f=fiip, domain, solution, alg, abstol, reltol) + do_tests(; f = fiip, domain, solution, alg, abstol, reltol) end # BatchIntegralFunction{false} @@ -171,7 +171,7 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob @info "Batched, oop infinity test" alg=nameof(typeof(alg)) problem=j bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) - do_tests(; f=bf, domain, solution, alg, abstol, reltol) + do_tests(; f = bf, domain, solution, alg, abstol, reltol) end # BatchIntegralFunction{true} @@ -184,5 +184,5 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob @info "Batched, iip infinity test" alg=nameof(typeof(alg)) problem=j bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0)) - do_tests(; f=bfiip, domain, solution, alg, abstol, reltol) + do_tests(; f = bfiip, domain, solution, alg, abstol, reltol) end diff --git a/test/interface_tests.jl b/test/interface_tests.jl index 4eea513e..272aed54 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -10,7 +10,8 @@ abstol = 1e-3 alg_req = Dict( QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1, allows_iip = true), - QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, + QuadratureRule(gausslegendre, n = 50) => ( + nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, allows_iip = false), GaussLegendre() => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, allows_iip = false), @@ -18,7 +19,7 @@ alg_req = Dict( max_dim = Inf, allows_iip = true), # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, # allows_iip = true), - VEGASMC(seed=42) => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = Inf, + VEGASMC(seed = 42) => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = Inf, allows_iip = true), CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, allows_iip = true), @@ -32,27 +33,28 @@ alg_req = Dict( max_dim = Inf, allows_iip = true), CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, allows_iip = true), - ArblibJL() => (nout=1, allows_batch=false, min_dim=1, max_dim=1, allows_iip=true), + ArblibJL() => ( + nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, allows_iip = true) ) integrands = [ (x, p) -> 1.0, - (x, p) -> x isa Number ? cos(x) : prod(cos.(x)), + (x, p) -> x isa Number ? cos(x) : prod(cos.(x)) ] iip_integrands = [(dx, x, p) -> (dx .= f(x, p)) for f in integrands] integrands_v = [(x, p, nout) -> collect(1.0:nout) - (x, p, nout) -> integrands[2](x, p) * collect(1.0:nout)] + (x, p, nout) -> integrands[2](x, p) * collect(1.0:nout)] iip_integrands_v = [(dx, x, p, nout) -> (dx .= f(x, p, nout)) for f in integrands_v] exact_sol = [ (ndim, nout, lb, ub) -> prod(ub - lb), - (ndim, nout, lb, ub) -> prod(sin.(ub) - sin.(lb)), + (ndim, nout, lb, ub) -> prod(sin.(ub) - sin.(lb)) ] exact_sol_v = [ (ndim, nout, lb, ub) -> prod(ub - lb) * collect(1.0:nout), - (ndim, nout, lb, ub) -> exact_sol[2](ndim, nout, lb, ub) * collect(1:nout), + (ndim, nout, lb, ub) -> exact_sol[2](ndim, nout, lb, ub) * collect(1:nout) ] batch_f(f) = (pts, p) -> begin @@ -67,7 +69,7 @@ end batch_iip_f(f) = (fevals, pts, p) -> begin ax = axes(pts) for i in ax[end] - x = pts[ax[begin:(end-1)]..., i] + x = pts[ax[begin:(end - 1)]..., i] fevals[i] = f(x, p) end nothing @@ -270,10 +272,11 @@ end for dim in 1:max_dim_test lb, ub = (ones(dim), 3ones(dim)) for nout in 1:max_nout_test - prob = IntegralProblem((dx, x, p) -> iip_integrands_v[i](dx, x, p, nout), + prob = IntegralProblem( + (dx, x, p) -> iip_integrands_v[i](dx, x, p, nout), lb, ub, nout = nout) if dim > req.max_dim || dim < req.min_dim || req.nout < nout || - !req.allows_iip + !req.allows_iip continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" @@ -354,7 +357,8 @@ end @testset "Allowed keyword test" begin f(u, p) = sum(sin.(u)) prob = IntegralProblem(f, ones(3), 3ones(3)) - @test_throws Integrals.CommonKwargError((:relztol => 1e-3, :abstol => 1e-3)) solve(prob, + @test_throws Integrals.CommonKwargError((:relztol => 1e-3, :abstol => 1e-3)) solve( + prob, HCubatureJL(); relztol = 1e-3, abstol = 1e-3) diff --git a/test/nested_ad_tests.jl b/test/nested_ad_tests.jl index 079d8fae..d4506aaa 100644 --- a/test/nested_ad_tests.jl +++ b/test/nested_ad_tests.jl @@ -11,7 +11,7 @@ my_solution = my_integration(my_parameters) @test ForwardDiff.jacobian(my_integration, my_parameters) == [0.0 8.0] @test ForwardDiff.jacobian(x -> ForwardDiff.jacobian(my_integration, x), my_parameters) == [0.0 0.0 - 0.0 4.0] + 0.0 4.0] ff(x, p) = sum(sin.(x .* p)) lb = ones(2) diff --git a/test/runtests.jl b/test/runtests.jl index 8d15e5b9..090b8fd5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,6 @@ using Pkg using SafeTestsets using Test - @time @safetestset "Quality Assurance" include("qa.jl") @time @safetestset "Interface Tests" include("interface_tests.jl") @time @safetestset "Derivative Tests" include("derivative_tests.jl") diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index 19a67871..9c7398b8 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -8,7 +8,7 @@ using Integrals, Test grid2 = rand(npoints) .* (ub - lb) .+ lb grid2 = [lb; sort(grid2); ub] - grid3 = rand(npoints+1).*(ub-lb) .+ lb # also test odd number of points + grid3 = rand(npoints + 1) .* (ub - lb) .+ lb # also test odd number of points grid3 = [lb; sort(grid3); ub] exact_sols = [1 / 6 * (ub^6 - lb^6), sin(ub) - sin(lb)]