Skip to content

Commit

Permalink
Merge pull request #240 from ArnoStrouwen/format
Browse files Browse the repository at this point in the history
reapply formatting
  • Loading branch information
ChrisRackauckas authored Feb 22, 2024
2 parents acb4339 + f2e5b5c commit 12c635f
Show file tree
Hide file tree
Showing 24 changed files with 314 additions and 239 deletions.
3 changes: 2 additions & 1 deletion .JuliaFormatter.toml
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
style = "sciml"
format_markdown = true
format_markdown = true
format_docstrings = true
2 changes: 1 addition & 1 deletion docs/pages.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,5 +6,5 @@ pages = ["index.md",
"basics/SampledIntegralProblem.md",
"basics/solve.md",
"basics/FAQ.md"],
"Solvers" => Any["solvers/IntegralSolvers.md"],
"Solvers" => Any["solvers/IntegralSolvers.md"]
]
7 changes: 4 additions & 3 deletions docs/src/basics/FAQ.md
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,9 @@ For badly-behaved integrands, such as (nearly) singular and highly oscillatory
functions, most algorithms will fail to converge and either throw an error or
silently return the incorrect result. In some cases Integrals.jl can provide an
error code when things go wrong, but otherwise you can always check if the error
estimate for the integral is less than the requested tolerance, e.g. `sol.resid
< max(abstol, reltol*norm(sol.u))`. Sometimes using a larger tolerance or higher
estimate for the integral is less than the requested tolerance, e.g.
`sol.resid < max(abstol, reltol*norm(sol.u))`.
Sometimes using a larger tolerance or higher
precision arithmetic may help.

## How can I integrate arbitrarily-spaced data?
Expand All @@ -60,4 +61,4 @@ Otherwise, feel free to open an issue or pull request.

## Can I take derivatives with respect to the limits of integration?

Currently this is not implemented.
Currently this is not implemented.
2 changes: 1 addition & 1 deletion docs/src/basics/IntegralFunction.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,4 @@
```@docs
SciMLBase.IntegralFunction
SciMLBase.BatchIntegralFunction
```
```
6 changes: 3 additions & 3 deletions docs/src/basics/SampledIntegralProblem.md
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ Say, by some means we have generated a dataset `x` and `y`:
```@example 1
using Integrals # hide
f = x -> x^2
x = range(0, 1, length=20)
x = range(0, 1, length = 20)
y = f.(x)
```

Expand Down Expand Up @@ -63,9 +63,9 @@ using Integrals # hide
f1 = x -> x^2
f2 = x -> x^3
f3 = x -> x^4
x = range(0, 1, length=20)
x = range(0, 1, length = 20)
y = [f1.(x) f2.(x) f3.(x)]
problem = SampledIntegralProblem(y, x; dim=1)
problem = SampledIntegralProblem(y, x; dim = 1)
method = SimpsonsRule()
solve(problem, method)
```
Expand Down
2 changes: 1 addition & 1 deletion docs/src/tutorials/differentiating_integrals.md
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ function testf(p)
sin(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1])
end
testf(p)
dp1 = Zygote.gradient(testf,p)
dp1 = Zygote.gradient(testf, p)
dp2 = FiniteDiff.finite_difference_gradient(testf, p)
dp3 = ForwardDiff.gradient(testf, p)
dp1[1] ≈ dp2 ≈ dp3
Expand Down
24 changes: 13 additions & 11 deletions ext/IntegralsArblibExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@ module IntegralsArblibExt
using Arblib
using Integrals

function Integrals.__solvebp_call(prob::IntegralProblem, alg::ArblibJL, sensealg, domain, p;
reltol = 1e-8, abstol = 1e-8, maxiters = nothing)

