From fc9fefc3873afdb9db4fa15c7b08051d304f8c80 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:40:04 -0500 Subject: [PATCH 01/26] add type stability tests to frontend --- test/interface_tests.jl | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/test/interface_tests.jl b/test/interface_tests.jl index 272aed5..3187ecf 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -34,7 +34,7 @@ alg_req = Dict( 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) + nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, allows_iip = false) ) integrands = [ @@ -44,7 +44,7 @@ integrands = [ 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)] + let f=integrands[2]; (x, p, nout) -> f(x, p) * collect(1.0:nout); end] iip_integrands_v = [(dx, x, p, nout) -> (dx .= f(x, p, nout)) for f in integrands_v] exact_sol = [ @@ -103,7 +103,7 @@ end for i in 1:length(integrands) prob = IntegralProblem(integrands[i], lb, ub) @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 end @@ -121,7 +121,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 end @@ -140,7 +140,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg if sol.u isa Number @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 @@ -162,7 +162,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg @test sol.u[1]≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 end @@ -180,7 +180,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg if sol.u isa Number @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 @@ -204,7 +204,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg if sol.u isa Number @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 @@ -223,13 +223,13 @@ end for (alg, req) in pairs(alg_req) for i in 1:length(integrands_v) for nout in 1:max_nout_test - prob = IntegralProblem((x, p) -> integrands_v[i](x, p, nout), lb, ub, + prob = IntegralProblem(let f=integrands_v[i], nout=nout; (x, p) -> f(x, p, nout); end, lb, ub, nout = nout) if req.min_dim > 1 || req.nout < nout continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg if nout == 1 @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2 @@ -250,10 +250,10 @@ end if dim > req.max_dim || dim < req.min_dim || req.nout < nout continue end - prob = IntegralProblem((x, p) -> integrands_v[i](x, p, nout), lb, ub, + prob = IntegralProblem(let f=integrands_v[i], nout=nout; (x, p) -> f(x, p, nout); end, lb, ub, nout = nout) @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg if nout == 1 @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2 @@ -272,15 +272,15 @@ 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(let f=iip_integrands_v[i]; + (dx, x, p) -> f(dx, x, p, nout); end, lb, ub, nout = nout) if dim > req.max_dim || dim < req.min_dim || req.nout < nout || !req.allows_iip continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg if nout == 1 @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2 @@ -304,7 +304,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 end @@ -325,7 +325,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 end @@ -346,7 +346,7 @@ end continue end @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob, alg, reltol = reltol, abstol = abstol) + sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @test sol.alg == alg @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 end From a8679a452a2c5d7d8c8df9e66a707a3d1d189523 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:41:38 -0500 Subject: [PATCH 02/26] add caching tests in infinity tests --- test/inf_integral_tests.jl | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl index c0763e1..fadfaad 100644 --- a/test/inf_integral_tests.jl +++ b/test/inf_integral_tests.jl @@ -148,11 +148,11 @@ 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)) + sol = solve(prob, alg; reltol, abstol) @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) - cache = @test_nowarn @inferred init(prob, alg; do_inf_transformation = Val(true)) + cache = @test_nowarn @inferred init(prob, alg) @test_nowarn @inferred solve!(cache) - @test_nowarn @inferred solve(prob, alg; do_inf_transformation = Val(true)) + @test_nowarn @inferred solve(prob, alg) end # IntegralFunction{false} @@ -201,3 +201,17 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0)) do_tests(; f = bfiip, domain, solution, alg, abstol, reltol) end + + +@testset "Caching interface" begin + # two distinct semi-infinite transformations should still work as expected + (; f, domain, solution) = problems[8] + prob = IntegralProblem(f, domain) + cache = init(prob, QuadGKJL(); abstol, reltol) + sol = solve!(cache).u + @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) + (; domain, solution) = problems[9] + cache.domain = domain + sol = solve!(cache).u + @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) +end From a1fbb71c35f13d05891736df09b7082a28c9d39d Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:42:08 -0500 Subject: [PATCH 03/26] arblib require correct inplace type --- ext/IntegralsArblibExt.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/IntegralsArblibExt.jl b/ext/IntegralsArblibExt.jl index e0c408a..6e9d5ed 100644 --- a/ext/IntegralsArblibExt.jl +++ b/ext/IntegralsArblibExt.jl @@ -15,7 +15,8 @@ function Integrals.__solvebp_call( if isinplace(prob) res = Acb(0) - y_ = similar(prob.f.integrand_prototype, typeof(res)) + @assert res isa eltype(prob.f.integrand_prototype) "Arblib require inplace prototypes to store Acb elements" + y_ = similar(prob.f.integrand_prototype) 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, From ded012874630ae9a370adbf283c229956ba4560f Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:43:15 -0500 Subject: [PATCH 04/26] add ChangeOfVariables algorithm and deprecate do_inf_transformation --- src/Integrals.jl | 165 ++++++++++++++++++++----------- src/algorithms_meta.jl | 15 +++ src/common.jl | 12 ++- src/infinity_handling.jl | 206 ++++++++++++++++++--------------------- 4 files changed, 226 insertions(+), 172 deletions(-) create mode 100644 src/algorithms_meta.jl diff --git a/src/Integrals.jl b/src/Integrals.jl index 91d4772..476217a 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -8,6 +8,7 @@ using Reexport, MonteCarloIntegration, QuadGK, HCubature @reexport using SciMLBase using LinearAlgebra +include("algorithms_meta.jl") include("common.jl") include("algorithms.jl") include("algorithms_sampled.jl") @@ -63,9 +64,33 @@ function checkkwargs(kwargs...) return nothing end +# Give a layer for meta algorithms above AD +__solve(args...; kwargs...) = __solvebp(args...; kwargs...) # Give a layer to intercept with AD __solvebp(args...; kwargs...) = __solvebp_call(args...; kwargs...) +function init_cacheval(alg::ChangeOfVariables, prob::IntegralProblem) + f, domain = alg.fu2gv(prob.f, prob.domain) + cache_alg = init_cacheval(alg.alg, remake(prob, f = f, domain = domain)) + return (alg = cache_alg,) +end + +function __solve(cache::IntegralCache, alg::ChangeOfVariables, sensealg, udomain, p; + kwargs...) + cacheval = cache.cacheval.alg + g, vdomain = alg.fu2gv(cache.f, udomain) + _cache = IntegralCache(Val(isinplace(g)), + g, + vdomain, + p, + cache.prob_kwargs, + alg.alg, + sensealg, + cache.kwargs, + cacheval) + return __solve(_cache, alg.alg, sensealg, vdomain, p; kwargs...) +end + 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 @@ -89,132 +114,156 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p reltol = 1e-8, abstol = 1e-8, maxiters = typemax(Int)) prob = build_problem(cache) - lb_, ub_ = domain lb, ub = map(first, domain) - if !isone(length(lb_)) || !isone(length(ub_)) + if any(!isone ∘ length, domain) error("QuadGKJL only accepts one-dimensional quadrature problems.") end - mid = ((lb + ub) / 2) - if prob.f isa BatchIntegralFunction - if isinplace(prob) + f = cache.f + mid = (lb + ub) / 2 + if f isa BatchIntegralFunction + if isinplace(f) # quadgk only works with vector buffers. If the buffer is an array, we have to # turn it into a vector of arrays - bu = prob.f.integrand_prototype - f = if bu isa AbstractVector - BatchIntegrand((y, x) -> prob.f(y, x, p), similar(bu)) + prototype = f.integrand_prototype + _f = if prototype isa AbstractVector + BatchIntegrand((y, u) -> f(y, u, p), similar(prototype)) else - fsize = size(bu)[begin:(end - 1)] - 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))) - return nothing + fsize = size(prototype)[begin:(end - 1)] + BatchIntegrand{Array{eltype(prototype), length(fsize)}}() do v, u + let y = similar(v, eltype(eltype(v)), fsize..., length(v)) + f(y, u, p) + map!(collect, v, eachslice(y; dims = ndims(y))) + end + return end end - val, err = quadgk(f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters, + val, err = quadgk(_f, 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 - u = prob.f(typeof(mid)[], p) - f = if u isa AbstractVector - BatchIntegrand((y, x) -> y .= prob.f(x, p), u) + prototype = f(typeof(mid)[], p) + _f = if prototype isa AbstractVector + BatchIntegrand((y, u) -> y .= f(u, p), prototype) else - BatchIntegrand{Array{eltype(u), ndims(u) - 1}}() do y, x - map!(collect, y, eachslice(prob.f(x, p); dims = ndims(u))) - return nothing + BatchIntegrand{Array{eltype(prototype), ndims(prototype) - 1}}() do v, u + y = f(u, p) + map!(collect, v, eachslice(y; dims = ndims(y))) + return end end - val, err = quadgk(f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters, + val, err = quadgk(_f, 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) end else - 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, + if isinplace(f) + result = f.integrand_prototype * mid # result may have different units than prototype + _f = (y, u) -> f(y, u, p) + 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 - f = x -> prob.f(x, p) - val, err = quadgk(f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters, + _f = u -> f(u, p) + val, err = quadgk(_f, 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) end end end -function __solvebp_call(prob::IntegralProblem, alg::HCubatureJL, sensealg, domain, p; +function init_cacheval(alg::HCubatureJL, prob::IntegralProblem) + lb, ub = prob.domain + lb isa Number && return + if isinplace(prob) + f = u -> prob.f.integrand_prototype + else + f = u -> prob.f(u, prob.p) + end + return HCubature.hcubature_buffer(f, lb, ub; norm = alg.norm) +end +function __solvebp_call(cache::IntegralCache, alg::HCubatureJL, sensealg, domain, p; reltol = 1e-8, abstol = 1e-8, maxiters = typemax(Int)) + prob = build_problem(cache) lb, ub = domain + f = cache.f - @assert prob.f isa IntegralFunction - if isinplace(prob) + @assert f isa IntegralFunction + if isinplace(f) # allocate a new output array at each evaluation since HCubature.jl doesn't support # inplace ops - f = x -> (dx = similar(prob.f.integrand_prototype); prob.f(dx, x, prob.p); dx) + prototype = f.integrand_prototype + _f = let f = f.f + u -> (y = similar(prototype); f(y, u, p); y) + end else - f = x -> prob.f(x, prob.p) + _f = let f = f.f + u -> f(u, p) + end end - if lb isa Number - val, err = hquadrature(f, lb, ub; - rtol = reltol, atol = abstol, + val, err = if lb isa Number + hquadrature(_f, lb, ub; + rtol = reltol, atol = abstol, buffer = cache.cacheval, maxevals = maxiters, norm = alg.norm, initdiv = alg.initdiv) else - val, err = hcubature(f, lb, ub; - rtol = reltol, atol = abstol, - maxevals = maxiters, norm = alg.norm, initdiv = alg.initdiv) + ret = _f((lb + ub) / 2) * (prod(ub - lb) / 2) # this calculation for type stability with vector endpoints + hcubature(_f, lb, ub; + rtol = reltol, atol = abstol, buffer = cache.cacheval, + maxevals = maxiters, norm = alg.norm, + initdiv = alg.initdiv)::Tuple{typeof(ret), typeof(alg.norm(ret))} end SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) end function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p; - reltol = 1e-8, abstol = 1e-8, + reltol = 1e-4, abstol = 1e-4, maxiters = 1000) lb, ub = domain mid = (lb + ub) / 2 - if prob.f isa BatchIntegralFunction + f = prob.f + if f isa BatchIntegralFunction + # MonteCarloIntegration v0.0.x passes points as rows of a matrix + # MonteCarloIntegration v0.1 passes batches as a vector of views of + # a matrix with points as columns of a matrix + # see https://github.com/ranjanan/MonteCarloIntegration.jl/issues/16 + # This is an ugly hack that is compatible with both + wrangle = x -> begin + xx = eltype(x) <: SubArray ? parent(first(x)) : x' + mid isa Number ? vec(xx) : xx + end if isinplace(prob) y = similar(prob.f.integrand_prototype, size(prob.f.integrand_prototype)[begin:(end - 1)]..., prob.f.max_batch) - # MonteCarloIntegration v0.0.x passes points as rows of a matrix - # MonteCarloIntegration v0.1 passes batches as a vector of views of - # a matrix with points as columns of a matrix - # see https://github.com/ranjanan/MonteCarloIntegration.jl/issues/16 - # 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)) + _f = x -> (f(y, wrangle(x), p); vec(y)) else - 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) + y = mid isa Number ? f(typeof(mid)[], p) : + f(Matrix{eltype(mid)}(undef, length(mid), 0), p) + _f = x -> vec(f(wrangle(x), p)) end else if isinplace(prob) @assert prob.f.integrand_prototype isa AbstractArray{<:Real}&&length(prob.f.integrand_prototype) == 1 "VEGAS only supports Float64-valued integrands" y = similar(prob.f.integrand_prototype) - f = x -> (prob.f(y, x, p); only(y)) + _f = x -> (prob.f(y, mid isa Number ? only(x) : x, p); only(y)) else y = prob.f(mid, p) - f = x -> only(prob.f(x, prob.p)) + _f = x -> only(prob.f(mid isa Number ? only(x) : x, prob.p)) end end - if prob.f isa BatchIntegralFunction + if f isa BatchIntegralFunction @assert prod(size(y)[begin:(end - 1)]) == 1&&eltype(y) <: Real "VEGAS only supports Float64-valued scalar integrands" else @assert length(y) == 1&&eltype(y) <: Real "VEGAS only supports Float64-valued scalar integrands" end - ncalls = prob.f isa BatchIntegralFunction ? prob.f.max_batch : alg.ncalls - out = vegas(f, lb, ub, rtol = reltol, atol = abstol, + ncalls = alg.ncalls + out = vegas(_f, lb, ub, rtol = reltol, atol = abstol, maxiter = maxiters, nbins = alg.nbins, debug = alg.debug, ncalls = ncalls, batch = prob.f isa BatchIntegralFunction) val, err, chi = out isa Tuple ? out : diff --git a/src/algorithms_meta.jl b/src/algorithms_meta.jl new file mode 100644 index 0000000..8e9fc30 --- /dev/null +++ b/src/algorithms_meta.jl @@ -0,0 +1,15 @@ +abstract type AbstractIntegralMetaAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end + +""" + ChangeOfVariables(fu2gv, alg) + +Apply a change of variables from `∫ f(u,p) du` to an equivalent integral `∫ g(v,p) dv` using +a helper function `fu2gv(f, u_domain) -> (g, v_domain)` where `f` and `g` should be +integral functions. Acts as a wrapper to algorithm `alg` +""" +# internal algorithm +struct ChangeOfVariables{T, A <: SciMLBase.AbstractIntegralAlgorithm} <: + AbstractIntegralMetaAlgorithm + fu2gv::T + alg::A +end diff --git a/src/common.jl b/src/common.jl index cd31327..4c82fea 100644 --- a/src/common.jl +++ b/src/common.jl @@ -19,15 +19,17 @@ function SciMLBase.init(prob::IntegralProblem{iip}, sensealg = ReCallVJP(ZygoteVJP()), do_inf_transformation = nothing, kwargs...) where {iip} checkkwargs(kwargs...) - prob = transformation_if_inf(prob, do_inf_transformation) - cacheval = init_cacheval(alg, prob) + do_inf_transformation === nothing || + @warn "do_inf_transformation is deprecated. All integral problems are transformed" + _alg = alg isa ChangeOfVariables ? alg : ChangeOfVariables(transformation_if_inf, alg) + cacheval = init_cacheval(_alg, prob) IntegralCache{iip, typeof(prob.f), typeof(prob.domain), typeof(prob.p), typeof(prob.kwargs), - typeof(alg), + typeof(_alg), typeof(sensealg), typeof(kwargs), typeof(cacheval)}(Val(iip), @@ -35,7 +37,7 @@ function SciMLBase.init(prob::IntegralProblem{iip}, prob.domain, prob.p, prob.kwargs, - alg, + _alg, sensealg, kwargs, cacheval) @@ -102,7 +104,7 @@ function SciMLBase.solve(prob::IntegralProblem, end function SciMLBase.solve!(cache::IntegralCache) - __solvebp(cache, cache.alg, cache.sensealg, cache.domain, cache.p; + __solve(cache, cache.alg, cache.sensealg, cache.domain, cache.p; cache.kwargs...) end diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl index 25cfa69..93cf60a 100644 --- a/src/infinity_handling.jl +++ b/src/infinity_handling.jl @@ -1,7 +1,92 @@ -_oftype(x, y) = oftype(x, y) -_oftype(::SubArray, y) = y +# generic change of variables code -function substitute_bounds(lb, ub) +function substitute_u(u2v, lb, ub) + bounds = map(u2v, lb, ub) + lb_sub = lb isa Number ? first(bounds) : map(first, bounds) + ub_sub = ub isa Number ? last(bounds) : map(last, bounds) + return (lb_sub, ub_sub) +end + +# without batching point container type should match the inputs +substitute_v(v2ujac, v, lb::Number, ub::Number) = v2ujac(only(v), lb, ub) +function substitute_v(v2ujac, v, lb::AbstractVector, ub::AbstractVector) + xjac = map((l, u, v) -> v2ujac(v, l, u), lb, ub, v) # ordering may influence container type + x = map(first, xjac) + jac = prod(last, xjac) + return x, jac +end + +function substitute_f(f::IntegralFunction{false}, v2ujac, lb, ub) + _f = f.f + IntegralFunction{false}(f.integrand_prototype) do v, p + u, jac = substitute_v(v2ujac, v, lb, ub) + return _f(u, p) * jac + end +end +function substitute_f(f::IntegralFunction{true}, v2ujac, lb, ub) + _f = f.f + prototype = similar(f.integrand_prototype) + vol = prod((ub - lb) / 2) # just to get the type of the jacobian determinant + IntegralFunction{true}(prototype * vol) do y, v, p + u, jac = substitute_v(v2ujac, v, lb, ub) + _f(prototype, u, p) + y .= prototype .* jac + return + end +end + +# with batching the point container type is assumed to be mutable +function substitute_bv(v2ujac, v::AbstractArray, lb::Number, ub::Number) + x = similar(v, typeof(one(eltype(v)) * (first(lb) + first(ub)))) + jac = similar(x) + for i in axes(v, 1) + x[i], jac[i] = v2ujac(v[i], lb, ub) + end + return x, jac +end +function substitute_bv(v2ujac, v::AbstractArray, lb::AbstractVector, ub::AbstractVector) + x = similar(v, typeof(one(eltype(v)) * (first(lb) + first(ub)))) + jac = similar(v, typeof(zero(eltype(v))*prod(lb)), size(v)[end]) + idx = CartesianIndices(axes(v)[begin:end-1]) + for i in axes(v)[end] + _jac = one(eltype(jac)) + for (ii, l, u) in zip(idx, lb, ub) + x[ii, i], j = v2ujac(v[ii, i] , l, u) + _jac *= j + end + jac[i] = _jac + end + return x, jac +end + +function substitute_f(f::BatchIntegralFunction{false}, v2ujac::F, lb, ub) where {F} + _f = f.f + BatchIntegralFunction{false}(f.integrand_prototype, max_batch = f.max_batch) do v, p + u, jac = substitute_bv(v2ujac, v, lb, ub) + y = _f(u, p) + return y .* reshape(jac, ntuple(d -> d == ndims(y) ? length(jac) : 1, ndims(y))) + end +end +function substitute_f(f::BatchIntegralFunction{true}, v2ujac::F, lb, ub) where {F} + _f = f.f + prototype = similar(f.integrand_prototype) + vol = prod((ub - lb) / 2) # just to get the type of the jacobian determinant + BatchIntegralFunction{true}(prototype * vol, max_batch = f.max_batch) do y, v, p + u, jac = substitute_bv(v2ujac, v, lb, ub) + _prototype = similar(prototype, size(y)) + _f(_prototype, u, p) + for (i, j) in zip(axes(y)[end], jac) + for iy in CartesianIndices(axes(y)[begin:(end - 1)]) + y[iy, i] = j * _prototype[iy, i] + end + end + return + end +end + +# specific changes of variables + +function u2t(lb, ub) mid = (lb + ub) / 2 # floating-point promotion if isinf(lb) && isinf(ub) lb_sub = flipsign(one(mid), lb) @@ -19,7 +104,7 @@ function substitute_bounds(lb, ub) return lb_sub, ub_sub # unitless end -function substitute_t(t::Number, lb::Number, ub::Number) +function t2ujac(t::Number, lb::Number, ub::Number) u = oneunit(eltype(lb)) # apply correct units if isinf(lb) && isinf(ub) @@ -32,114 +117,17 @@ function substitute_t(t::Number, lb::Number, ub::Number) den = inv(1 - flipsign(t, ub)) lb + t * den * u, den^2 * u else - den = (ub - lb) * oftype(t, 0.5) + den = (ub - lb) * oftype(float(one(u)), 0.5) lb + (1 + t) * den, den end end -function substitute_t(t::AbstractVector, lb::AbstractVector, ub::AbstractVector) - 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]) - jac *= dj - end - return _oftype(t, x), jac -end -function substitute_f(f::F, t, p, lb, ub) where {F} - x, jac = substitute_t(t, lb, ub) - return f(x, p) * jac -end -function substitute_f(f::F, dt, t, p, lb, ub) where {F} - x, jac = substitute_t(t, lb, ub) - f(dt, x, p) - dt .*= jac - return +function transformation_if_inf(f, domain) + lb, ub = promote(domain...) + tdomain = substitute_u(u2t, lb, ub) + g = substitute_f(f, t2ujac, lb, ub) + return g, tdomain end -function substitute_t(t::AbstractVector, lb::Number, ub::Number) - 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) - end - return x, jac -end -function substitute_t(t::AbstractArray, lb::AbstractVector, ub::AbstractVector) - 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) - end - return x, jac -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))) -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)]) - dt[idt, i] *= j - end - end - return -end - -function transformation_if_inf(prob, ::Val{true}) - lb, ub = promote(prob.domain...) - 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) - 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, - 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, - 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, - f.integrand_prototype) - else - IntegralFunction{false}(let f = f.f - (t, p) -> substitute_f(f, t, p, lb, ub) - end, - f.integrand_prototype) - end - end - return remake(prob, f = f_sub, domain = (lb_sub, ub_sub)) -end - -function transformation_if_inf(prob, ::Nothing) - lb, ub = prob.domain - if any(isinf, lb) || any(isinf, ub) - return transformation_if_inf(prob, Val(true)) - else - return transformation_if_inf(prob, Val(false)) - end -end - -function transformation_if_inf(prob, ::Val{false}) - return prob -end - -function transformation_if_inf(prob, do_inf_transformation = nothing) - transformation_if_inf(prob, do_inf_transformation) -end +# to implement more transformations, define functions u2v and v2ujac and then a wrapper +# similar to transformation_if_inf to pass to ChangeOfVariables From 20cd125050bb627122d17ae359c660a4879b5d9c Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:43:49 -0500 Subject: [PATCH 05/26] make MCIntegration wrapper type stable --- ext/IntegralsMCIntegrationExt.jl | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/ext/IntegralsMCIntegrationExt.jl b/ext/IntegralsMCIntegrationExt.jl index e4d6779..994b366 100644 --- a/ext/IntegralsMCIntegrationExt.jl +++ b/ext/IntegralsMCIntegrationExt.jl @@ -2,46 +2,53 @@ module IntegralsMCIntegrationExt using MCIntegration, Integrals +_oftype(::Number, x) = only(x) +_oftype(y, x) = oftype(y, x) + function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, domain, p; reltol = nothing, abstol = nothing, maxiters = 10) lb, ub = domain - mid = vec(collect((lb + ub) / 2)) + mid = (lb + ub) / 2 + tmp = vec(collect(mid)) 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") else - if isinplace(prob) - f0 = similar(prob.f.integrand_prototype) + f0 = if isinplace(prob) + _f0 = similar(prob.f.integrand_prototype) f_ = (x, f, c) -> begin n = 0 for v in x - mid[n += 1] = first(v) + tmp[n += 1] = first(v) end - prob.f(f0, mid, p) - f .= vec(f0) + prob.f(_f0, _oftype(mid, tmp), p) + f .= vec(_f0) end + _f0 else - f0 = prob.f(mid, p) f_ = (x, c) -> begin n = 0 for v in x - mid[n += 1] = first(v) + tmp[n += 1] = first(v) end - fx = prob.f(mid, p) + fx = prob.f(_oftype(mid, tmp), p) fx isa AbstractArray ? vec(fx) : fx end + prob.f(mid, p) 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...) - out, err, chi = if f0 isa Number + # the package itself is not type-stable + out::typeof(f0), err::typeof(f0), chi2 = 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) + chi::typeof(f0) = map(sqrt, chi2) + SciMLBase.build_solution(prob, alg, out, err, chi = chi, + retcode = ReturnCode.Success) end end From 85d749c8d364aa6616fa3296b166f58de3a14d10 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:45:35 -0500 Subject: [PATCH 06/26] Cubature add let bindings for type stable closures --- ext/IntegralsCubatureExt.jl | 64 ++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 30 deletions(-) diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl index 9033291..ed6ae76 100644 --- a/ext/IntegralsCubatureExt.jl +++ b/ext/IntegralsCubatureExt.jl @@ -13,7 +13,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, mid = (lb + ub) / 2 # we get to pick fdim or not based on the IntegralFunction and its output dimensions - y = if prob.f isa BatchIntegralFunction + prototype = if prob.f isa BatchIntegralFunction isinplace(prob.f) ? prob.f.integrand_prototype : mid isa Number ? prob.f(eltype(mid)[], p) : prob.f(Matrix{eltype(mid)}(undef, length(mid), 0), p) @@ -22,21 +22,22 @@ function Integrals.__solvebp_call(prob::IntegralProblem, isinplace(prob.f) ? prob.f.integrand_prototype : prob.f(mid, p) end - @assert eltype(y)<:Real "Cubature.jl is only compatible with real-valued integrands" + @assert eltype(prototype)<:Real "Cubature.jl is only compatible with real-valued integrands" if prob.f isa BatchIntegralFunction - if y isa AbstractVector # this branch could be omitted since the following one should work similarly + if prototype isa AbstractVector # this branch could be omitted since the following one should work similarly if isinplace(prob) # dx is a Vector, but we provide the integrand a vector of the same type as # y, which needs to be resized since the number of batch points changes. - dy = similar(y) - f = (x, dx) -> begin - resize!(dy, length(dx)) - prob.f(dy, x, p) - dx .= dy + f = let y = similar(prototype) + (u, v) -> begin + resize!(y, length(v)) + prob.f(y, u, p) + v .= y + end end else - f = (x, dx) -> (dx .= prob.f(x, p)) + f = (u, v) -> (v .= prob.f(u, p)) end if mid isa Number if alg isa CubatureJLh @@ -59,50 +60,52 @@ function Integrals.__solvebp_call(prob::IntegralProblem, maxevals = maxiters) end end - elseif y isa AbstractArray - bfsize = size(y)[begin:(end - 1)] - bfdim = prod(bfsize) + elseif prototype isa AbstractArray + fsize = size(prototype)[begin:(end - 1)] + fdim = prod(fsize) if isinplace(prob) # dx is a Matrix, but to provide a buffer of the same type as y, we make # would like to make views of a larger buffer, but CubatureJL doesn't set # a hard limit for max_batch, so we allocate a new buffer with the needed size - f = (x, dx) -> begin - dy = similar(y, bfsize..., size(dx, 2)) - prob.f(dy, x, p) - dx .= reshape(dy, bfdim, size(dx, 2)) + f = let fsize = fsize + (u, v) -> begin + y = similar(prototype, fsize..., size(v, 2)) + prob.f(y, u, p) + v .= reshape(y, fdim, size(v, 2)) + end end else - f = (x, dx) -> (dx .= reshape(prob.f(x, p), bfdim, size(dx, 2))) + f = (u, v) -> (v .= reshape(prob.f(u, p), fdim, size(v, 2))) end if mid isa Number if alg isa CubatureJLh - val_, err = Cubature.hquadrature_v(bfdim, f, lb, ub; + val_, err = Cubature.hquadrature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) else - val_, err = Cubature.pquadrature_v(bfdim, f, lb, ub; + val_, err = Cubature.pquadrature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) end else if alg isa CubatureJLh - val_, err = Cubature.hcubature_v(bfdim, f, lb, ub; + val_, err = Cubature.hcubature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) else - val_, err = Cubature.pcubature_v(bfdim, f, lb, ub; + val_, err = Cubature.pcubature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) end end - val = reshape(val_, bfsize...) + val = reshape(val_, fsize...) else error("BatchIntegralFunction integrands must be arrays for Cubature.jl") end else - if y isa Real + if prototype isa Real # no inplace in this case, since the integrand_prototype would be mutable - f = x -> prob.f(x, p) + f = u -> prob.f(u, p) if lb isa Number if alg isa CubatureJLh val, err = Cubature.hquadrature(f, lb, ub; @@ -124,14 +127,15 @@ function Integrals.__solvebp_call(prob::IntegralProblem, maxevals = maxiters) end end - elseif y isa AbstractArray - fsize = size(y) - fdim = length(y) + elseif prototype isa AbstractArray + fsize = size(prototype) + fdim = length(prototype) if isinplace(prob) - dy = similar(y) - f = (x, v) -> (prob.f(dy, x, p); v .= vec(dy)) + f = let y = similar(prototype) + (u, v) -> (prob.f(y, u, p); v .= vec(y)) + end else - f = (x, v) -> (v .= vec(prob.f(x, p))) + f = (u, v) -> (v .= vec(prob.f(u, p))) end if mid isa Number if alg isa CubatureJLh From 3ff7c5e80989014fd9e1790a29ee8bafbe24b2e4 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:48:52 -0500 Subject: [PATCH 07/26] make Cuba type stable --- ext/IntegralsCubaExt.jl | 51 +++++++++++++++++++++++------------------ 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl index f5c1697..40c709c 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, @@ -34,17 +34,19 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori end if isinplace(prob) - fsize = size(prob.f.integrand_prototype)[begin:(end - 1)] - y = similar(prob.f.integrand_prototype, fsize..., nvec) + y = prob.f.integrand_prototype + fsize = size(y)[begin:(end - 1)] ax = map(_ -> (:), fsize) - f = function (x, dx) - dy = @view(y[ax..., begin:(begin + size(dx, 2) - 1)]) - prob.f(dy, scale(x), p) - dx .= reshape(dy, :, size(dx, 2)) .* vol + f = let y = similar(y, fsize..., nvec) + function (x, dx) + dy = @view(y[ax..., begin:(begin + size(dx, 2) - 1)]) + prob.f(dy, scale(x), p) + dx .= reshape(dy, :, size(dx, 2)) .* vol + end end else y = mid isa Number ? prob.f(typeof(mid)[], p) : - prob.f(Matrix{typeof(mid)}(undef, length(mid), 0), p) + prob.f(Matrix{eltype(mid)}(undef, length(mid), 0), p) fsize = size(y)[begin:(end - 1)] f = (x, dx) -> dx .= reshape(prob.f(scale(x), p), :, size(dx, 2)) .* vol end @@ -60,8 +62,13 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori end if isinplace(prob) - y = similar(prob.f.integrand_prototype) - f = (x, dx) -> dx .= vec(prob.f(y, scale(x), p)) .* vol + y = prob.f.integrand_prototype + f = let y = similar(y) + (x, dx) -> begin + prob.f(y, scale(x), p) + dx .= vec(y) .* vol + end + end else y = prob.f(mid, p) f = (x, dx) -> dx .= Iterators.flatten(prob.f(scale(x), p)) .* vol @@ -69,21 +76,21 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori ncomp = length(y) end - if alg isa CubaVegas - out = Cuba.vegas(f, ndim, ncomp; rtol = reltol, + out = if alg isa CubaVegas + Cuba.vegas(f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, seed = alg.seed, minevals = alg.minevals, nstart = alg.nstart, nincrease = alg.nincrease, gridno = alg.gridno) elseif alg isa CubaSUAVE - out = Cuba.suave(f, ndim, ncomp; rtol = reltol, + Cuba.suave(f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, seed = alg.seed, minevals = alg.minevals, nnew = alg.nnew, nmin = alg.nmin, flatness = alg.flatness) elseif alg isa CubaDivonne - out = Cuba.divonne(f, ndim, ncomp; rtol = reltol, + Cuba.divonne(f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, seed = alg.seed, minevals = alg.minevals, @@ -91,26 +98,26 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori maxpass = alg.maxpass, border = alg.border, maxchisq = alg.maxchisq, mindeviation = alg.mindeviation) elseif alg isa CubaCuhre - out = Cuba.cuhre(f, ndim, ncomp; rtol = reltol, + Cuba.cuhre(f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, minevals = alg.minevals, key = alg.key) end # out.integral is a Vector{Float64}, but we want to return it to the shape of the integrand - if prob.f isa BatchIntegralFunction + val = if prob.f isa BatchIntegralFunction if y isa AbstractVector - val = out.integral[1] + out.integral[1] else - val = reshape(out.integral, fsize) + reshape(out.integral, fsize) end else if y isa Real - val = out.integral[1] + out.integral[1] elseif y isa AbstractVector - val = out.integral + out.integral else - val = reshape(out.integral, size(y)) + reshape(out.integral, size(y)) end end From d40d099dc5f26802855a0281daae976258ae96aa Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:49:05 -0500 Subject: [PATCH 08/26] update docs --- docs/src/basics/solve.md | 4 ---- 1 file changed, 4 deletions(-) diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md index d513a1e..c572b00 100644 --- a/docs/src/basics/solve.md +++ b/docs/src/basics/solve.md @@ -3,7 +3,3 @@ ```@docs solve(prob::IntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm) ``` - -Additionally, the extra keyword arguments are splatted to the library calls, so -see the documentation of the integrator library for all the extra details. -These extra keyword arguments are not guaranteed to act uniformly. From c1b0b1a27b616361a2a0ce90c786acfabb7aba6c Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 09:51:34 -0500 Subject: [PATCH 09/26] apply format --- ext/IntegralsCubaExt.jl | 4 ++-- ext/IntegralsForwardDiffExt.jl | 2 +- ext/IntegralsZygoteExt.jl | 1 - src/infinity_handling.jl | 6 +++--- test/inf_integral_tests.jl | 3 +-- test/interface_tests.jl | 17 ++++++++++++----- 6 files changed, 19 insertions(+), 14 deletions(-) diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl index 40c709c..c67be23 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, diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index ee840cd..c6cec18 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -21,7 +21,7 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, prob = Integrals.build_problem(cache) dprob = remake(prob, f = df) dcache = init( - dprob, alg; sensealg = sensealg, do_inf_transformation = Val(false), kwargs...) + dprob, alg; sensealg = sensealg, kwargs...) Integrals.__solvebp_call(dcache, alg, sensealg, domain, p; kwargs...) else Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...) diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index 1b38956..c04e811 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -89,7 +89,6 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal dp_cache = init(dp_prob, alg; sensealg = sensealg, - do_inf_transformation = Val(false), cache.kwargs...) project_p = ProjectTo(p) diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl index 93cf60a..f552179 100644 --- a/src/infinity_handling.jl +++ b/src/infinity_handling.jl @@ -46,12 +46,12 @@ function substitute_bv(v2ujac, v::AbstractArray, lb::Number, ub::Number) end function substitute_bv(v2ujac, v::AbstractArray, lb::AbstractVector, ub::AbstractVector) x = similar(v, typeof(one(eltype(v)) * (first(lb) + first(ub)))) - jac = similar(v, typeof(zero(eltype(v))*prod(lb)), size(v)[end]) - idx = CartesianIndices(axes(v)[begin:end-1]) + jac = similar(v, typeof(zero(eltype(v)) * prod(lb)), size(v)[end]) + idx = CartesianIndices(axes(v)[begin:(end - 1)]) for i in axes(v)[end] _jac = one(eltype(jac)) for (ii, l, u) in zip(idx, lb, ub) - x[ii, i], j = v2ujac(v[ii, i] , l, u) + x[ii, i], j = v2ujac(v[ii, i], l, u) _jac *= j end jac[i] = _jac diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl index fadfaad..2d28a9e 100644 --- a/test/inf_integral_tests.jl +++ b/test/inf_integral_tests.jl @@ -202,7 +202,6 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob do_tests(; f = bfiip, domain, solution, alg, abstol, reltol) end - @testset "Caching interface" begin # two distinct semi-infinite transformations should still work as expected (; f, domain, solution) = problems[8] @@ -210,7 +209,7 @@ end cache = init(prob, QuadGKJL(); abstol, reltol) sol = solve!(cache).u @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) - (; domain, solution) = problems[9] + (; domain, solution) = problems[9] cache.domain = domain sol = solve!(cache).u @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) diff --git a/test/interface_tests.jl b/test/interface_tests.jl index 3187ecf..e58ed99 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -44,7 +44,9 @@ integrands = [ iip_integrands = [(dx, x, p) -> (dx .= f(x, p)) for f in integrands] integrands_v = [(x, p, nout) -> collect(1.0:nout) - let f=integrands[2]; (x, p, nout) -> f(x, p) * collect(1.0:nout); end] + let f = integrands[2] + (x, p, nout) -> f(x, p) * collect(1.0:nout) + end] iip_integrands_v = [(dx, x, p, nout) -> (dx .= f(x, p, nout)) for f in integrands_v] exact_sol = [ @@ -223,7 +225,9 @@ end for (alg, req) in pairs(alg_req) for i in 1:length(integrands_v) for nout in 1:max_nout_test - prob = IntegralProblem(let f=integrands_v[i], nout=nout; (x, p) -> f(x, p, nout); end, lb, ub, + prob = IntegralProblem(let f = integrands_v[i], nout = nout + (x, p) -> f(x, p, nout) + end, lb, ub, nout = nout) if req.min_dim > 1 || req.nout < nout continue @@ -250,7 +254,9 @@ end if dim > req.max_dim || dim < req.min_dim || req.nout < nout continue end - prob = IntegralProblem(let f=integrands_v[i], nout=nout; (x, p) -> f(x, p, nout); end, lb, ub, + prob = IntegralProblem(let f = integrands_v[i], nout = nout + (x, p) -> f(x, p, nout) + end, lb, ub, nout = nout) @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout" sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol) @@ -272,8 +278,9 @@ end for dim in 1:max_dim_test lb, ub = (ones(dim), 3ones(dim)) for nout in 1:max_nout_test - prob = IntegralProblem(let f=iip_integrands_v[i]; - (dx, x, p) -> f(dx, x, p, nout); end, + prob = IntegralProblem(let f = iip_integrands_v[i] + (dx, x, p) -> f(dx, x, p, nout) + end, lb, ub, nout = nout) if dim > req.max_dim || dim < req.min_dim || req.nout < nout || !req.allows_iip From d8e8f532f4497edd5a08139e6b05ea90f474c343 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 10:16:50 -0500 Subject: [PATCH 10/26] add VEGAS to tests with seed --- Project.toml | 2 ++ src/Integrals.jl | 2 ++ src/algorithms.jl | 7 ++++--- test/interface_tests.jl | 4 ++-- 4 files changed, 10 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 34abb67..02f55d0 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MonteCarloIntegration = "4886b29c-78c9-11e9-0a6e-41e1f4161f7b" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" @@ -48,6 +49,7 @@ MCIntegration = "0.4.2" MonteCarloIntegration = "0.2" Pkg = "1.10" QuadGK = "2.9" +Random = "1.10" Reexport = "1.0" SafeTestsets = "0.1" SciMLBase = "2.6" diff --git a/src/Integrals.jl b/src/Integrals.jl index 476217a..51680bb 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -7,6 +7,7 @@ end using Reexport, MonteCarloIntegration, QuadGK, HCubature @reexport using SciMLBase using LinearAlgebra +using Random include("algorithms_meta.jl") include("common.jl") @@ -224,6 +225,7 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p; lb, ub = domain mid = (lb + ub) / 2 f = prob.f + alg.seed !== nothing && Random.seed!(alg.seed) if f isa BatchIntegralFunction # MonteCarloIntegration v0.0.x passes points as rows of a matrix # MonteCarloIntegration v0.1 passes batches as a vector of views of diff --git a/src/algorithms.jl b/src/algorithms.jl index 4222b99..8b32fe8 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -53,7 +53,7 @@ end HCubatureJL(; initdiv = 1, norm = norm) = HCubatureJL(initdiv, norm) """ - VEGAS(; nbins = 100, ncalls = 1000, debug=false) + VEGAS(; nbins = 100, ncalls = 1000, debug=false, seed = nothing) Multidimensional adaptive Monte Carlo integration from MonteCarloIntegration.jl. Importance sampling is used to reduce variance. @@ -80,12 +80,13 @@ year={1978}, publisher={Elsevier} } """ -struct VEGAS <: SciMLBase.AbstractIntegralAlgorithm +struct VEGAS{S} <: SciMLBase.AbstractIntegralAlgorithm nbins::Int ncalls::Int debug::Bool + seed::S end -VEGAS(; nbins = 100, ncalls = 1000, debug = false) = VEGAS(nbins, ncalls, debug) +VEGAS(; nbins = 100, ncalls = 1000, debug = false, seed = nothing) = VEGAS(nbins, ncalls, debug, seed) """ GaussLegendre{C, N, W} diff --git a/test/interface_tests.jl b/test/interface_tests.jl index e58ed99..920681e 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -17,8 +17,8 @@ alg_req = Dict( allows_iip = false), 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), + VEGAS(seed = 109) => (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, allows_iip = true), CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, From c38b1f8f8de4af409f6d96f6676f1b6de499346e Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 17:52:56 -0500 Subject: [PATCH 11/26] fix bugs --- ext/IntegralsForwardDiffExt.jl | 10 ++++++--- ext/IntegralsZygoteExt.jl | 37 ++++++++++++++++++++++++++++------ src/infinity_handling.jl | 22 +++++++++++--------- 3 files changed, 50 insertions(+), 19 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index c6cec18..3f6548b 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -3,6 +3,12 @@ using Integrals isdefined(Base, :get_extension) ? (using ForwardDiff) : (using ..ForwardDiff) ### Forward-Mode AD Intercepts +function Integrals._evaluate!(f, y, u, p::Union{D, AbstractArray{<:D}}) where {T, V, P, D <: ForwardDiff.Dual{T, V, P}} + dy = similar(y, replace_dualvaltype(eltype(p), eltype(y))) + f(dy, u, p) + return dy +end + # Default to direct AD on solvers function Integrals.__solvebp(cache, alg, sensealg, domain, p::Union{D, AbstractArray{<:D}}; @@ -90,13 +96,11 @@ function Integrals.__solvebp( prob = Integrals.build_problem(cache) dp_prob = remake(prob, f = dfdp, p = rawp) - # the infinity transformation was already applied to f so we don't apply it to dfdp dp_cache = init(dp_prob, alg; sensealg = sensealg, - do_inf_transformation = Val(false), cache.kwargs...) - dual = Integrals.__solvebp_call(dp_cache, alg, sensealg, domain, rawp; kwargs...) + dual = solve!(dp_cache) res = reinterpret(reshape, DT, dual.u) # unwrap the dual when the primal would return a scalar diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index c04e811..dad0874 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -14,6 +14,36 @@ ChainRulesCore.@non_differentiable Integrals.checkkwargs(kwargs...) ChainRulesCore.@non_differentiable Integrals.isinplace(f, args...) # fixes #99 ChainRulesCore.@non_differentiable Integrals.init_cacheval(alg, prob) +function ChainRulesCore.rrule(::typeof(Integrals.transformation_if_inf), f, domain) + function transformation_if_inf_pullback(Δ) + return NoTangent(), Δ... + end + return Integrals.transformation_if_inf(f, domain), transformation_if_inf_pullback +end + +# we will need to implement the following adjoints when we compute ∂f/∂u +function ChainRulesCore.rrule(::typeof(Integrals.substitute_v), args...) + function substitute_v_pullback(_) + return NoTangent(), ntuple(_ -> NoTangent(), length(args))... + end + return Integrals.substitute_v(args...), substitute_v_pullback +end +function ChainRulesCore.rrule(::typeof(Integrals.substitute_bv), args...) + function substitute_bv_pullback(_) + return NoTangent(), ntuple(_ -> NoTangent(), length(args))... + end + return Integrals.substitute_bv(args...), substitute_bv_pullback +end +function ChainRulesCore.rrule(::typeof(Integrals._evaluate!), f, y, u, p) + out, back = Zygote.pullback(y, u, p) do y, u, p + b = Zygote.Buffer(y) + f(b, u, p) + return copy(b) + end + out, Δ -> (NoTangent(), NoTangent(), back(Δ)...) +end + + function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, sensealg, domain, p; kwargs...) @@ -92,12 +122,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal cache.kwargs...) project_p = ProjectTo(p) - dp = project_p(Integrals.__solvebp_call(dp_cache, - alg, - sensealg, - domain, - p; - kwargs...).u) + dp = project_p(solve!(dp_cache).u) lb, ub = domain if lb isa Number diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl index f552179..b822362 100644 --- a/src/infinity_handling.jl +++ b/src/infinity_handling.jl @@ -29,8 +29,8 @@ function substitute_f(f::IntegralFunction{true}, v2ujac, lb, ub) vol = prod((ub - lb) / 2) # just to get the type of the jacobian determinant IntegralFunction{true}(prototype * vol) do y, v, p u, jac = substitute_v(v2ujac, v, lb, ub) - _f(prototype, u, p) - y .= prototype .* jac + _y = _evaluate!(f, prototype, u, p) + y .= _y .* jac return end end @@ -59,7 +59,7 @@ function substitute_bv(v2ujac, v::AbstractArray, lb::AbstractVector, ub::Abstrac return x, jac end -function substitute_f(f::BatchIntegralFunction{false}, v2ujac::F, lb, ub) where {F} +function substitute_f(f::BatchIntegralFunction{false}, v2ujac, lb, ub) _f = f.f BatchIntegralFunction{false}(f.integrand_prototype, max_batch = f.max_batch) do v, p u, jac = substitute_bv(v2ujac, v, lb, ub) @@ -67,23 +67,25 @@ function substitute_f(f::BatchIntegralFunction{false}, v2ujac::F, lb, ub) where return y .* reshape(jac, ntuple(d -> d == ndims(y) ? length(jac) : 1, ndims(y))) end end -function substitute_f(f::BatchIntegralFunction{true}, v2ujac::F, lb, ub) where {F} +function substitute_f(f::BatchIntegralFunction{true}, v2ujac, lb, ub) _f = f.f prototype = similar(f.integrand_prototype) vol = prod((ub - lb) / 2) # just to get the type of the jacobian determinant BatchIntegralFunction{true}(prototype * vol, max_batch = f.max_batch) do y, v, p u, jac = substitute_bv(v2ujac, v, lb, ub) _prototype = similar(prototype, size(y)) - _f(_prototype, u, p) - for (i, j) in zip(axes(y)[end], jac) - for iy in CartesianIndices(axes(y)[begin:(end - 1)]) - y[iy, i] = j * _prototype[iy, i] - end - end + _y = _evaluate!(_f, _prototype, u, p) + y .= _y .* reshape(jac, ntuple(d -> d == ndims(y) ? length(jac) : 1, ndims(y))) return end end +# we need this function for autodiff compatibility where internal buffers need special types +function _evaluate!(f, y, u, p) + f(y, u, p) + return y +end + # specific changes of variables function u2t(lb, ub) From ce8be8366070bc5782b9e0bc2a954e2e456449e6 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 18:06:42 -0500 Subject: [PATCH 12/26] apply format --- ext/IntegralsForwardDiffExt.jl | 3 ++- ext/IntegralsZygoteExt.jl | 1 - test/derivative_tests.jl | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 3f6548b..1396bda 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -3,7 +3,8 @@ using Integrals isdefined(Base, :get_extension) ? (using ForwardDiff) : (using ..ForwardDiff) ### Forward-Mode AD Intercepts -function Integrals._evaluate!(f, y, u, p::Union{D, AbstractArray{<:D}}) where {T, V, P, D <: ForwardDiff.Dual{T, V, P}} +function Integrals._evaluate!(f, y, u, + p::Union{D, AbstractArray{<:D}}) where {T, V, P, D <: ForwardDiff.Dual{T, V, P}} dy = similar(y, replace_dualvaltype(eltype(p), eltype(y))) f(dy, u, p) return dy diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index dad0874..fd48dc9 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -43,7 +43,6 @@ function ChainRulesCore.rrule(::typeof(Integrals._evaluate!), f, y, u, p) out, Δ -> (NoTangent(), NoTangent(), back(Δ)...) end - function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, sensealg, domain, p; kwargs...) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index ba22e8f..9f0cd35 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -18,8 +18,8 @@ alg_req = Dict( 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), 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)) From eaaa7d29ce863f50fc08d75359186f30f7126cea Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 21:18:24 -0500 Subject: [PATCH 13/26] format --- src/algorithms.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 8b32fe8..e66d06e 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -86,7 +86,9 @@ struct VEGAS{S} <: SciMLBase.AbstractIntegralAlgorithm debug::Bool seed::S end -VEGAS(; nbins = 100, ncalls = 1000, debug = false, seed = nothing) = VEGAS(nbins, ncalls, debug, seed) +function VEGAS(; nbins = 100, ncalls = 1000, debug = false, seed = nothing) + VEGAS(nbins, ncalls, debug, seed) +end """ GaussLegendre{C, N, W} From 3a2edd9aebb4ef82c8b80efb01191ecf932ed3f6 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 22:07:44 -0500 Subject: [PATCH 14/26] add a buffer argument to QuadGKJL HCubatureJL --- src/Integrals.jl | 47 ++++++++++++++++++++++++----------------------- src/algorithms.jl | 22 +++++++++++++++------- 2 files changed, 39 insertions(+), 30 deletions(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index 51680bb..bbb71bc 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -92,23 +92,30 @@ function __solve(cache::IntegralCache, alg::ChangeOfVariables, sensealg, udomain return __solve(_cache, alg.alg, sensealg, vdomain, p; kwargs...) end -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 - NT = Base.promote_op(nrm, RT) - return DT, RT, NT +function get_prototype(prob::IntegralProblem) + f = prob.f + f.integrand_prototype !== nothing && return f.integrand_prototype + isinplace(f) && throw(ArgumentError("in-place integral functions require a prototype")) + lb, ub = prob.domain + mid = (lb + ub) / 2 + p = prob.p + if f isa BatchIntegralFunction + mid isa Number ? f(eltype(mid)[], p) : + f(Matrix{eltype(mid)}(undef, length(mid), 0), p) + else + f(mid, p) + end end + function init_cacheval(alg::QuadGKJL, prob::IntegralProblem) - lb, ub = map(first, prob.domain) - DT, RT, NT = quadgk_prob_types(prob.f, lb, ub, prob.p, alg.norm) - return (isconcretetype(RT) ? QuadGK.alloc_segbuf(DT, RT, NT) : nothing) -end -function refresh_cacheval(cacheval, alg::QuadGKJL, prob) - lb, ub = map(first, prob.domain) - DT, RT, NT = quadgk_prob_types(prob.f, lb, ub, prob.p, alg.norm) - isconcretetype(RT) || return nothing - T = QuadGK.Segment{DT, RT, NT} - return (cacheval isa Vector{T} ? cacheval : QuadGK.alloc_segbuf(DT, RT, NT)) + alg.buffer === nothing && return + lb, ub = prob.domain + mid = (lb + ub) / 2 + prototype = get_prototype(prob) * mid + TX = typeof(mid) + TI = typeof(prototype) + TE = typeof(alg.norm(prototype)) + QuadGK.alloc_segbuf(TX, TI, TE) end function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p; @@ -175,14 +182,8 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p end function init_cacheval(alg::HCubatureJL, prob::IntegralProblem) - lb, ub = prob.domain - lb isa Number && return - if isinplace(prob) - f = u -> prob.f.integrand_prototype - else - f = u -> prob.f(u, prob.p) - end - return HCubature.hcubature_buffer(f, lb, ub; norm = alg.norm) + alg.buffer === nothing && return + HCubature.hcubature_buffer(x -> get_prototype(prob), lb, ub; norm = alg.norm) end function __solvebp_call(cache::IntegralCache, alg::HCubatureJL, sensealg, domain, p; reltol = 1e-8, abstol = 1e-8, diff --git a/src/algorithms.jl b/src/algorithms.jl index e66d06e..4841f0a 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -1,10 +1,12 @@ """ - QuadGKJL(; order = 7, norm=norm) + QuadGKJL(; order = 7, norm = norm, buffer = nothing) One-dimensional Gauss-Kronrod integration from QuadGK.jl. This method also takes the optional arguments `order` and `norm`. Which are the order of the integration rule -and the norm for calculating the error, respectively +and the norm for calculating the error, respectively. +Lastly, the `buffer` keyword, if set, will allocate a buffer to reuse +for multiple integrals, which may require evaluating the integrand. ## References @@ -18,20 +20,23 @@ pages={1133--1145}, year={1997} } """ -struct QuadGKJL{F} <: SciMLBase.AbstractIntegralAlgorithm where {F} +struct QuadGKJL{F, B} <: SciMLBase.AbstractIntegralAlgorithm order::Int norm::F + buffer::B end -QuadGKJL(; order = 7, norm = norm) = QuadGKJL(order, norm) +QuadGKJL(; order = 7, norm = norm, buffer = nothing) = QuadGKJL(order, norm, buffer) """ - HCubatureJL(; norm=norm, initdiv=1) + HCubatureJL(; norm=norm, initdiv=1, buffer = nothing) Multidimensional "h-adaptive" integration from HCubature.jl. This method also takes the optional arguments `initdiv` and `norm`. Which are the initial number of segments each dimension of the integration domain is divided into, and the norm for calculating the error, respectively. +Lastly, the `buffer` keyword, if set, will allocate a buffer to reuse +for multiple integrals, which may require evaluating the integrand. ## References @@ -46,11 +51,14 @@ year={1980}, publisher={Elsevier} } """ -struct HCubatureJL{F} <: SciMLBase.AbstractIntegralAlgorithm where {F} +struct HCubatureJL{F, B} <: SciMLBase.AbstractIntegralAlgorithm initdiv::Int norm::F + buffer::B +end +function HCubatureJL(; initdiv = 1, norm = norm, buffer = nothing) + HCubatureJL(initdiv, norm, buffer) end -HCubatureJL(; initdiv = 1, norm = norm) = HCubatureJL(initdiv, norm) """ VEGAS(; nbins = 100, ncalls = 1000, debug=false, seed = nothing) From e314c389d02196dde3aa946b0650220e4eaeddfc Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 22:47:12 -0500 Subject: [PATCH 15/26] use get_prototype utility --- ext/IntegralsCubaExt.jl | 64 ++++++++++++++---------------- ext/IntegralsCubatureExt.jl | 68 +++++++++++++++----------------- ext/IntegralsMCIntegrationExt.jl | 41 +++++++++---------- src/Integrals.jl | 3 +- 4 files changed, 84 insertions(+), 92 deletions(-) diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl index c67be23..82bad2a 100644 --- a/ext/IntegralsCubaExt.jl +++ b/ext/IntegralsCubaExt.jl @@ -18,8 +18,12 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori throw(ArgumentError("Cuba.jl only supports real-valued integrands")) # we could support other types by multiplying by the jacobian determinant at the end - if prob.f isa BatchIntegralFunction - nvec = min(maxiters, prob.f.max_batch) + f = prob.f + prototype = Integrals.get_prototype(prob) + if f isa BatchIntegralFunction + fsize = size(prototype)[begin:(end - 1)] + ncomp = prod(fsize) + nvec = min(maxiters, f.max_batch) # nvec == 1 in Cuba will change vectors to matrices, so we won't support it when # batching nvec > 1 || @@ -33,26 +37,21 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori scale = x -> scale_x!(view(_x, :, 1:size(x, 2)), ub, lb, x) end - if isinplace(prob) - y = prob.f.integrand_prototype - fsize = size(y)[begin:(end - 1)] - ax = map(_ -> (:), fsize) - f = let y = similar(y, fsize..., nvec) - function (x, dx) - dy = @view(y[ax..., begin:(begin + size(dx, 2) - 1)]) - prob.f(dy, scale(x), p) - dx .= reshape(dy, :, size(dx, 2)) .* vol + if isinplace(f) + ax = ntuple(_ -> (:), length(fsize)) + _f = let y_ = similar(prototype, fsize..., nvec) + function (u, _y) + y = @view(y_[ax..., begin:(begin + size(_y, 2) - 1)]) + f(y, scale(u), p) + _y .= reshape(y, size(_y)) .* vol end end else - y = mid isa Number ? prob.f(typeof(mid)[], p) : - prob.f(Matrix{eltype(mid)}(undef, length(mid), 0), p) - fsize = size(y)[begin:(end - 1)] - f = (x, dx) -> dx .= reshape(prob.f(scale(x), p), :, size(dx, 2)) .* vol + _f = (u, y) -> y .= reshape(f(scale(u), p), size(y)) .* vol end - ncomp = prod(fsize) else nvec = 1 + ncomp = length(prototype) if mid isa Real scale = x -> scale_x(ub, lb, only(x)) @@ -61,36 +60,33 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori scale = x -> scale_x!(_x, ub, lb, x) end - if isinplace(prob) - y = prob.f.integrand_prototype - f = let y = similar(y) - (x, dx) -> begin - prob.f(y, scale(x), p) - dx .= vec(y) .* vol + if isinplace(f) + _f = let y = similar(prototype) + (u, _y) -> begin + f(y, scale(u), p) + _y .= vec(y) .* vol end end else - y = prob.f(mid, p) - f = (x, dx) -> dx .= Iterators.flatten(prob.f(scale(x), p)) .* vol + _f = (u, y) -> y .= Iterators.flatten(f(scale(u), p)) .* vol end - ncomp = length(y) end out = if alg isa CubaVegas - Cuba.vegas(f, ndim, ncomp; rtol = reltol, + Cuba.vegas(_f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, seed = alg.seed, minevals = alg.minevals, nstart = alg.nstart, nincrease = alg.nincrease, gridno = alg.gridno) elseif alg isa CubaSUAVE - Cuba.suave(f, ndim, ncomp; rtol = reltol, + Cuba.suave(_f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, seed = alg.seed, minevals = alg.minevals, nnew = alg.nnew, nmin = alg.nmin, flatness = alg.flatness) elseif alg isa CubaDivonne - Cuba.divonne(f, ndim, ncomp; rtol = reltol, + Cuba.divonne(_f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, seed = alg.seed, minevals = alg.minevals, @@ -98,26 +94,26 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori maxpass = alg.maxpass, border = alg.border, maxchisq = alg.maxchisq, mindeviation = alg.mindeviation) elseif alg isa CubaCuhre - Cuba.cuhre(f, ndim, ncomp; rtol = reltol, + Cuba.cuhre(_f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, minevals = alg.minevals, key = alg.key) end # out.integral is a Vector{Float64}, but we want to return it to the shape of the integrand - val = if prob.f isa BatchIntegralFunction - if y isa AbstractVector + val = if f isa BatchIntegralFunction + if prototype isa AbstractVector out.integral[1] else reshape(out.integral, fsize) end else - if y isa Real + if prototype isa Real out.integral[1] - elseif y isa AbstractVector + elseif prototype isa AbstractVector out.integral else - reshape(out.integral, size(y)) + reshape(out.integral, size(prototype)) end end diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl index ed6ae76..57327cc 100644 --- a/ext/IntegralsCubatureExt.jl +++ b/ext/IntegralsCubatureExt.jl @@ -13,49 +13,43 @@ function Integrals.__solvebp_call(prob::IntegralProblem, mid = (lb + ub) / 2 # we get to pick fdim or not based on the IntegralFunction and its output dimensions - prototype = if prob.f isa BatchIntegralFunction - isinplace(prob.f) ? prob.f.integrand_prototype : - mid isa Number ? prob.f(eltype(mid)[], p) : - prob.f(Matrix{eltype(mid)}(undef, length(mid), 0), p) - else - # we evaluate the oop function to decide whether the output should be vectorized - isinplace(prob.f) ? prob.f.integrand_prototype : prob.f(mid, p) - end + f = prob.f + prototype = Integrals.get_prototype(prob) @assert eltype(prototype)<:Real "Cubature.jl is only compatible with real-valued integrands" - if prob.f isa BatchIntegralFunction + if f isa BatchIntegralFunction if prototype isa AbstractVector # this branch could be omitted since the following one should work similarly - if isinplace(prob) + if isinplace(f) # dx is a Vector, but we provide the integrand a vector of the same type as # y, which needs to be resized since the number of batch points changes. - f = let y = similar(prototype) + _f = let y = similar(prototype) (u, v) -> begin resize!(y, length(v)) - prob.f(y, u, p) + f(y, u, p) v .= y end end else - f = (u, v) -> (v .= prob.f(u, p)) + _f = (u, v) -> (v .= f(u, p)) end if mid isa Number if alg isa CubatureJLh - val, err = Cubature.hquadrature_v(f, lb, ub; + val, err = Cubature.hquadrature_v(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) else - val, err = Cubature.pquadrature_v(f, lb, ub; + val, err = Cubature.pquadrature_v(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) end else if alg isa CubatureJLh - val, err = Cubature.hcubature_v(f, lb, ub; + val, err = Cubature.hcubature_v(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) else - val, err = Cubature.pcubature_v(f, lb, ub; + val, err = Cubature.pcubature_v(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) end @@ -63,37 +57,37 @@ function Integrals.__solvebp_call(prob::IntegralProblem, elseif prototype isa AbstractArray fsize = size(prototype)[begin:(end - 1)] fdim = prod(fsize) - if isinplace(prob) + if isinplace(f) # dx is a Matrix, but to provide a buffer of the same type as y, we make # would like to make views of a larger buffer, but CubatureJL doesn't set # a hard limit for max_batch, so we allocate a new buffer with the needed size - f = let fsize = fsize + _f = let fsize = fsize (u, v) -> begin y = similar(prototype, fsize..., size(v, 2)) - prob.f(y, u, p) + f(y, u, p) v .= reshape(y, fdim, size(v, 2)) end end else - f = (u, v) -> (v .= reshape(prob.f(u, p), fdim, size(v, 2))) + _f = (u, v) -> (v .= reshape(f(u, p), fdim, size(v, 2))) end if mid isa Number if alg isa CubatureJLh - val_, err = Cubature.hquadrature_v(fdim, f, lb, ub; + val_, err = Cubature.hquadrature_v(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) else - val_, err = Cubature.pquadrature_v(fdim, f, lb, ub; + val_, err = Cubature.pquadrature_v(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) end else if alg isa CubatureJLh - val_, err = Cubature.hcubature_v(fdim, f, lb, ub; + val_, err = Cubature.hcubature_v(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) else - val_, err = Cubature.pcubature_v(fdim, f, lb, ub; + val_, err = Cubature.pcubature_v(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) end @@ -105,24 +99,24 @@ function Integrals.__solvebp_call(prob::IntegralProblem, else if prototype isa Real # no inplace in this case, since the integrand_prototype would be mutable - f = u -> prob.f(u, p) + _f = u -> f(u, p) if lb isa Number if alg isa CubatureJLh - val, err = Cubature.hquadrature(f, lb, ub; + val, err = Cubature.hquadrature(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) else - val, err = Cubature.pquadrature(f, lb, ub; + val, err = Cubature.pquadrature(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) end else if alg isa CubatureJLh - val, err = Cubature.hcubature(f, lb, ub; + val, err = Cubature.hcubature(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) else - val, err = Cubature.pcubature(f, lb, ub; + val, err = Cubature.pcubature(_f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) end @@ -131,29 +125,29 @@ function Integrals.__solvebp_call(prob::IntegralProblem, fsize = size(prototype) fdim = length(prototype) if isinplace(prob) - f = let y = similar(prototype) - (u, v) -> (prob.f(y, u, p); v .= vec(y)) + _f = let y = similar(prototype) + (u, v) -> (f(y, u, p); v .= vec(y)) end else - f = (u, v) -> (v .= vec(prob.f(u, p))) + _f = (u, v) -> (v .= vec(f(u, p))) end if mid isa Number if alg isa CubatureJLh - val_, err = Cubature.hquadrature(fdim, f, lb, ub; + val_, err = Cubature.hquadrature(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) else - val_, err = Cubature.pquadrature(fdim, f, lb, ub; + val_, err = Cubature.pquadrature(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) end else if alg isa CubatureJLh - val_, err = Cubature.hcubature(fdim, f, lb, ub; + val_, err = Cubature.hcubature(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) else - val_, err = Cubature.pcubature(fdim, f, lb, ub; + val_, err = Cubature.pcubature(fdim, _f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters, error_norm = alg.error_norm) end diff --git a/ext/IntegralsMCIntegrationExt.jl b/ext/IntegralsMCIntegrationExt.jl index 994b366..4255bac 100644 --- a/ext/IntegralsMCIntegrationExt.jl +++ b/ext/IntegralsMCIntegrationExt.jl @@ -12,41 +12,42 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, tmp = vec(collect(mid)) var = Continuous(vec([tuple(a, b) for (a, b) in zip(lb, ub)])) - if prob.f isa BatchIntegralFunction + f = prob.f + if f isa BatchIntegralFunction error("VEGASMC doesn't support batching. See https://github.com/numericalEFT/MCIntegration.jl/issues/29") else - f0 = if isinplace(prob) - _f0 = similar(prob.f.integrand_prototype) - f_ = (x, f, c) -> begin - n = 0 - for v in x - tmp[n += 1] = first(v) + prototype = Integrals.get_prototype(prob) + if isinplace(prob) + _f = let y = similar(prototype) + (u, _y, c) -> begin + n = 0 + for v in u + tmp[n += 1] = first(v) + end + f(y, _oftype(mid, tmp), p) + _y .= vec(y) end - prob.f(_f0, _oftype(mid, tmp), p) - f .= vec(_f0) end - _f0 else - f_ = (x, c) -> begin + _f = (u, c) -> begin n = 0 - for v in x + for v in u tmp[n += 1] = first(v) end - fx = prob.f(_oftype(mid, tmp), p) - fx isa AbstractArray ? vec(fx) : fx + y = f(_oftype(mid, tmp), p) + y isa AbstractArray ? vec(y) : y end - prob.f(mid, p) end - dof = ones(Int, length(f0)) # each composite Continuous var gets 1 dof - res = integrate(f_; var, dof, inplace = isinplace(prob), type = eltype(f0), + dof = ones(Int, length(prototype)) # each composite Continuous var gets 1 dof + res = integrate(_f; var, dof, inplace = isinplace(prob), type = eltype(prototype), solver = :vegasmc, niter = maxiters, verbose = -2, print = -2, alg.kws...) # the package itself is not type-stable - out::typeof(f0), err::typeof(f0), chi2 = if f0 isa Number + out::typeof(prototype), err::typeof(prototype), chi2 = if prototype isa Number map(only, (res.mean, res.stdev, res.chi2)) else - map(a -> reshape(a, size(f0)), (res.mean, res.stdev, res.chi2)) + map(a -> reshape(a, size(prototype)), (res.mean, res.stdev, res.chi2)) end - chi::typeof(f0) = map(sqrt, chi2) + chi::typeof(prototype) = map(sqrt, chi2) SciMLBase.build_solution(prob, alg, out, err, chi = chi, retcode = ReturnCode.Success) end diff --git a/src/Integrals.jl b/src/Integrals.jl index bbb71bc..eacd8d8 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -183,7 +183,8 @@ end function init_cacheval(alg::HCubatureJL, prob::IntegralProblem) alg.buffer === nothing && return - HCubature.hcubature_buffer(x -> get_prototype(prob), lb, ub; norm = alg.norm) + prototype = get_prototype(prob) + HCubature.hcubature_buffer(x -> prototype, lb, ub; norm = alg.norm) end function __solvebp_call(cache::IntegralCache, alg::HCubatureJL, sensealg, domain, p; reltol = 1e-8, abstol = 1e-8, From 3311201eefbf633bdc22552758c2635c9abfae62 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 25 Feb 2024 22:47:40 -0500 Subject: [PATCH 16/26] fix issue with forwarddiff cache type --- ext/IntegralsForwardDiffExt.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 1396bda..f09309c 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -26,9 +26,16 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, IntegralFunction{true}(cache.f.f, dprototype) end prob = Integrals.build_problem(cache) - dprob = remake(prob, f = df) - dcache = init( - dprob, alg; sensealg = sensealg, kwargs...) + dcache = Integrals.IntegralCache(cache.iip, + df, + domain, + p, + cache.prob_kwargs, + alg, + sensealg, + cache.kwargs, + cache.cacheval + ) Integrals.__solvebp_call(dcache, alg, sensealg, domain, p; kwargs...) else Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...) From 20b1f708845769c5a2db5f14bd8ce14842bf0151 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 26 Feb 2024 07:05:52 -0500 Subject: [PATCH 17/26] pass tests --- test/gaussian_quadrature_tests.jl | 7 ++++--- test/nested_ad_tests.jl | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/test/gaussian_quadrature_tests.jl b/test/gaussian_quadrature_tests.jl index d3a2387..cd243e2 100644 --- a/test/gaussian_quadrature_tests.jl +++ b/test/gaussian_quadrature_tests.jl @@ -50,8 +50,9 @@ prob = IntegralProblem(f, -5, 3, 3.3) alg = GaussLegendre() sol = solve(prob, alg) @test isnothing(sol.chi) -@test sol.alg === alg -@test sol.prob === prob +# These tests don't pass with ChangeOfVariables +# @test sol.alg === alg +# @test sol.prob === prob @test isnothing(sol.resid) @test SciMLBase.successful_retcode(sol) @test sol.u ≈ -exp(3) * 3.3 + 3.3 / exp(5) - 40 + cos(5) - cos(3) @@ -87,7 +88,7 @@ alg = GaussLegendre(subintervals = 1) alg = GaussLegendre(subintervals = 17) @test sol.u ≈ sqrt(π) / 2 -# Make sure broadcasting correctly handles the argument p +# Make sure broadcasting correctly handles the argument p f = (x, p) -> 1 + x + x^p[1] - cos(x * p[2]) + exp(x) * p[3] p = [0.3, 1.3, -0.5] prob = IntegralProblem(f, 2.0, 6.3, p) diff --git a/test/nested_ad_tests.jl b/test/nested_ad_tests.jl index d4506aa..65840ea 100644 --- a/test/nested_ad_tests.jl +++ b/test/nested_ad_tests.jl @@ -8,8 +8,8 @@ function my_integration(p) return solve(my_problem, CubatureJLh(), reltol = 1e-3, abstol = 1e-3) # Errors end 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) == +@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] From 7a9f3a90a72b772f5ff93ffc84749cdd6f6eac71 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 16:39:11 -0500 Subject: [PATCH 18/26] fix inference failure in GaussLegendre --- ext/IntegralsFastGaussQuadratureExt.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/ext/IntegralsFastGaussQuadratureExt.jl b/ext/IntegralsFastGaussQuadratureExt.jl index 1dc5c39..704145d 100644 --- a/ext/IntegralsFastGaussQuadratureExt.jl +++ b/ext/IntegralsFastGaussQuadratureExt.jl @@ -15,7 +15,13 @@ Integrals.gausslegendre(n) = FastGaussQuadrature.gausslegendre(n) function gauss_legendre(f::F, p, lb, ub, nodes, weights) where {F} scale = (ub - lb) / 2 shift = (lb + ub) / 2 - I = mapreduce((w, x) -> w * f(scale * x + shift, p), +, weights, nodes) + # TODO reuse Integrals.evalrule instead + x0, xs = Iterators.peel(nodes) + w0, ws = Iterators.peel(weights) + I = w0 * f(scale * x0 + shift, p) + for (w, x) in zip(ws, xs) + I += w * f(scale * x + shift, p) + end return scale * I end function composite_gauss_legendre(f::F, p, lb, ub, nodes, weights, subintervals) where {F} From c39fba5473ea81be8795096095a7b8b0204c4f3d Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 16:40:53 -0500 Subject: [PATCH 19/26] update test deps --- Project.toml | 4 +--- test/qa.jl | 1 + test/runtests.jl | 1 - 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 02f55d0..552f89b 100644 --- a/Project.toml +++ b/Project.toml @@ -47,7 +47,6 @@ HCubature = "1.5.2" LinearAlgebra = "1.10" MCIntegration = "0.4.2" MonteCarloIntegration = "0.2" -Pkg = "1.10" QuadGK = "2.9" Random = "1.10" Reexport = "1.0" @@ -69,11 +68,10 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "Arblib", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature", "MCIntegration"] +test = ["Aqua", "Arblib", "StaticArrays", "FiniteDiff", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature", "MCIntegration"] diff --git a/test/qa.jl b/test/qa.jl index ccf1183..9992475 100644 --- a/test/qa.jl +++ b/test/qa.jl @@ -1,4 +1,5 @@ using Integrals, Aqua +using Test @testset "Aqua" begin Aqua.find_persistent_tasks_deps(Integrals) Aqua.test_ambiguities(Integrals, recursive = false) diff --git a/test/runtests.jl b/test/runtests.jl index 090b8fd..bfdbfe9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,3 @@ -using Pkg using SafeTestsets using Test From 519b8c826207d67c6ba3e459923a388c01f7adf4 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 16:41:20 -0500 Subject: [PATCH 20/26] restore gaussian quad tests --- test/gaussian_quadrature_tests.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/gaussian_quadrature_tests.jl b/test/gaussian_quadrature_tests.jl index cd243e2..ad8eb41 100644 --- a/test/gaussian_quadrature_tests.jl +++ b/test/gaussian_quadrature_tests.jl @@ -50,9 +50,8 @@ prob = IntegralProblem(f, -5, 3, 3.3) alg = GaussLegendre() sol = solve(prob, alg) @test isnothing(sol.chi) -# These tests don't pass with ChangeOfVariables -# @test sol.alg === alg -# @test sol.prob === prob +@test sol.alg === alg +@test sol.prob === prob @test isnothing(sol.resid) @test SciMLBase.successful_retcode(sol) @test sol.u ≈ -exp(3) * 3.3 + 3.3 / exp(5) - 40 + cos(5) - cos(3) From 5d1bebd19ed21175f09788540cd244c747e4819f Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 16:53:10 -0500 Subject: [PATCH 21/26] test that solution contains original problem --- test/inf_integral_tests.jl | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl index 2d28a9e..5de996f 100644 --- a/test/inf_integral_tests.jl +++ b/test/inf_integral_tests.jl @@ -7,12 +7,12 @@ abstol = 1e-3 alg_req = Dict( 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 = 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), CubatureJLh() => ( - nout = Inf, allows_batch = true, min_dim = 1, + 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, @@ -206,11 +206,16 @@ end # two distinct semi-infinite transformations should still work as expected (; f, domain, solution) = problems[8] prob = IntegralProblem(f, domain) - cache = init(prob, QuadGKJL(); abstol, reltol) - sol = solve!(cache).u - @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) + alg = QuadGKJL() + cache = init(prob, alg; abstol, reltol) + sol = solve!(cache) + @test abs(only(sol.u) - solution) < max(abstol, reltol * abs(solution)) + @test sol.prob == IntegralProblem(f, domain) + @test sol.alg == alg (; domain, solution) = problems[9] cache.domain = domain - sol = solve!(cache).u - @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution)) + sol = solve!(cache) + @test abs(only(sol.u) - solution) < max(abstol, reltol * abs(solution)) + @test sol.prob == IntegralProblem(f, domain) + @test sol.alg == alg end From 95e27e23908458e37753cef3d8b07c11a47dd1d4 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 17:53:21 -0500 Subject: [PATCH 22/26] attempt to fix derivative tests --- ext/IntegralsZygoteExt.jl | 16 ++++++++++++---- src/Integrals.jl | 10 +++++++++- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index fd48dc9..bb5d457 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -14,11 +14,19 @@ ChainRulesCore.@non_differentiable Integrals.checkkwargs(kwargs...) ChainRulesCore.@non_differentiable Integrals.isinplace(f, args...) # fixes #99 ChainRulesCore.@non_differentiable Integrals.init_cacheval(alg, prob) -function ChainRulesCore.rrule(::typeof(Integrals.transformation_if_inf), f, domain) - function transformation_if_inf_pullback(Δ) - return NoTangent(), Δ... +function ChainRulesCore.rrule(::typeof(Integrals.__solve), cache::Integrals.IntegralCache, + alg::Integrals.ChangeOfVariables, sensealg, udomain, p; + kwargs...) + _cache, vdomain = Integrals._change_variables(cache, alg, sensealg, udomain, p) + sol, back = Zygote.pullback((args...) -> Integrals.__solve(args...; kwargs...), + _cache, alg.alg, sensealg, vdomain, p) + function change_of_variables_pullback(Δ) + return (NoTangent(), back(Δ)...) end - return Integrals.transformation_if_inf(f, domain), transformation_if_inf_pullback + prob = Integrals.build_problem(cache) + _sol = SciMLBase.build_solution( + prob, alg.alg, sol.u, sol.resid, chi = sol.chi, retcode = sol.retcode, stats = sol.stats) + return _sol, change_of_variables_pullback end # we will need to implement the following adjoints when we compute ∂f/∂u diff --git a/src/Integrals.jl b/src/Integrals.jl index eacd8d8..844eaf1 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -78,6 +78,14 @@ end function __solve(cache::IntegralCache, alg::ChangeOfVariables, sensealg, udomain, p; kwargs...) + _cache, vdomain = _change_variables(cache, alg, sensealg, udomain, p) + sol = __solve(_cache, alg.alg, sensealg, vdomain, p; kwargs...) + prob = build_problem(cache) + return SciMLBase.build_solution( + prob, alg.alg, sol.u, sol.resid, chi = sol.chi, retcode = sol.retcode, stats = sol.stats) +end + +function _change_variables(cache, alg, sensealg, udomain, p) cacheval = cache.cacheval.alg g, vdomain = alg.fu2gv(cache.f, udomain) _cache = IntegralCache(Val(isinplace(g)), @@ -89,7 +97,7 @@ function __solve(cache::IntegralCache, alg::ChangeOfVariables, sensealg, udomain sensealg, cache.kwargs, cacheval) - return __solve(_cache, alg.alg, sensealg, vdomain, p; kwargs...) + return _cache, vdomain end function get_prototype(prob::IntegralProblem) From 9a39c80ff576dfe24d562c0f7ebc3ff5c8b8bd0a Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 17:58:32 -0500 Subject: [PATCH 23/26] avoid repetition of build_solution --- ext/IntegralsArblibExt.jl | 5 +---- src/Integrals.jl | 5 +---- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/ext/IntegralsArblibExt.jl b/ext/IntegralsArblibExt.jl index 6e9d5ed..0eac48a 100644 --- a/ext/IntegralsArblibExt.jl +++ b/ext/IntegralsArblibExt.jl @@ -21,16 +21,13 @@ function Integrals.__solvebp_call( 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) end + SciMLBase.build_solution(prob, alg, val, get_radius(val), retcode = ReturnCode.Success) end function get_radius(ball) diff --git a/src/Integrals.jl b/src/Integrals.jl index 844eaf1..1627304 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -156,7 +156,6 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p end val, err = quadgk(_f, 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 prototype = f(typeof(mid)[], p) _f = if prototype isa AbstractVector @@ -170,7 +169,6 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p end val, err = quadgk(_f, 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) end else if isinplace(f) @@ -179,14 +177,13 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p 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 _f = u -> f(u, p) val, err = quadgk(_f, 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) end end + SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) end function init_cacheval(alg::HCubatureJL, prob::IntegralProblem) From 74e4defdb37e235d34c88f50a99bb380cbefcb46 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 18:06:52 -0500 Subject: [PATCH 24/26] get_prototype for HCubatureJL --- src/Integrals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index 1627304..e7d53ac 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -217,7 +217,7 @@ function __solvebp_call(cache::IntegralCache, alg::HCubatureJL, sensealg, domain rtol = reltol, atol = abstol, buffer = cache.cacheval, maxevals = maxiters, norm = alg.norm, initdiv = alg.initdiv) else - ret = _f((lb + ub) / 2) * (prod(ub - lb) / 2) # this calculation for type stability with vector endpoints + ret = get_prototype(prob) * (prod(ub - lb) / 2) # this calculation for type stability with vector endpoints hcubature(_f, lb, ub; rtol = reltol, atol = abstol, buffer = cache.cacheval, maxevals = maxiters, norm = alg.norm, From ded5ee090ea715f622efd10e65401704c2e41bd8 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 2 Mar 2024 18:28:20 -0500 Subject: [PATCH 25/26] fix test after rebase --- test/inf_integral_tests.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl index 5de996f..ee41083 100644 --- a/test/inf_integral_tests.jl +++ b/test/inf_integral_tests.jl @@ -204,7 +204,9 @@ end @testset "Caching interface" begin # two distinct semi-infinite transformations should still work as expected - (; f, domain, solution) = problems[8] + f = (x, p) -> pdf(Normal(0.00, 1.00), x) + domain = (0.0, -Inf) + solution = -0.5 prob = IntegralProblem(f, domain) alg = QuadGKJL() cache = init(prob, alg; abstol, reltol) @@ -212,7 +214,8 @@ end @test abs(only(sol.u) - solution) < max(abstol, reltol * abs(solution)) @test sol.prob == IntegralProblem(f, domain) @test sol.alg == alg - (; domain, solution) = problems[9] + domain = (-Inf, 0.0) + solution = 0.5 cache.domain = domain sol = solve!(cache) @test abs(only(sol.u) - solution) < max(abstol, reltol * abs(solution)) From e6e9a5b29cd42b56659f971fd8bc35c334f52d6d Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 3 Mar 2024 12:23:21 -0500 Subject: [PATCH 26/26] fix buffer bugs --- src/Integrals.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index e7d53ac..5edee6d 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -117,7 +117,7 @@ end function init_cacheval(alg::QuadGKJL, prob::IntegralProblem) alg.buffer === nothing && return - lb, ub = prob.domain + lb, ub = map(first, prob.domain) mid = (lb + ub) / 2 prototype = get_prototype(prob) * mid TX = typeof(mid) @@ -189,6 +189,7 @@ end function init_cacheval(alg::HCubatureJL, prob::IntegralProblem) alg.buffer === nothing && return prototype = get_prototype(prob) + lb, ub = map(x -> x isa Number ? tuple(x) : x, prob.domain) HCubature.hcubature_buffer(x -> prototype, lb, ub; norm = alg.norm) end function __solvebp_call(cache::IntegralCache, alg::HCubatureJL, sensealg, domain, p;