diff (#95)

* Add diff

* tests pass

* show for adjtrans

* add more tests

* Update test_calculus.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "QuasiArrays"
 uuid = "c4ea9172-b204-11e9-377d-29865faadc5c"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.10"
+version = "0.11"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/QuasiArrays.jl

@@ -88,6 +88,7 @@ include("quasibroadcast.jl")
 include("abstractquasiarraymath.jl")
 include("quasireducedim.jl")
 
+
 include("quasiarray.jl")
 include("quasiarraymath.jl")
 
@@ -104,6 +105,8 @@ include("quasidiagonal.jl")
 include("quasifill.jl")
 include("dense.jl")
 
+include("calculus.jl")
+
 include("quasikron.jl")
 
 promote_leaf_eltypes(x::AbstractQuasiArray{T}) where {T<:Number} = T

src/calculus.jl

@@ -0,0 +1,94 @@
+###
+# sum/cumsum
+###
+
+# support overloading sum by MemoryLayout
+_sum(V::AbstractQuasiArray, dims) = sum_layout(MemoryLayout(V), V, dims)
+_sum(V::AbstractQuasiArray, ::Colon) = sum_layout(MemoryLayout(V), V, :)
+
+_cumsum(A, dims) = cumsum_layout(MemoryLayout(A), A, dims)
+cumsum(A::AbstractQuasiArray; dims::Integer) = _cumsum(A, dims)
+cumsum(x::AbstractQuasiVector) = cumsum(x, dims=1)
+
+# sum is equivalent to hitting by ones(n) on the left or right
+
+cumsum_layout(::QuasiArrayLayout, A, d::Int) = QuasiArray(cumsum(parent(A),dims=d), axes(A))
+
+for Sum in (:sum, :cumsum)
+    Sum_Lay = Symbol(Sum, "_layout")
+    @eval function $Sum_Lay(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, d::Int)
+        a = arguments(LAY, V)
+        if d == 1
+            *($Sum(first(a); dims=1), tail(a)...)
+        else
+            @assert d == 2
+            *(Base.front(a)..., $Sum(last(a); dims=2))
+        end
+    end
+end
+
+function cumsum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims)
+    a = arguments(LAY, V)
+    apply(*, cumsum(a[1]; dims=dims), tail(a)...)
+end
+
+function sum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
+    a = arguments(LAY, V)
+    first(apply(*, sum(a[1]; dims=1), tail(a)...))
+end
+
+sum_layout(::MemoryLayout, A, dims) = sum_size(size(A), A, dims)
+sum_size(::NTuple{N,Integer}, A, dims) where N = _sum(identity, A, dims)
+cumsum_layout(::MemoryLayout, A, dims) = cumsum_size(size(A), A, dims)
+cumsum_size(::NTuple{N,Integer}, A, dims) where N = error("Not implemented")
+
+
+####
+# diff
+####
+
+@inline diff(a::AbstractQuasiArray; dims::Integer=1) = diff_layout(MemoryLayout(a), a, dims)
+function diff_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims...)
+    a = arguments(LAY, V)
+    *(diff(a[1]), tail(a)...)
+end
+
+function diff_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, dims=1)
+    a = arguments(LAY, V)
+    @assert dims == 1 #for type stability, for now
+    # if dims == 1
+        *(diff(a[1]), tail(a)...)
+    # else
+    #     *(front(a)..., diff(a[end]; dims=dims))
+    # end
+end
+
+diff_layout(::MemoryLayout, A, dims...) = diff_size(size(A), A, dims...)
+diff_size(sz, a, dims...) = error("diff not implemented for $(typeof(a))")
+
+diff(x::Inclusion; dims::Integer=1) = ones(eltype(x), diffaxes(x))
+diff(c::AbstractQuasiFill{<:Any,1}; dims::Integer=1) =  zeros(eltype(c), diffaxes(axes(c,1)))
+function diff(c::AbstractQuasiFill{<:Any,2}; dims::Integer=1)
+    a,b = axes(c)
+    if dims == 1
+        zeros(eltype(c), diffaxes(a), b)
+    else
+        zeros(eltype(c), a, diffaxes(b))
+    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
+
+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)
+    D = diff(A.parent; dims=dims)
+    a,b = axes(A)
+    if dims ==  1
+        QuasiMatrix(D ./ diff(a.domain), (diffaxes(a), b))
+    else
+        QuasiMatrix(D ./ permutedims(diff(b.domain)), (a, diffaxes(b)))
+    end
+end

src/quasiadjtrans.jl