function Integrals.__solvebp_call(
prob::IntegralProblem, alg::ArblibJL, sensealg, domain, p;
reltol = 1e-8, abstol = 1e-8, maxiters = nothing)
lb_, ub_ = domain
lb, ub = map(first, domain)
if !isone(length(lb_)) || !isone(length(ub_))
Expand All @@ -17,16 +17,18 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::ArblibJL, sensealg
res = Acb(0)
y_ = similar(prob.f.integrand_prototype, typeof(res))
f_ = (y, x; kws...) -> (prob.f(y_, x, p; kws...); Arblib.set!(y, only(y_)))
val = Arblib.integrate!(f_, res, lb, ub, atol=abstol, rtol=reltol,
check_analytic=alg.check_analytic, take_prec=alg.take_prec,
warn_on_no_convergence=alg.warn_on_no_convergence, opts=alg.opts)
SciMLBase.build_solution(prob, alg, val, get_radius(val), retcode = ReturnCode.Success)
val = Arblib.integrate!(f_, res, lb, ub, atol = abstol, rtol = reltol,
check_analytic = alg.check_analytic, take_prec = alg.take_prec,
warn_on_no_convergence = alg.warn_on_no_convergence, opts = alg.opts)
SciMLBase.build_solution(
prob, alg, val, get_radius(val), retcode = ReturnCode.Success)
else
f_ = (x; kws...) -> only(prob.f(x, p; kws...))
val = Arblib.integrate(f_, lb, ub, atol=abstol, rtol=reltol,
check_analytic=alg.check_analytic, take_prec=alg.take_prec,
warn_on_no_convergence=alg.warn_on_no_convergence, opts=alg.opts)
SciMLBase.build_solution(prob, alg, val, get_radius(val), retcode = ReturnCode.Success)
val = Arblib.integrate(f_, lb, ub, atol = abstol, rtol = reltol,
check_analytic = alg.check_analytic, take_prec = alg.take_prec,
warn_on_no_convergence = alg.warn_on_no_convergence, opts = alg.opts)
SciMLBase.build_solution(
prob, alg, val, get_radius(val), retcode = ReturnCode.Success)
end
end

Expand Down
7 changes: 4 additions & 3 deletions ext/IntegralsCubaExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@ module IntegralsCubaExt

using Integrals, Cuba
import Integrals: transformation_if_inf,
scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm,
CubaSUAVE, CubaDivonne, CubaCuhre
scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm,
CubaSUAVE, CubaDivonne, CubaCuhre

function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm,
sensealg,
Expand All @@ -22,7 +22,8 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
nvec = min(maxiters, prob.f.max_batch)
# nvec == 1 in Cuba will change vectors to matrices, so we won't support it when
# batching
nvec > 1 || throw(ArgumentError("BatchIntegralFunction must take multiple batch points"))
nvec > 1 ||
throw(ArgumentError("BatchIntegralFunction must take multiple batch points"))

if mid isa Real
_x = zeros(typeof(mid), nvec)
Expand Down
2 changes: 1 addition & 1 deletion ext/IntegralsFastGaussQuadratureExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::Integrals.GaussLeg
sensealg, domain, p;
reltol = nothing, abstol = nothing,
maxiters = nothing) where {C}
if !all(isonelength, domain)
if !all(isone length, domain)
error("GaussLegendre only accepts one-dimensional quadrature problems.")
end
@assert prob.f isa IntegralFunction
Expand Down
4 changes: 2 additions & 2 deletions ext/IntegralsForwardDiffExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@ end

# Manually split for the pushforward
function Integrals.__solvebp(cache, alg, sensealg, domain,
p::Union{D,AbstractArray{<:D}};
kwargs...) where {T, V, P, D<:ForwardDiff.Dual{T, V, P}}
p::Union{D, AbstractArray{<:D}};
kwargs...) where {T, V, P, D <: ForwardDiff.Dual{T, V, P}}

# we need the output type to avoid perturbation confusion while unwrapping nested duals
# We compute a vector-valued integral of the primal and dual simultaneously
Expand Down
15 changes: 8 additions & 7 deletions ext/IntegralsMCIntegrationExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,10 @@ module IntegralsMCIntegrationExt
using MCIntegration, Integrals

