Support SubBasisLayout

src/ContinuumArrays.jl

@@ -2,10 +2,10 @@ module ContinuumArrays
 using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays
 import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
-                first, last, show, isempty, findfirst, findlast, findall, Slice
+                first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted
-import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, 
-                        adjointlayout, LdivApplyStyle, arguments, broadcastlayout, lazy_getindex,
+import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport,
+                        adjointlayout, LdivApplyStyle, arguments, _arguments, call, broadcastlayout, lazy_getindex,
                         sublayout, ApplyLayout, BroadcastLayout, combine_mul_styles
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
@@ -16,7 +16,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform
 
 ####
 # Interval indexing support
@@ -48,7 +48,7 @@ for find in (:findfirst, :findlast)
     @eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
 end
 
-function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion) 
+function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion)
     r = findfirst(f,d)
     r === nothing ? eltype(d)[] : [r]
 end
@@ -79,13 +79,13 @@ end
 AffineQuasiVector(A::AA, x::X, b::B) where {AA,X,B} =
     AffineQuasiVector{promote_type(eltype(AA), eltype(X), eltype(B)),AA,X,B}(A,x,b)
 
-AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))    
+AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))
 AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)
 
 AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)
 
 axes(A::AffineQuasiVector) = axes(A.x)
-getindex(A::AffineQuasiVector, k::Number) = A.A*A.x[k] .+ A.b    
+getindex(A::AffineQuasiVector, k::Number) = A.A*A.x[k] .+ A.b
 inbounds_getindex(A::AffineQuasiVector{<:Any,<:Any,<:Inclusion}, k::Number) = A.A*k .+ A.b
 isempty(A::AffineQuasiVector) = isempty(A.x)
 ==(a::AffineQuasiVector, b::AffineQuasiVector) = a.A == b.A && a.x == b.x && a.b == b.b
@@ -116,6 +116,11 @@ for find in (:findfirst, :findlast, :findall)
     @eval $find(f::Base.Fix2{typeof(isequal)}, d::AffineQuasiVector) = $find(isequal(d.A \ (f.x .- d.b)), d.x)
 end
 
+minimum(d::AffineQuasiVector{<:Real,<:Real,<:Inclusion}) = signbit(d.A) ? last(d) : first(d)
+maximum(d::AffineQuasiVector{<:Real,<:Real,<:Inclusion}) = signbit(d.A) ? first(d) : last(d)
+
+union(d::AffineQuasiVector{<:Real,<:Real,<:Inclusion}) = Inclusion(minimum(d)..maximum(d))
+
 
 
 

src/bases/bases.jl

@@ -5,18 +5,27 @@ abstract type Weight{T} <: LazyQuasiVector{T} end
 const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}
 
 struct WeightLayout <: MemoryLayout end
-struct BasisLayout <: MemoryLayout end
-struct AdjointBasisLayout <: MemoryLayout end
+abstract type AbstractBasisLayout <: MemoryLayout end
+struct BasisLayout <: AbstractBasisLayout end
+struct SubBasisLayout <: AbstractBasisLayout end
+struct MappedBasisLayout <: AbstractBasisLayout end
+
+abstract type AbstractAdjointBasisLayout <: MemoryLayout end
+struct AdjointBasisLayout <: AbstractAdjointBasisLayout end
+struct AdjointSubBasisLayout <: AbstractAdjointBasisLayout end
+struct AdjointMappedBasisLayout <: AbstractAdjointBasisLayout end
 
 MemoryLayout(::Type{<:Basis}) = BasisLayout()
 MemoryLayout(::Type{<:Weight}) = WeightLayout()
 
 adjointlayout(::Type, ::BasisLayout) = AdjointBasisLayout()
-transposelayout(::Type{<:Real}, ::BasisLayout) = AdjointBasisLayout()
+adjointlayout(::Type, ::SubBasisLayout) = AdjointSubBasisLayout()
+adjointlayout(::Type, ::MappedBasisLayout) = AdjointMappedBasisLayout()
 broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = BasisLayout()
+broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = SubBasisLayout()
 
