QuasiArrays v0.5 (#85)

* tests pass

* update tests

* v0.6.5

* Hcat improvements, transform matrix

* Update Project.toml

* add tests

* Test factorization

Author: dlfivefifty

Project.toml

@@ -1,11 +1,12 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.6.4"
+version = "0.7.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
@@ -16,15 +17,15 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-ArrayLayouts = "0.5, 0.6"
+ArrayLayouts = "0.6"
 BandedMatrices = "0.16"
 BlockArrays = "0.15.1"
 FillArrays = "0.11"
-InfiniteArrays = "0.9, 0.10"
-Infinities = "0.0.1, 0.0.2, 0.1"
+InfiniteArrays = "0.10"
+Infinities = "0.1"
 IntervalSets = "0.5"
-LazyArrays = "0.20, 0.21"
-QuasiArrays = "0.4.1"
+LazyArrays = "0.21"
+QuasiArrays = "0.5.1"
 StaticArrays = "0.12, 1"
 julia = "1.5"
 

src/ContinuumArrays.jl

@@ -8,7 +8,7 @@ import Base.Broadcast: materialize, BroadcastStyle, broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, most, combine_mul_styles, AbstractArrayApplyStyle,
                         adjointlayout, arguments, _mul_arguments, call, broadcastlayout, layout_getindex, UnknownLayout,
                         sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles, applylayout,
-                        simplifiable, _simplify
+                        simplifiable, _simplify, AbstractLazyLayout
 import LinearAlgebra: pinv, dot, norm2
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 import BlockArrays: block, blockindex, unblock, blockedrange, _BlockedUnitRange, _BlockArray
@@ -17,7 +17,8 @@ import ArrayLayouts: mul
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
-                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize
+                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize,
+                    AbstractQuasiFill, UnionDomain
 import InfiniteArrays: Infinity, InfAxes
 
 export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine, ..
@@ -78,6 +79,7 @@ end
     (inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)    
 
 checkpoints(d::AbstractInterval{T}) where T = width(d) .* SVector{3,float(T)}(0.823972,0.01,0.3273484) .+ leftendpoint(d)
+checkpoints(d::UnionDomain) = mapreduce(checkpoints,union,d.domains)
 checkpoints(x::Inclusion) = checkpoints(x.domain)
 checkpoints(A::AbstractQuasiMatrix) = checkpoints(axes(A,1))
 

src/bases/bases.jl

@@ -139,10 +139,21 @@ grid(T::TransformFactorization) = T.grid
 
 \(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiVector{T}) where T = a.iplan \  convert(Array{T}, b[a.grid])
 \(a::TransformFactorization, b::AbstractQuasiVector) = a.plan * convert(Array, b[a.grid])
-
 \(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVector) = a.iplan \  b
 \(a::TransformFactorization, b::AbstractVector) = a.plan * b
 
+\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiMatrix{T}) where T = a \  convert(Array{T}, b[a.grid,:])
+\(a::TransformFactorization, b::AbstractQuasiMatrix) = a \ convert(Array, b[a.grid,:])
+function \(a::TransformFactorization, b::AbstractMatrix)
+    c = a \ b[:,1]
+    ret = Array{eltype(c)}(undef, length(c), size(b,2))
+    ret[:,1] = c
+    for k = 2:size(ret,2)
+        ret[:,k] = a \ b[:,k]
+    end
+    ret
+end
+
 function _factorize(::AbstractBasisLayout, L)
     p = grid(L)
     TransformFactorization(p, nothing, factorize(L[p,:]))
@@ -177,11 +188,11 @@ end
 transform_ldiv(A, B, _) = factorize(A) \ B
 transform_ldiv(A, B) = transform_ldiv(A, B, size(A))
 
-copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector}) =
-    transform_ldiv(L.A, L.B)
-
-copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
-    transform_ldiv(L.A, L.B)
+copy(L::Ldiv{<:AbstractBasisLayout}) = transform_ldiv(L.A, L.B)
+# TODO: redesign to use simplifiable(\, A, B)
+copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = transform_ldiv(L.A, L.B)
+copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
+copy(L::Ldiv{<:AbstractBasisLayout,<:AbstractLazyLayout}) = transform_ldiv(L.A, L.B)
 
 struct WeightedFactorization{T, WW, FAC<:Factorization{T}} <: Factorization{T}
     w::WW

src/basisconcat.jl

