From e1ccf3fa6c64f5da1148c9364970f63d372461f9 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Tue, 7 Jul 2026 18:41:02 -0400 Subject: [PATCH 1/2] Fix sampled matrix integration over dimensional axes Co-Authored-By: Chris Rackauckas --- Project.toml | 4 +++- src/sampled.jl | 8 ++++++-- test/sampled_tests.jl | 12 +++++++++++- 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index bf5b4d26..d0e713d6 100644 --- a/Project.toml +++ b/Project.toml @@ -71,6 +71,7 @@ SciMLLogging = "1.10.1, 2" SciMLTesting = "1" StaticArrays = "1.9.8" Test = "1.10" +Unitful = "1" Zygote = "0.7.10" julia = "1.10" @@ -90,6 +91,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["ADTypes", "Arblib", "StaticArrays", "SafeTestsets", "SciMLTesting", "Test", "Distributions", "ChainRulesCore", "FastGaussQuadrature", "FastTanhSinhQuadrature", "Cuba", "Cubature", "MCIntegration", "DifferentiationInterface", "HAdaptiveIntegration"] +test = ["ADTypes", "Arblib", "StaticArrays", "SafeTestsets", "SciMLTesting", "Test", "Distributions", "ChainRulesCore", "FastGaussQuadrature", "FastTanhSinhQuadrature", "Cuba", "Cubature", "MCIntegration", "DifferentiationInterface", "HAdaptiveIntegration", "Unitful"] diff --git a/src/sampled.jl b/src/sampled.jl index d56e8bd9..8641a91f 100644 --- a/src/sampled.jl +++ b/src/sampled.jl @@ -72,8 +72,12 @@ function evalrule(data::AbstractMatrix{T}, weights, dim) where {T} n == 0 && throw(ArgumentError("No points to integrate")) if dim == 2 || dim == ndims(data) # Integration over columns (last dimension) - the common case - out = zeros(T, m) - @inbounds for i in 1:n + w1 = weights[1] + out = Vector{Base.promote_op(*, typeof(w1), T)}(undef, m) + @inbounds @simd for j in 1:m + out[j] = w1 * data[j, 1] + end + @inbounds for i in 2:n w = weights[i] @simd for j in 1:m out[j] += w * data[j, i] diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index 9c7398b8..1f8e26ee 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -1,4 +1,4 @@ -using Integrals, Test +using Integrals, Test, Unitful @testset "Sampled Integration" begin lb = 0.4 ub = 1.1 @@ -32,6 +32,16 @@ using Integrals, Test end end +@testset "Sampled integration with dimensional axes" begin + data = ones(2, 3) + axis = [0.0, 1.0, 2.0]u"m" + + prob = SampledIntegralProblem(data, axis; dim = 2) + sol = solve(prob, TrapezoidalRule()) + + @test sol.u == [2.0, 2.0]u"m" +end + @testset "Caching interface" begin x = 0.0:0.1:1.0 y = sin.(x) From e4198cfd4687b1ccf51317e55eff4a9e10aea63f Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Tue, 7 Jul 2026 21:00:26 -0400 Subject: [PATCH 2/2] Avoid non-public promote_op in sampled integration Co-Authored-By: Chris Rackauckas --- src/sampled.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/sampled.jl b/src/sampled.jl index 8641a91f..b2e85c70 100644 --- a/src/sampled.jl +++ b/src/sampled.jl @@ -73,7 +73,8 @@ function evalrule(data::AbstractMatrix{T}, weights, dim) where {T} if dim == 2 || dim == ndims(data) # Integration over columns (last dimension) - the common case w1 = weights[1] - out = Vector{Base.promote_op(*, typeof(w1), T)}(undef, m) + out_type = m == 0 ? typeof(w1 * zero(T)) : typeof(w1 * data[1, 1]) + out = Vector{out_type}(undef, m) @inbounds @simd for j in 1:m out[j] = w1 * data[j, 1] end