Skip to content

Commit 706e116

Browse files
dlfivefiftydlfivefifty
authored
diff (#95)
* Add diff * tests pass * show for adjtrans * add more tests * Update test_calculus.jl
1 parent 631e150 commit 706e116

File tree

9 files changed

+197
-64
lines changed

9 files changed

+197
-64
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "QuasiArrays"
22
uuid = "c4ea9172-b204-11e9-377d-29865faadc5c"
33
authors = ["Sheehan Olver <solver@mac.com>"]
4-
version = "0.10"
4+
version = "0.11"
55

66
[deps]
77
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

‎src/QuasiArrays.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,7 @@ include("quasibroadcast.jl")
8888
include("abstractquasiarraymath.jl")
8989
include("quasireducedim.jl")
9090

91+
9192
include("quasiarray.jl")
9293
include("quasiarraymath.jl")
9394

@@ -104,6 +105,8 @@ include("quasidiagonal.jl")
104105
include("quasifill.jl")
105106
include("dense.jl")
106107

108+
include("calculus.jl")
109+
107110
include("quasikron.jl")
108111

109112
promote_leaf_eltypes(x::AbstractQuasiArray{T}) where {T<:Number} = T

‎src/calculus.jl‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
###
2+
# sum/cumsum
3+
###
4+
5+
# support overloading sum by MemoryLayout
6+
_sum(V::AbstractQuasiArray, dims) = sum_layout(MemoryLayout(V), V, dims)
7+
_sum(V::AbstractQuasiArray, ::Colon) = sum_layout(MemoryLayout(V), V, :)
8+
9+
_cumsum(A, dims) = cumsum_layout(MemoryLayout(A), A, dims)
10+
cumsum(A::AbstractQuasiArray; dims::Integer) = _cumsum(A, dims)
11+
cumsum(x::AbstractQuasiVector) = cumsum(x, dims=1)
12+
13+
# sum is equivalent to hitting by ones(n) on the left or right
14+
15+
cumsum_layout(::QuasiArrayLayout, A, d::Int) = QuasiArray(cumsum(parent(A),dims=d), axes(A))
16+
17+
for Sum in (:sum, :cumsum)
18+
Sum_Lay = Symbol(Sum, "_layout")
19+
@eval function $Sum_Lay(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, d::Int)
20+
a = arguments(LAY, V)
21+
if d == 1
22+
*($Sum(first(a); dims=1), tail(a)...)
23+
else
24+
@assert d == 2
25+
*(Base.front(a)..., $Sum(last(a); dims=2))
26+
end
27+
end
28+
end
29+
30+
function cumsum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims)
31+
a = arguments(LAY, V)
32+
apply(*, cumsum(a[1]; dims=dims), tail(a)...)
33+
end
34+
35+
function sum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
36+
a = arguments(LAY, V)
37+
first(apply(*, sum(a[1]; dims=1), tail(a)...))
38+
end
39+
40+
sum_layout(::MemoryLayout, A, dims) = sum_size(size(A), A, dims)
41+
sum_size(::NTuple{N,Integer}, A, dims) where N = _sum(identity, A, dims)
42+
cumsum_layout(::MemoryLayout, A, dims) = cumsum_size(size(A), A, dims)
43+
cumsum_size(::NTuple{N,Integer}, A, dims) where N = error("Not implemented")
44+
45+
46+
####
47+
# diff
48+
####
49+
50+
@inline diff(a::AbstractQuasiArray; dims::Integer=1) = diff_layout(MemoryLayout(a), a, dims)
51+
function diff_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims...)
52+
a = arguments(LAY, V)
53+
*(diff(a[1]), tail(a)...)
54+
end
55+
56+
function diff_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, dims=1)
57+
a = arguments(LAY, V)
58+
@assert dims == 1 #for type stability, for now
59+
# if dims == 1
60+
*(diff(a[1]), tail(a)...)
61+
# else
62+
# *(front(a)..., diff(a[end]; dims=dims))
63+
# end
64+
end
65+
66+
diff_layout(::MemoryLayout, A, dims...) = diff_size(size(A), A, dims...)
67+
diff_size(sz, a, dims...) = error("diff not implemented for $(typeof(a))")
68+
69+
diff(x::Inclusion; dims::Integer=1) = ones(eltype(x), diffaxes(x))
70+
diff(c::AbstractQuasiFill{<:Any,1}; dims::Integer=1) = zeros(eltype(c), diffaxes(axes(c,1)))
71+
function diff(c::AbstractQuasiFill{<:Any,2}; dims::Integer=1)
72+
a,b = axes(c)
73+
if dims == 1
74+
zeros(eltype(c), diffaxes(a), b)
75+
else
76+
zeros(eltype(c), a, diffaxes(b))
77+
end
78+
end
79+
80+
81+
diffaxes(a::Inclusion{<:Any,<:AbstractVector}) = Inclusion(a.domain[1:end-1])
82+
diffaxes(a::OneTo) = oneto(length(a)-1)
83+
diffaxes(a) = a # default is differentiation does not change axes
84+
85+
diff(b::QuasiVector; dims::Integer=1) = QuasiVector(diff(b.parent) ./ diff(b.axes[1]), (diffaxes(axes(b,1)),))
86+
function diff(A::QuasiMatrix; dims::Integer=1)
87+
D = diff(A.parent; dims=dims)
88+
a,b = axes(A)
89+
if dims == 1
90+
QuasiMatrix(D ./ diff(a.domain), (diffaxes(a), b))
91+
else
92+
QuasiMatrix(D ./ permutedims(diff(b.domain)), (a, diffaxes(b)))
93+
end
94+
end