@@ -108,3 +108,6 @@ function QuasiArrays._getindex(::Type{IND}, A::HvcatBasis{T}, (x,j)::IND) where
     J = Int(block(Jj))
     hvcat(A.n, zeros(J-1)..., A.args[J][x, blockindex(Jj)], zeros(length(A.args)-J)...)::T
 end
+
+
+@simplify *(D::Derivative, H::ApplyQuasiMatrix{<:Any,typeof(hcat)}) = hcat((Ref(D) .* H.args)...)

src/maps.jl

@@ -91,8 +91,9 @@ broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AbstractAffineQuasiVector,
 broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
 broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)
 
-Base.@propagate_inbounds checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval}) = QuasiArrays._affine_checkindex(inds, r)
-Base.@propagate_inbounds checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AbstractAffineQuasiVector) = QuasiArrays._affine_checkindex(inds, r)
+Base.@propagate_inbounds _affine_checkindex(inds, r) = isempty(r) | (checkindex(Bool, inds, Base.to_indices(inds, (first(r),))...) & checkindex(Bool, inds,  Base.to_indices(inds, (last(r),))...))
+Base.@propagate_inbounds checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval}) = _affine_checkindex(inds, r)
+Base.@propagate_inbounds checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AbstractAffineQuasiVector) = _affine_checkindex(inds, r)
 
 minimum(d::AbstractAffineQuasiVector) = signbit(d.A) ? last(d) : first(d)
 maximum(d::AbstractAffineQuasiVector) = signbit(d.A) ? first(d) : last(d)
@@ -112,7 +113,7 @@ end
 AffineMap(domain::AbstractQuasiVector{T}, range::AbstractQuasiVector{V}) where {T,V} =
     AffineMap{promote_type(T,V), typeof(domain),typeof(range)}(domain,range)
 
-measure(x::Inclusion) = last(x)-first(x)
+measure(x::Inclusion{<:Any,<:AbstractInterval}) = last(x)-first(x)
 
 function getproperty(A::AffineMap, d::Symbol)
     domain, range = getfield(A, :domain), getfield(A, :range)

src/operators.jl

@@ -68,7 +68,6 @@ end
 DiracDelta{T}(x, axis::A) where {T,A<:AbstractQuasiVector} = DiracDelta{T,A}(x, axis)
 DiracDelta{T}(x, domain) where T = DiracDelta{T}(x, Inclusion(domain))
 DiracDelta(x, domain) = DiracDelta{float(eltype(x))}(x, Inclusion(domain))
-DiracDelta(axis::AbstractQuasiVector) = DiracDelta(zero(float(eltype(axis))), axis)
 DiracDelta(x) = DiracDelta(x, Inclusion(x))
 DiracDelta{T}() where T = DiracDelta(zero(T))
 DiracDelta() = DiracDelta{Float64}()
@@ -129,3 +128,6 @@ end
 const Identity{T,D} = QuasiDiagonal{T,Inclusion{T,D}}
 
 Identity(d::Inclusion) = QuasiDiagonal(d)
+
+@simplify *(D::Derivative, x::Inclusion) = ones(promote_type(eltype(D),eltype(x)), x)
+@simplify *(D::Derivative, c::AbstractQuasiFill) = zeros(promote_type(eltype(D),eltype(c)), axes(c,1))

test/runtests.jl

@@ -7,25 +7,37 @@ import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, App
 
 
 @testset "Inclusion" begin
