Support ExpansionLayout in ContinuumArrays.jl (#80)

* Support ExpansionLayout in ContinuumArrays.jl

* Update quasifill.jl

* InclusionLayout, simplify layout_broadcasted

* PolynomialLayout

* Increase coverage

* fix tests

Author: dlfivefifty

src/QuasiArrays.jl

@@ -34,7 +34,7 @@ import Base: ones, zeros, one, zero, fill
 
 import Base.Broadcast: materialize, materialize!, BroadcastStyle, AbstractArrayStyle, Style, broadcasted, Broadcasted, Unknown,
                         newindex, broadcastable, preprocess, _eachindex, _broadcast_getindex, broadcast_shape,
-                        DefaultArrayStyle, axistype, throwdm, instantiate, combine_eltypes, eltypes
+                        DefaultArrayStyle, axistype, throwdm, instantiate, combine_eltypes, eltypes, combine_styles
 
 import LinearAlgebra: transpose, adjoint, checkeltype_adjoint, checkeltype_transpose, Diagonal,
                         AbstractTriangular, pinv, inv, promote_leaf_eltypes, power_by_squaring,

src/abstractquasiarray.jl

@@ -504,7 +504,7 @@ function isequal(A::AbstractQuasiArray, B::AbstractQuasiArray)
     return true
 end
 
-function (==)(A::AbstractQuasiArray, B::AbstractQuasiArray)
+function _equals(_, _, A, B)
     if axes(A) != axes(B)
         return false
     end
@@ -521,8 +521,13 @@ function (==)(A::AbstractQuasiArray, B::AbstractQuasiArray)
 end
 
 
+(==)(A::AbstractQuasiArray, B::AbstractQuasiArray) = _equals(MemoryLayout(A), MemoryLayout(B), A, B)
+
 ##
 # show
 ##
 
-show(io::IO, A::AbstractQuasiArray) = summary(io, A)
+show(io::IO, A::AbstractQuasiArray) = summary(io, A)
+
+struct QuasiArrayLayout <: MemoryLayout end
+MemoryLayout(::Type{<:AbstractQuasiArray}) = QuasiArrayLayout()

src/indices.jl

@@ -142,6 +142,10 @@ Inclusion(domain) = Inclusion{eltype(domain)}(domain)
 Inclusion(S::Inclusion) = S
 Inclusion(S::Slice) = Inclusion(S.indices)
 
+struct PolynomialLayout <: MemoryLayout end
+
+MemoryLayout(::Type{<:Inclusion}) = PolynomialLayout()
+
 convert(::Type{Inclusion}, d::Inclusion) = d
 convert(::Type{Inclusion{T}}, d::Inclusion) where T = Inclusion{T}(d)
 convert(::Type{AbstractVector}, d::Inclusion{<:Any,<:AbstractVector}) =

src/lazyquasiarrays.jl

@@ -102,6 +102,10 @@ subbroadcaststyle(::LazyQuasiArrayStyle{2}, ::Type{<:Tuple{Number,Number}}) = De
 subbroadcaststyle(::LazyQuasiArrayStyle{2}, ::Type{<:Tuple{Number,JR}}) where JR = Base.BroadcastStyle(JR)
 subbroadcaststyle(::LazyQuasiArrayStyle{2}, ::Type{<:Tuple{KR,Number}}) where KR = Base.BroadcastStyle(KR)
 
+layout_broadcasted(_, f, args...) = Broadcasted{typeof(combine_styles(args...))}(f, args)
+
+broadcasted(S::LazyQuasiArrayStyle, f, args...) = layout_broadcasted(map(MemoryLayout,args), f, args...)
+
 
 struct BroadcastQuasiArray{T, N, F, Args} <: LazyQuasiArray{T, N}
     f::F
@@ -137,7 +141,7 @@ broadcasted(A::BroadcastQuasiArray) = instantiate(broadcasted(A.f, A.args...))
 axes(A::BroadcastQuasiArray) = axes(broadcasted(A))
 size(A::BroadcastQuasiArray) = map(length, axes(A))
 
-function ==(A::BroadcastQuasiArray, B::BroadcastQuasiArray)
+function _equals(::BroadcastLayout, ::BroadcastLayout, A, B)
     A.f == B.f && all(A.args .== B.args) && return true
     error("Not implemented")
 end
@@ -170,6 +174,9 @@ call(b::BroadcastLayout, a::SubQuasiArray) = call(b, parent(a))
 summary(io::IO, A::BroadcastQuasiArray) = _broadcastarray_summary(io, A)
 summary(io::IO, A::ApplyQuasiArray) = _applyarray_summary(io, A)
 
+_mul_summary(_, io, A) = _applyarray_summary(io, A)
+summary(io::IO, A::ApplyQuasiArray{<:Any,N,typeof(*)}) where N = _mul_summary(MemoryLayout(A), io, A)
+
 for op in (:+, :-, :*, :\, :/)
     @eval begin
         function summary(io::IO, A::BroadcastQuasiArray{<:Any,N,typeof($op)}) where N

src/matmul.jl

@@ -86,9 +86,6 @@ function copyto!(dest::MulQuasiArray, src::MulQuasiArray)
     dest
 end
 
-
-struct QuasiArrayLayout <: MemoryLayout end
-MemoryLayout(::Type{<:AbstractQuasiArray}) = QuasiArrayLayout()
 Array(M::Mul{<:Any,<:Any,<:Any,<:AbstractQuasiMatrix}) = eltype(M)[M[k,j] for k in axes(M)[1], j in axes(M)[2]]
 Array(M::Mul{<:Any,<:Any,<:Any,<:AbstractQuasiVector}) = eltype(M)[M[k] for k in axes(M)[1]]
 _quasi_mul(M, _) = QuasiArray(M)
@@ -132,8 +129,12 @@ _rdiv_scal_reduce(x::Number, Z) = (Z / x,)
 _rdiv_scal_reduce(x::Number, Z::AbstractArray, Y...) = (Y..., Z/x)
 _rdiv_scal_reduce(x::Number, Z, Y...) = (_rdiv_scal_reduce(x, Y...)..., Z)
 
-*(x::Number, A::MulQuasiArray) = ApplyQuasiArray(*, _lmul_scal_reduce(x, arguments(A)...)...)
-*(A::MulQuasiArray, x::Number) = ApplyQuasiArray(*, _rmul_scal_reduce(x, reverse(arguments(A))...)...)
+broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(*), x::Number, A::MulQuasiArray{<:Any,N}) where N = 
+    ApplyQuasiArray(*, _lmul_scal_reduce(x, arguments(A)...)...)
+broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(*), A::MulQuasiArray{<:Any,N}, x::Number) where N = 
+    ApplyQuasiArray(*, _rmul_scal_reduce(x, reverse(arguments(A))...)...)
 
