Simplify Ldiv

Project.toml

@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -13,18 +13,19 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 
 [compat]
-ArrayLayouts = "0.4.3"
+ArrayLayouts = "0.4.7"
 BandedMatrices = "0.15.17"
 FillArrays = "0.9.3"
 InfiniteArrays = "0.8"
-IntervalSets = "0.3.2, 0.4, 0.5"
-LazyArrays = "0.17.1"
-QuasiArrays = "0.3.1"
+IntervalSets = "0.4, 0.5"
+LazyArrays = "0.18"
+QuasiArrays = "0.3.4"
 julia = "1.5"
 
 [extras]
+Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["FastTransforms", "Test"]
+test = ["Base64", "FastTransforms", "Test"]

src/ContinuumArrays.jl

@@ -88,7 +88,7 @@ function dot(x::Inclusion{T,<:AbstractInterval}, y::Inclusion{V,<:AbstractInterv
     a,b = endpoints(x.domain)
     convert(TV, b^3 - a^3)/3
 end
-    
+
 
 for find in (:findfirst, :findlast)
     @eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
@@ -106,18 +106,15 @@ function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::
 end
 
 
-BroadcastStyle(::Type{<:Inclusion}) = LazyQuasiArrayStyle{1}()
-BroadcastStyle(::Type{<:QuasiAdjoint{<:Any,<:Inclusion}}) = LazyQuasiArrayStyle{2}()
-BroadcastStyle(::Type{<:QuasiTranspose{<:Any,<:Inclusion}}) = LazyQuasiArrayStyle{2}()
-
-
 ###
 # Maps
 ###
 
 # Affine map represents A*x .+ b
 abstract type AbstractAffineQuasiVector{T,AA,X,B} <: AbstractQuasiVector{T} end
 
+show(io::IO, ::MIME"text/plain", a::AbstractAffineQuasiVector) = print(io, "$(a.A) * $(a.x) .+ ($(a.b))")
+
 struct AffineQuasiVector{T,AA,X,B} <: AbstractAffineQuasiVector{T,AA,X,B}
     A::AA
     x::X
@@ -135,7 +132,7 @@ AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*
 axes(A::AbstractAffineQuasiVector) = axes(A.x)
 affine_getindex(A, k) = A.A*A.x[k] .+ A.b
 getindex(A::AbstractAffineQuasiVector, k::Number) = affine_getindex(A, k)
-function getindex(A::AbstractAffineQuasiVector, k::Inclusion) 
+function getindex(A::AbstractAffineQuasiVector, k::Inclusion)
     @boundscheck A.x[k] # throws bounds error if k ≠ x
     A
 end
@@ -188,7 +185,7 @@ struct AffineMap{T,D,R} <: AbstractAffineQuasiVector{T,T,D,T}
     range::R
 end
 
-AffineMap(domain::AbstractQuasiVector{T}, range::AbstractQuasiVector{V}) where {T,V} = 
+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)
@@ -223,6 +220,19 @@ affine(a::AbstractQuasiVector, b) = affine(a, Inclusion(b))
 affine(a, b) = affine(Inclusion(a), Inclusion(b))
 
 
+# mapped vectors
+const AffineMappedQuasiVector = SubQuasiArray{<:Any, 1, <:Any, <:Tuple{AbstractAffineQuasiVector}}
+const AffineMappedQuasiMatrix = SubQuasiArray{<:Any, 2, <:Any, <:Tuple{AbstractAffineQuasiVector,Slice}}
+
+==(a::AffineMappedQuasiVector, b::AffineMappedQuasiVector) = parentindices(a) == parentindices(b) && parent(a) == parent(b)
+
+_sum(V::AffineMappedQuasiVector, ::Colon) = parentindices(V)[1].A \ sum(parent(V))
+
+# pretty print for bases
+show(io::IO, P::AffineMappedQuasiMatrix) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
+show(io::IO, P::AffineMappedQuasiVector) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
+show(io::IO, ::MIME"text/plain", P::AffineMappedQuasiMatrix) = show(io, P)
+show(io::IO, ::MIME"text/plain", P::AffineMappedQuasiVector) = show(io, P)
 
 const QInfAxes = Union{Inclusion,AbstractAffineQuasiVector}
 

src/bases/bases.jl

@@ -10,6 +10,12 @@ struct BasisLayout <: AbstractBasisLayout end
 struct SubBasisLayout <: AbstractBasisLayout end
 struct MappedBasisLayout <: AbstractBasisLayout end
 struct WeightedBasisLayout <: AbstractBasisLayout end
+struct SubWeightedBasisLayout <: AbstractBasisLayout end
+struct MappedWeightedBasisLayout <: AbstractBasisLayout end
+
+SubBasisLayouts = Union{SubBasisLayout,SubWeightedBasisLayout}
+WeightedBasisLayouts = Union{WeightedBasisLayout,SubWeightedBasisLayout,MappedWeightedBasisLayout}
+MappedBasisLayouts = Union{MappedBasisLayout,MappedWeightedBasisLayout}
 
 abstract type AbstractAdjointBasisLayout <: AbstractQuasiLazyLayout end
 struct AdjointBasisLayout <: AbstractAdjointBasisLayout end
@@ -19,11 +25,21 @@ struct AdjointMappedBasisLayout <: AbstractAdjointBasisLayout end
 MemoryLayout(::Type{<:Basis}) = BasisLayout()
 MemoryLayout(::Type{<:Weight}) = WeightLayout()
 
-adjointlayout(::Type, ::BasisLayout) = AdjointBasisLayout()
+adjointlayout(::Type, ::AbstractBasisLayout) = AdjointBasisLayout()
 adjointlayout(::Type, ::SubBasisLayout) = AdjointSubBasisLayout()
-adjointlayout(::Type, ::MappedBasisLayout) = AdjointMappedBasisLayout()
+adjointlayout(::Type, ::MappedBasisLayouts) = AdjointMappedBasisLayout()
 broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = WeightedBasisLayout()
 broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = WeightedBasisLayout()
+broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::MappedBasisLayouts) = MappedWeightedBasisLayout()
+
+# A sub of a weight is still a weight
+sublayout(::WeightLayout, _) = WeightLayout()
+
+## Weighted basis interface
+unweightedbasis(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = last(P.args)
+unweightedbasis(V::SubQuasiArray) = view(unweightedbasis(parent(V)), parentindices(V)...)
+
+
 
 # Default is lazy
 ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
@@ -48,25 +64,25 @@ end
 @inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)},<:Any,<:AbstractQuasiVector}) =
     transform_ldiv(L.A, L.B)
 