-    x = Inclusion(-1..1)
-    @test eltype(x) == Float64
-    @test x[0.1] ≡ 0.1
-    @test x[0] ≡ x[0.0] ≡ 0.0
-
-    x = Inclusion(-1.0..1)
-    X = QuasiDiagonal(x)
-    @test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
-    @test Base.BroadcastStyle(typeof(x)) == LazyQuasiArrayStyle{1}()
-    @test x .* x isa BroadcastQuasiArray
-    @test (x.*x)[0.1] == 0.1^2
-
-    @test exp.(x)[1.0] == exp(1.0)
-    @test exp.(-x)[1.0] == exp(-1.0)
-    @test exp.(x.^2)[1.0] == exp(1.0)
-    @test exp.(-x.^2/(20/2))[1.0] == exp(-1.0^2/(20/2))
-
-    @test dot(x,x) ≈ 2/3
-    @test norm(x) ≈ sqrt(2/3)
+    @testset "basics" begin
+        x = Inclusion(-1..1)
+        @test eltype(x) == Float64
+        @test x[0.1] ≡ 0.1
+        @test x[0] ≡ x[0.0] ≡ 0.0
+    end
+
+    @testset "broadcast" begin
+        x = Inclusion(-1.0..1)
+        X = QuasiDiagonal(x)
+        @test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
+        @test Base.BroadcastStyle(typeof(x)) == LazyQuasiArrayStyle{1}()
+        @test x .* x isa BroadcastQuasiArray
+        @test (x.*x)[0.1] == 0.1^2
+
+        @test exp.(x)[1.0] == exp(1.0)
+        @test exp.(-x)[1.0] == exp(-1.0)
+        @test exp.(x.^2)[1.0] == exp(1.0)
+        @test exp.(-x.^2/(20/2))[1.0] == exp(-1.0^2/(20/2))
+
+        @test dot(x,x) ≈ 2/3
+        @test norm(x) ≈ sqrt(2/3)
+    end
+
+    @testset "Derivative" begin
+        x = Inclusion(-1..1)
+        D = Derivative(x)
+        @test D*x ≡ QuasiOnes(x)
+        @test D^2 * x ≡ QuasiZeros(x)
+        @test D*[x D*x] == [D*x D^2*x]
+    end
 end
 
 include("test_maps.jl")
@@ -76,7 +88,7 @@ end
         @test_throws BoundsError H[[0.1,2.1], 1]
         @test MemoryLayout(typeof(H)) == BasisLayout()
         @test ApplyStyle(*, typeof(H), typeof([1,2])) isa MulStyle
-        
+
         f = H*[1,2]
         @test f isa ApplyQuasiArray
         @test axes(f) == (Inclusion(1.0..3.0),)
@@ -189,6 +201,8 @@ end
         f = L*[1,2,4]
 
         D = Derivative(axes(L,1))
+        @test copy(D) == D
+
         @test D*L isa MulQuasiMatrix
         @test length((D*L).args) == 2
         @test eltype(D*L) == Float64
@@ -316,6 +330,8 @@ end
             @test (L \ x) == [1,2,3]
             @test factorize(L[:,2:end-1]) isa ContinuumArrays.ProjectionFactorization
             @test L[:,1:2] \ x == [1,2]
+            @test L \ [x one(x)] ≈ [L\x L\one(x)]
+            @test factorize(L) \ QuasiOnes(x, Base.OneTo(3)) ≈ L \ QuasiOnes(x, Base.OneTo(3)) ≈ ones(3,3)
 
             L = LinearSpline(range(0,1; length=10_000))
             x = axes(L,1)
@@ -494,6 +510,8 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
         g = grid(F)
         @test T \ exp.(x) == F \ exp.(x) == F \ exp.(g) == chebyshevtransform(exp.(g), Val(1))
         @test all(checkpoints(T) .∈ Ref(axes(T,1)))
+
+        @test F \ [exp.(x) cos.(x)] ≈ [F \ exp.(x) F \ cos.(x)]
     end
 
     @testset "Weighted" begin
@@ -533,7 +551,7 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
             @test 2 ∉ m
             @test 0.1 ∈ mi
             @test -0.1 ∉ mi
-            
+
             @test m[findfirst(isequal(0.1), m)] ≈ 0.1
             @test m[findlast(isequal(0.1), m)] ≈ 0.1
             @test m[findall(isequal(0.1), m)] ≈ [0.1]

test/test_basisconcat.jl

@@ -1,5 +1,5 @@
 using ContinuumArrays, BlockArrays, Test
-import ContinuumArrays: PiecewiseBasis, VcatBasis, HvcatBasis, arguments, ApplyLayout
+import ContinuumArrays: PiecewiseBasis, VcatBasis, HvcatBasis, arguments, ApplyLayout, checkpoints
 
 @testset "ConcatBasis" begin
     @testset "PiecewiseBasis" begin
@@ -8,6 +8,7 @@ import ContinuumArrays: PiecewiseBasis, VcatBasis, HvcatBasis, arguments, ApplyL
         S = PiecewiseBasis(S1, S2)
 
         @test S == S
+        @test checkpoints(S) == union(checkpoints(S1), checkpoints(S2))
 
         @test S[0.5,1:4] == [S1[0.5,1:2]; zeros(2)]
         @test S[2.5,1:4] == [zeros(2); S2[2.5,1:2]]