add tests

JuliaApproximation · Jan 26, 2025 · 97dd159 · 97dd159
1 parent c4d78f8
commit 97dd159
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 13 deletions.
diff --git a/src/QuasiArrays.jl b/src/QuasiArrays.jl
@@ -3,7 +3,7 @@ using Base, LinearAlgebra, LazyArrays, ArrayLayouts, DomainSets, FillArrays, Sta
 import Base: getindex, size, axes, axes1, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
                 Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, copy, copyto!, zero,
                 map, eachindex, eltype, first, last, firstindex, lastindex, in, reshape, permutedims, all,
-                isreal, iszero, isempty, empty, isapprox, fill!, getproperty, showarg
+                isreal, iszero, isone, isempty, empty, isapprox, fill!, getproperty, showarg
 import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinline_meta,
                 DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type,
                 safe_tail, front, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,

diff --git a/src/calculus.jl b/src/calculus.jl
@@ -69,24 +69,29 @@ function diff_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, order.
 end
 
 diff_layout(::MemoryLayout, A, order...; dims...) = diff_size(size(A), A, order...; dims...)
-diff_size(sz, a, order...; dims...) = error("diff not implemented for $(typeof(a))")
+diff_size(sz, a; dims...) = error("diff not implemented for $(typeof(a))")
+function diff_size(sz, a, order::Int; dims...)
+    order < 0 && throw(ArgumentError("order must be non-negative"))
+    order == 0 && return a
+    isone(order) ? diff(a) : diff(diff(a), order-1)
+end
 
 diff(x::Inclusion; dims::Integer=1) = ones(eltype(x), diffaxes(x))
-diff(x::Inclusion, order::Int; dims::Integer=1) = fill(ifelse(isone(order), one(eltype(x)), zero(eltype(x))), diffaxes(x))
-diff(c::AbstractQuasiFill{<:Any,1}, order...; dims::Integer=1) =  zeros(eltype(c), diffaxes(axes(c,1)))
+diff(x::Inclusion, order::Int; dims::Integer=1) = fill(ifelse(isone(order), one(eltype(x)), zero(eltype(x))), diffaxes(x,order))
+diff(c::AbstractQuasiFill{<:Any,1}, order...; dims::Integer=1) =  zeros(eltype(c), diffaxes(axes(c,1),order...))
 function diff(c::AbstractQuasiFill{<:Any,2}, order...; dims::Integer=1)
     a,b = axes(c)
     if dims == 1
-        zeros(eltype(c), diffaxes(a), b)
+        zeros(eltype(c), diffaxes(a, order...), b)
     else
-        zeros(eltype(c), a, diffaxes(b))
+        zeros(eltype(c), a, diffaxes(b, order...))
     end
 end
 
 
-diffaxes(a::Inclusion{<:Any,<:AbstractVector}) = Inclusion(a.domain[1:end-1])
-diffaxes(a::OneTo) = oneto(length(a)-1)
-diffaxes(a) = a # default is differentiation does not change axes
+diffaxes(a::Inclusion{<:Any,<:AbstractVector}, order=1) = Inclusion(a.domain[1:end-order])
+diffaxes(a::OneTo, order=1) = oneto(length(a)-order)
+diffaxes(a, order...) = a # default is differentiation does not change axes
 
 diff(b::QuasiVector; dims::Integer=1) = QuasiVector(diff(b.parent) ./ diff(b.axes[1]), (diffaxes(axes(b,1)),))
 function diff(A::QuasiMatrix; dims::Integer=1)

diff --git a/test/test_calculus.jl b/test/test_calculus.jl
@@ -30,19 +30,25 @@ using QuasiArrays, IntervalSets, Test
 
     @testset "Diff" begin
         x = range(0, 1; length=10_000)
-        @test diff(Inclusion(x)) == ones(Inclusion(x[1:end-1]))
-        @test diff(ones(Inclusion(x))) == zeros(Inclusion(x[1:end-1]))
+        @test diff(Inclusion(x)) == diff(Inclusion(x),1) == ones(Inclusion(x[1:end-1]))
+        @test diff(Inclusion(x),2) == diff(diff(Inclusion(x))) == zeros(Inclusion(x[1:end-2]))
+        @test diff(ones(Inclusion(x))) == diff(ones(Inclusion(x)),1) == zeros(Inclusion(x[1:end-1]))
+        @test diff(ones(Inclusion(x)),2) == diff(diff(ones(Inclusion(x)))) == zeros(Inclusion(x[1:end-2]))
 
         @test diff(ones(Inclusion(x), Inclusion(x))) == zeros(Inclusion(x[1:end-1]), Inclusion(x))
+        @test diff(ones(Inclusion(x), Inclusion(x)), 2) == zeros(Inclusion(x[1:end-2]), Inclusion(x))
         @test diff(ones(Inclusion(x), Inclusion(x)); dims=2) == zeros(Inclusion(x), Inclusion(x[1:end-1]))
+        @test diff(ones(Inclusion(x), Inclusion(x)), 2; dims=2) == zeros(Inclusion(x), Inclusion(x[1:end-2]))
 
         b = QuasiVector(exp.(x), x)
 
         @test diff(b) ≈ b[Inclusion(x[1:end-1])] atol=1E-2
+        @test diff(b,2) ≈ b[Inclusion(x[1:end-2])] atol=1E-1
 
 
         A = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
         @test diff(A; dims=1)[:,0] == diff(A[:,0])
+        @test diff(A,2; dims=1)[:,0] == diff(diff(A[:,0]))
         @test diff(A; dims=2)[1,:] == diff(A[1,:])
 
         @testset "* diff" begin
@@ -57,8 +63,10 @@ using QuasiArrays, IntervalSets, Test
 
     @testset "Interval" begin
         @test diff(Inclusion(0.0..1)) ≡ ones(Inclusion(0.0..1))
-        @test diff(ones(Inclusion(0.0..1))) ≡ zeros(Inclusion(0.0..1))
-        @test diff(ones(Inclusion(0.0..1), Base.OneTo(3))) ≡ zeros(Inclusion(0.0..1), Base.OneTo(3))
+        @test diff(Inclusion(0.0..1),1) ≡ fill(1.0,Inclusion(0.0..1))
+        @test diff(Inclusion(0.0..1),2) ≡ fill(0.0,Inclusion(0.0..1))
+        @test diff(ones(Inclusion(0.0..1))) ≡ diff(ones(Inclusion(0.0..1)),1) ≡ diff(ones(Inclusion(0.0..1)),2) ≡ zeros(Inclusion(0.0..1))
+        @test diff(ones(Inclusion(0.0..1), Base.OneTo(3))) ≡ diff(ones(Inclusion(0.0..1), Base.OneTo(3)),2) ≡ zeros(Inclusion(0.0..1), Base.OneTo(3))
         @test diff(ones(Inclusion(0.0..1), Base.OneTo(3)); dims=2) ≡ zeros(Inclusion(0.0..1), Base.OneTo(2))
         @test diff(ones(Base.OneTo(3), Inclusion(0.0..1))) ≡ zeros(Base.OneTo(2), Inclusion(0.0..1))
         @test diff(ones(Base.OneTo(3), Inclusion(0.0..1)); dims=2) ≡ zeros(Base.OneTo(3), Inclusion(0.0..1))