function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, domain, p;
reltol = nothing, abstol = nothing, maxiters = 10)
reltol = nothing, abstol = nothing, maxiters = 10)
lb, ub = domain
mid = vec(collect((lb + ub) / 2))
var = Continuous(vec([tuple(a,b) for (a,b) in zip(lb, ub)]))
var = Continuous(vec([tuple(a, b) for (a, b) in zip(lb, ub)]))

if prob.f isa BatchIntegralFunction
error("VEGASMC doesn't support batching. See https://github.com/numericalEFT/MCIntegration.jl/issues/29")
Expand All @@ -16,7 +16,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg,
f_ = (x, f, c) -> begin
n = 0
for v in x
mid[n+=1] = first(v)
mid[n += 1] = first(v)
end
prob.f(f0, mid, p)
f .= vec(f0)
Expand All @@ -26,21 +26,22 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg,
f_ = (x, c) -> begin
n = 0
for v in x
mid[n+=1] = first(v)
mid[n += 1] = first(v)
end
fx = prob.f(mid, p)
fx isa AbstractArray ? vec(fx) : fx
end
end
dof = ones(Int, length(f0)) # each composite Continuous var gets 1 dof
res = integrate(f_; var, dof, inplace=isinplace(prob), type=eltype(f0),
solver=:vegasmc, niter=maxiters, verbose=-2, print=-2, alg.kws...)
res = integrate(f_; var, dof, inplace = isinplace(prob), type = eltype(f0),
solver = :vegasmc, niter = maxiters, verbose = -2, print = -2, alg.kws...)
out, err, chi = if f0 isa Number
map(only, (res.mean, res.stdev, res.chi2))
else
map(a -> reshape(a, size(f0)), (res.mean, res.stdev, res.chi2))
end
SciMLBase.build_solution(prob, alg, out, err, chi=chi, retcode = ReturnCode.Success)
SciMLBase.build_solution(
prob, alg, out, err, chi = chi, retcode = ReturnCode.Success)
end
end

