From 32e438e6742388da9eb4188554e6f28c9f97a406 Mon Sep 17 00:00:00 2001 From: Pihlajamaa Date: Thu, 7 Sep 2023 13:23:21 +0200 Subject: [PATCH 01/16] add Trapezoidal rule + tests --- Project.toml | 26 +++++------ src/Integrals.jl | 98 ++++++++++++++++++++++++++++++++++++++++- src/algorithms.jl | 64 +++++++++++++++++++++++++++ test/interface_tests.jl | 6 ++- test/runtests.jl | 4 ++ test/sampled_tests.jl | 45 +++++++++++++++++++ 6 files changed, 227 insertions(+), 16 deletions(-) create mode 100644 test/sampled_tests.jl diff --git a/Project.toml b/Project.toml index 80a892db..989da2f8 100644 --- a/Project.toml +++ b/Project.toml @@ -13,10 +13,22 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +[weakdeps] +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[extensions] +IntegralsFastGaussQuadratureExt = "FastGaussQuadrature" +IntegralsForwardDiffExt = "ForwardDiff" +IntegralsZygoteExt = ["Zygote", "ChainRulesCore"] + [compat] ChainRulesCore = "0.10.7, 1" CommonSolve = "0.2" 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" @@ -26,16 +38,11 @@ Requires = "1" SciMLBase = "1.70" Zygote = "0.4.22, 0.5, 0.6" julia = "1.6" -FastGaussQuadrature = "0.5" - -[extensions] -IntegralsForwardDiffExt = "ForwardDiff" -IntegralsZygoteExt = ["Zygote", "ChainRulesCore"] -IntegralsFastGaussQuadratureExt = "FastGaussQuadrature" [extras] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" @@ -44,13 +51,6 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" -FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" [targets] test = ["SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature"] - -[weakdeps] -ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" -FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" diff --git a/src/Integrals.jl b/src/Integrals.jl index 07bb525e..f3a93bbb 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -147,5 +147,101 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, lb, ub, p; SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success) end -export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre +is_sampled_problem(prob::IntegralProblem) = prob.f isa AbstractArray +import SciMLBase.IntegralProblem # this is type piracy, and belongs in SciMLBase +function IntegralProblem(y::AbstractArray, lb, ub, args...; kwargs...) + IntegralProblem{false}(y, lb, ub, args...; kwargs...) +end +function construct_grid(prob, alg, lb, ub, dim) + x = alg.spec + if x isa Integer + @assert length(ub) == length(lb) == 1 "Multidimensional integration is not supported with the Trapezoidal method" + grid = range(lb[1], ub[1], length=x) + else + grid = x + @assert ndims(grid) == 1 "Multidimensional integration is not supported with the Trapezoidal method" + end + + @assert lb[1] ≈ grid[begin] "Lower bound in `IntegralProblem` must coincide with that of the grid" + @assert ub[1] ≈ grid[end] "Upper bound in `IntegralProblem` must coincide with that of the grid" + if is_sampled_problem(prob) + @assert size(prob.f, dim) == length(grid) "Integrand and grid must be of equal length along the integrated dimension" + @assert axes(prob.f, dim) == axes(grid,1) "Grid and integrand array must use same indexing along integrated dimension" + end + return grid +end + +@inline dimension(::Val{D}) where D = D +function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, lb, ub, p; kwargs...) where {S,D} + # since all AbstractRange types are equidistant by design, we can rely on that + @assert prob.batch == 0 + # using `Val`s for dimensionality is required to make `selectdim` not allocate + dim = dimension(D) + p = p + if is_sampled_problem(prob) + @assert alg.spec isa AbstractArray "For pre-sampled problems where the integrand is an array, the integration grid must also be an array." + end + + grid = construct_grid(prob, alg, lb, ub, dim) + + err = Inf64 + if is_sampled_problem(prob) + # inlining is required in order to not allocate + @inline function integrand(i) + # integrate along dimension `dim` + selectdim(prob.f, dim, i) + end + else + if isinplace(prob) + y = zeros(eltype(lb), prob.nout) + integrand = i -> @inbounds (prob.f(y, grid[i], p); y) + else + integrand = i -> @inbounds prob.f(grid[i], p) + end + end + + firstidx, lastidx = firstindex(grid), lastindex(grid) + + out = integrand(firstidx) + if isbits(out) + # fast path for equidistant grids + if alg.spec isa Integer + dx = grid[begin+1] - grid[begin] + for i in (firstidx+1):(lastidx-1) + out += 2*integrand(i) + end + out += integrand(lastidx) + out *= dx/2 + # irregular grids: + elseif alg.spec isa AbstractVector + out *= (grid[firstidx + 1] - grid[firstidx]) + for i in (firstidx+1):(lastidx-1) + @inbounds out += integrand(i) * (grid[i + 1] - grid[i-1]) + end + out += integrand(lastidx) * (grid[lastidx] - grid[lastidx-1]) + out /= 2 + end + else # same, but inplace, broadcasted + out = zeros(eltype(out), size(out)...) # to prevent aliasing + if grid isa AbstractRange + dx = grid[begin+1] - grid[begin] + for i in (firstidx+1):(lastidx-1) + out .+= 2.0 .* integrand(i) + end + out .+= integrand(lastidx) + out .*= dx/2 + else + out .*= (grid[firstidx + 1] - grid[firstidx]) + for i in (firstidx+1):(lastidx-1) + @inbounds out .+= integrand(i) .* (grid[i + 1] - grid[i-1]) + end + out .+= integrand(lastidx) .* (grid[lastidx] - grid[lastidx-1]) + out ./= 2 + end + end + + return SciMLBase.build_solution(prob, alg, out, err, retcode = ReturnCode.Success) +end + +export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, Trapezoidal end # module diff --git a/src/algorithms.jl b/src/algorithms.jl index e49122a4..725058c9 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -122,3 +122,67 @@ function GaussLegendre(; n = 250, subintervals = 1, nodes = nothing, weights = n end return GaussLegendre(nodes, weights, subintervals) end + +""" + Trapezoidal{S, DIM} + +Struct for evaluating an integral via Trapezoidal rule. +The field `spec` contains either the number of gridpoints or an array of specified gridpoints + +The Trapezoidal rule supports integration of pre-sampled data, stored in an array, as well as integration of +functions. It does not support batching or integration over multidimensional spaces. + +To use the Trapezoidal rule to integrate a function on a regular grid with `n` points: + +```@example trapz1 +using Integrals +f = (x, p) -> x^9 +n = 1000 +method = Trapezoidal(n) +problem = IntegralProblem(f, 0.0, 1.0) +solve(problem, method) +``` + +To use the Trapezoidal rule to integrate a function on an predefined irregular grid, see the following example. +Note that the lower and upper bound of integration must coincide with the first and last element of the grid. + +```@example trapz2 +using Integrals +f = (x, p) -> x^9 +x = sort(rand(1000)) +x = [0.0; x; 1.0] +method = Trapezoidal(x) +problem = IntegralProblem(f, 0.0, 1.0) +solve(problem, method) +``` + +To use the Trapezoidal rule to integrate a set of sampled data, see the following example. +By default, the integration occurs over the first dimension of the input array. +```@example trapz3 +using Integrals +x = sort(rand(1000)) +x = [0.0; x; 1.0] +y1 = x' .^ 4 +y2 = x' .^ 9 +y = [y1; y2] +method = Trapezoidal(x; dim=2) +problem = IntegralProblem(y, 0.0, 1.0) +solve(problem, method) +``` +""" +struct Trapezoidal{S, DIM} <: SciMLBase.AbstractIntegralAlgorithm + spec::S + function Trapezoidal(npoints::I; dim=1) where I<:Integer + @assert npoints > 1 + return new{I, Val(dim)}(npoints) + end + function Trapezoidal(grid::V; dim=1) where V<:AbstractVector + npoints = length(grid) + @assert npoints > 1 + @assert isfinite(first(grid)) + @assert isfinite(last(grid)) + @assert issorted(grid) "The gridpoints must be sorted from low to high." + return new{V, Val(dim)}(grid) + end +end + diff --git a/test/interface_tests.jl b/test/interface_tests.jl index 6789f653..ad3ff377 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -8,7 +8,7 @@ reltol = 1e-3 abstol = 1e-3 algs = [QuadGKJL(), HCubatureJL(), CubatureJLh(), CubatureJLp(), #VEGAS(), #CubaVegas(), - CubaSUAVE(), CubaDivonne(), CubaCuhre()] + CubaSUAVE(), CubaDivonne(), CubaCuhre(), Trapezoidal(1000)] alg_req = Dict(QuadGKJL() => (nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, allows_iip = false), @@ -27,7 +27,9 @@ alg_req = Dict(QuadGKJL() => (nout = 1, allows_batch = false, min_dim = 1, max_d CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, allows_iip = true), CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, - allows_iip = true)) + allows_iip = true), + Trapezoidal(1000) => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = 1, + allows_iip = true)) integrands = [ (x, p) -> 1.0, diff --git a/test/runtests.jl b/test/runtests.jl index ef990ac6..f3570386 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,3 +22,7 @@ end @time @safetestset "Gaussian Quadrature Tests" begin include("gaussian_quadrature_tests.jl") end + +@time @safetestset "Sampled Integration Tests" begin + include("sampled_tests.jl") +end diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl new file mode 100644 index 00000000..fa94d3a3 --- /dev/null +++ b/test/sampled_tests.jl @@ -0,0 +1,45 @@ +using Integrals, Test +@testset "Sampled Integration" begin + lb = 0.0 + ub = 1.0 + npoints = 1000 + for method = [Trapezoidal] # Simpson's later + nouts = [1,2,1,2] + for (i,f) = enumerate([(x,p) -> x^5, (x,p) -> [x^5, x^5], (out, x,p) -> (out[1] = x^5; out), (out, x, p) -> (out[1] = x^5; out[2] = x^5; out)]) + + exact = 1/6 + prob = IntegralProblem(f, lb, ub, nout=nouts[i]) + + # AbstractRange + error1 = solve(prob, method(npoints)).u .- exact + @test all(error1 .< 10^-4) + + # AbstractVector equidistant + error2 = solve(prob, method(collect(range(lb, ub, length=npoints)))).u .- exact + @test all(error2 .< 10^-4) + + # AbstractVector irregular + grid = rand(npoints) + grid = [lb; sort(grid); ub] + error3 = solve(prob, method(grid)).u .- exact + @test all(error3 .< 10^-4) + + + end + exact = 1/6 + + grid = rand(npoints) + grid = [lb; sort(grid); ub] + # single dimensional y + y = grid .^ 5 + prob = IntegralProblem(y, lb, ub) + error4 = solve(prob, method(grid, dim=1)).u .- exact + @test all(error4 .< 10^-4) + + # along dim=2 + y = ([grid grid]') .^ 5 + prob = IntegralProblem(y, lb, ub) + error5 = solve(prob, method(grid, dim=2)).u .- exact + @test all(error5 .< 10^-4) + end +end \ No newline at end of file From 881804b4b4dfb840787e70da068105a4d7cd6b06 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Fri, 8 Sep 2023 11:48:58 +0200 Subject: [PATCH 02/16] fix correctness bug --- src/Integrals.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index f3a93bbb..50dbc5b2 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -203,9 +203,10 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, firstidx, lastidx = firstindex(grid), lastindex(grid) out = integrand(firstidx) + if isbits(out) # fast path for equidistant grids - if alg.spec isa Integer + if grid isa AbstractRange dx = grid[begin+1] - grid[begin] for i in (firstidx+1):(lastidx-1) out += 2*integrand(i) @@ -222,7 +223,7 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, out /= 2 end else # same, but inplace, broadcasted - out = zeros(eltype(out), size(out)...) # to prevent aliasing + out = copy(out) # to prevent aliasing if grid isa AbstractRange dx = grid[begin+1] - grid[begin] for i in (firstidx+1):(lastidx-1) From afc0f12cdeaa9ddda4d5ef4dbf4d533710fd0ec0 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Fri, 8 Sep 2023 11:51:20 +0200 Subject: [PATCH 03/16] small change for consistency --- src/Integrals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index 50dbc5b2..11991a3c 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -214,7 +214,7 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, out += integrand(lastidx) out *= dx/2 # irregular grids: - elseif alg.spec isa AbstractVector + else out *= (grid[firstidx + 1] - grid[firstidx]) for i in (firstidx+1):(lastidx-1) @inbounds out += integrand(i) * (grid[i + 1] - grid[i-1]) From ade9f6bcbdfe1b30e2faa39117497656565ec7bc Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Fri, 8 Sep 2023 12:24:21 +0200 Subject: [PATCH 04/16] remove issorted check, and return scalar instead of 0-dim arr --- src/Integrals.jl | 8 ++++++-- src/algorithms.jl | 3 +-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index 11991a3c..beba7d11 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -171,6 +171,9 @@ function construct_grid(prob, alg, lb, ub, dim) return grid end +@inline myselectdim(y::AbstractArray{T,dims}, d, i) where {T,dims} = selectdim(y, d, i) +@inline myselectdim(y::AbstractArray{T,1}, _, i) where {T} = @inbounds y[i] + @inline dimension(::Val{D}) where D = D function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, lb, ub, p; kwargs...) where {S,D} # since all AbstractRange types are equidistant by design, we can rely on that @@ -186,10 +189,11 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, err = Inf64 if is_sampled_problem(prob) + data = prob.f # inlining is required in order to not allocate @inline function integrand(i) - # integrate along dimension `dim` - selectdim(prob.f, dim, i) + # integrate along dimension `dim`, returning a n-1 dimensional array, or scalar if n=1 + myselectdim(data, dim, i) end else if isinplace(prob) diff --git a/src/algorithms.jl b/src/algorithms.jl index 725058c9..d083cfbe 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -144,7 +144,7 @@ solve(problem, method) ``` To use the Trapezoidal rule to integrate a function on an predefined irregular grid, see the following example. -Note that the lower and upper bound of integration must coincide with the first and last element of the grid. +Note that the lower and upper bound of integration must coincide with the first and last element of the grid. ```@example trapz2 using Integrals @@ -181,7 +181,6 @@ struct Trapezoidal{S, DIM} <: SciMLBase.AbstractIntegralAlgorithm @assert npoints > 1 @assert isfinite(first(grid)) @assert isfinite(last(grid)) - @assert issorted(grid) "The gridpoints must be sorted from low to high." return new{V, Val(dim)}(grid) end end From 7265f246d721ab1f3ed7a1c07e23ac04f1030e65 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Fri, 8 Sep 2023 12:35:19 +0200 Subject: [PATCH 05/16] small efficiency improvement --- src/Integrals.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index beba7d11..98eed99e 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -212,11 +212,12 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, # fast path for equidistant grids if grid isa AbstractRange dx = grid[begin+1] - grid[begin] + out /= 2 for i in (firstidx+1):(lastidx-1) - out += 2*integrand(i) + out += integrand(i) end - out += integrand(lastidx) - out *= dx/2 + out += integrand(lastidx)/2 + out *= dx # irregular grids: else out *= (grid[firstidx + 1] - grid[firstidx]) @@ -230,11 +231,12 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, out = copy(out) # to prevent aliasing if grid isa AbstractRange dx = grid[begin+1] - grid[begin] + out ./= 2 for i in (firstidx+1):(lastidx-1) - out .+= 2.0 .* integrand(i) + out .+= integrand(i) end - out .+= integrand(lastidx) - out .*= dx/2 + out .+= integrand(lastidx) ./ 2 + out .*= dx else out .*= (grid[firstidx + 1] - grid[firstidx]) for i in (firstidx+1):(lastidx-1) @@ -244,7 +246,6 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, out ./= 2 end end - return SciMLBase.build_solution(prob, alg, out, err, retcode = ReturnCode.Success) end From c4a9208082b5f746d81ab39e3219e6031229f775 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Fri, 8 Sep 2023 13:31:59 +0200 Subject: [PATCH 06/16] clarify error msgs --- src/Integrals.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Integrals.jl b/src/Integrals.jl index 98eed99e..708842d7 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -154,8 +154,8 @@ function IntegralProblem(y::AbstractArray, lb, ub, args...; kwargs...) end function construct_grid(prob, alg, lb, ub, dim) x = alg.spec + @assert length(ub) == length(lb) == 1 "Multidimensional integration is not supported with the Trapezoidal method" if x isa Integer - @assert length(ub) == length(lb) == 1 "Multidimensional integration is not supported with the Trapezoidal method" grid = range(lb[1], ub[1], length=x) else grid = x @@ -182,7 +182,7 @@ function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, dim = dimension(D) p = p if is_sampled_problem(prob) - @assert alg.spec isa AbstractArray "For pre-sampled problems where the integrand is an array, the integration grid must also be an array." + @assert alg.spec isa AbstractArray "For pre-sampled problems where the integrand is an array, the integration grid must also be specified by an array." end grid = construct_grid(prob, alg, lb, ub, dim) From efaf2e1e23d37b153dba09d218d3f69cc152be11 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Sat, 16 Sep 2023 23:03:00 +0200 Subject: [PATCH 07/16] implement trapezoidal rule for sampled data --- src/Integrals.jl | 123 +++++++--------------------------------- src/algorithms.jl | 63 ++------------------ src/common.jl | 10 ++++ src/trapezoidal.jl | 55 ++++++++++++++++++ test/interface_tests.jl | 6 +- test/sampled_tests.jl | 61 ++++++++------------ 6 files changed, 114 insertions(+), 204 deletions(-) create mode 100644 src/trapezoidal.jl diff --git a/src/Integrals.jl b/src/Integrals.jl index 708842d7..3c8c9ce6 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -8,10 +8,29 @@ using Reexport, MonteCarloIntegration, QuadGK, HCubature @reexport using SciMLBase using LinearAlgebra +############ should be removed once PR in SciMLBase is merged +struct SampledIntegralProblem{Y, X, D, K} <: SciMLBase.AbstractIntegralProblem{false} + y::Y + x::X + dim::D + kwargs::K + function SampledIntegralProblem(y::AbstractArray, x::AbstractVector; + dim = 1, + kwargs...) + @assert dim<=ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`" + @assert length(x)==size(y, dim) "The integrand `y` must have the same length as the sampling points `x` along the integrated dimension." + @assert axes(x, 1)==axes(y, dim) "The integrand `y` must obey the same indexing as the sampling points `x` along the integrated dimension." + new{typeof(y), typeof(x), Val{dim}, typeof(kwargs)}(y, x, Val(dim), kwargs) + end +end +export SampledIntegralProblem +############# + include("common.jl") include("init.jl") include("algorithms.jl") include("infinity_handling.jl") +include("trapezoidal.jl") abstract type QuadSensitivityAlg end struct ReCallVJP{V} @@ -147,107 +166,5 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, lb, ub, p; SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success) end -is_sampled_problem(prob::IntegralProblem) = prob.f isa AbstractArray -import SciMLBase.IntegralProblem # this is type piracy, and belongs in SciMLBase -function IntegralProblem(y::AbstractArray, lb, ub, args...; kwargs...) - IntegralProblem{false}(y, lb, ub, args...; kwargs...) -end -function construct_grid(prob, alg, lb, ub, dim) - x = alg.spec - @assert length(ub) == length(lb) == 1 "Multidimensional integration is not supported with the Trapezoidal method" - if x isa Integer - grid = range(lb[1], ub[1], length=x) - else - grid = x - @assert ndims(grid) == 1 "Multidimensional integration is not supported with the Trapezoidal method" - end - - @assert lb[1] ≈ grid[begin] "Lower bound in `IntegralProblem` must coincide with that of the grid" - @assert ub[1] ≈ grid[end] "Upper bound in `IntegralProblem` must coincide with that of the grid" - if is_sampled_problem(prob) - @assert size(prob.f, dim) == length(grid) "Integrand and grid must be of equal length along the integrated dimension" - @assert axes(prob.f, dim) == axes(grid,1) "Grid and integrand array must use same indexing along integrated dimension" - end - return grid -end - -@inline myselectdim(y::AbstractArray{T,dims}, d, i) where {T,dims} = selectdim(y, d, i) -@inline myselectdim(y::AbstractArray{T,1}, _, i) where {T} = @inbounds y[i] - -@inline dimension(::Val{D}) where D = D -function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, lb, ub, p; kwargs...) where {S,D} - # since all AbstractRange types are equidistant by design, we can rely on that - @assert prob.batch == 0 - # using `Val`s for dimensionality is required to make `selectdim` not allocate - dim = dimension(D) - p = p - if is_sampled_problem(prob) - @assert alg.spec isa AbstractArray "For pre-sampled problems where the integrand is an array, the integration grid must also be specified by an array." - end - - grid = construct_grid(prob, alg, lb, ub, dim) - - err = Inf64 - if is_sampled_problem(prob) - data = prob.f - # inlining is required in order to not allocate - @inline function integrand(i) - # integrate along dimension `dim`, returning a n-1 dimensional array, or scalar if n=1 - myselectdim(data, dim, i) - end - else - if isinplace(prob) - y = zeros(eltype(lb), prob.nout) - integrand = i -> @inbounds (prob.f(y, grid[i], p); y) - else - integrand = i -> @inbounds prob.f(grid[i], p) - end - end - - firstidx, lastidx = firstindex(grid), lastindex(grid) - - out = integrand(firstidx) - - if isbits(out) - # fast path for equidistant grids - if grid isa AbstractRange - dx = grid[begin+1] - grid[begin] - out /= 2 - for i in (firstidx+1):(lastidx-1) - out += integrand(i) - end - out += integrand(lastidx)/2 - out *= dx - # irregular grids: - else - out *= (grid[firstidx + 1] - grid[firstidx]) - for i in (firstidx+1):(lastidx-1) - @inbounds out += integrand(i) * (grid[i + 1] - grid[i-1]) - end - out += integrand(lastidx) * (grid[lastidx] - grid[lastidx-1]) - out /= 2 - end - else # same, but inplace, broadcasted - out = copy(out) # to prevent aliasing - if grid isa AbstractRange - dx = grid[begin+1] - grid[begin] - out ./= 2 - for i in (firstidx+1):(lastidx-1) - out .+= integrand(i) - end - out .+= integrand(lastidx) ./ 2 - out .*= dx - else - out .*= (grid[firstidx + 1] - grid[firstidx]) - for i in (firstidx+1):(lastidx-1) - @inbounds out .+= integrand(i) .* (grid[i + 1] - grid[i-1]) - end - out .+= integrand(lastidx) .* (grid[lastidx] - grid[lastidx-1]) - out ./= 2 - end - end - return SciMLBase.build_solution(prob, alg, out, err, retcode = ReturnCode.Success) -end - -export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, Trapezoidal +export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, TrapezoidalRule end # module diff --git a/src/algorithms.jl b/src/algorithms.jl index d083cfbe..f4b9040d 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -124,64 +124,9 @@ function GaussLegendre(; n = 250, subintervals = 1, nodes = nothing, weights = n end """ - Trapezoidal{S, DIM} - -Struct for evaluating an integral via Trapezoidal rule. -The field `spec` contains either the number of gridpoints or an array of specified gridpoints - -The Trapezoidal rule supports integration of pre-sampled data, stored in an array, as well as integration of -functions. It does not support batching or integration over multidimensional spaces. - -To use the Trapezoidal rule to integrate a function on a regular grid with `n` points: - -```@example trapz1 -using Integrals -f = (x, p) -> x^9 -n = 1000 -method = Trapezoidal(n) -problem = IntegralProblem(f, 0.0, 1.0) -solve(problem, method) -``` - -To use the Trapezoidal rule to integrate a function on an predefined irregular grid, see the following example. -Note that the lower and upper bound of integration must coincide with the first and last element of the grid. - -```@example trapz2 -using Integrals -f = (x, p) -> x^9 -x = sort(rand(1000)) -x = [0.0; x; 1.0] -method = Trapezoidal(x) -problem = IntegralProblem(f, 0.0, 1.0) -solve(problem, method) -``` - -To use the Trapezoidal rule to integrate a set of sampled data, see the following example. -By default, the integration occurs over the first dimension of the input array. -```@example trapz3 -using Integrals -x = sort(rand(1000)) -x = [0.0; x; 1.0] -y1 = x' .^ 4 -y2 = x' .^ 9 -y = [y1; y2] -method = Trapezoidal(x; dim=2) -problem = IntegralProblem(y, 0.0, 1.0) -solve(problem, method) -``` + TrapezoidalRule + +Struct for evaluating an integral via the trapezoidal rule. """ -struct Trapezoidal{S, DIM} <: SciMLBase.AbstractIntegralAlgorithm - spec::S - function Trapezoidal(npoints::I; dim=1) where I<:Integer - @assert npoints > 1 - return new{I, Val(dim)}(npoints) - end - function Trapezoidal(grid::V; dim=1) where V<:AbstractVector - npoints = length(grid) - @assert npoints > 1 - @assert isfinite(first(grid)) - @assert isfinite(last(grid)) - return new{V, Val(dim)}(grid) - end +struct TrapezoidalRule <: SciMLBase.AbstractIntegralAlgorithm end - diff --git a/src/common.jl b/src/common.jl index d13a0165..ea41f152 100644 --- a/src/common.jl +++ b/src/common.jl @@ -80,6 +80,12 @@ function SciMLBase.solve(prob::IntegralProblem, solve!(init(prob, alg; kwargs...)) end +function SciMLBase.solve(prob::SampledIntegralProblem, + alg::SciMLBase.AbstractIntegralAlgorithm; + kwargs...) + __solvebp(prob, alg; kwargs...) +end + function SciMLBase.solve!(cache::IntegralCache) __solvebp(cache, cache.alg, cache.sensealg, cache.lb, cache.ub, cache.p; cache.kwargs...) @@ -94,3 +100,7 @@ end function __solvebp_call(cache::IntegralCache, args...; kwargs...) __solvebp_call(build_problem(cache), args...; kwargs...) end + +@inline _selectdim(y::AbstractArray{T, dims}, d, i) where {T, dims} = selectdim(y, d, i) +@inline _selectdim(y::AbstractArray{T, 1}, _, i) where {T} = @inbounds y[i] +@inline dimension(::Val{D}) where {D} = D diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl new file mode 100644 index 00000000..a58a4fb6 --- /dev/null +++ b/src/trapezoidal.jl @@ -0,0 +1,55 @@ +function __solvebp_call(prob::SampledIntegralProblem, alg::TrapezoidalRule; kwargs...) + dim = dimension(prob.dim) + err = Inf64 + data = prob.y + grid = prob.x + # inlining is required in order to not allocate + integrand = @inline function (i) + # integrate along dimension `dim`, returning a n-1 dimensional array, or scalar if n=1 + _selectdim(data, dim, i) + end + + firstidx, lastidx = firstindex(grid), lastindex(grid) + + out = integrand(firstidx) + + if isbits(out) + # fast path for equidistant grids + if grid isa AbstractRange + dx = step(grid) + out /= 2 + for i in (firstidx + 1):(lastidx - 1) + out += integrand(i) + end + out += integrand(lastidx) / 2 + out *= dx + # irregular grids: + else + out *= (grid[firstidx + 1] - grid[firstidx]) + for i in (firstidx + 1):(lastidx - 1) + @inbounds out += integrand(i) * (grid[i + 1] - grid[i - 1]) + end + out += integrand(lastidx) * (grid[lastidx] - grid[lastidx - 1]) + out /= 2 + end + else # same, but inplace, broadcasted + out = copy(out) # to prevent aliasing + if grid isa AbstractRange + dx = grid[begin + 1] - grid[begin] + out ./= 2 + for i in (firstidx + 1):(lastidx - 1) + out .+= integrand(i) + end + out .+= integrand(lastidx) ./ 2 + out .*= dx + else + out .*= (grid[firstidx + 1] - grid[firstidx]) + for i in (firstidx + 1):(lastidx - 1) + @inbounds out .+= integrand(i) .* (grid[i + 1] - grid[i - 1]) + end + out .+= integrand(lastidx) .* (grid[lastidx] - grid[lastidx - 1]) + out ./= 2 + end + end + return SciMLBase.build_solution(prob, alg, out, err, retcode = ReturnCode.Success) +end diff --git a/test/interface_tests.jl b/test/interface_tests.jl index ad3ff377..6789f653 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -8,7 +8,7 @@ reltol = 1e-3 abstol = 1e-3 algs = [QuadGKJL(), HCubatureJL(), CubatureJLh(), CubatureJLp(), #VEGAS(), #CubaVegas(), - CubaSUAVE(), CubaDivonne(), CubaCuhre(), Trapezoidal(1000)] + CubaSUAVE(), CubaDivonne(), CubaCuhre()] alg_req = Dict(QuadGKJL() => (nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, allows_iip = false), @@ -27,9 +27,7 @@ alg_req = Dict(QuadGKJL() => (nout = 1, allows_batch = false, min_dim = 1, max_d CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, allows_iip = true), CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, - allows_iip = true), - Trapezoidal(1000) => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = 1, - allows_iip = true)) + allows_iip = true)) integrands = [ (x, p) -> 1.0, diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index fa94d3a3..03d838cf 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -1,45 +1,30 @@ using Integrals, Test @testset "Sampled Integration" begin - lb = 0.0 - ub = 1.0 + lb = 0.4 + ub = 1.1 npoints = 1000 - for method = [Trapezoidal] # Simpson's later - nouts = [1,2,1,2] - for (i,f) = enumerate([(x,p) -> x^5, (x,p) -> [x^5, x^5], (out, x,p) -> (out[1] = x^5; out), (out, x, p) -> (out[1] = x^5; out[2] = x^5; out)]) - - exact = 1/6 - prob = IntegralProblem(f, lb, ub, nout=nouts[i]) - - # AbstractRange - error1 = solve(prob, method(npoints)).u .- exact - @test all(error1 .< 10^-4) - - # AbstractVector equidistant - error2 = solve(prob, method(collect(range(lb, ub, length=npoints)))).u .- exact - @test all(error2 .< 10^-4) - - # AbstractVector irregular - grid = rand(npoints) - grid = [lb; sort(grid); ub] - error3 = solve(prob, method(grid)).u .- exact - @test all(error3 .< 10^-4) - + grid1 = range(lb, ub, length = npoints) + grid2 = rand(npoints).*(ub-lb) .+ lb + grid2 = [lb; sort(grid2); ub] + + exact_sols = [1 / 6 * (ub^6 - lb^6), sin(ub) - sin(lb)] + for method in [TrapezoidalRule] # Simpson's later + for grid in [grid1, grid2] + for (i, f) in enumerate([x -> x^5, x -> cos(x)]) + exact = exact_sols[i] + # single dimensional y + y = f.(grid) + prob = SampledIntegralProblem(y, grid) + error = solve(prob, method()).u .- exact + @test all(error .< 10^-4) + + # along dim=2 + y = f.([grid grid]') + prob = SampledIntegralProblem(y, grid; dim=2) + error = solve(prob, method()).u .- exact + @test all(error .< 10^-4) + end end - exact = 1/6 - - grid = rand(npoints) - grid = [lb; sort(grid); ub] - # single dimensional y - y = grid .^ 5 - prob = IntegralProblem(y, lb, ub) - error4 = solve(prob, method(grid, dim=1)).u .- exact - @test all(error4 .< 10^-4) - - # along dim=2 - y = ([grid grid]') .^ 5 - prob = IntegralProblem(y, lb, ub) - error5 = solve(prob, method(grid, dim=2)).u .- exact - @test all(error5 .< 10^-4) end end \ No newline at end of file From aa677c4df3f38045c5fea3a4bb2ff6517aa0b531 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Tue, 19 Sep 2023 20:57:15 +0200 Subject: [PATCH 08/16] add trapezoidal rule for sampled data, attempt 2 --- docs/pages.jl | 1 + docs/src/basics/SampledIntegralProblem.md | 73 +++++++++++++++++++++++ src/Integrals.jl | 1 + src/algorithms.jl | 15 ++++- src/common.jl | 3 - src/sampled.jl | 57 ++++++++++++++++++ src/trapezoidal.jl | 70 ++++++---------------- test/sampled_tests.jl | 2 + 8 files changed, 167 insertions(+), 55 deletions(-) create mode 100644 docs/src/basics/SampledIntegralProblem.md create mode 100644 src/sampled.jl diff --git a/docs/pages.jl b/docs/pages.jl index 76c058fc..99ef1719 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,6 +2,7 @@ pages = ["index.md", "Tutorials" => Any["tutorials/numerical_integrals.md", "tutorials/differentiating_integrals.md"], "Basics" => Any["basics/IntegralProblem.md", + "basics/SampledIntegralProblem.md", "basics/solve.md", "basics/FAQ.md"], "Solvers" => Any["solvers/IntegralSolvers.md"], diff --git a/docs/src/basics/SampledIntegralProblem.md b/docs/src/basics/SampledIntegralProblem.md new file mode 100644 index 00000000..02317cd4 --- /dev/null +++ b/docs/src/basics/SampledIntegralProblem.md @@ -0,0 +1,73 @@ +# Integrating pre-sampled data + +In some cases, instead of a function that acts as integrand, +one only possesses a list of data points `y` at a set of sampling +locations `x`, that must be integrated. This package contains functionality +for doing that. + +## Example + +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) +y = f.(x) +``` + +Now, we can integrate this data set as follows: + +```example 1 +problem = SampledIntegralProblem(y, x) +method = TrapezoidalRule() +solve(problem, method) +``` + +The exact aswer is of course \$ 1/3 \$. + +## Details + +### Non-equidistant grids + +If the sampling points `x` are provided as an `AbstractRange` +(constructed with the `range` function for example), faster methods are used that take advantage of +the fact that the points are equidistantly spaced. Otherwise, general methods are used for +non-uniform grids. + +Example: + +```example 2 +using Integrals # hide +f = x -> x^7 +x = [0.0; sort(rand(1000)); 1.0] +y = f.(x) +problem = SampledIntegralProblem(y, x) +method = TrapezoidalRule() +solve(problem, method) +``` + +### Evaluating multiple integrals at once + +If the provided data set `y` is a multidimensional array, the integrals are evaluated across only one +of its axes. For performance reasons, the last axis of the array `y` is chosen by default, but this can be modified with the `dim` +keyword argument to the problem definition. + +```example 3 +using Integrals # hide +f1 = x -> x^2 +f2 = x -> x^3 +f3 = x -> x^4 +x = range(0, 1, length=20) +y = [f1.(x) f2.(x) f3.(x)] +problem = SampledIntegralProblem(y, x; dim=1) +method = TrapezoidalRule() +solve(problem, method) +``` + +### Supported methods + +Right now, only the `TrapezoidalRule` is supported, [see wikipedia](https://en.wikipedia.org/wiki/Trapezoidal_rule). + +```@docs +TrapezoidalRule +``` \ No newline at end of file diff --git a/src/Integrals.jl b/src/Integrals.jl index f407ae53..02b5bc5d 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -13,6 +13,7 @@ include("init.jl") include("algorithms.jl") include("infinity_handling.jl") include("quadrules.jl") +include("sampled.jl") include("trapezoidal.jl") abstract type QuadSensitivityAlg end diff --git a/src/algorithms.jl b/src/algorithms.jl index 37f474d0..0adefdbd 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -127,10 +127,23 @@ end TrapezoidalRule Struct for evaluating an integral via the trapezoidal rule. + + +Example with sampled data: + +``` +using Integrals +f = x -> x^2 +x = range(0, 1, length=20) +y = f.(x) +problem = SampledIntegralProblem(y, x) +method = TrapezoidalRul() +solve(problem, method) +``` """ struct TrapezoidalRule <: SciMLBase.AbstractIntegralAlgorithm end - + """ QuadratureRule(q; n=250) diff --git a/src/common.jl b/src/common.jl index ea41f152..39876049 100644 --- a/src/common.jl +++ b/src/common.jl @@ -101,6 +101,3 @@ function __solvebp_call(cache::IntegralCache, args...; kwargs...) __solvebp_call(build_problem(cache), args...; kwargs...) end -@inline _selectdim(y::AbstractArray{T, dims}, d, i) where {T, dims} = selectdim(y, d, i) -@inline _selectdim(y::AbstractArray{T, 1}, _, i) where {T} = @inbounds y[i] -@inline dimension(::Val{D}) where {D} = D diff --git a/src/sampled.jl b/src/sampled.jl new file mode 100644 index 00000000..9cdcbe21 --- /dev/null +++ b/src/sampled.jl @@ -0,0 +1,57 @@ +abstract type AbstractWeights end + +# must have field `n` for length, and a field `h` for stepsize +abstract type UniformWeights <: AbstractWeights end +@inline Base.iterate(w::UniformWeights) = (0 == w.n) ? nothing : (w[1], 1) +@inline Base.iterate(w::UniformWeights, i) = (i == w.n) ? nothing : (w[i+1], i+1) +Base.length(w::UniformWeights) = w.n +Base.eltype(w::UniformWeights) = typeof(w.h) +Base.size(w::UniformWeights) = (length(w), ) + +# must contain field `x` which are the sampling points +abstract type NonuniformWeights <: AbstractWeights end +@inline Base.iterate(w::NonuniformWeights) = (0 == length(w.x)) ? nothing : (w[firstindex(w.x)], firstindex(w.x)) +@inline Base.iterate(w::NonuniformWeights, i) = (i == lastindex(w.x)) ? nothing : (w[i+1], i+1) +Base.length(w::NonuniformWeights) = length(w.x) +Base.eltype(w::NonuniformWeights) = eltype(w.x) +Base.size(w::NonuniformWeights) = (length(w), ) + +_eachslice(data::AbstractArray; dims=ndims(data)) = eachslice(data; dims=dims) +_eachslice(data::AbstractArray{T, 1}; dims=ndims(data)) where T = data + + +# these can be removed when the Val(dim) is removed from SciMLBase +dimension(::Val{D}) where {D} = D +dimension(D::Int) = D + + +function evalrule(data::AbstractArray, weights, dim) + f = _eachslice(data, dims=dim) + firstidx, lastidx = firstindex(f), lastindex(f) + out = f[firstidx]*weights[firstidx] + if isbits(out) + for i in firstidx+1:lastidx + @inbounds out += f[i]*weights[i] + end + else + for i in firstidx+1:lastidx + @inbounds out .+= f[i] .* weights[i] + end + end + return out + +end + + +# can be reused for other sampled rules +function __solvebp_call(prob::SampledIntegralProblem, alg::TrapezoidalRule; kwargs...) + dim = dimension(prob.dim) + err = nothing + data = prob.y + grid = prob.x + weights = find_weights(grid, alg) + I = evalrule(data, weights, dim) + return SciMLBase.build_solution(prob, alg, I, err, retcode = ReturnCode.Success) +end + + diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl index a58a4fb6..1056fca3 100644 --- a/src/trapezoidal.jl +++ b/src/trapezoidal.jl @@ -1,55 +1,23 @@ -function __solvebp_call(prob::SampledIntegralProblem, alg::TrapezoidalRule; kwargs...) - dim = dimension(prob.dim) - err = Inf64 - data = prob.y - grid = prob.x - # inlining is required in order to not allocate - integrand = @inline function (i) - # integrate along dimension `dim`, returning a n-1 dimensional array, or scalar if n=1 - _selectdim(data, dim, i) - end +struct TrapezoidalUniformWeights <: UniformWeights + n::Int + h::Float64 +end + +@inline Base.getindex(w::TrapezoidalUniformWeights, i) = ifelse((i == 1) || (i == w.n), w.h*0.5 , w.h) - firstidx, lastidx = firstindex(grid), lastindex(grid) - out = integrand(firstidx) +struct TrapezoidalNonuniformWeights{X<:AbstractArray} <: NonuniformWeights + x::X +end - if isbits(out) - # fast path for equidistant grids - if grid isa AbstractRange - dx = step(grid) - out /= 2 - for i in (firstidx + 1):(lastidx - 1) - out += integrand(i) - end - out += integrand(lastidx) / 2 - out *= dx - # irregular grids: - else - out *= (grid[firstidx + 1] - grid[firstidx]) - for i in (firstidx + 1):(lastidx - 1) - @inbounds out += integrand(i) * (grid[i + 1] - grid[i - 1]) - end - out += integrand(lastidx) * (grid[lastidx] - grid[lastidx - 1]) - out /= 2 - end - else # same, but inplace, broadcasted - out = copy(out) # to prevent aliasing - if grid isa AbstractRange - dx = grid[begin + 1] - grid[begin] - out ./= 2 - for i in (firstidx + 1):(lastidx - 1) - out .+= integrand(i) - end - out .+= integrand(lastidx) ./ 2 - out .*= dx - else - out .*= (grid[firstidx + 1] - grid[firstidx]) - for i in (firstidx + 1):(lastidx - 1) - @inbounds out .+= integrand(i) .* (grid[i + 1] - grid[i - 1]) - end - out .+= integrand(lastidx) .* (grid[lastidx] - grid[lastidx - 1]) - out ./= 2 - end - end - return SciMLBase.build_solution(prob, alg, out, err, retcode = ReturnCode.Success) +@inline function Base.getindex(w::TrapezoidalNonuniformWeights, i) + x = w.x + (i == firstindex(x)) && return (x[i + 1] - x[i])*0.5 + (i == lastindex(x)) && return (x[i] - x[i - 1])*0.5 + return (x[i + 1] - x[i - 1])*0.5 end + +function find_weights(x::AbstractVector, ::TrapezoidalRule) + x isa AbstractRange && return TrapezoidalUniformWeights(length(x), step(x)) + return TrapezoidalNonuniformWeights(x) +end \ No newline at end of file diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index 03d838cf..c394f6a8 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -16,12 +16,14 @@ using Integrals, Test # single dimensional y y = f.(grid) prob = SampledIntegralProblem(y, grid) + @show solve(prob, method()).u, exact error = solve(prob, method()).u .- exact @test all(error .< 10^-4) # along dim=2 y = f.([grid grid]') prob = SampledIntegralProblem(y, grid; dim=2) + @show solve(prob, method()).u, exact error = solve(prob, method()).u .- exact @test all(error .< 10^-4) end From e7aefc9586eeabf7de5cca210dd2554820c5405c Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Tue, 19 Sep 2023 23:25:44 +0200 Subject: [PATCH 09/16] now also works for Julia<= 1.9 Apparently the `eachslice behaviour changed in 1.9 causing tests to fail... --- src/sampled.jl | 28 +++++++++++++++++++--------- test/sampled_tests.jl | 2 -- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/sampled.jl b/src/sampled.jl index 9cdcbe21..992dae45 100644 --- a/src/sampled.jl +++ b/src/sampled.jl @@ -27,19 +27,29 @@ dimension(D::Int) = D function evalrule(data::AbstractArray, weights, dim) f = _eachslice(data, dims=dim) - firstidx, lastidx = firstindex(f), lastindex(f) - out = f[firstidx]*weights[firstidx] + f1, statef = iterate(f) + w1, statew = iterate(weights) + out = w1 * f1 + nextf = iterate(f, statef) + nextw = iterate(weights, statew) if isbits(out) - for i in firstidx+1:lastidx - @inbounds out += f[i]*weights[i] + while nextf !== nothing + fi, statef = nextf + wi, statew = nextw + out += wi * fi + nextf = iterate(f, statef) + nextw = iterate(weights, statew) end - else - for i in firstidx+1:lastidx - @inbounds out .+= f[i] .* weights[i] + else + while nextf !== nothing + fi, statef = nextf + wi, statew = nextw + out .+= wi .* fi + nextf = iterate(f, statef) + nextw = iterate(weights, statew) end end - return out - + return out end diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index c394f6a8..03d838cf 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -16,14 +16,12 @@ using Integrals, Test # single dimensional y y = f.(grid) prob = SampledIntegralProblem(y, grid) - @show solve(prob, method()).u, exact error = solve(prob, method()).u .- exact @test all(error .< 10^-4) # along dim=2 y = f.([grid grid]') prob = SampledIntegralProblem(y, grid; dim=2) - @show solve(prob, method()).u, exact error = solve(prob, method()).u .- exact @test all(error .< 10^-4) end From 608af768b907993f5986dfd4915bdcd31c5251dd Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 20 Sep 2023 16:44:34 -0400 Subject: [PATCH 10/16] add init interface for SampledIntegralProblem --- src/common.jl | 67 +++++++++++++++++++++++++++++++++++++++++----- src/sampled.jl | 32 +++++++++++++--------- src/trapezoidal.jl | 4 +-- 3 files changed, 83 insertions(+), 20 deletions(-) diff --git a/src/common.jl b/src/common.jl index 39876049..e00ae2c9 100644 --- a/src/common.jl +++ b/src/common.jl @@ -80,12 +80,6 @@ function SciMLBase.solve(prob::IntegralProblem, solve!(init(prob, alg; kwargs...)) end -function SciMLBase.solve(prob::SampledIntegralProblem, - alg::SciMLBase.AbstractIntegralAlgorithm; - kwargs...) - __solvebp(prob, alg; kwargs...) -end - function SciMLBase.solve!(cache::IntegralCache) __solvebp(cache, cache.alg, cache.sensealg, cache.lb, cache.ub, cache.p; cache.kwargs...) @@ -101,3 +95,64 @@ function __solvebp_call(cache::IntegralCache, args...; kwargs...) __solvebp_call(build_problem(cache), args...; kwargs...) end + +mutable struct SampledIntegralCache{Y, X, D, PK, A, K, Tc} + y::Y + x::X + dim::D + prob_kwargs::PK + alg::A + kwargs::K + isfresh::Bool # state of whether weights have been calculated + cacheval::Tc # store alg weights here +end + +function Base.setproperty!(cache::SampledIntegralCache, name::Symbol, x) + if name === :x + setfield!(cache, :isfresh, true) + end + setfield!(cache, name, x) +end + +function SciMLBase.init(prob::SampledIntegralProblem, + alg::SciMLBase.AbstractIntegralAlgorithm; + kwargs...) + NamedTuple(kwargs) == NamedTuple() || throw(ArgumentError("There are no keyword arguments allowed to `solve`")) + + cacheval = init_cacheval(alg, prob) + isfresh = true + + SampledIntegralCache( + prob.y, + prob.x, + prob.dim, + prob.kwargs, + alg, + kwargs, + isfresh, + cacheval) +end + + +""" +```julia +solve(prob::SampledIntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...) +``` + +## Keyword Arguments + +There are no keyword arguments used to solve `SampledIntegralProblem`s +""" +function SciMLBase.solve(prob::SampledIntegralProblem, + alg::SciMLBase.AbstractIntegralAlgorithm; + kwargs...) + solve!(init(prob, alg; kwargs...)) +end + +function SciMLBase.solve!(cache::SampledIntegralCache) + __solvebp(cache, cache.alg; cache.kwargs...) +end + +function build_problem(cache::SampledIntegralCache) + SampledIntegralProblem(cache.y, cache.x; dim = dimension(cache.dim), cache.prob_kwargs...) +end diff --git a/src/sampled.jl b/src/sampled.jl index 992dae45..d132a12f 100644 --- a/src/sampled.jl +++ b/src/sampled.jl @@ -9,7 +9,7 @@ Base.eltype(w::UniformWeights) = typeof(w.h) Base.size(w::UniformWeights) = (length(w), ) # must contain field `x` which are the sampling points -abstract type NonuniformWeights <: AbstractWeights end +abstract type NonuniformWeights <: AbstractWeights end @inline Base.iterate(w::NonuniformWeights) = (0 == length(w.x)) ? nothing : (w[firstindex(w.x)], firstindex(w.x)) @inline Base.iterate(w::NonuniformWeights, i) = (i == lastindex(w.x)) ? nothing : (w[i+1], i+1) Base.length(w::NonuniformWeights) = length(w.x) @@ -22,7 +22,7 @@ _eachslice(data::AbstractArray{T, 1}; dims=ndims(data)) where T = data # these can be removed when the Val(dim) is removed from SciMLBase dimension(::Val{D}) where {D} = D -dimension(D::Int) = D +dimension(D::Int) = D function evalrule(data::AbstractArray, weights, dim) @@ -40,7 +40,7 @@ function evalrule(data::AbstractArray, weights, dim) nextf = iterate(f, statef) nextw = iterate(weights, statew) end - else + else while nextf !== nothing fi, statef = nextf wi, statew = nextw @@ -49,19 +49,27 @@ function evalrule(data::AbstractArray, weights, dim) nextw = iterate(weights, statew) end end - return out + return out end -# can be reused for other sampled rules -function __solvebp_call(prob::SampledIntegralProblem, alg::TrapezoidalRule; kwargs...) - dim = dimension(prob.dim) +# can be reused for other sampled rules, which should implement find_weights(x, alg) + +function init_cacheval(alg::SciMLBase.AbstractIntegralAlgorithm, prob::SampledIntegralProblem) + find_weights(prob.x, alg) +end + +function __solvebp_call(cache::SampledIntegralCache, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...) + dim = dimension(cache.dim) err = nothing - data = prob.y - grid = prob.x - weights = find_weights(grid, alg) + data = cache.y + grid = cache.x + if cache.isfresh + cache.cacheval = find_weights(grid, alg) + cache.isfresh = false + end + weights = cache.cacheval I = evalrule(data, weights, dim) + prob = build_problem(cache) return SciMLBase.build_solution(prob, alg, I, err, retcode = ReturnCode.Success) end - - diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl index 1056fca3..d1365e1c 100644 --- a/src/trapezoidal.jl +++ b/src/trapezoidal.jl @@ -10,7 +10,7 @@ struct TrapezoidalNonuniformWeights{X<:AbstractArray} <: NonuniformWeights x::X end -@inline function Base.getindex(w::TrapezoidalNonuniformWeights, i) +@inline function Base.getindex(w::TrapezoidalNonuniformWeights, i) x = w.x (i == firstindex(x)) && return (x[i + 1] - x[i])*0.5 (i == lastindex(x)) && return (x[i] - x[i - 1])*0.5 @@ -20,4 +20,4 @@ end function find_weights(x::AbstractVector, ::TrapezoidalRule) x isa AbstractRange && return TrapezoidalUniformWeights(length(x), step(x)) return TrapezoidalNonuniformWeights(x) -end \ No newline at end of file +end From 2eb0090264c4b492ac1d576bc83e02d2fe49981a Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 20 Sep 2023 18:05:33 -0400 Subject: [PATCH 11/16] bump sciml version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 989da2f8..5f66712c 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,7 @@ MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3" QuadGK = "2.5" Reexport = "0.2, 1.0" Requires = "1" -SciMLBase = "1.70" +SciMLBase = "1.98" Zygote = "0.4.22, 0.5, 0.6" julia = "1.6" From 346c0ac9121b2ef33200b33b8049deacd546dd39 Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 20 Sep 2023 18:10:53 -0400 Subject: [PATCH 12/16] use zip for evalrule iterator --- src/sampled.jl | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/src/sampled.jl b/src/sampled.jl index d132a12f..16453b54 100644 --- a/src/sampled.jl +++ b/src/sampled.jl @@ -26,27 +26,23 @@ dimension(D::Int) = D function evalrule(data::AbstractArray, weights, dim) - f = _eachslice(data, dims=dim) - f1, statef = iterate(f) - w1, statew = iterate(weights) + fw = zip(_eachslice(data, dims=dim), weights) + next = iterate(fw) + next === nothing && throw(ArgumentError("No points to integrate")) + (f1, w1), state = next out = w1 * f1 - nextf = iterate(f, statef) - nextw = iterate(weights, statew) + next = iterate(fw, state) if isbits(out) - while nextf !== nothing - fi, statef = nextf - wi, statew = nextw + while next !== nothing + (fi, wi), state = next out += wi * fi - nextf = iterate(f, statef) - nextw = iterate(weights, statew) + next = iterate(fw, state) end else - while nextf !== nothing - fi, statef = nextf - wi, statew = nextw + while next !== nothing + (fi, wi), state = next out .+= wi .* fi - nextf = iterate(f, statef) - nextw = iterate(weights, statew) + next = iterate(fw, state) end end return out From 9ab8ab2d3809b7c1e5ea32e6813f6a66e6e219d8 Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 20 Sep 2023 18:42:12 -0400 Subject: [PATCH 13/16] add docs and tests --- docs/src/tutorials/caching_interface.md | 47 +++++++++++++++++++++++++ test/sampled_tests.jl | 43 +++++++++++++++++++++- 2 files changed, 89 insertions(+), 1 deletion(-) diff --git a/docs/src/tutorials/caching_interface.md b/docs/src/tutorials/caching_interface.md index 13fe41de..02020524 100644 --- a/docs/src/tutorials/caching_interface.md +++ b/docs/src/tutorials/caching_interface.md @@ -50,3 +50,50 @@ Note that the types of these variables is not allowed to change. If it is necessary to change the integrand `f` instead of defining a new `IntegralProblem`, consider using [FunctionWrappers.jl](https://github.com/yuyichao/FunctionWrappers.jl). + +## Caching for sampled integral problems + +For sampled integral problems, it is possible to cache the weights and reuse +them for multiple data sets. +```@example cache2 +using Integrals + +x = 0.0:0.1:1.0 +y = sin.(x) + +prob = SampledIntegralProblem(y, x) +alg = TrapezoidalRule() + +cache = init(prob, alg) +sol1 = solve!(cache) +``` + +```@example cache2 +cache.y = cos.(x) # use .= to update in-place +sol2 = solve!(cache) +``` +If the grid is modified, the weights are recomputed. +```@example cache2 +cache.x = 0.0:0.2:2.0 +cache.y = sin.(cache.x) +sol3 = solve!(cache) +``` + +For multi-dimensional datasets, the integration dimension can also be changed +```@example cache3 +using Integrals + +x = 0.0:0.1:1.0 +y = sin.(x) .* cos.(x') + +prob = SampledIntegralProblem(y, x) +alg = TrapezoidalRule() + +cache = init(prob, alg) +sol1 = solve!(cache) +``` + +```@example cache3 +cache.dim = 1 +sol2 = solve!(cache) +``` \ No newline at end of file diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index 03d838cf..ef60ded8 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -27,4 +27,45 @@ using Integrals, Test end end end -end \ No newline at end of file +end + +@testset "Caching interface" begin + + x = 0.0:0.1:1.0 + y = sin.(x) + + prob = SampledIntegralProblem(y, x) + alg = TrapezoidalRule() + + cache = init(prob, alg) + sol1 = solve!(cache) + + @test sol1 == solve(prob, alg) + + cache.y = cos.(x) # use .= to update in-place + sol2 = solve!(cache) + + @test sol2 == solve(SampledIntegralProblem(cache.y, cache.x), alg) + + cache.x = 0.0:0.2:2.0 + cache.y = sin.(cache.x) + sol3 = solve!(cache) + + @test sol3 == solve(SampledIntegralProblem(cache.y, cache.x), alg) + + x = 0.0:0.1:1.0 + y = sin.(x) .* cos.(x') + + prob = SampledIntegralProblem(y, x) + alg = TrapezoidalRule() + + cache = init(prob, alg) + sol1 = solve!(cache) + + @test sol1 == solve(prob, alg) + + cache.dim = 1 + sol2 = solve!(cache) + + @test sol2 == solve(SampledIntegralProblem(y, x, dim=1), alg) +end From 813a4601e7bc62314f3fb15827e0923f00bbedcf Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Thu, 21 Sep 2023 18:39:14 +0200 Subject: [PATCH 14/16] Change type parameter in TrapezoidalUniformWeights Co-authored-by: Christopher Rackauckas --- src/trapezoidal.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl index d1365e1c..1ddd78b7 100644 --- a/src/trapezoidal.jl +++ b/src/trapezoidal.jl @@ -1,6 +1,6 @@ -struct TrapezoidalUniformWeights <: UniformWeights +struct TrapezoidalUniformWeights{T} <: UniformWeights n::Int - h::Float64 + h::T end @inline Base.getindex(w::TrapezoidalUniformWeights, i) = ifelse((i == 1) || (i == w.n), w.h*0.5 , w.h) From 7c3d71043b2bf52ebba86f8d3b939d5734063163 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Thu, 21 Sep 2023 18:40:01 +0200 Subject: [PATCH 15/16] Improve type safety Co-authored-by: Christopher Rackauckas --- src/trapezoidal.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl index 1ddd78b7..b51cdd98 100644 --- a/src/trapezoidal.jl +++ b/src/trapezoidal.jl @@ -12,9 +12,9 @@ end @inline function Base.getindex(w::TrapezoidalNonuniformWeights, i) x = w.x - (i == firstindex(x)) && return (x[i + 1] - x[i])*0.5 - (i == lastindex(x)) && return (x[i] - x[i - 1])*0.5 - return (x[i + 1] - x[i - 1])*0.5 + (i == firstindex(x)) && return (x[i + 1] - x[i])/2 + (i == lastindex(x)) && return (x[i] - x[i - 1])/2 + return (x[i + 1] - x[i - 1])/2 end function find_weights(x::AbstractVector, ::TrapezoidalRule) From 65cf88d8c8d8656b6ac2ba760cba4e8ae726e91f Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Thu, 21 Sep 2023 18:40:31 +0200 Subject: [PATCH 16/16] Improve type safety Co-authored-by: Christopher Rackauckas --- src/trapezoidal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl index b51cdd98..ff6777de 100644 --- a/src/trapezoidal.jl +++ b/src/trapezoidal.jl @@ -3,7 +3,7 @@ struct TrapezoidalUniformWeights{T} <: UniformWeights h::T end -@inline Base.getindex(w::TrapezoidalUniformWeights, i) = ifelse((i == 1) || (i == w.n), w.h*0.5 , w.h) +@inline Base.getindex(w::TrapezoidalUniformWeights, i) = ifelse((i == 1) || (i == w.n), w.h/2 , w.h) struct TrapezoidalNonuniformWeights{X<:AbstractArray} <: NonuniformWeights