-function copy(P::Ldiv{BasisLayout,BasisLayout})
+function copy(P::Ldiv{<:AbstractBasisLayout,<:AbstractBasisLayout})
     A, B = P.A, P.B
     A == B || throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
     SquareEye{eltype(P)}((axes(A,2),))
 end
-function copy(P::Ldiv{SubBasisLayout,SubBasisLayout})
+function copy(P::Ldiv{<:SubBasisLayouts,<:SubBasisLayouts})
     A, B = P.A, P.B
     parent(A) == parent(B) ||
         throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
     Eye{eltype(P)}((axes(A,2),axes(B,2)))
 end
 
-@inline function copy(P::Ldiv{MappedBasisLayout,MappedBasisLayout})
+@inline function copy(P::Ldiv{<:MappedBasisLayouts,<:MappedBasisLayouts})
     A, B = P.A, P.B
     demap(A)\demap(B)
 end
 
-@inline copy(L::Ldiv{BasisLayout,SubBasisLayout}) = apply(\, L.A, ApplyQuasiArray(L.B))
-@inline function copy(L::Ldiv{SubBasisLayout,BasisLayout}) 
+@inline copy(L::Ldiv{<:AbstractBasisLayout,<:SubBasisLayouts}) = apply(\, L.A, ApplyQuasiArray(L.B))
+@inline function copy(L::Ldiv{<:SubBasisLayouts,<:AbstractBasisLayout}) 
     P = parent(L.A)
     kr, jr = parentindices(L.A)
     layout_getindex(apply(\, P, L.B), jr, :) # avoid sparse arrays