@@ -168,8 +168,8 @@ map(f, tvs::QuasiTransposeAbsVec...) = transpose(map((xs...) -> transpose(f(tran
 
 # QuasiAdjoint/QuasiTranspose-vector * vector
 dot(a::AbstractQuasiArray, b::AbstractQuasiArray) = ArrayLayouts.dot(a, b)
-@inline copy(d::Dot{<:Any,<:Any,<:AbstractQuasiArray,<:AbstractQuasiArray}) = _dot(size(d.A,1), d.A, d.B)
-_dot(sz, a, b) = Base.invoke(dot, NTuple{2,Any}, a, b)
+@inline copy(d::Dot{<:Any,<:Any,<:AbstractQuasiArray,<:AbstractQuasiArray}) = dot_size(size(d.A,1), d.A, d.B)
+dot_size(sz, a, b) = Base.invoke(dot, NTuple{2,Any}, a, b)
 
 
 *(u::QuasiAdjointAbsVec, v::AbstractQuasiVector) = dot(u.parent, v)
@@ -247,4 +247,24 @@ arguments(::ApplyLayout{typeof(vcat)}, A::QuasiTranspose) = map(transpose, argum
 arguments(::ApplyLayout{typeof(hcat)}, A::QuasiTranspose) = map(transpose, arguments(ApplyLayout{typeof(vcat)}(), parent(A)))
 
 
-copy(M::Mul{ApplyLayout{typeof(vcat)},QuasiArrayLayout}) = vcat((arguments(vcat, M.A) .* Ref(M.B))...)
+copy(M::Mul{ApplyLayout{typeof(vcat)},QuasiArrayLayout}) = vcat((arguments(vcat, M.A) .* Ref(M.B))...)
+
+
+###
+# show/summary
+###
+
+for shw in (:show, :summary)
+    @eval begin
+        function $shw(io::IO, Ac::QuasiAdjoint)
+            print(io, "adjoint(")
+            $shw(io, parent(Ac))
+            print(io, ")")
+        end
+        function $shw(io::IO, Ac::QuasiTranspose)
+            print(io, "transpose(")
+            $shw(io, parent(Ac))
+            print(io, ")")
+        end
+    end
+end

src/quasireducedim.jl

@@ -262,42 +262,4 @@ end
 # argmax(A::AbstractQuasiArray; dims=:) = findmax(A; dims=dims)[2]
 
 
-# support overloading sum by MemoryLayout
-_sum(V::AbstractQuasiArray, dims) = sum_layout(MemoryLayout(V), V, dims)
-_sum(V::AbstractQuasiArray, ::Colon) = sum_layout(MemoryLayout(V), V, :)
-
-_cumsum(A, dims) = cumsum_layout(MemoryLayout(A), A, dims)
-cumsum(A::AbstractQuasiArray; dims::Integer) = _cumsum(A, dims)
-cumsum(x::AbstractQuasiVector) = cumsum(x, dims=1)
-
-# sum is equivalent to hitting by ones(n) on the left or right
-
-cumsum_layout(::QuasiArrayLayout, A, ::Colon) = QuasiArray(cumsum(parent(A)), axes(A))
-cumsum_layout(::QuasiArrayLayout, A, d::Int) = QuasiArray(cumsum(parent(A),dims=d), axes(A))
-
-for Sum in (:sum, :cumsum)
-    Sum_Lay = Symbol(Sum, "_layout")
-    @eval function $Sum_Lay(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, d::Int)
-        a = arguments(LAY, V)
-        if d == 1
-            *($Sum(first(a); dims=1), tail(a)...)
-        else
-            @assert d == 2
-            *(Base.front(a)..., $Sum(last(a); dims=2))
-        end
-    end
-end
-
-function cumsum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims)
-    a = arguments(LAY, V)
-    apply(*, cumsum(a[1]; dims=dims), tail(a)...)
-end
-
-function sum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
-    a = arguments(LAY, V)
-    first(apply(*, sum(a[1]; dims=1), tail(a)...))
-end
-
-sum_layout(::MemoryLayout, A, dims) = sum_size(size(A), A, dims)
-sum_size(::NTuple{N,Int}, A, dims) where N = _sum(identity, A, dims)
 

test/runtests.jl

@@ -13,6 +13,8 @@ include("test_continuous.jl")
 include("test_matmul.jl")
 include("test_quasifill.jl")
 
+include("test_calculus.jl")
+
 include("test_quasiconcat.jl")
 include("test_quasikron.jl")
 

test/test_calculus.jl

@@ -0,0 +1,65 @@
+using QuasiArrays, IntervalSets, Test
+
+@testset "Calculus" begin
+    @testset "sum" begin
+        A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
+        @test sum(A) ≈ sum(A.parent)
+        @test sum(A; dims=1) ≈ QuasiArray(sum(A.parent; dims=1), (1:1, 1:0.5:2))
+        @test sum(A; dims=2) ≈ QuasiArray(sum(A.parent; dims=2), (0:0.5:0.5, 1:1))
+    end
+
+    @testset "* sum" begin
+        b = QuasiVector(randn(3), 1:0.5:2)
+        A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
+        B = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
+
+        @test sum(ApplyQuasiArray(*, A, b)) ≈ sum(A*b)
+        @test sum(ApplyQuasiArray(*, A, B)) ≈ sum(A*B)
+        @test sum(ApplyQuasiArray(*, A, B); dims=1) ≈ sum(A*B; dims=1)
+        @test sum(ApplyQuasiArray(*, A, B); dims=2) ≈ sum(A*B; dims=2)
+
+        @test sum(b) ≈ last(cumsum(b)) ≈ cumsum(b)[2]
+        @test cumsum(B; dims=1)[2:2,:] ≈ sum(B; dims=1)
+        @test cumsum(B; dims=2)[:,0.5:0.5] ≈ sum(B; dims=2)
+
+        @test cumsum(ApplyQuasiArray(*, A, b)) ≈ cumsum(A*b)
+        @test cumsum(ApplyQuasiArray(*, A, B); dims=1) ≈ cumsum(A*B; dims=1)
+        @test cumsum(ApplyQuasiArray(*, A, B); dims=2) ≈ cumsum(A*B; dims=2)
+    end
+
+    @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(ones(Inclusion(x), Inclusion(x))) == zeros(Inclusion(x[1:end-1]), Inclusion(x))
+        @test diff(ones(Inclusion(x), Inclusion(x)); dims=2) == zeros(Inclusion(x), Inclusion(x[1:end-1]))
+
+        b = QuasiVector(exp.(x), x)
+
+        @test diff(b) ≈ b[Inclusion(x[1:end-1])] atol=1E-2
+
+
+        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; dims=2)[1,:] == diff(A[1,:])
+
+        @testset "* diff" begin
+            b = QuasiVector(randn(3), 1:0.5:2)
+            A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
+            B = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
+    
+            @test @inferred(diff(ApplyQuasiArray(*, A, b))) ≈ diff(A*b)
+            @test @inferred(diff(ApplyQuasiArray(*, A, B))) ≈ diff(A*B)
+        end
+    end
+
+    @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(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))
+    end
+end

test/test_quasiadjtrans.jl

@@ -1,7 +1,7 @@
 # This file is based on a part of Julia. License is MIT: https://julialang.org/license
 
 
-using QuasiArrays, Test, LinearAlgebra
+using QuasiArrays, Test, Base64, LinearAlgebra
 import QuasiArrays: MemoryLayout
 
 
@@ -250,4 +250,12 @@ import QuasiArrays: MemoryLayout
         @test MemoryLayout(typeof(A)) == MemoryLayout(typeof(A')) == MemoryLayout(typeof(transpose(A)))
         @test QuasiArray(A') == QuasiArray(A.parent',(Base.OneTo(2),[0,0.5]))
     end
+
+    @testset "show" begin
+        x = QuasiArray([1,2],[0,0.5])
+        @test stringmime("text/plain", x') == "adjoint(QuasiVector([1, 2], [0.0, 0.5]))"
+        @test stringmime("text/plain", transpose(x)) == "transpose(QuasiVector([1, 2], [0.0, 0.5]))"
+        @test summary(x') == "adjoint(QuasiVector{Int64, Tuple{Vector{Float64}}})"
+        @test summary(transpose(x)) == "transpose(QuasiVector{Int64, Tuple{Vector{Float64}}})"
+    end
 end

test/test_quasireducedim.jl

@@ -2,9 +2,6 @@ using QuasiArrays, Test
 
 @testset "reducedim" begin
     A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
-    @test sum(A) ≈ sum(A.parent)
-    @test sum(A; dims=1) ≈ QuasiArray(sum(A.parent; dims=1), (1:1, 1:0.5:2))
-    @test sum(A; dims=2) ≈ QuasiArray(sum(A.parent; dims=2), (0:0.5:0.5, 1:1))
 
     @test prod(A) ≈ prod(A.parent)
     @test prod(A; dims=1) ≈ QuasiArray(prod(A.parent; dims=1), (1:1, 1:0.5:2))
@@ -21,22 +18,4 @@ using QuasiArrays, Test
     @test minimum(A) ≈ minimum(A.parent)
     @test_broken minimum(A; dims=1) ≈ QuasiArray(minimum(A.parent; dims=1), (1:1, 1:0.5:2))
     @test_broken minimum(A; dims=2) ≈ QuasiArray(minimum(A.parent; dims=2), (0:0.5:0.5, 1:1))
-
-    @testset "* sum" begin
-        b = QuasiVector(randn(3), 1:0.5:2)
-        B = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
-
-        @test sum(ApplyQuasiArray(*, A, b)) ≈ sum(A*b)
-        @test sum(ApplyQuasiArray(*, A, B)) ≈ sum(A*B)
-        @test sum(ApplyQuasiArray(*, A, B); dims=1) ≈ sum(A*B; dims=1)
-        @test sum(ApplyQuasiArray(*, A, B); dims=2) ≈ sum(A*B; dims=2)
-
-        @test sum(b) ≈ last(cumsum(b)) ≈ cumsum(b)[2]
-        @test cumsum(B; dims=1)[2:2,:] ≈ sum(B; dims=1)
-        @test cumsum(B; dims=2)[:,0.5:0.5] ≈ sum(B; dims=2)
-
-        @test cumsum(ApplyQuasiArray(*, A, b)) ≈ cumsum(A*b)
-        @test cumsum(ApplyQuasiArray(*, A, B); dims=1) ≈ cumsum(A*B; dims=1)
-        @test cumsum(ApplyQuasiArray(*, A, B); dims=2) ≈ cumsum(A*B; dims=2)
-    end
 end