From 393b1b39690d7519c31d9bc5ceed763b21b0b9b2 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 18 Feb 2024 16:43:05 -0500 Subject: [PATCH 1/4] add CExtensionAlgorithm type --- src/algorithms_extension.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/algorithms_extension.jl b/src/algorithms_extension.jl index b79eaa8..811a4ae 100644 --- a/src/algorithms_extension.jl +++ b/src/algorithms_extension.jl @@ -1,8 +1,9 @@ ## Extension Algorithms abstract type AbstractIntegralExtensionAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end +abstract type AbstractIntegralCExtensionAlgorithm <: AbstractIntegralExtensionAlgorithm end -abstract type AbstractCubaAlgorithm <: AbstractIntegralExtensionAlgorithm end +abstract type AbstractCubaAlgorithm <: AbstractIntegralCExtensionAlgorithm end """ CubaVegas() @@ -148,7 +149,7 @@ function CubaCuhre(; flags = 0, minevals = 0, key = 0) return CubaCuhre(flags, minevals, key) end -abstract type AbstractCubatureJLAlgorithm <: AbstractIntegralExtensionAlgorithm end +abstract type AbstractCubatureJLAlgorithm <: AbstractIntegralCExtensionAlgorithm end """ CubatureJLh(; error_norm=Cubature.INDIVIDUAL) @@ -215,7 +216,7 @@ documentation for additional details the algorithm arguments and on implementing high-precision integrands. Additionally, the error estimate is included in the return value of the integral, representing a ball. """ -struct ArblibJL{O} <: AbstractIntegralExtensionAlgorithm +struct ArblibJL{O} <: AbstractIntegralCExtensionAlgorithm check_analytic::Bool take_prec::Bool warn_on_no_convergence::Bool From 5eed7b1f3c20fcef7f79fcfd49307ebbbe157174 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 18 Feb 2024 16:43:54 -0500 Subject: [PATCH 2/4] default to direct AD for ForwardDiff of non-C algorithms --- ext/IntegralsForwardDiffExt.jl | 37 ++++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 300b606..ec19992 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -3,25 +3,35 @@ using Integrals isdefined(Base, :get_extension) ? (using ForwardDiff) : (using ..ForwardDiff) ### Forward-Mode AD Intercepts -#= Direct AD on solvers with QuadGK and HCubature -# incompatible with iip since types must change -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 +# Default to direct AD on solvers +function Integrals.__solvebp(cache, alg, sensealg, domain, + p::Union{D,AbstractArray{<:D}}; + kwargs...) where {T, V, P, D<:ForwardDiff.Dual{T, V, P}} -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...) + if isinplace(cache.f) + prototype = cache.f.integrand_prototype + elt = eltype(prototype) + ForwardDiff.can_dual(elt) || throw(ArgumentError("ForwardDiff of in-place integrands only supports prototypes with real elements")) + dprototype = similar(prototype, replace_dualvaltype(D, elt)) + df = if cache.f isa BatchIntegralFunction + BatchIntegralFunction{true}(cache.f.f, dprototype) + else + IntegralFunction{true}(cache.f.f, dprototype) + end + prob = Integrals.build_problem(cache) + dprob = remake(prob, f = df) + dcache = init(dprob, alg; sensealg = sensealg, do_inf_transformation=Val(false), kwargs...) + Integrals.__solvebp_call(dcache, alg, sensealg, domain, p; kwargs...) + else + Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...) + end end -=# + # TODO: add the pushforward for derivative w.r.t lb, and ub (and then combinations?) # Manually split for the pushforward -function Integrals.__solvebp(cache, alg, sensealg, domain, +function Integrals.__solvebp(cache, alg::Integrals.AbstractIntegralCExtensionAlgorithm, sensealg, domain, p::Union{D,AbstractArray{<:D}}; kwargs...) where {T, V, P, D<:ForwardDiff.Dual{T, V, P}} @@ -73,6 +83,7 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, end end + DT <: Real || throw(ArgumentError("differentiating algorithms in C")) ForwardDiff.can_dual(elt) || ForwardDiff.throw_cannot_dual(elt) rawp = p isa D ? reinterpret(V, [p]) : copy(reinterpret(V, vec(p))) From 1a3919417055f94c305632e57d3b7bd0351092bb Mon Sep 17 00:00:00 2001 From: lxvm Date: Sun, 18 Feb 2024 16:44:54 -0500 Subject: [PATCH 3/4] expand FAQ section for C algorithms --- docs/src/basics/FAQ.md | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index baaff81..82345c5 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -6,7 +6,7 @@ The in-place interface allows evaluating vector-valued integrands without allocating an output array. This can be beneficial for reducing allocations when integrating many functions simultaneously or to make use of existing in-place code. However, note that not all algorithms use in-place operations under the -hood, i.e. `HCubatureJL()`, and may still allocate. +hood, i.e. [`HCubatureJL`](@ref), and may still allocate. You can construct an `IntegralFunction(f, prototype)`, where `f` is of the form `f(y, u, p)` where `prototype` is of the desired type and shape of `y`. @@ -22,16 +22,17 @@ different points, which maximizes the parallelism for a given algorithm. You can construct an out-of-place `BatchIntegralFunction(bf)` where `bf` is of the form `bf(u, p) = stack(x -> f(x, p), eachslice(u; dims=ndims(u)))`, where `f` is the (unbatched) integrand. +For interoperability with as many algorithms as possible, it is important that your out-of-place batch integrand accept an **empty** array of quadrature points and still return an output with a size and type consistent with the non-empty case. You can construct an in-place `BatchIntegralFunction(bf, prototype)`, where `bf` is of the form `bf(y, u, p) = foreach((y,x) -> f(y,x,p), eachslice(y, dims=ndims(y)), eachslice(x, dims=ndims(x)))`. Note that not all algorithms use in-place batched operations under the hood, -i.e. `QuadGKJL()`. +i.e. [`QuadGKJL`](@ref). ## What should I do if my solution is not converged? -Certain algorithms, such as `QuadratureRule` used a fixed number of points to +Certain algorithms, such as [`QuadratureRule`](@ref) used a fixed number of points to calculate an integral and cannot provide an error estimate. In this case, you have to increase the number of points and check the convergence yourself, which will depend on the accuracy of the rule you choose. @@ -46,7 +47,7 @@ precision arithmetic may help. ## How can I integrate arbitrarily-spaced data? -See `SampledIntegralProblem`. +See [`SampledIntegralProblem`](@ref). ## How can I integrate on arbitrary geometries? @@ -58,6 +59,13 @@ because that is what lower-level packages implement. Fixed quadrature rules from other packages can be used with `QuadratureRule`. Otherwise, feel free to open an issue or pull request. +## My integrand works with algorithm X but fails on algorithm Y + +While bugs are not out of the question, certain algorithms, especially those implemented in C, are not compatible with arbitrary Julia types and have to return specific numeric types or arrays thereof. +In some cases, such as [`ArblibJL`](@ref), it is also expected that the integrand work with a custom quadrature point type. +Moreover, some algorithms, such as [`VEGAS`](@ref), only support scalar integrands. +For more details see the [solver page](solvers). + ## Can I take derivatives with respect to the limits of integration? Currently this is not implemented. \ No newline at end of file From 2787daf7eff77af6c3ed7f9fcf6c665e27f78bcd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 22 Feb 2024 00:09:09 -0500 Subject: [PATCH 4/4] Update docs/src/basics/FAQ.md Co-authored-by: Morten Piibeleht --- docs/src/basics/FAQ.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index 82345c5..1ec079d 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -64,7 +64,7 @@ Otherwise, feel free to open an issue or pull request. While bugs are not out of the question, certain algorithms, especially those implemented in C, are not compatible with arbitrary Julia types and have to return specific numeric types or arrays thereof. In some cases, such as [`ArblibJL`](@ref), it is also expected that the integrand work with a custom quadrature point type. Moreover, some algorithms, such as [`VEGAS`](@ref), only support scalar integrands. -For more details see the [solver page](solvers). +For more details see the [solver page](@ref solvers). ## Can I take derivatives with respect to the limits of integration?