@@ -83,7 +99,7 @@ end
 _grid(_, P) = error("Overload Grid")
 _grid(::MappedBasisLayout, P) = igetindex.(Ref(parentindices(P)[1]), grid(demap(P)))
 _grid(::SubBasisLayout, P) = grid(parent(P))
-_grid(::WeightedBasisLayout, P) = grid(last(P.args))
+_grid(::WeightedBasisLayouts, P) = grid(unweightedbasis(P))
 grid(P) = _grid(MemoryLayout(typeof(P)), P)
 
 struct TransformFactorization{T,Grid,Plan,IPlan} <: Factorization{T}
@@ -121,6 +137,11 @@ end
 \(a::ProjectionFactorization, b::AbstractVector) = (a.F \ b)[a.inds]
 
 _factorize(::SubBasisLayout, L) = ProjectionFactorization(factorize(parent(L)), parentindices(L)[2])
+# function _factorize(::MappedBasisLayout, L)
+#     kr, jr = parentindices(L)
+#     P = parent(L)
+#     ProjectionFactorization(factorize(view(P,:,jr)), parentindices(L)[2])
+# end
 
 transform_ldiv(A, B, _) = factorize(A) \ B
 transform_ldiv(A, B) = transform_ldiv(A, B, axes(A))
@@ -131,18 +152,6 @@ copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector}) =
 copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
     transform_ldiv(L.A, L.B)
 
-function copy(L::Ldiv{ApplyLayout{typeof(*)},<:AbstractBasisLayout})
-    args = arguments(ApplyLayout{typeof(*)}(), L.A)
-    @assert length(args) == 2 # temporary
-    apply(\, last(args), apply(\, first(args), L.B))
-end
-
-
-function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
-    p,T = factorize(L.A)
-    T \ L.B[p]
-end
-
 
 ##
 # Algebra
@@ -151,9 +160,10 @@ end
 # struct ExpansionLayout <: MemoryLayout end
 # applylayout(::Type{typeof(*)}, ::BasisLayout, _) = ExpansionLayout()
 
-const Expansion{T,Space<:Basis,Coeffs<:AbstractVector} = ApplyQuasiVector{T,typeof(*),<:Tuple{Space,Coeffs}}
+const Expansion{T,Space<:AbstractQuasiMatrix,Coeffs<:AbstractVector} = ApplyQuasiVector{T,typeof(*),<:Tuple{Space,Coeffs}}
 
-basis(v::AbstractQuasiVector) = v.args[1]
+
+basis(v::Expansion) = v.args[1]
 
 for op in (:*, :\)
     @eval function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Number, f::Expansion)
@@ -169,11 +179,19 @@ for op in (:*, :/)
 end
 
 
-function broadcastbasis(::typeof(+), a, b)
-    a ≠ b && error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
+function _broadcastbasis(::typeof(+), _, _, a, b)
+    try
+        a ≠ b && error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
+    catch
+        error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
+    end
     a
 end
 
+_broadcastbasis(::typeof(+), ::MappedBasisLayouts, ::MappedBasisLayouts, a, b) = broadcastbasis(+, demap(a), demap(b))[basismap(a), :]
+
+broadcastbasis(::typeof(+), a, b) = _broadcastbasis(+, MemoryLayout(a), MemoryLayout(b), a, b)
+
 broadcastbasis(::typeof(-), a, b) = broadcastbasis(+, a, b)
 
 for op in (:+, :-)
@@ -228,48 +246,43 @@ end
     (Derivative(axes(P,1))*P*kr.A)[kr,jr]
 end
 