‎src/quasiadjtrans.jl‎

Lines changed: 23 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -168,8 +168,8 @@ map(f, tvs::QuasiTransposeAbsVec...) = transpose(map((xs...) -> transpose(f(tran
168168

169169
# QuasiAdjoint/QuasiTranspose-vector * vector
170170
dot(a::AbstractQuasiArray, b::AbstractQuasiArray) = ArrayLayouts.dot(a, b)
171-
@inline copy(d::Dot{<:Any,<:Any,<:AbstractQuasiArray,<:AbstractQuasiArray}) = _dot(size(d.A,1), d.A, d.B)
172-
_dot(sz, a, b) = Base.invoke(dot, NTuple{2,Any}, a, b)
171+
@inline copy(d::Dot{<:Any,<:Any,<:AbstractQuasiArray,<:AbstractQuasiArray}) = dot_size(size(d.A,1), d.A, d.B)
172+
dot_size(sz, a, b) = Base.invoke(dot, NTuple{2,Any}, a, b)
173173

174174

175175
*(u::QuasiAdjointAbsVec, v::AbstractQuasiVector) = dot(u.parent, v)
@@ -247,4 +247,24 @@ arguments(::ApplyLayout{typeof(vcat)}, A::QuasiTranspose) = map(transpose, argum
247247
arguments(::ApplyLayout{typeof(hcat)}, A::QuasiTranspose) = map(transpose, arguments(ApplyLayout{typeof(vcat)}(), parent(A)))
248248

249249

250-
copy(M::Mul{ApplyLayout{typeof(vcat)},QuasiArrayLayout}) = vcat((arguments(vcat, M.A) .* Ref(M.B))...)
250+
copy(M::Mul{ApplyLayout{typeof(vcat)},QuasiArrayLayout}) = vcat((arguments(vcat, M.A) .* Ref(M.B))...)
251+
252+
253+
###
254+
# show/summary
255+
###
256+
257+
for shw in (:show, :summary)
258+
@eval begin
259+
function $shw(io::IO, Ac::QuasiAdjoint)
260+
print(io, "adjoint(")
261+
$shw(io, parent(Ac))
262+
print(io, ")")
263+
end
264+
function $shw(io::IO, Ac::QuasiTranspose)
265+
print(io, "transpose(")
266+
$shw(io, parent(Ac))
267+
print(io, ")")
268+
end
269+
end
270+
end

‎src/quasireducedim.jl‎

Lines changed: 0 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -262,42 +262,4 @@ end
262262
# argmax(A::AbstractQuasiArray; dims=:) = findmax(A; dims=dims)[2]
263263

264264

265-
# support overloading sum by MemoryLayout
266-
_sum(V::AbstractQuasiArray, dims) = sum_layout(MemoryLayout(V), V, dims)
267-
_sum(V::AbstractQuasiArray, ::Colon) = sum_layout(MemoryLayout(V), V, :)
268-
269-
_cumsum(A, dims) = cumsum_layout(MemoryLayout(A), A, dims)
270-
cumsum(A::AbstractQuasiArray; dims::Integer) = _cumsum(A, dims)
271-
cumsum(x::AbstractQuasiVector) = cumsum(x, dims=1)
272-
273-
# sum is equivalent to hitting by ones(n) on the left or right
274-
275-
cumsum_layout(::QuasiArrayLayout, A, ::Colon) = QuasiArray(cumsum(parent(A)), axes(A))
276-
cumsum_layout(::QuasiArrayLayout, A, d::Int) = QuasiArray(cumsum(parent(A),dims=d), axes(A))
277-
278-
for Sum in (:sum, :cumsum)
279-
Sum_Lay = Symbol(Sum, "_layout")
280-
@eval function $Sum_Lay(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, d::Int)
281-
a = arguments(LAY, V)
282-
if d == 1
283-
*($Sum(first(a); dims=1), tail(a)...)
284-
else
285-
@assert d == 2
286-
*(Base.front(a)..., $Sum(last(a); dims=2))
287-
end
288-
end
289-
end
290-
291-
function cumsum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims)
292-
a = arguments(LAY, V)
293-
apply(*, cumsum(a[1]; dims=dims), tail(a)...)
294-
end
295-
296-
function sum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
297-
a = arguments(LAY, V)
298-
first(apply(*, sum(a[1]; dims=1), tail(a)...))
299-
end
300-
301-
sum_layout(::MemoryLayout, A, dims) = sum_size(size(A), A, dims)
302-
sum_size(::NTuple{N,Int}, A, dims) where N = _sum(identity, A, dims)
303265

‎test/runtests.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,8 @@ include("test_continuous.jl")
1313
include("test_matmul.jl")
1414
include("test_quasifill.jl")
1515

16+
include("test_calculus.jl")
17+
1618
include("test_quasiconcat.jl")
1719
include("test_quasikron.jl")
1820

‎test/test_calculus.jl‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
using QuasiArrays, IntervalSets, Test
2+
3+
@testset "Calculus" begin
4+
@testset "sum" begin
5+
A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
6+
@test sum(A) ≈ sum(A.parent)
7+
@test sum(A; dims=1) ≈ QuasiArray(sum(A.parent; dims=1), (1:1, 1:0.5:2))
8+
@test sum(A; dims=2) ≈ QuasiArray(sum(A.parent; dims=2), (0:0.5:0.5, 1:1))
9+
end
10+
11+
@testset "* sum" begin
12+
b = QuasiVector(randn(3), 1:0.5:2)
13+
A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
14+
B = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
15+
16+
@test sum(ApplyQuasiArray(*, A, b)) ≈ sum(A*b)
17+
@test sum(ApplyQuasiArray(*, A, B)) ≈ sum(A*B)
18+
@test sum(ApplyQuasiArray(*, A, B); dims=1) ≈ sum(A*B; dims=1)
19+
@test sum(ApplyQuasiArray(*, A, B); dims=2) ≈ sum(A*B; dims=2)
20+
21+
@test sum(b) ≈ last(cumsum(b)) ≈ cumsum(b)[2]
22+
@test cumsum(B; dims=1)[2:2,:] ≈ sum(B; dims=1)
23+
@test cumsum(B; dims=2)[:,0.5:0.5] ≈ sum(B; dims=2)
24+
25+
@test cumsum(ApplyQuasiArray(*, A, b)) ≈ cumsum(A*b)
26+
@test cumsum(ApplyQuasiArray(*, A, B); dims=1) ≈ cumsum(A*B; dims=1)
27+
@test cumsum(ApplyQuasiArray(*, A, B); dims=2) ≈ cumsum(A*B; dims=2)
28+
end
29+
30+
@testset "Diff" begin
31+
x = range(0, 1; length=10_000)
32+
@test diff(Inclusion(x)) == ones(Inclusion(x[1:end-1]))
33+
@test diff(ones(Inclusion(x))) == zeros(Inclusion(x[1:end-1]))
34+
35+
@test diff(ones(Inclusion(x), Inclusion(x))) == zeros(Inclusion(x[1:end-1]), Inclusion(x))
36+
@test diff(ones(Inclusion(x), Inclusion(x)); dims=2) == zeros(Inclusion(x), Inclusion(x[1:end-1]))
37+
38+
b = QuasiVector(exp.(x), x)
39+
40+
@test diff(b) ≈ b[Inclusion(x[1:end-1])] atol=1E-2
41+
42+
43+
A = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
44+
@test diff(A; dims=1)[:,0] == diff(A[:,0])
45+
@test diff(A; dims=2)[1,:] == diff(A[1,:])
46+
47+
@testset "* diff" begin
48+
b = QuasiVector(randn(3), 1:0.5:2)
49+
A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
50+
B = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
51+
52+
@test @inferred(diff(ApplyQuasiArray(*, A, b))) ≈ diff(A*b)
53+
@test @inferred(diff(ApplyQuasiArray(*, A, B))) ≈ diff(A*B)
54+
end
55+
end
56+
57+
@testset "Interval" begin
58+
@test diff(Inclusion(0.0..1)) ≡ ones(Inclusion(0.0..1))
59+
@test diff(ones(Inclusion(0.0..1))) ≡ zeros(Inclusion(0.0..1))
60+
@test diff(ones(Inclusion(0.0..1), Base.OneTo(3))) ≡ zeros(Inclusion(0.0..1), Base.OneTo(3))
61+
@test diff(ones(Inclusion(0.0..1), Base.OneTo(3)); dims=2) ≡ zeros(Inclusion(0.0..1), Base.OneTo(2))
62+
@test diff(ones(Base.OneTo(3), Inclusion(0.0..1))) ≡ zeros(Base.OneTo(2), Inclusion(0.0..1))
63+
@test diff(ones(Base.OneTo(3), Inclusion(0.0..1)); dims=2) ≡ zeros(Base.OneTo(3), Inclusion(0.0..1))
64+
end
65+
end

‎test/test_quasiadjtrans.jl‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
# This file is based on a part of Julia. License is MIT: https://julialang.org/license
22

33

4-
using QuasiArrays, Test, LinearAlgebra
4+
using QuasiArrays, Test, Base64, LinearAlgebra
55
import QuasiArrays: MemoryLayout
66

77

@@ -250,4 +250,12 @@ import QuasiArrays: MemoryLayout
250250
@test MemoryLayout(typeof(A)) == MemoryLayout(typeof(A')) == MemoryLayout(typeof(transpose(A)))
251251
@test QuasiArray(A') == QuasiArray(A.parent',(Base.OneTo(2),[0,0.5]))
252252
end
253+
254+
@testset "show" begin
255+
x = QuasiArray([1,2],[0,0.5])
256+
@test stringmime("text/plain", x') == "adjoint(QuasiVector([1, 2], [0.0, 0.5]))"
257+
@test stringmime("text/plain", transpose(x)) == "transpose(QuasiVector([1, 2], [0.0, 0.5]))"
258+
@test summary(x') == "adjoint(QuasiVector{Int64, Tuple{Vector{Float64}}})"
259+
@test summary(transpose(x)) == "transpose(QuasiVector{Int64, Tuple{Vector{Float64}}})"
260+
end
253261
end

‎test/test_quasireducedim.jl‎

Lines changed: 0 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,6 @@ using QuasiArrays, Test
22

33
@testset "reducedim" begin
44
A = QuasiArray(randn(2,3), (0:0.5:0.5, 1:0.5:2))
5-
@test sum(A) ≈ sum(A.parent)
6-
@test sum(A; dims=1) ≈ QuasiArray(sum(A.parent; dims=1), (1:1, 1:0.5:2))
7-
@test sum(A; dims=2) ≈ QuasiArray(sum(A.parent; dims=2), (0:0.5:0.5, 1:1))
85

96
@test prod(A) ≈ prod(A.parent)
107
@test prod(A; dims=1) ≈ QuasiArray(prod(A.parent; dims=1), (1:1, 1:0.5:2))
@@ -21,22 +18,4 @@ using QuasiArrays, Test
2118
@test minimum(A) ≈ minimum(A.parent)
2219
@test_broken minimum(A; dims=1) ≈ QuasiArray(minimum(A.parent; dims=1), (1:1, 1:0.5:2))
2320
@test_broken minimum(A; dims=2) ≈ QuasiArray(minimum(A.parent; dims=2), (0:0.5:0.5, 1:1))
24-
25-
@testset "* sum" begin
26-
b = QuasiVector(randn(3), 1:0.5:2)
27-
B = QuasiArray(randn(3,2), (1:0.5:2,0:0.5:0.5))
28-
29-
@test sum(ApplyQuasiArray(*, A, b)) ≈ sum(A*b)
30-
@test sum(ApplyQuasiArray(*, A, B)) ≈ sum(A*B)
31-
@test sum(ApplyQuasiArray(*, A, B); dims=1) ≈ sum(A*B; dims=1)
32-
@test sum(ApplyQuasiArray(*, A, B); dims=2) ≈ sum(A*B; dims=2)
33-
34-
@test sum(b) ≈ last(cumsum(b)) ≈ cumsum(b)[2]
35-
@test cumsum(B; dims=1)[2:2,:] ≈ sum(B; dims=1)
36-
@test cumsum(B; dims=2)[:,0.5:0.5] ≈ sum(B; dims=2)
37-
38-
@test cumsum(ApplyQuasiArray(*, A, b)) ≈ cumsum(A*b)
39-
@test cumsum(ApplyQuasiArray(*, A, B); dims=1) ≈ cumsum(A*B; dims=1)
40-
@test cumsum(ApplyQuasiArray(*, A, B); dims=2) ≈ cumsum(A*B; dims=2)
41-
end
4221
end

0 commit comments

Comments
 (0)