-\(x::Number, A::MulQuasiArray) = ApplyQuasiArray(*, _ldiv_scal_reduce(x, arguments(A)...)...)
-/(A::MulQuasiArray, x::Number) = ApplyQuasiArray(*, _rdiv_scal_reduce(x, reverse(arguments(A))...)...)
+broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(\), x::Number, A::MulQuasiArray{<:Any,N}) where N = 
+    ApplyQuasiArray(*, _ldiv_scal_reduce(x, arguments(A)...)...)
+broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(/), A::MulQuasiArray{<:Any,N}, x::Number) where N = 
+    ApplyQuasiArray(*, _ldiv_scal_reduce(x, arguments(A)...)...)

src/quasifill.jl

@@ -98,7 +98,6 @@ copy(F::QuasiFill) = QuasiFill(F.value, F.axes)
 convert(::Type{T}, F::T) where T<:QuasiFill = F
 
 
-
 getindex(F::QuasiFill{<:Any,0}) = getindex_value(F)
 
 
@@ -191,6 +190,40 @@ for op in (:+, :-, :*, :/, :\)
 end
 
 
+########
+# show
+#######
+
+function summary(io::IO, F::QuasiOnes)
+    print(io, "ones(")
+    summary(io, F.axes[1])
+    for a in tail(F.axes)
+        print(io, ", ")
+        summary(io, a)
+    end
+    print(io, ")")
+end
+
+function summary(io::IO, F::QuasiZeros)
+    print(io, "zeros(")
+    summary(io, F.axes[1])
+    for a in tail(F.axes)
+        print(io, ", ")
+        summary(io, a)
+    end
+    print(io, ")")
+end
+
+function summary(io::IO, F::QuasiFill)
+    print(io, "fill($(F.value), ")
+    summary(io, F.axes[1])
+    for a in tail(F.axes)
+        print(io, ", ")
+        summary(io, a)
+    end
+    print(io, ")")
+end
+
 #########
 #  Special matrix types
 #########