-function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix})
-    args = arguments(L.B)
-    # this is a temporary hack
-    if args isa Tuple{AbstractQuasiMatrix,Number}
-        (L.A \  first(args))*last(args)
-    elseif args isa Tuple{Number,AbstractQuasiMatrix}
-        first(args)*(L.A \ last(args))
-    else
-        error("Not implemented")
-    end
-end
-
-
 # we represent as a Mul with a banded matrix
 sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
 sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()
+sublayout(::WeightedBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedWeightedBasisLayout()
+sublayout(::WeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubWeightedBasisLayout()
 
 @inline sub_materialize(::AbstractBasisLayout, V::AbstractQuasiArray) = V
 @inline sub_materialize(::AbstractBasisLayout, V::AbstractArray) = V
 
 demap(x) = x
-demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Any,<:Slice}}) = parent(V)
+demap(x::BroadcastQuasiArray) = BroadcastQuasiArray(x.f, map(demap, arguments(x))...)
+demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) = parent(V)
+demap(V::SubQuasiArray{<:Any,1,<:Any,<:Tuple{<:AbstractAffineQuasiVector}}) = parent(V)
 function demap(V::SubQuasiArray{<:Any,2}) 
     kr, jr = parentindices(V)
     demap(parent(V)[kr,:])[:,jr]
 end
 
+basismap(x::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) = parentindices(x)[1]
+basismap(x::SubQuasiArray{<:Any,1,<:Any,<:Tuple{<:AbstractAffineQuasiVector}}) = parentindices(x)[1]
+basismap(x::BroadcastQuasiArray) = basismap(x.args[1])
+
 
 ##
 # SubLayout behaves like ApplyLayout{typeof(*)}
 
-combine_mul_styles(::SubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
-_mul_arguments(::SubBasisLayout, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
-arguments(::SubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
-call(::SubBasisLayout, ::SubQuasiArray) = *
+combine_mul_styles(::SubBasisLayouts) = combine_mul_styles(ApplyLayout{typeof(*)}())
+_mul_arguments(::SubBasisLayouts, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
+arguments(::SubBasisLayouts, A) = arguments(ApplyLayout{typeof(*)}(), A)
+call(::SubBasisLayouts, ::SubQuasiArray) = *
 
 combine_mul_styles(::AdjointSubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
 _mul_arguments(::AdjointSubBasisLayout, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
 arguments(::AdjointSubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
 call(::AdjointSubBasisLayout, ::SubQuasiArray) = *
 
-copy(M::Mul{AdjointSubBasisLayout,SubBasisLayout}) = apply(*, arguments(M.A)..., arguments(M.B)...)
+copy(M::Mul{AdjointSubBasisLayout,<:SubBasisLayouts}) = apply(*, arguments(M.A)..., arguments(M.B)...)
 
 function arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
     A = parent(V)
@@ -297,8 +310,8 @@ function __sum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
     first(apply(*, sum(a[1]; dims=1), tail(a)...))
 end
 
-function __sum(::MappedBasisLayout, V::AbstractQuasiArray, dims)
-    kr, jr = parentindices(V)
+function __sum(::MappedBasisLayouts, V::AbstractQuasiArray, dims)
+    kr = basismap(V)
     @assert kr isa AbstractAffineQuasiVector
     sum(demap(V); dims=dims)/kr.A
 end

src/operators.jl

@@ -85,6 +85,9 @@ end
 @simplify *(A::QuasiAdjoint{<:Any,<:DiracDelta}, B::AbstractQuasiVector) = B[parent(A).x]
 @simplify *(A::QuasiAdjoint{<:Any,<:DiracDelta}, B::AbstractQuasiMatrix) = B[parent(A).x,:]
 
+show(io::IO, δ::DiracDelta) = print(io, "δ at $(δ.x) over $(axes(δ,1))")
+show(io::IO, ::MIME"text/plain", δ::DiracDelta) = show(io, δ)
+
 #########
 # Derivative
 #########