Work on QuasiKron (#84)

* Work on QuasiKron

* v0.9.3

* Update ci.yml

Author: dlfivefifty

.github/workflows/ci.yml

@@ -12,7 +12,6 @@ jobs:
         version:
           - '1.6'
           - '1'
-          - '^1.8.0-0'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

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

src/quasikron.jl

@@ -18,14 +18,39 @@ QuasiKron{T}(a::AbstractQuasiOrMatrix...) where T = QuasiKron{mapreduce(eltype,p
 QuasiKron(a...) = QuasiKron{mapreduce(eltype,promote_type,a)}(a...)
 
 
-quasikron(a::Inclusion...) = Inclusion(ProductDomain(map(domain,a)...))
+quasikron(a::Union{Inclusion,OneTo}...) = Inclusion(ProductDomain(map(domain,a)...))
 quasikron(a...) = QuasiKron(a...)
 
 _kronaxes(::Tuple{}...) = ()
 _kronaxes(a::Tuple...) = (quasikron(map(first,a)...), _kronaxes(map(tail,a)...)...)
 axes(K::QuasiKron) = _kronaxes(map(axes,K.args)...)
 
-function _getindex(::Type{IND}, K::QuasiKron{<:Any,1}, i::IND) where IND
+function _getindex(::Type{IND}, K::QuasiKron, i::IND) where IND
     @boundscheck checkbounds(K, i...)
     prod(getindex.(K.args,i...))
+end
+
+
+MemoryLayout(::Type{<:QuasiKron{<:Any, <:Any, <:Args}}) where Args = kronlayout(tuple_type_memorylayouts(Args)...)
+kronlayout(_...) = QuasiLazyLayout()
+cardinality(K::QuasiKron) = prod(cardinality, K.args)
+
+struct ArrayQuasiVector{T, Arr} <: AbstractQuasiVector{T}
+    array::Arr
+end
+
+ArrayQuasiVector(A::AbstractArray{T}) where T = ArrayQuasiVector{T, typeof(A)}(A)
+
+axes(a::ArrayQuasiVector) = (quasikron(axes(a.array)...),)
+
+function _getindex(::Type{IND}, K::ArrayQuasiVector, i::IND) where IND
+    @boundscheck checkbounds(K, i...)
+    K.array[first(i)...]
+end
+
+function _getindex(::Type{IND}, f::ApplyQuasiVector{<:Any, typeof(*), <:Tuple{QuasiKron,ArrayQuasiVector}}, ij::IND) where IND
+    K,c = f.args
+    A,B = K.args
+    i,j = ij[1]
+    transpose(A[i,:]) * c.array * B[j,:]
 end

test/test_quasikron.jl

@@ -1,4 +1,5 @@
 using QuasiArrays, DomainSets, StaticArrays, Test
+import QuasiArrays: ArrayQuasiVector
 
 @testset "Kron" begin
     @testset "InclusionKron" begin
@@ -12,12 +13,53 @@ using QuasiArrays, DomainSets, StaticArrays, Test
     end
 
     @testset "QuasiKron" begin
-        a = QuasiVector(randn(5), 0:0.5:2)
-        b = QuasiVector([5,6,8], [1,3,4])
-        K = quasikron(a, b)
-        @test axes(K) == (quasikron(axes(a,1),axes(b,1)),)
-        @test K[SVector(0.5,3)] == a[0.5]b[3]
-        @test_throws BoundsError K[SVector(0.6,3)]
-        @test_throws BoundsError K[SVector(0.5,2)]
+        @testset "vec" begin
+            a = QuasiVector(randn(5), 0:0.5:2)
+            b = QuasiVector([5,6,8], [1,3,4])
+            K = quasikron(a, b)
+            @test axes(K) == (quasikron(axes(a,1),axes(b,1)),)
+            @test K[SVector(0.5,3)] == a[0.5]b[3]
+            @test_throws BoundsError K[SVector(0.6,3)]
+            @test_throws BoundsError K[SVector(0.5,2)]
+        end
+        @testset "mat" begin
+            a = QuasiMatrix(randn(5,3), 0:0.5:2, Base.OneTo(3))
+            b = QuasiMatrix(randn(3,2), [1,3,4], Base.OneTo(2))
+            K = QuasiKron(a,b)
+            @test SVector(0.5,3) in axes(K,1)
+            @test SVector(3,2) in axes(K,2)
+            @test K[SVector(0.5,3), SVector(3,2)] ≈ a[0.5,3]*b[3,2]
+        end
+    end
+
+    @testset "ArrayQuasiVector" begin
+        A = randn(2,3)
+        a = ArrayQuasiVector(A)
+        for k in axes(A,1), j in axes(A,2)
+            @test a[(k,j)] == A[k,j]
+        end
+        @test_throws BoundsError a[(3,1)]
+        @test_throws BoundsError a[(1,4)]
+
+        A = randn(2,3,4)
+        a = ArrayQuasiVector(A)
+        for k in axes(A,1), j in axes(A,2), l in axes(A,3)
+            @test a[(k,j,l)] == A[k,j,l]
+        end
+        @test_throws BoundsError a[(3,1,1)]
+        @test_throws BoundsError a[(1,4,1)]
+        @test_throws BoundsError a[(1,1,5)]
+
+        @testset "mul" begin
+            a = QuasiMatrix(randn(5,3), 0:0.5:2, Base.OneTo(3))
+            b = QuasiMatrix(randn(3,2), [1,3,4], Base.OneTo(2))
+            K = QuasiKron(a,b)
+
+            C = randn(3,2)
+            c = ArrayQuasiVector(C)
+            f = K*c
+            @test axes(f,1) == axes(K,1)
+            @test f[SVector(0.5,3)] ≈ a[0.5,:]'*C*b[3,:]
+        end
     end
 end