src/quasireducedim.jl

@@ -266,19 +266,37 @@ end
 _sum(V::AbstractQuasiArray, dims) = __sum(MemoryLayout(V), V, dims)
 _sum(V::AbstractQuasiArray, ::Colon) = __sum(MemoryLayout(V), V, :)
 
-# sum is equivalent to hitting by ones(n) on the left or rifght
-function __sum(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
-        *(most(a)..., sum(last(a); dims=2))
+_cumsum(A, dims) = __cumsum(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(::QuasiArrayLayout, A, ::Colon) = QuasiArray(cumsum(parent(A)), axes(A))
+__cumsum(::QuasiArrayLayout, A, d::Int) = QuasiArray(cumsum(parent(A),dims=d), axes(A))
+
+for Sum in (:sum, :cumsum)
+    __Sum = Symbol("__", Sum)
+    @eval function $__Sum(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
+            *(most(a)..., $Sum(last(a); dims=2))
+        end
     end
 end
+
+function __cumsum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims)
+    a = arguments(LAY, V)
+    apply(*, cumsum(a[1]; dims=dims), tail(a)...)
+end
+
 function __sum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
     a = arguments(LAY, V)
     first(apply(*, sum(a[1]; dims=1), tail(a)...))
 end
 
-__sum(_, A, dims) = _sum(identity, A, dims)
+__sum(_, A, dims) = _sum(identity, A, dims)
+

test/test_matmul.jl

@@ -7,7 +7,6 @@ import QuasiArrays: apply
         @test rowsupport(A,0.1) == colsupport(A,0.1) == 0.1
         b = Inclusion(0:0.1:1)
         Ab = A*b
-        @test Ab isa QuasiArray
         @test Ab[0.1] ≈ 0.1^2
         @test_throws DimensionMismatch A*Inclusion(1:2)
 

test/test_quasifill.jl

@@ -704,4 +704,10 @@ import QuasiArrays: AbstractQuasiFill
         @test_throws DimensionMismatch I+B
         @test B*I ≈ I*B ≈ B
     end
+
+    @testset "show" begin
+        @test stringmime("text/plain",ones(Inclusion([1,2,3]))) == "ones(Inclusion([1, 2, 3]))"
+        @test stringmime("text/plain",zeros(Inclusion([1,2,3]))) == "zeros(Inclusion([1, 2, 3]))"
+        @test stringmime("text/plain",fill(2,Inclusion([1,2,3]))) == "fill(2, Inclusion([1, 2, 3]))"
+    end
 end

test/test_quasilazy.jl

@@ -1,4 +1,4 @@
-using QuasiArrays, LazyArrays, ArrayLayouts, Test
+using QuasiArrays, LazyArrays, ArrayLayouts, Base64, Test
 import QuasiArrays: QuasiLazyLayout, QuasiArrayApplyStyle, LazyQuasiMatrix, LazyQuasiArrayStyle
 import LazyArrays: MulStyle, ApplyStyle
 
@@ -78,6 +78,11 @@ Base.getindex(A::MyQuasiLazyMatrix, x::Float64, y::Float64) = A.A[x,y]
 
             @test BroadcastQuasiArray(*, x, ApplyQuasiArray(^, A, 2)) * y ≈ (x .* A^2) * y
         end
+
+        @testset "summary" begin
+            A = ApplyQuasiArray(*, ones(Inclusion([1,2,3]), Inclusion([4,5])), fill(2,Inclusion([4,5])))
+            @test stringmime("text/plain", A) == "(ones(Inclusion([1, 2, 3]), Inclusion([4, 5]))) * (fill(2, Inclusion([4, 5])))"
+        end
     end
     @testset "\\" begin
         A = QuasiArray(rand(3,3),(0:0.5:1,0:0.5:1))