From 68ae9d6ea68b9804b4787e636a189288aa41b886 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:24:24 -0400 Subject: [PATCH 01/27] bump dependencies --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 5f66712c..68db8e53 100644 --- a/Project.toml +++ b/Project.toml @@ -31,11 +31,11 @@ Distributions = "0.23, 0.24, 0.25" FastGaussQuadrature = "0.5" ForwardDiff = "0.10" HCubature = "1.4" -MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3" +MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3, 0.1" QuadGK = "2.5" Reexport = "0.2, 1.0" Requires = "1" -SciMLBase = "1.98" +SciMLBase = "2.5" Zygote = "0.4.22, 0.5, 0.6" julia = "1.6" From bb98991e867b81c5dd0fbe983f4798e7310c9bf3 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:25:38 -0400 Subject: [PATCH 02/27] update Integrals to SciMLBasev2 --- src/Integrals.jl | 75 ++++++++++++++++++++++++---------------- src/algorithms.jl | 6 +++- src/common.jl | 62 +++++++++++++++++++++------------ src/infinity_handling.jl | 54 ++++++++++++++++++++++------- src/quadrules.jl | 9 ++--- 5 files changed, 138 insertions(+), 68 deletions(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index 02b5bc5d..4d8a28da 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -71,45 +71,47 @@ function quadgk_prob_types(f, lb::T, ub::T, p, nrm) where {T} return DT, RT, NT end function init_cacheval(alg::QuadGKJL, prob::IntegralProblem) - DT, RT, NT = quadgk_prob_types(prob.f, prob.lb, prob.ub, prob.p, alg.norm) + lb, ub = 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) - DT, RT, NT = quadgk_prob_types(prob.f, prob.lb, prob.ub, prob.p, alg.norm) + lb, ub = 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)) end -function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, lb, ub, p; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int)) +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 if isinplace(prob) || lb isa AbstractArray || ub isa AbstractArray error("QuadGKJL only accepts one-dimensional quadrature problems.") end - @assert prob.batch == 0 - @assert prob.nout == 1 + @assert prob.f isa IntegralFunction - p = p f = x -> prob.f(x, 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, QuadGKJL(), val, err, retcode = ReturnCode.Success) end -function __solvebp_call(prob::IntegralProblem, alg::HCubatureJL, sensealg, lb, ub, p; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int)) - p = p +function __solvebp_call(prob::IntegralProblem, alg::HCubatureJL, sensealg, domain, p; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int)) + lb, ub = domain + @assert prob.f isa IntegralFunction if isinplace(prob) - dx = zeros(eltype(lb), prob.nout) - f = x -> (prob.f(dx, x, prob.p); dx) + # 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) else f = x -> prob.f(x, prob.p) end - @assert prob.batch == 0 if lb isa Number val, err = hquadrature(f, lb, ub; @@ -123,30 +125,45 @@ function __solvebp_call(prob::IntegralProblem, alg::HCubatureJL, sensealg, lb, u SciMLBase.build_solution(prob, HCubatureJL(), val, err, retcode = ReturnCode.Success) end -function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, lb, ub, p; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int)) - p = p - @assert prob.nout == 1 - if prob.batch == 0 +function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int)) + lb, ub = domain + mid = (lb + ub) / 2 + if prob.f isa BatchIntegralFunction if isinplace(prob) - dx = zeros(eltype(lb), prob.nout) - f = x -> (prob.f(dx, x, p); dx[1]) + y = similar(prob.f.integrand_prototype, + size(prob.f.integrand_prototype)[begin:(end - 1)]..., + prob.f.max_batch) + f = x -> (prob.f(y, x', p); vec(y)) else - f = x -> prob.f(x, prob.p) + y = prob.f(mid isa Number ? typeof(mid)[] : + Matrix{eltype(mid)}(undef, length(mid), 0), + p) + f = x -> prob.f(x', p) end else if isinplace(prob) - dx = zeros(eltype(lb), prob.batch) - f = x -> (prob.f(dx, x', p); dx) + @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)) else - f = x -> prob.f(x', p) + y = prob.f(mid, p) + f = x -> prob.f(x, prob.p) end end - ncalls = prob.batch == 0 ? alg.ncalls : prob.batch + + if prob.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 val, err, chi = vegas(f, lb, ub, rtol = reltol, atol = abstol, maxiter = maxiters, nbins = alg.nbins, debug = alg.debug, - ncalls = ncalls, batch = prob.batch != 0) + ncalls = ncalls, batch = prob.f isa BatchIntegralFunction) SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success) end diff --git a/src/algorithms.jl b/src/algorithms.jl index 0adefdbd..ee694d58 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -63,6 +63,10 @@ each dimension of the integration domain is divided into, the number of function calls per iteration of the algorithm, and whether debug info should be printed, respectively. +## Limitations + +This algorithm can only integrate `Float64`-valued functions + ## References @article{lepage1978new, @@ -132,7 +136,7 @@ Struct for evaluating an integral via the trapezoidal rule. Example with sampled data: ``` -using Integrals +using Integrals f = x -> x^2 x = range(0, 1, length=20) y = f.(x) diff --git a/src/common.jl b/src/common.jl index 01e6ca6a..c1d4a03c 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1,11 +1,8 @@ -mutable struct IntegralCache{iip, F, B, P, PK, A, S, K, Tc} +mutable struct IntegralCache{iip, F, D, P, PK, A, S, K, Tc} iip::Val{iip} f::F - lb::B - ub::B - nout::Int + domain::D p::P - batch::Int prob_kwargs::PK alg::A sensealg::S @@ -18,16 +15,16 @@ SciMLBase.isinplace(::IntegralCache{iip}) where {iip} = iip init_cacheval(::SciMLBase.AbstractIntegralAlgorithm, args...) = nothing function SciMLBase.init(prob::IntegralProblem{iip}, - alg::SciMLBase.AbstractIntegralAlgorithm; - sensealg = ReCallVJP(ZygoteVJP()), - do_inf_transformation = nothing, kwargs...) where {iip} + alg::SciMLBase.AbstractIntegralAlgorithm; + 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) IntegralCache{iip, typeof(prob.f), - typeof(prob.lb), + typeof(prob.domain), typeof(prob.p), typeof(prob.kwargs), typeof(alg), @@ -35,11 +32,8 @@ function SciMLBase.init(prob::IntegralProblem{iip}, typeof(kwargs), typeof(cacheval)}(Val(iip), prob.f, - prob.lb, - prob.ub, - prob.nout, + prob.domain, prob.p, - prob.batch, prob.kwargs, alg, sensealg, @@ -47,6 +41,31 @@ function SciMLBase.init(prob::IntegralProblem{iip}, cacheval) end +function Base.getproperty(cache::IntegralCache, name::Symbol) + if name === :lb + domain = getfield(cache, :domain) + lb, ub = domain + return lb + elseif name === :ub + domain = getfield(cache, :domain) + lb, ub = domain + return ub + end + return getfield(cache, name) +end +function Base.setproperty!(cache::IntegralCache, name::Symbol, x) + if name === :lb + lb, ub = cache.domain + setfield!(cache, :domain, (oftype(lb, x), ub)) + return x + elseif name === :ub + lb, ub = cache.domain + setfield!(cache, :domain, (lb, oftype(ub, x))) + return x + end + return setfield!(cache, name, x) +end + # Throw error if alg is not provided, as defaults are not implemented. function SciMLBase.solve(::IntegralProblem; kwargs...) checkkwargs(kwargs...) @@ -75,19 +94,18 @@ These common arguments are: - `reltol` (relative tolerance in changes of the objective value) """ function SciMLBase.solve(prob::IntegralProblem, - alg::SciMLBase.AbstractIntegralAlgorithm; - kwargs...) + alg::SciMLBase.AbstractIntegralAlgorithm; + kwargs...) solve!(init(prob, alg; kwargs...)) end function SciMLBase.solve!(cache::IntegralCache) - __solvebp(cache, cache.alg, cache.sensealg, cache.lb, cache.ub, cache.p; + __solvebp(cache, cache.alg, cache.sensealg, cache.domain, cache.p; cache.kwargs...) end function build_problem(cache::IntegralCache{iip}) where {iip} - IntegralProblem{iip}(cache.f, cache.lb, cache.ub, cache.p; - nout = cache.nout, batch = cache.batch, cache.prob_kwargs...) + IntegralProblem{iip}(cache.f, cache.domain, cache.p; cache.prob_kwargs...) end # fallback method for existing algorithms which use no cache @@ -114,8 +132,8 @@ function Base.setproperty!(cache::SampledIntegralCache, name::Symbol, x) end function SciMLBase.init(prob::SampledIntegralProblem, - alg::SciMLBase.AbstractIntegralAlgorithm; - kwargs...) + alg::SciMLBase.AbstractIntegralAlgorithm; + kwargs...) NamedTuple(kwargs) == NamedTuple() || throw(ArgumentError("There are no keyword arguments allowed to `solve`")) @@ -142,8 +160,8 @@ solve(prob::SampledIntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kw There are no keyword arguments used to solve `SampledIntegralProblem`s """ function SciMLBase.solve(prob::SampledIntegralProblem, - alg::SciMLBase.AbstractIntegralAlgorithm; - kwargs...) + alg::SciMLBase.AbstractIntegralAlgorithm; + kwargs...) solve!(init(prob, alg; kwargs...)) end diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl index 5913d835..0cf5638c 100644 --- a/src/infinity_handling.jl +++ b/src/infinity_handling.jl @@ -18,7 +18,7 @@ function substitute_bounds(lb, ub) end return lb_sub, ub_sub end -function substitute_f_scalar(t, p, f, lb, ub) +function substitute_f(t, p, f, lb::Number, ub::Number) if isinf(lb) && isinf(ub) return f(t / (1 - t^2), p) * (1 + t^2) / (1 - t^2)^2 elseif isinf(lb) @@ -29,7 +29,19 @@ function substitute_f_scalar(t, p, f, lb, ub) return f(t, p) end end -function substitute_f_vector(t, p, f, lb, ub) +function substitute_f_iip(dt, dy, t, p, f, lb::Number, ub::Number) + if isinf(lb) && isinf(ub) + f(dt, t / (1 - t^2), p) + dt .= dy .* ((1 + t^2) / (1 - t^2)^2) + elseif isinf(lb) + return f(ub + (t / (1 + t)), p) * 1 / ((1 + t)^2) + elseif isinf(ub) + return f(lb + (t / (1 - t)), p) * 1 / ((1 - t)^2) + else + return f(t, p) + end +end +function substitute_f(t, p, f, lb::AbstractVector, ub::AbstractVector) x = similar(t) jac_diag = similar(t) for i in eachindex(lb) @@ -49,7 +61,7 @@ function substitute_f_vector(t, p, f, lb, ub) end f(x, p) * prod(jac_diag) end -function substitute_f_vector_iip(dt, t, p, f, lb, ub) +function substitute_f_iip(dt, t, p, f, lb, ub) x = similar(t) jac_diag = similar(t) for i in eachindex(lb) @@ -72,28 +84,46 @@ function substitute_f_vector_iip(dt, t, p, f, lb, ub) end function transformation_if_inf(prob, ::Val{true}) - lb = prob.lb - ub = prob.ub + lb, ub = prob.domain f = prob.f if lb isa Number lb_sub, ub_sub = substitute_bounds(lb, ub) - f_sub = (t, p) -> substitute_f_scalar(t, p, f, lb, ub) + # f_sub = (t, p) -> substitute_f_scalar(t, p, f, lb, ub) else bounds = substitute_bounds.(lb, ub) lb_sub = first.(bounds) ub_sub = last.(bounds) - if isinplace(prob) - f_sub = (dt, t, p) -> substitute_f_vector_iip(dt, t, p, f, lb, ub) + # if isinplace(prob) + # f_sub = (dt, t, p) -> substitute_f_vector_iip(dt, t, p, f, lb, ub) + # else + # f_sub = (t, p) -> substitute_f_vector(t, p, f, lb, ub) + # end + end + f_sub = if isinplace(prob) + if f isa BatchIntegralFunction + BatchIntegralFunction{true}((dt, t, p) -> substitute_f_iip(dt, t, p, f, lb, ub), + f.integrand_prototype, + max_batch = f.max_batch) + else + IntegralFunction{true}((dt, t, p) -> substitute_f_iip(dt, t, p, f, lb, ub), + f.integrand_prototype) + end + else + if f isa BatchIntegralFunction + BatchIntegralFunction{false}((t, p) -> substitute_f(t, p, f, lb, ub), + f.integrand_prototype) else - f_sub = (t, p) -> substitute_f_vector(t, p, f, lb, ub) + IntegralFunction{false}((t, p) -> substitute_f(t, p, f, lb, ub), + f.integrand_prototype) end end - return remake(prob, f = f_sub, lb = lb_sub, ub = ub_sub) + return remake(prob, f = f_sub, domain = (lb_sub, ub_sub)) end function transformation_if_inf(prob, ::Nothing) - if (prob.lb isa Number && prob.ub isa Number && (prob.ub == Inf || prob.lb == -Inf)) || - -Inf in prob.lb || Inf in prob.ub + lb, ub = prob.domain + if (lb isa Number && ub isa Number && (ub == Inf || lb == -Inf)) || + -Inf in lb || Inf in ub return transformation_if_inf(prob, Val(true)) else return transformation_if_inf(prob, Val(false)) diff --git a/src/quadrules.jl b/src/quadrules.jl index 956d4377..3d3c4531 100644 --- a/src/quadrules.jl +++ b/src/quadrules.jl @@ -23,15 +23,16 @@ function init_cacheval(alg::QuadratureRule, ::IntegralProblem) end function Integrals.__solvebp_call(cache::IntegralCache, alg::QuadratureRule, - sensealg, lb, ub, p; - reltol = nothing, abstol = nothing, - maxiters = nothing) + sensealg, domain, p; + reltol = nothing, abstol = nothing, + maxiters = nothing) prob = build_problem(cache) if isinplace(prob) error("QuadratureRule does not support inplace integrands.") end - @assert prob.batch == 0 + @assert prob.f isa IntegralFunction + lb, ub = domain val = evalrule(cache.f, cache.p, lb, ub, cache.cacheval...) err = nothing From b2798ca760db10eb08c6e824516e24e04f2d28ce Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:25:54 -0400 Subject: [PATCH 03/27] format --- src/sampled.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/sampled.jl b/src/sampled.jl index 947ad962..118f26cc 100644 --- a/src/sampled.jl +++ b/src/sampled.jl @@ -51,13 +51,13 @@ end # can be reused for other sampled rules, which should implement find_weights(x, alg) function init_cacheval(alg::SciMLBase.AbstractIntegralAlgorithm, - prob::SampledIntegralProblem) + prob::SampledIntegralProblem) find_weights(prob.x, alg) end function __solvebp_call(cache::SampledIntegralCache, - alg::SciMLBase.AbstractIntegralAlgorithm; - kwargs...) + alg::SciMLBase.AbstractIntegralAlgorithm; + kwargs...) dim = dimension(cache.dim) err = nothing data = cache.y From 187954eba97e80673b9ad089ac5f14d39535b1b8 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:26:13 -0400 Subject: [PATCH 04/27] update FastGauss... to SciMLBasev2 --- ext/IntegralsFastGaussQuadratureExt.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/ext/IntegralsFastGaussQuadratureExt.jl b/ext/IntegralsFastGaussQuadratureExt.jl index 7dd1a55a..f9a004ff 100644 --- a/ext/IntegralsFastGaussQuadratureExt.jl +++ b/ext/IntegralsFastGaussQuadratureExt.jl @@ -30,14 +30,16 @@ function composite_gauss_legendre(f, p, lb, ub, nodes, weights, subintervals) end function Integrals.__solvebp_call(prob::IntegralProblem, alg::Integrals.GaussLegendre{C}, - sensealg, lb, ub, p; - reltol = nothing, abstol = nothing, - maxiters = nothing) where {C} + sensealg, domain, p; + reltol = nothing, abstol = nothing, + maxiters = nothing) where {C} + lb, ub = domain + mid = (lb + ub) / 2 if isinplace(prob) || lb isa AbstractArray || ub isa AbstractArray error("GaussLegendre only accepts one-dimensional quadrature problems.") end - @assert prob.batch == 0 - @assert prob.nout == 1 + @assert prob.f isa IntegralFunction + @assert !isinplace(prob) if C val = composite_gauss_legendre(prob.f, prob.p, lb, ub, alg.nodes, alg.weights, alg.subintervals) From 5389a23265cee062557f65f0b60cf7907f17ad80 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:26:31 -0400 Subject: [PATCH 05/27] update IntegralsCuba to SciMLBasev2 --- lib/IntegralsCuba/src/IntegralsCuba.jl | 143 ++++++++++++------------- 1 file changed, 69 insertions(+), 74 deletions(-) diff --git a/lib/IntegralsCuba/src/IntegralsCuba.jl b/lib/IntegralsCuba/src/IntegralsCuba.jl index 26bfd7a5..d4294b43 100644 --- a/lib/IntegralsCuba/src/IntegralsCuba.jl +++ b/lib/IntegralsCuba/src/IntegralsCuba.jl @@ -116,111 +116,99 @@ struct CubaCuhre <: AbstractCubaAlgorithm end function CubaVegas(; flags = 0, seed = 0, minevals = 0, nstart = 1000, nincrease = 500, - gridno = 0) + gridno = 0) CubaVegas(flags, seed, minevals, nstart, nincrease, gridno) end function CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000, nmin = 2, - flatness = 25.0) + flatness = 25.0) CubaSUAVE(flags, seed, minevals, nnew, nmin, flatness) end function CubaDivonne(; flags = 0, seed = 0, minevals = 0, - key1 = 47, key2 = 1, key3 = 1, maxpass = 5, border = 0.0, - maxchisq = 10.0, mindeviation = 0.25) + key1 = 47, key2 = 1, key3 = 1, maxpass = 5, border = 0.0, + maxchisq = 10.0, mindeviation = 0.25) CubaDivonne(flags, seed, minevals, key1, key2, key3, maxpass, border, maxchisq, mindeviation) end CubaCuhre(; flags = 0, minevals = 0, key = 0) = CubaCuhre(flags, minevals, key) function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm, - sensealg, - lb, ub, p; - reltol = 1e-8, abstol = 1e-8, - maxiters = alg isa CubaSUAVE ? 1000000 : typemax(Int)) + sensealg, + domain, p; + reltol = 1e-8, abstol = 1e-8, + maxiters = alg isa CubaSUAVE ? 1000000 : typemax(Int)) @assert maxiters>=1000 "maxiters for $alg should be larger than 1000" - if lb isa Number && prob.batch == 0 - _x = Float64[lb] - elseif lb isa Number - _x = zeros(eltype(lb), length(lb), prob.batch) - elseif prob.batch == 0 - _x = zeros(eltype(lb), length(lb)) - else - _x = zeros(eltype(lb), length(lb), prob.batch) - end + lb, ub = domain + mid = (lb + ub) / 2 + ndim = length(mid) + (vol = prod(map(-, ub, lb))) isa Real || + 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 == 1 in Cuba will change vectors to matrices, so we won't support it when + # batching + (nvec = prob.f.max_batch) > 1 || + throw(ArgumentError("BatchIntegralFunction must take multiple batch points")) + + if mid isa Real + _x = zeros(typeof(mid), prob.f.max_batch) + scale = x -> scale_x!(resize!(_x, length(x)), ub, lb, vec(x)) + else + _x = zeros(eltype(mid), length(mid), prob.f.max_batch) + scale = x -> scale_x!(view(_x, :, 1:size(x, 2)), ub, lb, x) + end - if prob.batch == 0 if isinplace(prob) + fsize = size(prob.f.integrand_prototype)[begin:(end - 1)] + y = similar(prob.f.integrand_prototype, fsize..., nvec) + ax = map(_ -> (:), fsize) f = function (x, dx) - prob.f(dx, scale_x!(_x, ub, lb, x), p) - dx .*= prod((y) -> y[1] - y[2], zip(ub, lb)) + dy = @view(y[ax..., begin:(begin + size(dx, 2) - 1)]) + prob.f(dy, scale(x), p) + dx .= reshape(dy, :, size(dx, 2)) .* vol end else - f = function (x, dx) - dx .= prob.f(scale_x!(_x, ub, lb, x), p) .* - prod((y) -> y[1] - y[2], zip(ub, lb)) - end + y = mid isa Number ? prob.f(typeof(mid)[], p) : + prob.f(Matrix{typeof(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 + ncomp = prod(fsize) else - if lb isa Number - if isinplace(prob) - f = function (x, dx) - #todo check scale_x! - prob.f(dx', scale_x!(view(_x, 1:length(x)), ub, lb, x), p) - dx .*= prod((y) -> y[1] - y[2], zip(ub, lb)) - end - else - if prob.f([lb ub], p) isa Vector - f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p)' .* - prod((y) -> y[1] - y[2], zip(ub, lb)) - end - else - f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p) .* - prod((y) -> y[1] - y[2], zip(ub, lb)) - end - end - end + nvec = 1 + + if mid isa Real + scale = x -> scale_x(ub, lb, only(x)) else - if isinplace(prob) - f = function (x, dx) - prob.f(dx, scale_x(ub, lb, x), p) - dx .*= prod((y) -> y[1] - y[2], zip(ub, lb)) - end - else - if prob.f([lb ub], p) isa Vector - f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p)' .* - prod((y) -> y[1] - y[2], zip(ub, lb)) - end - else - f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p) .* - prod((y) -> y[1] - y[2], zip(ub, lb)) - end - end - end + _x = zeros(eltype(mid), length(mid)) + scale = x -> scale_x!(_x, ub, lb, x) end - end - ndim = length(lb) - - nvec = prob.batch == 0 ? 1 : prob.batch + if isinplace(prob) + y = similar(prob.f.integrand_prototype) + f = (x, dx) -> dx .= vec(prob.f(y, scale(x), p)) .* vol + else + y = prob.f(mid, p) + f = (x, dx) -> dx .= Iterators.flatten(prob.f(scale(x), p)) .* vol + end + ncomp = length(y) + end if alg isa CubaVegas - out = Cuba.vegas(f, ndim, prob.nout; rtol = reltol, + out = 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, prob.nout; rtol = reltol, + out = 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, prob.nout; rtol = reltol, + out = Cuba.divonne(f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, seed = alg.seed, minevals = alg.minevals, @@ -228,19 +216,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, prob.nout; rtol = reltol, + out = Cuba.cuhre(f, ndim, ncomp; rtol = reltol, atol = abstol, nvec = nvec, maxevals = maxiters, flags = alg.flags, minevals = alg.minevals, key = alg.key) end - if isinplace(prob) || prob.batch != 0 - val = out.integral - else - if prob.nout == 1 && prob.f(lb, p) isa Number + # out.integral is a Vector{Float64}, but we want to return it to the shape of the integrand + if prob.f isa BatchIntegralFunction + if y isa AbstractVector val = out.integral[1] else + val = reshape(out.integral, fsize) + end + else + if y isa Real + val = out.integral[1] + elseif y isa AbstractVector val = out.integral + else + val = reshape(out.integral, size(y)) end end From 6763300576b172cea9c6c4f2a1bd4b6d2dffcde8 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:26:48 -0400 Subject: [PATCH 06/27] update IntegralsCubature to SciMLBasev2 --- .../src/IntegralsCubature.jl | 162 +++++++++++++++++- 1 file changed, 158 insertions(+), 4 deletions(-) diff --git a/lib/IntegralsCubature/src/IntegralsCubature.jl b/lib/IntegralsCubature/src/IntegralsCubature.jl index 939bb76a..3a3f2ce7 100644 --- a/lib/IntegralsCubature/src/IntegralsCubature.jl +++ b/lib/IntegralsCubature/src/IntegralsCubature.jl @@ -49,16 +49,169 @@ end CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) function Integrals.__solvebp_call(prob::IntegralProblem, - alg::AbstractCubatureJLAlgorithm, - sensealg, lb, ub, p; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int)) + alg::AbstractCubatureJLAlgorithm, + sensealg, domain, p; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int)) + lb, ub = domain + 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 + 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 + + @assert eltype(y)<: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 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 + end + else + f = (x, dx) -> (dx .= prob.f(x, p)) + end + if mid isa Number + if alg isa CubatureJLh + val, err = Cubature.hquadrature_v(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + 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; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + val, err = Cubature.pcubature_v(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + end + end + elseif y isa AbstractArray + fsize = size(y)[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, fsize..., size(dx, 2)) + prob.f(dy, x, p) + dx .= reshape(dy, fdim, size(dx, 2)) + end + else + f = (x, dx) -> (dx .= reshape(prob.f(x, p), fdim, size(dx, 2))) + end + if mid isa Number + if alg isa CubatureJLh + 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; + 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; + reltol = reltol, abstol = abstol, + maxevals = maxiters, error_norm = alg.error_norm) + else + 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_, fsize...) + else + error("BatchIntegralFunction integrands must be arrays for Cubature.jl") + end + else + if y isa Real + # no inplace in this case, since the integrand_prototype would be mutable + f = x -> prob.f(x, p) + if lb isa Number + if alg isa CubatureJLh + val, err = Cubature.hquadrature(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + 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; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + val, err = Cubature.pcubature(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + end + end + elseif y isa AbstractArray + fsize = size(y) + fdim = length(y) + if isinplace(prob) + dy = similar(y) + f = (x, v) -> (prob.f(dy, x, p); v .= vec(dy)) + else + f = (x, v) -> (v .= vec(prob.f(x, p))) + end + if mid isa Number + if alg isa CubatureJLh + 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; + 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; + reltol = reltol, abstol = abstol, + maxevals = maxiters, error_norm = alg.error_norm) + else + val_, err = Cubature.pcubature(fdim, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, error_norm = alg.error_norm) + end + end + val = reshape(val_, fsize) + else + error("IntegralFunctions must be scalars or arrays for Cubature.jl") + end + end + + #= nout = prob.nout if nout == 1 # the output of prob.f could be either scalar or a vector of length 1, however # the behavior of the output of the integration routine is undefined (could differ # across algorithms) # Cubature will output a real number in when called without nout/fdim + if prob.batch == 0 if isinplace(prob) dx = zeros(eltype(lb), prob.nout) @@ -181,6 +334,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, end end end + =# SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) end From 6454e1cc2cb3a9a307b080705e4a5fc05dce217b Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:28:34 -0400 Subject: [PATCH 07/27] update AD ForwardDiff to SciMLBasev2 --- ext/IntegralsForwardDiffExt.jl | 135 +++++++++++++++++++-------------- 1 file changed, 78 insertions(+), 57 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 83d4a064..1e3dde7d 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -4,88 +4,109 @@ isdefined(Base, :get_extension) ? (using ForwardDiff) : (using ..ForwardDiff) ### Forward-Mode AD Intercepts # Direct AD on solvers with QuadGK and HCubature -function Integrals.__solvebp(cache, alg::QuadGKJL, sensealg, lb, ub, - p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; - kwargs...) where {T, V, P, N} - Integrals.__solvebp_call(cache, alg, sensealg, lb, ub, p; kwargs...) +function Integrals.__solvebp(cache, alg::QuadGKJL, sensealg, domain, + p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; + kwargs...) where {T, V, P, N} + Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...) end -function Integrals.__solvebp(cache, alg::HCubatureJL, sensealg, lb, ub, - p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; - kwargs...) where {T, V, P, N} - Integrals.__solvebp_call(cache, alg, sensealg, lb, ub, p; kwargs...) +function Integrals.__solvebp(cache, alg::HCubatureJL, sensealg, domain, + p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; + kwargs...) where {T, V, P, N} + Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...) end # Manually split for the pushforward -function Integrals.__solvebp(cache, alg, sensealg, lb, ub, - p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; - kwargs...) where {T, V, P, N} - primal = Integrals.__solvebp_call(cache, alg, sensealg, lb, ub, ForwardDiff.value.(p); - kwargs...) - - nout = cache.nout * P +function Integrals.__solvebp(cache, alg, sensealg, domain, + p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; + kwargs...) where {T, V, P, N} + # we need the output type to avoid perturbation confusion while unwrapping nested duals + # We compute a vector-valued integral of the primal and dual simultaneously if isinplace(cache) - dfdp = function (out, x, p) - dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p) - if cache.batch > 0 - dx = cache.nout == 1 ? similar(dualp, size(x, ndims(x))) : - similar(dualp, cache.nout, size(x, ndims(x))) - else - dx = similar(dualp, cache.nout) - end - cache.f(dx, x, dualp) + elt = eltype(cache.f.integrand_prototype) + DT = replace_dualvaltype(eltype(p), elt) + len = duallen(p) + dual_prototype = similar(cache.f.integrand_prototype, + len, + size(cache.f.integrand_prototype)...) - ys = reinterpret(ForwardDiff.Dual{T, V, P}, dx) - idx = 0 - for y in ys - for p in ForwardDiff.partials(y) - out[idx += 1] = p - end - end + dfdp_ = function (out, x, p) + dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p) + dout = reinterpret(reshape, DT, out) + cache.f(dout, x, dualp) return out end + dfdp = if cache.f isa BatchIntegralFunction + BatchIntegralFunction{true}(dfdp_, dual_prototype) + else + IntegralFunction{true}(dfdp_, dual_prototype) + end else - dfdp = function (x, p) + lb, ub = domain + mid = (lb + ub) / 2 + y = if cache.f isa BatchIntegralFunction + mid isa Number ? cache.f(eltype(mid)[], p) : + cache.f(Matrix{eltype(mid)}(undef, length(mid), 0), p) + else + cache.f(mid, p) + end + DT = y isa AbstractArray ? eltype(y) : typeof(y) + elt = unwrap_dualvaltype(DT) + + dfdp_ = function (x, p) dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p) ys = cache.f(x, dualp) - if cache.batch > 0 - out = similar(p, V, nout, size(x, ndims(x))) - else - out = similar(p, V, nout) - end - - idx = 0 - for y in ys - for p in ForwardDiff.partials(y) - out[idx += 1] = p - end - end - - return out + ys_ = ys isa AbstractArray ? ys : [ys] + # we need to reshape in order for batching to be consistent + return reinterpret(reshape, elt, ys_) + end + dfdp = if cache.f isa BatchIntegralFunction + BatchIntegralFunction{false}(dfdp_, nothing) + else + IntegralFunction{false}(dfdp_, nothing) end end + ForwardDiff.can_dual(elt) || ForwardDiff.throw_cannot_dual(elt) rawp = copy(reinterpret(V, p)) prob = Integrals.build_problem(cache) - dp_prob = remake(prob, f = dfdp, nout = nout, p = rawp) + 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, lb, ub, rawp; kwargs...) + dual = Integrals.__solvebp_call(dp_cache, alg, sensealg, domain, rawp; kwargs...) - res = similar(p, cache.nout) - partials = reinterpret(typeof(first(res).partials), dual.u) - for idx in eachindex(res) - res[idx] = ForwardDiff.Dual{T, V, P}(primal.u[idx], partials[idx]) - end - if primal.u isa Number - res = first(res) - end - SciMLBase.build_solution(prob, alg, res, primal.resid) + res = reinterpret(reshape, DT, dual.u) + # TODO: if y is a Number/Dual, then return a Number/dual, not an array + # if y isa AbstractArray + # reinterpret(reshape, ForwardDiff.Dual{T, DT, P}, dual.u) + # else + # ForwardDiff.Dual + # end + SciMLBase.build_solution(prob, alg, res, dual.resid) +end + +duallen(::Type{T}) where {T} = 1 +duallen(::T) where {T} = duallen(T) +duallen(::AbstractArray{T}) where {T} = duallen(T) +function duallen(::Type{ForwardDiff.Dual{T, V, P}}) where {T, V, P} + len = duallen(V) + return len * (P + 1) +end + +replace_dualvaltype(::Type{T}, ::Type{S}) where {T, S} = S +function replace_dualvaltype(::Type{ForwardDiff.Dual{T, V, P}}, + ::Type{S}) where {T, V, P, S} + return ForwardDiff.Dual{T, replace_dualvaltype(V, S), P} +end + +unwrap_dualvaltype(::Type{T}) where {T} = T +function unwrap_dualvaltype(::Type{ForwardDiff.Dual{T, V, P}}) where {T, V, P} + unwrap_dualvaltype(V) end end From edfa55442009d82d7e48047f6185ce0a69ac7995 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:28:46 -0400 Subject: [PATCH 08/27] update AD Zygote to SciMLBasev2 --- ext/IntegralsZygoteExt.jl | 101 +++++++++++++++++++++++++++++--------- 1 file changed, 79 insertions(+), 22 deletions(-) diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index 986c522d..bfecfab1 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -3,29 +3,81 @@ using Integrals if isdefined(Base, :get_extension) using Zygote import ChainRulesCore - import ChainRulesCore: NoTangent, ProjectTo + import ChainRulesCore: Tangent, NoTangent, ProjectTo else using ..Zygote import ..Zygote.ChainRulesCore - import ..Zygote.ChainRulesCore: NoTangent, ProjectTo + import ..Zygote.ChainRulesCore: Tangent, NoTangent, ProjectTo end ChainRulesCore.@non_differentiable Integrals.checkkwargs(kwargs...) -ChainRulesCore.@non_differentiable Integrals.isinplace(f, n) # fixes #99 +ChainRulesCore.@non_differentiable Integrals.isinplace(f, args...) # fixes #99 -function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, sensealg, lb, ub, - p; - kwargs...) - out = Integrals.__solvebp_call(cache, alg, sensealg, lb, ub, p; kwargs...) +function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, sensealg, domain, + p; + kwargs...) + # TODO: integrate the primal and dual in the same call to the quadrature library + out = Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...) # the adjoint will be the integral of the input sensitivities, so it maps the # sensitivity of the output to an object of the type of the parameters function quadrature_adjoint(Δ) # https://juliadiff.org/ChainRulesCore.jl/dev/design/many_tangents.html#manytypes - y = cache.nout == 1 ? Δ[1] : Δ # interpret the output as scalar - # this will not be type-stable, but I believe it is unavoidable due to two ambiguities: - # 1. Δ is the output of the algorithm, and when nout = 1 it is undefined whether the - # output of the algorithm must be a scalar or a vector of length 1 - # 2. when nout = 1 the integrand can either be a scalar or a vector of length 1 + if isinplace(cache) + # zygote doesn't support mutation, so we build an oop pullback + dx = similar(cache.f.integrand_prototype) + _f = x -> cache.f(dx, x, p) + if sensealg.vjp isa Integrals.ZygoteVJP + if cache.f isa BatchIntegralFunction + # TODO: let the user pass a batched jacobian so we can return a BatchIntegralFunction + dfdp_ = function (x, p) + x_ = x isa AbstractArray ? reshape(x, size(x)..., 1) : [x] + z, back = Zygote.pullback(p) do p + _dx = Zygote.Buffer(dx, size(dx)[begin:(end - 1)]..., 1) + cache.f(_dx, x_, p) + copy(_dx) + end + return back(z .= (Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) : + [Δ]))[1] + end + dfdp = IntegralFunction{false}(dfdp_, nothing) + else + dfdp_ = function (x, p) + _, back = Zygote.pullback(p) do p + _dx = Zygote.Buffer(dx) + cache.f(_dx, x, p) + copy(_dx) + end + back(Δ)[1] + end + dfdp = IntegralFunction{false}(dfdp_, nothing) + end + elseif sensealg.vjp isa Integrals.ReverseDiffVJP + error("TODO") + end + else + _f = x -> cache.f(x, p) + if sensealg.vjp isa Integrals.ZygoteVJP + if cache.f isa BatchIntegralFunction + # TODO: let the user pass a batched jacobian so we can return a BatchIntegralFunction + dfdp_ = function (x, p) + x_ = x isa AbstractArray ? reshape(x, size(x)..., 1) : [x] + z, back = Zygote.pullback(p -> cache.f(x_, p), p) + return back(z .= (Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) : + [Δ]))[1] + end + dfdp = IntegralFunction{false}(dfdp_, nothing) + else + dfdp_ = function (x, p) + z, back = Zygote.pullback(p -> cache.f(x, p), p) + back(z isa Number ? only(Δ) : Δ)[1] + end + dfdp = IntegralFunction{false}(dfdp_, nothing) + end + elseif sensealg.vjp isa Integrals.ReverseDiffVJP + error("TODO") + end + end + #= if isinplace(cache) dx = zeros(cache.nout) _f = x -> cache.f(dx, x, p) @@ -83,9 +135,10 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal error("TODO") end end - + =# prob = Integrals.build_problem(cache) - dp_prob = remake(prob, f = dfdp, nout = length(p)) + # dp_prob = remake(prob, f = dfdp) # fails because we change iip + dp_prob = IntegralProblem(dfdp, prob.domain, prob.p; prob.kwargs...) # the infinity transformation was already applied to f so we don't apply it to dfdp dp_cache = init(dp_prob, alg; @@ -97,15 +150,20 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal dp = project_p(Integrals.__solvebp_call(dp_cache, alg, sensealg, - lb, - ub, + domain, p; kwargs...).u) + lb, ub = domain if lb isa Number - dlb = cache.batch > 0 ? -_f([lb]) : -_f(lb) - dub = cache.batch > 0 ? _f([ub]) : _f(ub) - return (NoTangent(), NoTangent(), NoTangent(), NoTangent(), dlb, dub, dp) + dlb = cache.f isa BatchIntegralFunction ? -_f([lb]) : -_f(lb) + dub = cache.f isa BatchIntegralFunction ? _f([ub]) : _f(ub) + return (NoTangent(), + NoTangent(), + NoTangent(), + NoTangent(), + Tangent{typeof(domain)}(dlb, dub), + dp) else # we need to compute 2*length(lb) integrals on the faces of the hypercube, as we # can see from writing the multidimensional integral as an iterated integral @@ -117,15 +175,14 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal # other kinds of domains. The only question is to determine ω in terms of f and # the deformation of the surface (e.g. consider integral over an ellipse and # asking for the derivative of the result w.r.t. the semiaxes of the ellipse) - return (NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), - NoTangent(), dp) + return (NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), dp) end end out, quadrature_adjoint end Zygote.@adjoint function Zygote.literal_getproperty(sol::SciMLBase.IntegralSolution, - ::Val{:u}) + ::Val{:u}) sol.u, Δ -> (SciMLBase.build_solution(sol.prob, sol.alg, Δ, sol.resid),) end end From 4f60f0619f6fe738f5d223ae30e3ae83886bc8ca Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 30 Oct 2023 17:28:58 -0400 Subject: [PATCH 09/27] update test for SciMLBasev2 --- test/derivative_tests.jl | 2 +- test/inf_integral_tests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 400e1bbe..170dc52a 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -198,7 +198,7 @@ dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p) ## iip Batch multi dim function g(dx, x, p) - dx .= dropdims(sum(x * p[1] .+ p[2] * p[3], dims = 1), dims = 1) + dx .= sum(x * p[1] .+ p[2] * p[3], dims = 1) end lb = [1.0, 1.0] diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl index b51ea690..7c99c9f7 100644 --- a/test/inf_integral_tests.jl +++ b/test/inf_integral_tests.jl @@ -58,4 +58,4 @@ prob = IntegralProblem(m2, SVector(-Inf, -Inf), SVector(Inf, Inf)) prob = @test_nowarn @inferred Integrals.transformation_if_inf(prob, Val(true)) @test_nowarn @inferred Integrals.__solvebp_call(prob, HCubatureJL(), Integrals.ReCallVJP(Integrals.ZygoteVJP()), - prob.lb, prob.ub, prob.p) + prob.domain, prob.p) From 64f80c34b1a2e00cce222a33848d0254910c580d Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 1 Nov 2023 15:45:02 -0400 Subject: [PATCH 10/27] bump SciML version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 68db8e53..13a153f5 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,7 @@ MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3, 0.1" QuadGK = "2.5" Reexport = "0.2, 1.0" Requires = "1" -SciMLBase = "2.5" +SciMLBase = "2.6" Zygote = "0.4.22, 0.5, 0.6" julia = "1.6" From cf803f33f5e4023589a0e40bc7d9ed58f818b37f Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 1 Nov 2023 15:52:59 -0400 Subject: [PATCH 11/27] unwrap output array from Forward pass when ok --- ext/IntegralsForwardDiffExt.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 1e3dde7d..11d8e577 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -82,13 +82,13 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, dual = Integrals.__solvebp_call(dp_cache, alg, sensealg, domain, rawp; kwargs...) res = reinterpret(reshape, DT, dual.u) - # TODO: if y is a Number/Dual, then return a Number/dual, not an array - # if y isa AbstractArray - # reinterpret(reshape, ForwardDiff.Dual{T, DT, P}, dual.u) - # else - # ForwardDiff.Dual - # end - SciMLBase.build_solution(prob, alg, res, dual.resid) + out = if (cache.f isa BatchIntegralFunction && cache.f.integrand_prototype isa AbstractVector) || + (cache.f isa IntegralFunction && !(y isa AbstractArray)) + only(res) + else + res + end + SciMLBase.build_solution(prob, alg, out, dual.resid) end duallen(::Type{T}) where {T} = 1 From 02710ca65cfb7dfd59fe0e3c52888c54ef693051 Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 1 Nov 2023 15:53:14 -0400 Subject: [PATCH 12/27] delete old comments --- ext/IntegralsZygoteExt.jl | 58 -------- .../src/IntegralsCubature.jl | 131 ------------------ src/infinity_handling.jl | 6 - 3 files changed, 195 deletions(-) diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index bfecfab1..d5dcc034 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -77,65 +77,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal error("TODO") end end - #= - if isinplace(cache) - dx = zeros(cache.nout) - _f = x -> cache.f(dx, x, p) - if sensealg.vjp isa Integrals.ZygoteVJP - dfdp = function (dx, x, p) - z, back = Zygote.pullback(p) do p - _dx = cache.nout == 1 ? - Zygote.Buffer(dx, eltype(y), size(x, ndims(x))) : - Zygote.Buffer(dx, eltype(y), cache.nout, size(x, ndims(x))) - cache.f(_dx, x, p) - copy(_dx) - end - z .= zero(eltype(z)) - for idx in 1:size(x, ndims(x)) - z isa Vector ? (z[idx] = y) : (z[:, idx] .= y) - dx[:, idx] .= back(z)[1] - z isa Vector ? (z[idx] = zero(eltype(z))) : - (z[:, idx] .= zero(eltype(z))) - end - end - elseif sensealg.vjp isa Integrals.ReverseDiffVJP - error("TODO") - end - else - _f = x -> cache.f(x, p) - if sensealg.vjp isa Integrals.ZygoteVJP - if cache.batch > 0 - dfdp = function (x, p) - z, back = Zygote.pullback(p -> cache.f(x, p), p) - # messy, there are 4 cases, some better in forward mode than reverse - # 1: length(y) == 1 and length(p) == 1 - # 2: length(y) > 1 and length(p) == 1 - # 3: length(y) == 1 and length(p) > 1 - # 4: length(y) > 1 and length(p) > 1 - z .= zero(eltype(z)) - out = zeros(eltype(p), size(p)..., size(x, ndims(x))) - for idx in 1:size(x, ndims(x)) - z isa Vector ? (z[idx] = y) : (z[:, idx] .= y) - out isa Vector ? (out[idx] = back(z)[1]) : - (out[:, idx] .= back(z)[1]) - z isa Vector ? (z[idx] = zero(y)) : - (z[:, idx] .= zero(eltype(y))) - end - out - end - else - dfdp = function (x, p) - _, back = Zygote.pullback(p -> cache.f(x, p), p) - back(y)[1] - end - end - - elseif sensealg.vjp isa Integrals.ReverseDiffVJP - error("TODO") - end - end - =# prob = Integrals.build_problem(cache) # dp_prob = remake(prob, f = dfdp) # fails because we change iip dp_prob = IntegralProblem(dfdp, prob.domain, prob.p; prob.kwargs...) diff --git a/lib/IntegralsCubature/src/IntegralsCubature.jl b/lib/IntegralsCubature/src/IntegralsCubature.jl index 3a3f2ce7..97ca8766 100644 --- a/lib/IntegralsCubature/src/IntegralsCubature.jl +++ b/lib/IntegralsCubature/src/IntegralsCubature.jl @@ -204,137 +204,6 @@ function Integrals.__solvebp_call(prob::IntegralProblem, end end - #= - nout = prob.nout - if nout == 1 - # the output of prob.f could be either scalar or a vector of length 1, however - # the behavior of the output of the integration routine is undefined (could differ - # across algorithms) - # Cubature will output a real number in when called without nout/fdim - - if prob.batch == 0 - if isinplace(prob) - dx = zeros(eltype(lb), prob.nout) - f = (x) -> (prob.f(dx, x, p); dx[1]) - else - f = (x) -> prob.f(x, p)[1] - end - if lb isa Number - if alg isa CubatureJLh - val, err = Cubature.hquadrature(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - 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; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - val, err = Cubature.pcubature(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - end - end - else - if isinplace(prob) - f = (x, dx) -> prob.f(dx, x, p) - else - f = (x, dx) -> (dx .= prob.f(x, p)) - end - if lb isa Number - if alg isa CubatureJLh - val, err = Cubature.hquadrature_v(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - 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; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - val, err = Cubature.pcubature_v(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - end - end - end - else - if prob.batch == 0 - if isinplace(prob) - f = (x, dx) -> (prob.f(dx, x, p); dx) - else - f = (x, dx) -> (dx .= prob.f(x, p)) - end - if lb isa Number - if alg isa CubatureJLh - val, err = Cubature.hquadrature(nout, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) - else - val, err = Cubature.pquadrature(nout, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) - end - else - if alg isa CubatureJLh - val, err = Cubature.hcubature(nout, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) - else - val, err = Cubature.pcubature(nout, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) - end - end - else - if isinplace(prob) - f = (x, dx) -> (prob.f(dx, x, p); dx) - else - f = (x, dx) -> (dx .= prob.f(x, p)) - end - - if lb isa Number - if alg isa CubatureJLh - val, err = Cubature.hquadrature_v(nout, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) - else - val, err = Cubature.pquadrature_v(nout, 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(nout, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) - else - val, err = Cubature.pcubature_v(nout, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) - end - end - end - end - =# SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) end diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl index 0cf5638c..f2e6e8e0 100644 --- a/src/infinity_handling.jl +++ b/src/infinity_handling.jl @@ -88,16 +88,10 @@ function transformation_if_inf(prob, ::Val{true}) f = prob.f if lb isa Number lb_sub, ub_sub = substitute_bounds(lb, ub) - # f_sub = (t, p) -> substitute_f_scalar(t, p, f, lb, ub) else bounds = substitute_bounds.(lb, ub) lb_sub = first.(bounds) ub_sub = last.(bounds) - # if isinplace(prob) - # f_sub = (dt, t, p) -> substitute_f_vector_iip(dt, t, p, f, lb, ub) - # else - # f_sub = (t, p) -> substitute_f_vector(t, p, f, lb, ub) - # end end f_sub = if isinplace(prob) if f isa BatchIntegralFunction From 77f0e83c1521ab4601129387b4b29c552d142b7e Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 1 Nov 2023 17:39:47 -0400 Subject: [PATCH 13/27] fix dual unwrapping in forwarddiff --- ext/IntegralsForwardDiffExt.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 11d8e577..11e84d09 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -24,6 +24,7 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, # we need the output type to avoid perturbation confusion while unwrapping nested duals # We compute a vector-valued integral of the primal and dual simultaneously if isinplace(cache) + y = cache.f.integrand_prototype elt = eltype(cache.f.integrand_prototype) DT = replace_dualvaltype(eltype(p), elt) len = duallen(p) @@ -82,8 +83,9 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, dual = Integrals.__solvebp_call(dp_cache, alg, sensealg, domain, rawp; kwargs...) res = reinterpret(reshape, DT, dual.u) - out = if (cache.f isa BatchIntegralFunction && cache.f.integrand_prototype isa AbstractVector) || - (cache.f isa IntegralFunction && !(y isa AbstractArray)) + # unwrap the dual when the primal would return a scalar + out = if (cache.f isa BatchIntegralFunction && y isa AbstractVector) || + !(y isa AbstractArray) only(res) else res From 24fce57bd6f7095ed359b22fb55a0141c23693a8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 02:22:51 -0400 Subject: [PATCH 14/27] Change to using extension packages --- Project.toml | 10 +- .../IntegralsCubaExt.jl | 134 +---------- lib/IntegralsCuba/LICENSE | 19 -- lib/IntegralsCuba/Project.toml | 19 -- lib/IntegralsCuba/test/runtests.jl | 1 - lib/IntegralsCubature/LICENSE | 19 -- lib/IntegralsCubature/Project.toml | 19 -- .../src/IntegralsCubature.jl | 212 ------------------ lib/IntegralsCubature/test/runtests.jl | 1 - src/Integrals.jl | 3 + src/algorithms.jl | 173 ++++++++++++++ test/runtests.jl | 8 - 12 files changed, 187 insertions(+), 431 deletions(-) rename lib/IntegralsCuba/src/IntegralsCuba.jl => ext/IntegralsCubaExt.jl (57%) delete mode 100644 lib/IntegralsCuba/LICENSE delete mode 100644 lib/IntegralsCuba/Project.toml delete mode 100644 lib/IntegralsCuba/test/runtests.jl delete mode 100644 lib/IntegralsCubature/LICENSE delete mode 100644 lib/IntegralsCubature/Project.toml delete mode 100644 lib/IntegralsCubature/src/IntegralsCubature.jl delete mode 100644 lib/IntegralsCubature/test/runtests.jl diff --git a/Project.toml b/Project.toml index 13a153f5..84ec6d34 100644 --- a/Project.toml +++ b/Project.toml @@ -15,11 +15,15 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" [weakdeps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1" +Cubature = "667455a9-e2ce-5579-9412-b964f529a492" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [extensions] +IntegralsCubaExt = "Cuba" +IntegralsCubatureExt = "Cubature" IntegralsFastGaussQuadratureExt = "FastGaussQuadrature" IntegralsForwardDiffExt = "ForwardDiff" IntegralsZygoteExt = ["Zygote", "ChainRulesCore"] @@ -27,6 +31,8 @@ IntegralsZygoteExt = ["Zygote", "ChainRulesCore"] [compat] ChainRulesCore = "0.10.7, 1" CommonSolve = "0.2" +Cuba = "2" +Cubature = "1" Distributions = "0.23, 0.24, 0.25" FastGaussQuadrature = "0.5" ForwardDiff = "0.10" @@ -41,6 +47,8 @@ julia = "1.6" [extras] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1" +Cubature = "667455a9-e2ce-5579-9412-b964f529a492" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" @@ -53,4 +61,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature"] +test = ["SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"] diff --git a/lib/IntegralsCuba/src/IntegralsCuba.jl b/ext/IntegralsCubaExt.jl similarity index 57% rename from lib/IntegralsCuba/src/IntegralsCuba.jl rename to ext/IntegralsCubaExt.jl index d4294b43..d485d312 100644 --- a/lib/IntegralsCuba/src/IntegralsCuba.jl +++ b/ext/IntegralsCubaExt.jl @@ -1,136 +1,8 @@ -module IntegralsCuba +module IntegralsCubaExt using Integrals, Cuba import Integrals: transformation_if_inf, scale_x, scale_x! -abstract type AbstractCubaAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end -""" - CubaVegas() - -Multidimensional adaptive Monte Carlo integration from Cuba.jl. -Importance sampling is used to reduce variance. - -## References - -@article{lepage1978new, -title={A new algorithm for adaptive multidimensional integration}, -author={Lepage, G Peter}, -journal={Journal of Computational Physics}, -volume={27}, -number={2}, -pages={192--203}, -year={1978}, -publisher={Elsevier} -} -""" -struct CubaVegas <: AbstractCubaAlgorithm - flags::Int - seed::Int - minevals::Int - nstart::Int - nincrease::Int - gridno::Int -end -""" - CubaSUAVE() - -Multidimensional adaptive Monte Carlo integration from Cuba.jl. -Suave stands for subregion-adaptive VEGAS. -Importance sampling and subdivision are thus used to reduce variance. - -## References - -@article{hahn2005cuba, -title={Cuba—a library for multidimensional numerical integration}, -author={Hahn, Thomas}, -journal={Computer Physics Communications}, -volume={168}, -number={2}, -pages={78--95}, -year={2005}, -publisher={Elsevier} -} -""" -struct CubaSUAVE{R} <: AbstractCubaAlgorithm where {R <: Real} - flags::Int - seed::Int - minevals::Int - nnew::Int - nmin::Int - flatness::R -end -""" - CubaDivonne() - -Multidimensional adaptive Monte Carlo integration from Cuba.jl. -Stratified sampling is used to reduce variance. - -## References - -@article{friedman1981nested, -title={A nested partitioning procedure for numerical multiple integration}, -author={Friedman, Jerome H and Wright, Margaret H}, -journal={ACM Transactions on Mathematical Software (TOMS)}, -volume={7}, -number={1}, -pages={76--92}, -year={1981}, -publisher={ACM New York, NY, USA} -} -""" -struct CubaDivonne{R1, R2, R3} <: - AbstractCubaAlgorithm where {R1 <: Real, R2 <: Real, R3 <: Real} - flags::Int - seed::Int - minevals::Int - key1::Int - key2::Int - key3::Int - maxpass::Int - border::R1 - maxchisq::R2 - mindeviation::R3 -end -""" - CubaCuhre() - -Multidimensional h-adaptive integration from Cuba.jl. - -## References - -@article{berntsen1991adaptive, -title={An adaptive algorithm for the approximate calculation of multiple integrals}, -author={Berntsen, Jarle and Espelid, Terje O and Genz, Alan}, -journal={ACM Transactions on Mathematical Software (TOMS)}, -volume={17}, -number={4}, -pages={437--451}, -year={1991}, -publisher={ACM New York, NY, USA} -} -""" -struct CubaCuhre <: AbstractCubaAlgorithm - flags::Int - minevals::Int - key::Int -end - -function CubaVegas(; flags = 0, seed = 0, minevals = 0, nstart = 1000, nincrease = 500, - gridno = 0) - CubaVegas(flags, seed, minevals, nstart, nincrease, gridno) -end -function CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000, nmin = 2, - flatness = 25.0) - CubaSUAVE(flags, seed, minevals, nnew, nmin, flatness) -end -function CubaDivonne(; flags = 0, seed = 0, minevals = 0, - key1 = 47, key2 = 1, key3 = 1, maxpass = 5, border = 0.0, - maxchisq = 10.0, mindeviation = 0.25) - CubaDivonne(flags, seed, minevals, key1, key2, key3, maxpass, border, maxchisq, - mindeviation) -end -CubaCuhre(; flags = 0, minevals = 0, key = 0) = CubaCuhre(flags, minevals, key) - function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm, sensealg, domain, p; @@ -243,6 +115,4 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori chi = out.probability, retcode = ReturnCode.Success) end -export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre - -end +end \ No newline at end of file diff --git a/lib/IntegralsCuba/LICENSE b/lib/IntegralsCuba/LICENSE deleted file mode 100644 index aa799652..00000000 --- a/lib/IntegralsCuba/LICENSE +++ /dev/null @@ -1,19 +0,0 @@ -Copyright (c) 2019 Chris Rackauckas - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. diff --git a/lib/IntegralsCuba/Project.toml b/lib/IntegralsCuba/Project.toml deleted file mode 100644 index 2acf64ca..00000000 --- a/lib/IntegralsCuba/Project.toml +++ /dev/null @@ -1,19 +0,0 @@ -name = "IntegralsCuba" -uuid = "e00cd5f1-6337-4131-8b37-28b2fe4cd6cb" -authors = ["Chris Rackauckas "] -version = "0.3.0" - -[deps] -Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1" -Integrals = "de52edbc-65ea-441a-8357-d3a637375a31" - -[compat] -Cuba = "2" -Integrals = "3" -julia = "1.6" - -[extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["Test"] diff --git a/lib/IntegralsCuba/test/runtests.jl b/lib/IntegralsCuba/test/runtests.jl deleted file mode 100644 index 8b137891..00000000 --- a/lib/IntegralsCuba/test/runtests.jl +++ /dev/null @@ -1 +0,0 @@ - diff --git a/lib/IntegralsCubature/LICENSE b/lib/IntegralsCubature/LICENSE deleted file mode 100644 index aa799652..00000000 --- a/lib/IntegralsCubature/LICENSE +++ /dev/null @@ -1,19 +0,0 @@ -Copyright (c) 2019 Chris Rackauckas - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. diff --git a/lib/IntegralsCubature/Project.toml b/lib/IntegralsCubature/Project.toml deleted file mode 100644 index a4f3b7b5..00000000 --- a/lib/IntegralsCubature/Project.toml +++ /dev/null @@ -1,19 +0,0 @@ -name = "IntegralsCubature" -uuid = "c31f79ba-6e32-46d4-a52f-182a8ac42a54" -authors = ["Chris Rackauckas "] -version = "0.2.3" - -[deps] -Cubature = "667455a9-e2ce-5579-9412-b964f529a492" -Integrals = "de52edbc-65ea-441a-8357-d3a637375a31" - -[compat] -Cubature = "1" -Integrals = "3" -julia = "1.6" - -[extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["Test"] diff --git a/lib/IntegralsCubature/src/IntegralsCubature.jl b/lib/IntegralsCubature/src/IntegralsCubature.jl deleted file mode 100644 index 97ca8766..00000000 --- a/lib/IntegralsCubature/src/IntegralsCubature.jl +++ /dev/null @@ -1,212 +0,0 @@ -module IntegralsCubature - -using Integrals, Cubature - -import Integrals: transformation_if_inf, scale_x, scale_x! -import Cubature: INDIVIDUAL, PAIRED, L1, L2, LINF - -abstract type AbstractCubatureJLAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end -""" - CubatureJLh() - -Multidimensional h-adaptive integration from Cubature.jl. -`error_norm` specifies the convergence criterion for vector valued integrands. -Defaults to `IntegralsCubature.INDIVIDUAL`, other options are -`IntegralsCubature.PAIRED`, `IntegralsCubature.L1`, `IntegralsCubature.L2`, or `IntegralsCubature.LINF`. - -## References - -@article{genz1980remarks, -title={Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region}, -author={Genz, Alan C and Malik, Aftab Ahmad}, -journal={Journal of Computational and Applied mathematics}, -volume={6}, -number={4}, -pages={295--302}, -year={1980}, -publisher={Elsevier} -} -""" -struct CubatureJLh <: AbstractCubatureJLAlgorithm - error_norm::Int32 -end -CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) - -""" - CubatureJLp() - -Multidimensional p-adaptive integration from Cubature.jl. -This method is based on repeatedly doubling the degree of the cubature rules, -until convergence is achieved. -The used cubature rule is a tensor product of Clenshaw–Curtis quadrature rules. -`error_norm` specifies the convergence criterion for vector valued integrands. -Defaults to `IntegralsCubature.INDIVIDUAL`, other options are -`IntegralsCubature.PAIRED`, `IntegralsCubature.L1`, `IntegralsCubature.L2`, or `IntegralsCubature.LINF`. -""" -struct CubatureJLp <: AbstractCubatureJLAlgorithm - error_norm::Int32 -end -CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) - -function Integrals.__solvebp_call(prob::IntegralProblem, - alg::AbstractCubatureJLAlgorithm, - sensealg, domain, p; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int)) - lb, ub = domain - 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 - 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 - - @assert eltype(y)<: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 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 - end - else - f = (x, dx) -> (dx .= prob.f(x, p)) - end - if mid isa Number - if alg isa CubatureJLh - val, err = Cubature.hquadrature_v(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - 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; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - val, err = Cubature.pcubature_v(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - end - end - elseif y isa AbstractArray - fsize = size(y)[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, fsize..., size(dx, 2)) - prob.f(dy, x, p) - dx .= reshape(dy, fdim, size(dx, 2)) - end - else - f = (x, dx) -> (dx .= reshape(prob.f(x, p), fdim, size(dx, 2))) - end - if mid isa Number - if alg isa CubatureJLh - 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; - 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; - reltol = reltol, abstol = abstol, - maxevals = maxiters, error_norm = alg.error_norm) - else - 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_, fsize...) - else - error("BatchIntegralFunction integrands must be arrays for Cubature.jl") - end - else - if y isa Real - # no inplace in this case, since the integrand_prototype would be mutable - f = x -> prob.f(x, p) - if lb isa Number - if alg isa CubatureJLh - val, err = Cubature.hquadrature(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - 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; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - else - val, err = Cubature.pcubature(f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters) - end - end - elseif y isa AbstractArray - fsize = size(y) - fdim = length(y) - if isinplace(prob) - dy = similar(y) - f = (x, v) -> (prob.f(dy, x, p); v .= vec(dy)) - else - f = (x, v) -> (v .= vec(prob.f(x, p))) - end - if mid isa Number - if alg isa CubatureJLh - 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; - 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; - reltol = reltol, abstol = abstol, - maxevals = maxiters, error_norm = alg.error_norm) - else - val_, err = Cubature.pcubature(fdim, f, lb, ub; - reltol = reltol, abstol = abstol, - maxevals = maxiters, error_norm = alg.error_norm) - end - end - val = reshape(val_, fsize) - else - error("IntegralFunctions must be scalars or arrays for Cubature.jl") - end - end - - SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) -end - -export CubatureJLh, CubatureJLp - -end diff --git a/lib/IntegralsCubature/test/runtests.jl b/lib/IntegralsCubature/test/runtests.jl deleted file mode 100644 index 8b137891..00000000 --- a/lib/IntegralsCubature/test/runtests.jl +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/Integrals.jl b/src/Integrals.jl index 4d8a28da..f442971d 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -168,4 +168,7 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p; end export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, QuadratureRule, TrapezoidalRule +export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre +export CubatureJLh, CubatureJLp + end # module diff --git a/src/algorithms.jl b/src/algorithms.jl index ee694d58..550f1dbd 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -168,3 +168,176 @@ struct QuadratureRule{Q} <: SciMLBase.AbstractIntegralAlgorithm end end QuadratureRule(q; n = 250) = QuadratureRule(q, n) + +## Extension Algorithms + +abstract type AbstractCubaAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end +""" + CubaVegas() + +Multidimensional adaptive Monte Carlo integration from Cuba.jl. +Importance sampling is used to reduce variance. + +## References + +@article{lepage1978new, +title={A new algorithm for adaptive multidimensional integration}, +author={Lepage, G Peter}, +journal={Journal of Computational Physics}, +volume={27}, +number={2}, +pages={192--203}, +year={1978}, +publisher={Elsevier} +} +""" +struct CubaVegas <: AbstractCubaAlgorithm + flags::Int + seed::Int + minevals::Int + nstart::Int + nincrease::Int + gridno::Int +end +""" + CubaSUAVE() + +Multidimensional adaptive Monte Carlo integration from Cuba.jl. +Suave stands for subregion-adaptive VEGAS. +Importance sampling and subdivision are thus used to reduce variance. + +## References + +@article{hahn2005cuba, +title={Cuba—a library for multidimensional numerical integration}, +author={Hahn, Thomas}, +journal={Computer Physics Communications}, +volume={168}, +number={2}, +pages={78--95}, +year={2005}, +publisher={Elsevier} +} +""" +struct CubaSUAVE{R} <: AbstractCubaAlgorithm where {R <: Real} + flags::Int + seed::Int + minevals::Int + nnew::Int + nmin::Int + flatness::R +end +""" + CubaDivonne() + +Multidimensional adaptive Monte Carlo integration from Cuba.jl. +Stratified sampling is used to reduce variance. + +## References + +@article{friedman1981nested, +title={A nested partitioning procedure for numerical multiple integration}, +author={Friedman, Jerome H and Wright, Margaret H}, +journal={ACM Transactions on Mathematical Software (TOMS)}, +volume={7}, +number={1}, +pages={76--92}, +year={1981}, +publisher={ACM New York, NY, USA} +} +""" +struct CubaDivonne{R1, R2, R3} <: + AbstractCubaAlgorithm where {R1 <: Real, R2 <: Real, R3 <: Real} + flags::Int + seed::Int + minevals::Int + key1::Int + key2::Int + key3::Int + maxpass::Int + border::R1 + maxchisq::R2 + mindeviation::R3 +end +""" + CubaCuhre() + +Multidimensional h-adaptive integration from Cuba.jl. + +## References + +@article{berntsen1991adaptive, +title={An adaptive algorithm for the approximate calculation of multiple integrals}, +author={Berntsen, Jarle and Espelid, Terje O and Genz, Alan}, +journal={ACM Transactions on Mathematical Software (TOMS)}, +volume={17}, +number={4}, +pages={437--451}, +year={1991}, +publisher={ACM New York, NY, USA} +} +""" +struct CubaCuhre <: AbstractCubaAlgorithm + flags::Int + minevals::Int + key::Int +end + +function CubaVegas(; flags = 0, seed = 0, minevals = 0, nstart = 1000, nincrease = 500, + gridno = 0) + CubaVegas(flags, seed, minevals, nstart, nincrease, gridno) +end +function CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000, nmin = 2, + flatness = 25.0) + CubaSUAVE(flags, seed, minevals, nnew, nmin, flatness) +end +function CubaDivonne(; flags = 0, seed = 0, minevals = 0, + key1 = 47, key2 = 1, key3 = 1, maxpass = 5, border = 0.0, + maxchisq = 10.0, mindeviation = 0.25) + CubaDivonne(flags, seed, minevals, key1, key2, key3, maxpass, border, maxchisq, + mindeviation) +end +CubaCuhre(; flags = 0, minevals = 0, key = 0) = CubaCuhre(flags, minevals, key) + +abstract type AbstractCubatureJLAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end +""" + CubatureJLh() + +Multidimensional h-adaptive integration from Cubature.jl. +`error_norm` specifies the convergence criterion for vector valued integrands. +Defaults to `Cubature.INDIVIDUAL`, other options are +`Cubature.PAIRED`, `Cubature.L1`, `Cubature.L2`, or `Cubature.LINF`. + +## References + +@article{genz1980remarks, +title={Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region}, +author={Genz, Alan C and Malik, Aftab Ahmad}, +journal={Journal of Computational and Applied mathematics}, +volume={6}, +number={4}, +pages={295--302}, +year={1980}, +publisher={Elsevier} +} +""" +struct CubatureJLh <: AbstractCubatureJLAlgorithm + error_norm::Int32 +end +CubatureJLh() = CubatureJLh(Cubature.INDIVIDUAL) + +""" + CubatureJLp() + +Multidimensional p-adaptive integration from Cubature.jl. +This method is based on repeatedly doubling the degree of the cubature rules, +until convergence is achieved. +The used cubature rule is a tensor product of Clenshaw–Curtis quadrature rules. +`error_norm` specifies the convergence criterion for vector valued integrands. +Defaults to `Cubature.INDIVIDUAL`, other options are +`Cubature.PAIRED`, `Cubature.L1`, `Cubature.L2`, or `Cubature.LINF`. +""" +struct CubatureJLp <: AbstractCubatureJLAlgorithm + error_norm::Int32 +end +CubatureJLp() = CubatureJLp(Cubature.INDIVIDUAL) diff --git a/test/runtests.jl b/test/runtests.jl index c61ae8b0..00a11307 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,14 +2,6 @@ using Pkg using SafeTestsets using Test -function dev_subpkg(subpkg) - subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) - Pkg.develop(PackageSpec(path = subpkg_path)) -end - -dev_subpkg("IntegralsCuba") -dev_subpkg("IntegralsCubature") - @time @safetestset "Interface Tests" begin include("interface_tests.jl") end From 88b4ea7be336defd35d67d3bda9c2c4b74bd862f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 02:30:40 -0400 Subject: [PATCH 15/27] Incorporate https://github.com/SciML/Integrals.jl/pull/185 --- ext/IntegralsCubatureExt.jl | 144 ++++++++++++++++++++++++++++++++++++ src/algorithms.jl | 4 +- 2 files changed, 146 insertions(+), 2 deletions(-) create mode 100644 ext/IntegralsCubatureExt.jl diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl new file mode 100644 index 00000000..03401731 --- /dev/null +++ b/ext/IntegralsCubatureExt.jl @@ -0,0 +1,144 @@ +module IntegralsCubatureExt + +using Integrals, Cubature + +import Integrals: transformation_if_inf, scale_x, scale_x! +import Cubature: INDIVIDUAL, PAIRED, L1, L2, LINF + +function Integrals.__solvebp_call(prob::IntegralProblem, + alg::AbstractCubatureJLAlgorithm, + sensealg, lb, ub, p; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int)) + nout = prob.nout + if nout == 1 + # the output of prob.f could be either scalar or a vector of length 1, however + # the behavior of the output of the integration routine is undefined (could differ + # across algorithms) + # Cubature will output a real number in when called without nout/fdim + if prob.batch == 0 + if isinplace(prob) + dx = zeros(eltype(lb), prob.nout) + f = (x) -> (prob.f(dx, x, p); dx[1]) + else + f = (x) -> prob.f(x, p)[1] + end + if lb isa Number + if alg isa CubatureJLh + val, err = Cubature.hquadrature(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + 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; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + val, err = Cubature.pcubature(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + end + end + else + if isinplace(prob) + f = (x, dx) -> prob.f(dx, x, p) + else + f = (x, dx) -> (dx .= prob.f(x, p)) + end + if lb isa Number + if alg isa CubatureJLh + val, err = Cubature.hquadrature_v(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + 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; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + else + val, err = Cubature.pcubature_v(f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters) + end + end + end + else + if prob.batch == 0 + if isinplace(prob) + f = (x, dx) -> (prob.f(dx, x, p); dx) + else + f = (x, dx) -> (dx .= prob.f(x, p)) + end + if lb isa Number + if alg isa CubatureJLh + val, err = Cubature.hquadrature(nout, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, + error_norm = alg.error_norm) + else + val, err = Cubature.pquadrature(nout, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, + error_norm = alg.error_norm) + end + else + if alg isa CubatureJLh + val, err = Cubature.hcubature(nout, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, + error_norm = alg.error_norm) + else + val, err = Cubature.pcubature(nout, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, + error_norm = alg.error_norm) + end + end + else + if isinplace(prob) + f = (x, dx) -> (prob.f(dx, x, p); dx) + else + f = (x, dx) -> (dx .= prob.f(x, p)) + end + + if lb isa Number + if alg isa CubatureJLh + val, err = Cubature.hquadrature_v(nout, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, + error_norm = alg.error_norm) + else + val, err = Cubature.pquadrature_v(nout, 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(nout, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, + error_norm = alg.error_norm) + else + val, err = Cubature.pcubature_v(nout, f, lb, ub; + reltol = reltol, abstol = abstol, + maxevals = maxiters, + error_norm = alg.error_norm) + end + end + end + end + SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) +end + +end \ No newline at end of file diff --git a/src/algorithms.jl b/src/algorithms.jl index 550f1dbd..55f69445 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -324,7 +324,7 @@ publisher={Elsevier} struct CubatureJLh <: AbstractCubatureJLAlgorithm error_norm::Int32 end -CubatureJLh() = CubatureJLh(Cubature.INDIVIDUAL) +CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) """ CubatureJLp() @@ -340,4 +340,4 @@ Defaults to `Cubature.INDIVIDUAL`, other options are struct CubatureJLp <: AbstractCubatureJLAlgorithm error_norm::Int32 end -CubatureJLp() = CubatureJLp(Cubature.INDIVIDUAL) +CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) \ No newline at end of file From c9bc978820e53e44cda0cb815bc6c5481da9603b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 02:44:48 -0400 Subject: [PATCH 16/27] fix tests to use subpackages --- test/derivative_tests.jl | 2 +- test/interface_tests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 170dc52a..87b326f7 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -1,5 +1,5 @@ using Integrals, Zygote, FiniteDiff, ForwardDiff#, SciMLSensitivity -using IntegralsCuba, IntegralsCubature +using Cuba, Cubature using Test ### One Dimensional diff --git a/test/interface_tests.jl b/test/interface_tests.jl index 2d597adc..13d7b376 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -1,5 +1,5 @@ using Integrals -using IntegralsCuba, IntegralsCubature +using Cuba, Cubature using Test max_dim_test = 2 From f9720b0acbe242642b01ee2225aef1a6d123cd72 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 03:06:40 -0400 Subject: [PATCH 17/27] move alg dispatches --- ext/IntegralsCubatureExt.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl index 03401731..7bc97c8f 100644 --- a/ext/IntegralsCubatureExt.jl +++ b/ext/IntegralsCubatureExt.jl @@ -5,6 +5,9 @@ using Integrals, Cubature import Integrals: transformation_if_inf, scale_x, scale_x! import Cubature: INDIVIDUAL, PAIRED, L1, L2, LINF +CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) +CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) + function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubatureJLAlgorithm, sensealg, lb, ub, p; From 7979058166ff56febc6c88922a6bc582d386a06d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 03:18:15 -0400 Subject: [PATCH 18/27] fix docs --- docs/Project.toml | 8 ++++---- docs/src/solvers/IntegralSolvers.md | 12 ++++++------ docs/src/tutorials/differentiating_integrals.md | 2 +- docs/src/tutorials/numerical_integrals.md | 6 +++--- src/algorithms.jl | 3 +-- 5 files changed, 15 insertions(+), 16 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index b2b74b4c..8f6fc2a7 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,19 +1,19 @@ [deps] +Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1" +Cubature = "667455a9-e2ce-5579-9412-b964f529a492" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Integrals = "de52edbc-65ea-441a-8357-d3a637375a31" -IntegralsCuba = "e00cd5f1-6337-4131-8b37-28b2fe4cd6cb" -IntegralsCubature = "c31f79ba-6e32-46d4-a52f-182a8ac42a54" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] +Cuba = "2" +Cubature = "1" Distributions = "0.25" Documenter = "1" FiniteDiff = "2" ForwardDiff = "0.10" Integrals = "3" -IntegralsCuba = "0.3" -IntegralsCubature = "0.2" Zygote = "0.6" diff --git a/docs/src/solvers/IntegralSolvers.md b/docs/src/solvers/IntegralSolvers.md index f9fd3c61..10faab2e 100644 --- a/docs/src/solvers/IntegralSolvers.md +++ b/docs/src/solvers/IntegralSolvers.md @@ -5,12 +5,12 @@ The following algorithms are available: - `QuadGKJL`: Uses QuadGK.jl. Requires `nout=1` and `batch=0`, in-place is not allowed. - `HCubatureJL`: Uses HCubature.jl. Requires `batch=0`. - `VEGAS`: Uses MonteCarloIntegration.jl. Requires `nout=1`. Works only for `>1`-dimensional integrations. - - `CubatureJLh`: h-Cubature from Cubature.jl. Requires `using IntegralsCubature`. - - `CubatureJLp`: p-Cubature from Cubature.jl. Requires `using IntegralsCubature`. - - `CubaVegas`: Vegas from Cuba.jl. Requires `using IntegralsCuba`, `nout=1`. - - `CubaSUAVE`: SUAVE from Cuba.jl. Requires `using IntegralsCuba`. - - `CubaDivonne`: Divonne from Cuba.jl. Requires `using IntegralsCuba`. Works only for `>1`-dimensional integrations. - - `CubaCuhre`: Cuhre from Cuba.jl. Requires `using IntegralsCuba`. Works only for `>1`-dimensional integrations. + - `CubatureJLh`: h-Cubature from Cubature.jl. Requires `using Cubature`. + - `CubatureJLp`: p-Cubature from Cubature.jl. Requires `using Cubature`. + - `CubaVegas`: Vegas from Cuba.jl. Requires `using Cuba`, `nout=1`. + - `CubaSUAVE`: SUAVE from Cuba.jl. Requires `using Cuba`. + - `CubaDivonne`: Divonne from Cuba.jl. Requires `using Cuba`. Works only for `>1`-dimensional integrations. + - `CubaCuhre`: Cuhre from Cuba.jl. Requires `using Cuba`. Works only for `>1`-dimensional integrations. - `GaussLegendre`: Uses Gauss-Legendre quadrature with nodes and weights from FastGaussQuadrature.jl. - `QuadratureRule`: Accepts a user-defined function that returns nodes and weights. diff --git a/docs/src/tutorials/differentiating_integrals.md b/docs/src/tutorials/differentiating_integrals.md index a8b7da1c..413d7870 100644 --- a/docs/src/tutorials/differentiating_integrals.md +++ b/docs/src/tutorials/differentiating_integrals.md @@ -6,7 +6,7 @@ calls. It integrates with ForwardDiff.jl for forward-mode automatic differentiat and Zygote.jl for reverse-mode automatic differentiation. For example: ```@example AD -using Integrals, ForwardDiff, FiniteDiff, Zygote, IntegralsCuba +using Integrals, ForwardDiff, FiniteDiff, Zygote, Cuba f(x, p) = sum(sin.(x .* p)) lb = ones(2) ub = 3ones(2) diff --git a/docs/src/tutorials/numerical_integrals.md b/docs/src/tutorials/numerical_integrals.md index 6f08a087..740e79cd 100644 --- a/docs/src/tutorials/numerical_integrals.md +++ b/docs/src/tutorials/numerical_integrals.md @@ -53,7 +53,7 @@ This means that a new output vector is created every time the function `f` is ca If we do not want these allocations, we can also define `f` in-position. ```@example integrate3 -using Integrals, IntegralsCubature +using Integrals, Cubature function f(y, u, p) y[1] = sum(sin.(u)) y[2] = sum(cos.(u)) @@ -74,7 +74,7 @@ The batch interface allows us to compute multiple points at once. For example, here we do allocation-free multithreading with Cubature.jl: ```@example integrate4 -using Integrals, IntegralsCubature, Base.Threads +using Integrals, Cubature, Base.Threads function f(y, u, p) Threads.@threads for i in 1:size(u, 2) y[1, i] = sum(sin.(@view(u[:, i]))) @@ -97,7 +97,7 @@ the change is a one-argument change: ```@example integrate5 using Integrals -using IntegralsCuba +using Cuba f(u, p) = sum(sin.(u)) prob = IntegralProblem(f, ones(3), 3ones(3)) sol = solve(prob, CubaCuhre(); reltol = 1e-3, abstol = 1e-3) diff --git a/src/algorithms.jl b/src/algorithms.jl index 55f69445..2b09744c 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -324,7 +324,7 @@ publisher={Elsevier} struct CubatureJLh <: AbstractCubatureJLAlgorithm error_norm::Int32 end -CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) + """ CubatureJLp() @@ -340,4 +340,3 @@ Defaults to `Cubature.INDIVIDUAL`, other options are struct CubatureJLp <: AbstractCubatureJLAlgorithm error_norm::Int32 end -CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) \ No newline at end of file From 9cfcadbc3aea9ee1f683d514ed9b08fcb7748a3d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 03:27:04 -0400 Subject: [PATCH 19/27] fix alg dispatches --- ext/IntegralsCubatureExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl index 7bc97c8f..4ccfe996 100644 --- a/ext/IntegralsCubatureExt.jl +++ b/ext/IntegralsCubatureExt.jl @@ -5,8 +5,8 @@ using Integrals, Cubature import Integrals: transformation_if_inf, scale_x, scale_x! import Cubature: INDIVIDUAL, PAIRED, L1, L2, LINF -CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) -CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) +Integrals.CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) +Integrals.CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubatureJLAlgorithm, From ebd0331fd378ae79ac8a2a708e535d83ddf19e85 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 03:28:39 -0400 Subject: [PATCH 20/27] import algs --- ext/IntegralsCubaExt.jl | 3 ++- ext/IntegralsCubatureExt.jl | 7 ++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl index d485d312..e999eab1 100644 --- a/ext/IntegralsCubaExt.jl +++ b/ext/IntegralsCubaExt.jl @@ -1,7 +1,8 @@ module IntegralsCubaExt using Integrals, Cuba -import Integrals: transformation_if_inf, scale_x, scale_x! +import Integrals: transformation_if_inf, scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm, + CubaSUAVE, CubaDivonne, CubaCuhre function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm, sensealg, diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl index 4ccfe996..cf19234f 100644 --- a/ext/IntegralsCubatureExt.jl +++ b/ext/IntegralsCubatureExt.jl @@ -2,11 +2,12 @@ module IntegralsCubatureExt using Integrals, Cubature -import Integrals: transformation_if_inf, scale_x, scale_x! +import Integrals: transformation_if_inf, scale_x, scale_x!, CubatureJLh, CubatureJLp, + AbstractCubatureJLAlgorithm import Cubature: INDIVIDUAL, PAIRED, L1, L2, LINF -Integrals.CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) -Integrals.CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) +CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) +CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubatureJLAlgorithm, From aebfa98063377bbcb0aa2ae13e7e18f33122d4c2 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 03:47:15 -0400 Subject: [PATCH 21/27] get cubature updates back --- ext/IntegralsCubatureExt.jl | 165 ++++++++++++++++++++++-------------- 1 file changed, 102 insertions(+), 63 deletions(-) diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl index cf19234f..d17a0b29 100644 --- a/ext/IntegralsCubatureExt.jl +++ b/ext/IntegralsCubatureExt.jl @@ -11,137 +11,176 @@ CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubatureJLAlgorithm, - sensealg, lb, ub, p; + sensealg, domain, p; reltol = 1e-8, abstol = 1e-8, maxiters = typemax(Int)) - nout = prob.nout - if nout == 1 - # the output of prob.f could be either scalar or a vector of length 1, however - # the behavior of the output of the integration routine is undefined (could differ - # across algorithms) - # Cubature will output a real number in when called without nout/fdim - if prob.batch == 0 + + lb, ub = domain + 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 + 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 + + @assert eltype(y)<: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 isinplace(prob) - dx = zeros(eltype(lb), prob.nout) - f = (x) -> (prob.f(dx, x, p); dx[1]) + # 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 + end else - f = (x) -> prob.f(x, p)[1] + f = (x, dx) -> (dx .= prob.f(x, p)) end - if lb isa Number + if mid isa Number if alg isa CubatureJLh - val, err = Cubature.hquadrature(f, lb, ub; + val, err = Cubature.hquadrature_v(f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) else - val, err = Cubature.pquadrature(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(f, lb, ub; + val, err = Cubature.hcubature_v(f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) else - val, err = Cubature.pcubature(f, lb, ub; + val, err = Cubature.pcubature_v(f, lb, ub; reltol = reltol, abstol = abstol, maxevals = maxiters) end end - else + elseif y isa AbstractArray + fsize = size(y)[begin:(end - 1)] + fdim = prod(fsize) if isinplace(prob) - f = (x, dx) -> prob.f(dx, x, p) + # 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, fsize..., size(dx, 2)) + prob.f(dy, x, p) + dx .= reshape(dy, fdim, size(dx, 2)) + end else - f = (x, dx) -> (dx .= prob.f(x, p)) + f = (x, dx) -> (dx .= reshape(prob.f(x, p), fdim, size(dx, 2))) end - if lb isa Number + if mid isa Number if alg isa CubatureJLh - val, err = Cubature.hquadrature_v(f, lb, ub; + val_, err = Cubature.hquadrature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters) + maxevals = maxiters, error_norm = alg.error_norm) else - val, err = Cubature.pquadrature_v(f, lb, ub; + val_, err = Cubature.pquadrature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters) + maxevals = maxiters, error_norm = alg.error_norm) end else if alg isa CubatureJLh - val, err = Cubature.hcubature_v(f, lb, ub; + val_, err = Cubature.hcubature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters) + maxevals = maxiters, error_norm = alg.error_norm) else - val, err = Cubature.pcubature_v(f, lb, ub; + val_, err = Cubature.pcubature_v(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters) + maxevals = maxiters, error_norm = alg.error_norm) end end + val = reshape(val_, fsize...) + else + error("BatchIntegralFunction integrands must be arrays for Cubature.jl") end else - if prob.batch == 0 - if isinplace(prob) - f = (x, dx) -> (prob.f(dx, x, p); dx) - else - f = (x, dx) -> (dx .= prob.f(x, p)) - end + if y isa Real + # no inplace in this case, since the integrand_prototype would be mutable + f = x -> prob.f(x, p) if lb isa Number if alg isa CubatureJLh - val, err = Cubature.hquadrature(nout, f, lb, ub; + val, err = Cubature.hquadrature(f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters) else - val, err = Cubature.pquadrature(nout, f, lb, ub; + val, err = Cubature.pquadrature(f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters) end else if alg isa CubatureJLh - val, err = Cubature.hcubature(nout, f, lb, ub; + val, err = Cubature.hcubature(f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters) else - val, err = Cubature.pcubature(nout, f, lb, ub; + val, err = Cubature.pcubature(f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters) end end - else + elseif y isa AbstractArray + fsize = size(y) + fdim = length(y) if isinplace(prob) - f = (x, dx) -> (prob.f(dx, x, p); dx) + dy = similar(y) + f = (x, v) -> (prob.f(dy, x, p); v .= vec(dy)) else - f = (x, dx) -> (dx .= prob.f(x, p)) + f = (x, v) -> (v .= vec(prob.f(x, p))) end - - if lb isa Number + if mid isa Number if alg isa CubatureJLh - val, err = Cubature.hquadrature_v(nout, f, lb, ub; + val_, err = Cubature.hquadrature(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters, error_norm = alg.error_norm) else - val, err = Cubature.pquadrature_v(nout, f, lb, ub; + val_, err = Cubature.pquadrature(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters, error_norm = alg.error_norm) end else if alg isa CubatureJLh - val, err = Cubature.hcubature_v(nout, f, lb, ub; + val_, err = Cubature.hcubature(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters, error_norm = alg.error_norm) else - val, err = Cubature.pcubature_v(nout, f, lb, ub; + val_, err = Cubature.pcubature(fdim, f, lb, ub; reltol = reltol, abstol = abstol, - maxevals = maxiters, - error_norm = alg.error_norm) + maxevals = maxiters, error_norm = alg.error_norm) end end + val = reshape(val_, fsize) + else + error("IntegralFunctions must be scalars or arrays for Cubature.jl") + end + end + + #= + nout = prob.nout + if nout == 1 + # the output of prob.f could be either scalar or a vector of length 1, however + # the behavior of the output of the integration routine is undefined (could differ + # across algorithms) + # Cubature will output a real number in when called without nout/fdim + if prob.batch == 0 + if isinplace(prob) + dx = zeros(eltype(lb), prob.nout) +@@ -181,6 +334,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, + end end end + =# SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success) end From 73c6d90f13d3e142cf3671e94faae455f04c43d9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 04:07:52 -0400 Subject: [PATCH 22/27] make sure to extend dispatch --- ext/IntegralsCubatureExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl index d17a0b29..158866d1 100644 --- a/ext/IntegralsCubatureExt.jl +++ b/ext/IntegralsCubatureExt.jl @@ -6,8 +6,8 @@ import Integrals: transformation_if_inf, scale_x, scale_x!, CubatureJLh, Cubatur AbstractCubatureJLAlgorithm import Cubature: INDIVIDUAL, PAIRED, L1, L2, LINF -CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) -CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) +Integrals.CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm) +Integrals.CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm) function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubatureJLAlgorithm, From 2f3226e97ffe8fed87c0f266da9bcb09a14d03d4 Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 2 Nov 2023 10:59:17 -0400 Subject: [PATCH 23/27] update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 84ec6d34..3901e93f 100644 --- a/Project.toml +++ b/Project.toml @@ -37,6 +37,7 @@ Distributions = "0.23, 0.24, 0.25" FastGaussQuadrature = "0.5" ForwardDiff = "0.10" HCubature = "1.4" +LinearAlgebra = "1.6" MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3, 0.1" QuadGK = "2.5" Reexport = "0.2, 1.0" From df38e87a016a257d6a669cdf0b29397139b421bd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 11:39:17 -0400 Subject: [PATCH 24/27] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3901e93f..d70d8924 100644 --- a/Project.toml +++ b/Project.toml @@ -44,7 +44,7 @@ Reexport = "0.2, 1.0" Requires = "1" SciMLBase = "2.6" Zygote = "0.4.22, 0.5, 0.6" -julia = "1.6" +julia = "1.9" [extras] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From 53186d2d220fa022b872d420eb49e2e3a293a849 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 2 Nov 2023 11:40:15 -0400 Subject: [PATCH 25/27] Update CI.yml --- .github/workflows/CI.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 17996d26..1b8d306e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -19,7 +19,6 @@ jobs: - Core version: - '1' - - '1.6' steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 From dba856948fd509f4de973ff35876c0c9fd9f5165 Mon Sep 17 00:00:00 2001 From: Lorenzo Van Munoz <66997677+lxvm@users.noreply.github.com> Date: Thu, 2 Nov 2023 13:02:45 -0400 Subject: [PATCH 26/27] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d70d8924..2434f139 100644 --- a/Project.toml +++ b/Project.toml @@ -37,7 +37,7 @@ Distributions = "0.23, 0.24, 0.25" FastGaussQuadrature = "0.5" ForwardDiff = "0.10" HCubature = "1.4" -LinearAlgebra = "1.6" +LinearAlgebra = "1.9" MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3, 0.1" QuadGK = "2.5" Reexport = "0.2, 1.0" From 3cf503caf973ebdc60b6a0d63a4ce88b0140e295 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 3 Nov 2023 05:06:48 -0400 Subject: [PATCH 27/27] Fix CI scripts for no lib --- .github/workflows/CI.yml | 2 +- .github/workflows/CompatHelper.yml | 2 +- .github/workflows/Documentation.yml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 1b8d306e..750ce00a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -38,7 +38,7 @@ jobs: - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 with: - directories: src,ext,lib/IntegralsCuba/src,lib/IntegralsCubature/src + directories: src,ext - uses: codecov/codecov-action@v3 with: files: lcov.info diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index fa4f72a9..73494545 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -23,4 +23,4 @@ jobs: - name: CompatHelper.main() env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - run: julia -e 'using CompatHelper; CompatHelper.main(;subdirs=["", "docs", "lib/IntegralsCubature", "lib/IntegralsCuba"])' + run: julia -e 'using CompatHelper; CompatHelper.main(;subdirs=["", "docs"])' diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 347d50f3..cde17110 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -25,7 +25,7 @@ jobs: run: julia --project=docs/ --code-coverage=user docs/make.jl - uses: julia-actions/julia-processcoverage@v1 with: - directories: src,lib/IntegralsCuba/src,lib/IntegralsCubature/src + directories: src - uses: codecov/codecov-action@v3 with: files: lcov.info