-combine_mul_styles(::BasisLayout) = LazyQuasiArrayApplyStyle()
-combine_mul_styles(::AdjointBasisLayout) = LazyQuasiArrayApplyStyle()
+combine_mul_styles(::AbstractBasisLayout) = LazyQuasiArrayApplyStyle()
+combine_mul_styles(::AbstractAdjointBasisLayout) = LazyQuasiArrayApplyStyle()
 
 ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
 pinv(J::Basis) = apply(pinv,J)
@@ -25,7 +34,7 @@ _multup(a::Tuple) = Mul(a...)
 _multup(a) = a
 
 
-function ==(A::Basis, B::Basis) 
+function ==(A::Basis, B::Basis)
     axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
     false
 end
@@ -36,42 +45,45 @@ ApplyStyle(::typeof(\), ::Type{<:Basis}, ::Type{<:AbstractQuasiVector}) = LdivAp
 ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiMatrix}) = LdivApplyStyle()
 ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiVector}) = LdivApplyStyle()
 
-copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
-function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(-)}}) 
+copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
+copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) = 
+    +(broadcast(\,Ref(L.A),arguments(L.B))...)
+
+function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)}})
+    a,b = arguments(L.B)
+    (L.A\a)-(L.A\b)
+end
+
+function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)},<:Any,<:AbstractQuasiVector})
     a,b = arguments(L.B)
     (L.A\a)-(L.A\b)
 end
 
+function copy(P::Ldiv{BasisLayout,BasisLayout})
+    A, B = P.A, P.B
+    A == B || throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
+    Eye(size(A,2))
+end
+function copy(P::Ldiv{SubBasisLayout,SubBasisLayout})
+    A, B = P.A, P.B
+    (parent(A) == parent(B) && parentindices(A) == parentindices(B)) ||
+        throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
+    Eye(size(A,2))
+end
+
+function copy(P::Ldiv{MappedBasisLayout,MappedBasisLayout})
+    A, B = P.A, P.B
+    demap(A)\demap(B)
+end
+# function copy(P::Ldiv{MappedBasisLayout,SubBasisLayout})
+#     A, B = P.A, P.B
+#     # use lazy_getindex to avoid sparse arrays
+#     lazy_getindex(parent(A)\parent(B),:,parentindices(B)[2])
+# end
+
 for Bas1 in (:Basis, :WeightedBasis), Bas2 in (:Basis, :WeightedBasis)
-    @eval begin
-        function copy(P::Ldiv{<:Any,<:Any,<:$Bas1,<:$Bas2})
-            A, B = P.A, P.B
-            A == B || throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
-            Eye(size(A,2))
-        end
-        function copy(P::Ldiv{<:Any,<:Any,<:SubQuasiArray{<:Any,2,<:$Bas1},<:SubQuasiArray{<:Any,2,<:$Bas2}})
-            A, B = P.A, P.B
-            (parent(A) == parent(B) && parentindices(A) == parentindices(B)) || 
-                throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
-            Eye(size(A,2))
-        end
-
-        function copy(P::Ldiv{<:Any,<:Any,<:SubQuasiArray{<:Any,2,<:$Bas1,<:Tuple{<:AffineQuasiVector,<:Slice}},
-                                          <:SubQuasiArray{<:Any,2,<:$Bas2,<:Tuple{<:AffineQuasiVector,<:Slice}}})
-            A, B = P.A, P.B
-            parent(A)\parent(B)
-        end   
-        function copy(P::Ldiv{<:Any,<:Any,<:SubQuasiArray{<:Any,2,<:$Bas1,<:Tuple{<:AffineQuasiVector,<:Slice}},
-                                          <:SubQuasiArray{<:Any,2,<:$Bas2,<:Tuple{<:AffineQuasiVector,<:Any}}})
-            A, B = P.A, P.B
-            # use lazy_getindex to avoid sparse arrays
-            lazy_getindex(parent(A)\parent(B),:,parentindices(B)[2])
-        end        
-
-        function ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2})
-            all(parentindices(A) == parentindices(B)) && parent(A) == parent(B)
-        end
-    end
+    @eval ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2}) =
+        all(parentindices(A) == parentindices(B)) && parent(A) == parent(B)
 end
 
 