Expand Down
7 changes: 4 additions & 3 deletions ext/IntegralsZygoteExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
# zygote doesn't support mutation, so we build an oop pullback
if sensealg.vjp isa Integrals.ZygoteVJP
if cache.f isa BatchIntegralFunction
dx = similar(cache.f.integrand_prototype, size(cache.f.integrand_prototype)[begin:end-1]..., 1)
dx = similar(cache.f.integrand_prototype,
size(cache.f.integrand_prototype)[begin:(end - 1)]..., 1)
_f = x -> (cache.f(dx, x, p); dx)
# TODO: let the user pass a batched jacobian so we can return a BatchIntegralFunction
dfdp_ = function (x, p)
Expand Down Expand Up @@ -104,7 +105,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
# TODO replace evaluation at endpoint (which anyone can do without Integrals.jl)
# with integration of dfdx uing the same quadrature
dlb = cache.f isa BatchIntegralFunction ? -batch_unwrap(_f([lb])) : -_f(lb)
dub = cache.f isa BatchIntegralFunction ? batch_unwrap(_f([ub])) : _f(ub)
dub = cache.f isa BatchIntegralFunction ? batch_unwrap(_f([ub])) : _f(ub)
return (NoTangent(),
NoTangent(),
NoTangent(),
Expand All @@ -128,7 +129,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
out, quadrature_adjoint
end

batch_unwrap(x::AbstractArray) = dropdims(x; dims=ndims(x))
batch_unwrap(x::AbstractArray) = dropdims(x; dims = ndims(x))

Zygote.@adjoint function Zygote.literal_getproperty(sol::SciMLBase.IntegralSolution,
::Val{:u})
Expand Down
20 changes: 11 additions & 9 deletions src/Integrals.jl
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,6 @@ end
# Give a layer to intercept with AD
__solvebp(args...; kwargs...) = __solvebp_call(args...; kwargs...)


function quadgk_prob_types(f, lb::T, ub::T, p, nrm) where {T}
DT = float(T) # we need to be careful to infer the same result as `evalrule`
RT = Base.promote_op(*, DT, Base.promote_op(f, DT, typeof(p))) # kernel
Expand Down Expand Up @@ -106,10 +105,10 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p
BatchIntegrand((y, x) -> prob.f(y, x, p), similar(bu))
else
fsize = size(bu)[begin:(end - 1)]
BatchIntegrand{Array{eltype(bu),ndims(bu)-1}}() do y, x
BatchIntegrand{Array{eltype(bu), ndims(bu) - 1}}() do y, x
y_ = similar(bu, fsize..., length(y))
prob.f(y_, x, p)
map!(collect, y, eachslice(y_; dims=ndims(bu)))
map!(collect, y, eachslice(y_; dims = ndims(bu)))
return nothing
end
end
Expand All @@ -121,8 +120,8 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p
f = if u isa AbstractVector
BatchIntegrand((y, x) -> y .= prob.f(x, p), u)
else
BatchIntegrand{Array{eltype(u),ndims(u)-1}}() do y, x
map!(collect, y, eachslice(prob.f(x, p); dims=ndims(u)))
BatchIntegrand{Array{eltype(u), ndims(u) - 1}}() do y, x
map!(collect, y, eachslice(prob.f(x, p); dims = ndims(u)))
return nothing
end
end
Expand All @@ -134,7 +133,8 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p
if isinplace(prob)
result = prob.f.integrand_prototype * mid # result may have different units than prototype
f = (y, x) -> prob.f(y, x, p)
val, err = quadgk!(f, result, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
val, err = quadgk!(
f, result, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
rtol = reltol, atol = abstol, order = alg.order, norm = alg.norm)
SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
else
Expand Down Expand Up @@ -189,8 +189,9 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p;
# This is an ugly hack that is compatible with both
f = x -> (prob.f(y, eltype(x) <: SubArray ? parent(first(x)) : x', p); vec(y))
else
y = prob.f(mid isa Number ? typeof(mid)[] :
Matrix{eltype(mid)}(undef, length(mid), 0),
y = prob.f(
mid isa Number ? typeof(mid)[] :
Matrix{eltype(mid)}(undef, length(mid), 0),
p)
f = x -> prob.f(eltype(x) <: SubArray ? parent(first(x)) : x', p)
end
Expand Down Expand Up @@ -221,7 +222,8 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p;
SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success)
end

export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule, TrapezoidalRule, SimpsonsRule
export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule,
TrapezoidalRule, SimpsonsRule
export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre
export CubatureJLh, CubatureJLp
export ArblibJL
Expand Down
3 changes: 1 addition & 2 deletions src/algorithms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -127,9 +127,8 @@ function GaussLegendre(; n = 250, subintervals = 1, nodes = nothing, weights = n
return GaussLegendre(nodes, weights, subintervals)
end


"""
QuadratureRule(q; n=250)
QuadratureRule(q; n=250)
Algorithm to construct and evaluate a quadrature rule `q` of `n` points computed from the
inputs as `x, w = q(n)`. It assumes the nodes and weights are for the standard interval
Expand Down
37 changes: 22 additions & 15 deletions src/algorithms_extension.jl
Original file line number Diff line number Diff line change
Expand Up @@ -123,13 +123,15 @@ end

function CubaVegas(; flags = 0, seed = 0, minevals = 0, nstart = 1000, nincrease = 500,
gridno = 0)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaVegas requires `using Cuba`")
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
error("CubaVegas requires `using Cuba`")
return CubaVegas(flags, seed, minevals, nstart, nincrease, gridno)
end

function CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000, nmin = 2,
flatness = 25.0)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaSUAVE requires `using Cuba`")
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
error("CubaSUAVE requires `using Cuba`")
return CubaSUAVE(flags, seed, minevals, nnew, nmin, flatness)
end

Expand All @@ -138,13 +140,15 @@ function CubaDivonne(; flags = 0, seed = 0, minevals = 0,
maxchisq = 10.0, mindeviation = 0.25,
xgiven = zeros(Cdouble, 0, 0),
nextra = 0, peakfinder = C_NULL)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaDivonne requires `using Cuba`")
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
error("CubaDivonne requires `using Cuba`")
return CubaDivonne(flags, seed, minevals, key1, key2, key3, maxpass, border, maxchisq,
mindeviation, xgiven, nextra, peakfinder)
end

function CubaCuhre(; flags = 0, minevals = 0, key = 0)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaCuhre requires `using Cuba`")
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
error("CubaCuhre requires `using Cuba`")
return CubaCuhre(flags, minevals, key)
end

Expand Down Expand Up @@ -174,8 +178,9 @@ publisher={Elsevier}
struct CubatureJLh <: AbstractCubatureJLAlgorithm
error_norm::Int32
end
function CubatureJLh(; error_norm=0)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && error("CubatureJLh requires `using Cubature`")
function CubatureJLh(; error_norm = 0)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) &&
error("CubatureJLh requires `using Cubature`")
return CubatureJLh(error_norm)
end

Expand All @@ -193,13 +198,12 @@ Defaults to `Cubature.INDIVIDUAL`, other options are
struct CubatureJLp <: AbstractCubatureJLAlgorithm
error_norm::Int32
end
function CubatureJLp(; error_norm=0)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && error("CubatureJLp requires `using Cubature`")
function CubatureJLp(; error_norm = 0)
isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) &&
error("CubatureJLp requires `using Cubature`")
return CubatureJLp(error_norm)
end



"""
ArblibJL(; check_analytic=false, take_prec=false, warn_on_no_convergence=false, opts=C_NULL)
Expand All @@ -221,8 +225,10 @@ struct ArblibJL{O} <: AbstractIntegralExtensionAlgorithm
warn_on_no_convergence::Bool
opts::O
end
function ArblibJL(; check_analytic=false, take_prec=false, warn_on_no_convergence=false, opts=C_NULL)
isnothing(Base.get_extension(@__MODULE__, :IntegralsArblibExt)) && error("ArblibJL requires `using Arblib`")
function ArblibJL(; check_analytic = false, take_prec = false,
warn_on_no_convergence = false, opts = C_NULL)
isnothing(Base.get_extension(@__MODULE__, :IntegralsArblibExt)) &&
error("ArblibJL requires `using Arblib`")
return ArblibJL(check_analytic, take_prec, warn_on_no_convergence, opts)
end

Expand All @@ -232,14 +238,15 @@ end
Markov-chain based Vegas algorithm from MCIntegration.jl
Refer to
[`MCIntegration.integrate`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/#MCIntegration.integrate-Tuple{Function})
[`MCIntegration.integrate`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/#MCIntegration.integrate-Tuple%7BFunction%7D)
for documentation on the keywords, which are passed directly to the solver with a set of
defaults that works for conforming integrands.
"""
struct VEGASMC{K<:NamedTuple} <: AbstractIntegralExtensionAlgorithm
struct VEGASMC{K <: NamedTuple} <: AbstractIntegralExtensionAlgorithm
kws::K
end
function VEGASMC(; kws...)
isnothing(Base.get_extension(@__MODULE__, :IntegralsMCIntegrationExt)) && error("VEGASMC requires `using MCIntegration`")
isnothing(Base.get_extension(@__MODULE__, :IntegralsMCIntegrationExt)) &&
error("VEGASMC requires `using MCIntegration`")
return VEGASMC(NamedTuple(kws))
end
Loading

0 comments on commit 12c635f

Please sign in to comment.