Fix sampled evalrule silent wrong results on OffsetArrays

Project.toml

@@ -65,6 +65,7 @@ LinearAlgebra = "1.10"
 MCIntegration = "0.4.2"
 MonteCarloIntegration = "0.2, 0.3"
 Mooncake = "0.5.6"
+OffsetArrays = "1"
 QuadGK = "2.11"
 Random = "1.10"
 Reexport = "1.2"
@@ -91,11 +92,12 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastTanhSinhQuadrature = "b650e0df-f744-4436-b963-b44034668c57"
 HAdaptiveIntegration = "eaa5ad34-b243-48e9-b09c-54bc0655cecf"
 MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 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", "Unitful"]
+test = ["ADTypes", "Arblib", "StaticArrays", "SafeTestsets", "SciMLTesting", "Test", "Distributions", "ChainRulesCore", "FastGaussQuadrature", "FastTanhSinhQuadrature", "Cuba", "Cubature", "MCIntegration", "DifferentiationInterface", "HAdaptiveIntegration", "Unitful", "OffsetArrays"]

src/sampled.jl

@@ -50,7 +50,7 @@ dimension(D::Int) = D
 # Specialized evalrule for 1D arrays (Vector) - returns scalar
 function evalrule(data::AbstractVector, weights, dim)
     # Only use fast indexing path when scalar indexing is efficient (e.g., not on GPU arrays)
-    if ArrayInterface.fast_scalar_indexing(data)
+    if ArrayInterface.fast_scalar_indexing(data) && !Base.has_offset_axes(data)
         n = length(data)
         n == 0 && throw(ArgumentError("No points to integrate"))
         out = data[1] * weights[1]
@@ -67,7 +67,7 @@ end
 # This is the most common case for vector-valued integrands
 function evalrule(data::AbstractMatrix{T}, weights, dim) where {T}
     # Only use fast indexing path when scalar indexing is efficient (e.g., not on GPU arrays)
-    if ArrayInterface.fast_scalar_indexing(data)
+    if ArrayInterface.fast_scalar_indexing(data) && !Base.has_offset_axes(data)
         m, n = size(data)
         n == 0 && throw(ArgumentError("No points to integrate"))
         if dim == 2 || dim == ndims(data)

test/sampled_tests.jl

@@ -1,4 +1,4 @@
-using Integrals, Test, Unitful
+using Integrals, OffsetArrays, Test, Unitful
 @testset "Sampled Integration" begin
     lb = 0.4
     ub = 1.1
@@ -42,6 +42,18 @@ end
     @test sol.u == [2.0, 2.0]u"m"
 end
 
+@testset "Sampled integration with OffsetArrays" begin
+    x = range(0, 1, length = 5)
+    y = OffsetArray(ones(2, 5), 0:1, 1:5)
+    sol = solve(SampledIntegralProblem(y, x; dim = 2), TrapezoidalRule())
+    @test collect(sol.u) == [1.0, 1.0]
+
+    yv = OffsetArray(ones(5), 0:4)
+    xv = OffsetArray(collect(x), 0:4)
+    solv = solve(SampledIntegralProblem(yv, xv), TrapezoidalRule())
+    @test solv.u == 1.0
+end
+
 @testset "Caching interface" begin
     x = 0.0:0.1:1.0
     y = sin.(x)