@@ -82,15 +94,15 @@ function transform(L)
     p,L[p,:]
 end
 
-function copy(L::Ldiv{BasisLayout,<:Any,<:Any,<:AbstractQuasiVector})
+function copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector})
     p,T = transform(L.A)
     T \ L.B[p]
 end
 
-copy(L::Ldiv{BasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
+copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
     copy(Ldiv{LazyLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
 
-function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
+function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
     p,T = transform(L.A)
      T \ L.B[p]
 end
@@ -101,7 +113,7 @@ end
 #     *(arguments(S)...)
 
 
-# Differentiation of sub-arrays 
+# Differentiation of sub-arrays
 function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}})
     A, B = M.args
     P = parent(B)
@@ -115,7 +127,7 @@ function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatri
     (Derivative(axes(P,1))*P*kr.A)[kr,jr]
 end
 
-function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix}) 
+function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix})
     args = arguments(L.B)
     # this is a temporary hack
     if args isa Tuple{AbstractQuasiMatrix,Number}
@@ -129,8 +141,29 @@ end
 
 
 # we represent as a Mul with a banded matrix
-sublayout(::BasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = ApplyLayout{typeof(*)}()
-sublayout(::BasisLayout, ::Type{<:Tuple{<:AffineQuasiVector,<:AbstractUnitRange}}) = BasisLayout()
+sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
+sublayout(::BasisLayout, ::Type{<:Tuple{<:AffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()
+
+demap(x) = x
+demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Any,<:Slice}}) = parent(V)
+function demap(V::SubQuasiArray{<:Any,2}) 
+    kr, jr = parentindices(V)
+    demap(parent(V)[kr,:])[:,jr]
+end
+
+
+##
+# SubLayout behaves like ApplyLayout{typeof(*)}
+
+combine_mul_styles(::SubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
+_arguments(::SubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
+call(::SubBasisLayout, ::SubQuasiArray) = *
+
+combine_mul_styles(::AdjointSubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
+_arguments(::AdjointSubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
+arguments(::AdjointSubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
+call(::AdjointSubBasisLayout, ::SubQuasiArray) = *
+
 function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
     A = parent(V)
     _,jr = parentindices(V)
@@ -139,8 +172,8 @@ function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:Abstract
     A,P
 end
 
-
-include("splines.jl")
+copy(L::Ldiv{BasisLayout,SubBasisLayout}) = apply(\, L.A, ApplyQuasiArray(L.B))
 
 
 
+include("splines.jl")

test/runtests.jl

@@ -1,5 +1,5 @@
 using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, ForwardDiff, Test
-import ContinuumArrays: ℵ₁, materialize, SimplifyStyle, AffineQuasiVector, BasisLayout, AdjointBasisLayout
+import ContinuumArrays: ℵ₁, materialize, SimplifyStyle, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, MappedBasisLayout
 import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
 import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout
 import ForwardDiff: Dual
@@ -70,7 +70,7 @@ end
     @test L[2.1,2] ≈ 0.9
     @test L[2.1,3] == L'[3,2.1] == transpose(L)[3,2.1] ≈ 0.1
     @test_throws BoundsError L[3.1,2]
-    
+
     @test L[[1.1,2.1], 1] == L'[1,[1.1,2.1]] == transpose(L)[1,[1.1,2.1]] ≈ [0.9,0.0]
     @test L[1.1,1:2] ≈ [0.9,0.1]
     @test L[[1.1,2.1], 1:2] ≈ [0.9 0.1; 0.0 0.9]
@@ -87,6 +87,34 @@ end
     @test δ'L ≈ [0.8, 0.2, 0.0]
 
     @test L'L == SymTridiagonal([1/3,2/3,1/3], [1/6,1/6])
+
+    @testset "==" begin
+        L = LinearSpline([1,2,3])
+        H = HeavisideSpline([1,2,3])
+        @test L == L
+        @test L ≠ H
+        H = HeavisideSpline([1,1.5,2.5,3])
+        @test_throws ArgumentError L == H
+    end
+
+    @testset "Adjoint layout" begin
+        L = LinearSpline([1,2,3])
+        @test MemoryLayout(typeof(L')) == AdjointBasisLayout()
+        @test [3,4,5]'*L' isa ApplyQuasiArray
+    end
+
+    @testset "Broadcast layout" begin
+        L = LinearSpline([1,2,3])
+        b = BroadcastQuasiArray(+, L*[3,4,5], L*[1.,2,3])
+        @test (L\b) == [4,6,8]
+        B = BroadcastQuasiArray(+, L, L)
+        @test L\B == 2Eye(3)
+
+        b = BroadcastQuasiArray(-, L*[3,4,5], L*[1.,2,3])
+        @test (L\b) == [2,2,2]
+        B = BroadcastQuasiArray(-, L, L)
+        @test L\B == 0Eye(3)
+    end
 end
 
 @testset "Derivative" begin
@@ -168,7 +196,7 @@ end
     @test_throws BoundsError L[0.1,1]
     @test_throws BoundsError L[1.1,0]
 
-    @test MemoryLayout(typeof(L[:,2:3])) isa ApplyLayout{typeof(*)}
+    @test MemoryLayout(typeof(L[:,2:3])) isa SubBasisLayout
     @test L\L[:,2:3] isa BandedMatrix
     @test L\L[:,2:3] == [0 0; 1 0; 0 1.0; 0 0]
 
@@ -222,8 +250,8 @@ end
     @test u[0.1] ≈ 0.00012678835289369413
 end
 
-@testset "Change-of-variables" begin
-    x = Inclusion(0..1) 
+@testset "AffineQuasiVector" begin
+    x = Inclusion(0..1)
     @test 2x isa AffineQuasiVector
     @test (2x)[0.1] == 0.2
     @test_throws BoundsError (2x)[2]
@@ -248,14 +276,33 @@ end
     @test findall(isequal(0.2),y) == [0.6]
     @test findall(isequal(2),y) == Float64[]
 
+    @test AffineQuasiVector(x)[0.1] == 0.1
+    @test minimum(y) == -1
+    @test maximum(y) == 1
+    @test union(y) == Inclusion(-1..1)
+    @test ContinuumArrays.inbounds_getindex(y,0.1) == y[0.1]
+    @test ContinuumArrays.inbounds_getindex(y,2.1) == 2*2.1 - 1
+
+    z = 1 .- x
+    @test minimum(z) == 0.0
+    @test maximum(z) == 1.0
+    @test union(z) == Inclusion(0..1)
+
+    @test !isempty(z)
+    @test z == z
+end
+
+@testset "Change-of-variables" begin
+    x = Inclusion(0..1)
+    y = 2x .- 1
     L = LinearSpline(range(-1,stop=1,length=10))
     @test L[y,:][0.1,:] == L[2*0.1-1,:]
 
     D = Derivative(axes(L,1))
     H = HeavisideSpline(L.points)
     @test H\((D*L) * 2) ≈ (H\(D*L))*2 ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]
 
-    @test MemoryLayout(typeof(L[y,:])) isa BasisLayout
+    @test MemoryLayout(typeof(L[y,:])) isa MappedBasisLayout
     a,b = arguments((D*L)[y,:])
     @test H[y,:]\a == Eye(9)
     @test H[y,:] \ (D*L)[y,:] isa BandedMatrix
@@ -264,4 +311,10 @@ end
     @test (D*L[y,:])[0.1,1] ≈ -9
     @test H[y,:] \ (D*L[y,:]) isa BandedMatrix
     @test H[y,:] \ (D*L[y,:]) ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]
-end
+
+    B = L[y,2:end-1]
+    @test MemoryLayout(typeof(B)) isa MappedBasisLayout
+    @test B[0.1,1] == L[2*0.1-1,2]
+    @test B\B == Eye(8)
+    @test L[y,:] \ B == Eye(10)[:,2:end-1]
+end