Only use transform_ldiv when feasible (#114)

* Only use transform_ldiv when feasible

* Update bases.jl

* Update bases.jl

* WeightLayout transform_ldiv

* transform only if finite dimensional columns

* tests pass

* Increase coverage

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.8.7"
+version = "0.9.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -17,15 +17,15 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-ArrayLayouts = "0.7"
+ArrayLayouts = "0.7.2"
 BandedMatrices = "0.16"
-BlockArrays = "0.15.1, 0.16"
+BlockArrays = "0.16"
 FillArrays = "0.11, 0.12"
-InfiniteArrays = "0.10, 0.11"
+InfiniteArrays = "0.11"
 Infinities = "0.1"
 IntervalSets = "0.5"
-LazyArrays = "0.21"
-QuasiArrays = "0.6, 0.7"
+LazyArrays = "0.21.13"
+QuasiArrays = "0.7"
 RecipesBase = "1.0"
 StaticArrays = "1.0"
 julia = "1.6"

src/ContinuumArrays.jl

@@ -15,7 +15,7 @@ import BlockArrays: block, blockindex, unblock, blockedrange, _BlockedUnitRange,
 import FillArrays: AbstractFill, getindex_value, SquareEye
 import ArrayLayouts: mul, ZerosLayout, ScalarLayout
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
-                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
+                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, QuasiArrayLayout,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize,
                     AbstractQuasiFill, UnionDomain, __sum
@@ -78,6 +78,7 @@ end
 @inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{AbstractArray{<:BlockIndex{1}}, Vararg{Any}}) =
     (inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)    
 
+checkpoints(x::Number) = x
 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)

src/bases/bases.jl

@@ -57,7 +57,7 @@ end
 
 
 @inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
-@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) = 
+@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) =
     transform_ldiv(L.A, L.B)
 for op in (:+, :-)
     @eval @inline copy(L::Ldiv{Lay,BroadcastLayout{typeof($op)},<:Any,<:AbstractQuasiVector}) where Lay<:MappedBasisLayouts =
@@ -75,7 +75,7 @@ end
 @inline 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(eltype(P))}((axes(A,2),)) # use double eltype for array-valued 
+    SquareEye{eltype(eltype(P))}((axes(A,2),)) # use double eltype for array-valued
 end
 @inline function copy(P::Ldiv{<:SubBasisLayouts,<:SubBasisLayouts})
     A, B = P.A, P.B
@@ -97,12 +97,15 @@ copy(L::Ldiv{<:MappedBasisLayouts,ApplyLayout{typeof(*)}}) = copy(Ldiv{UnknownLa
 copy(L::Ldiv{<:MappedBasisLayouts,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = transform_ldiv(L.A, L.B)
 
 @inline copy(L::Ldiv{<:AbstractBasisLayout,<:SubBasisLayouts}) = apply(\, L.A, ApplyQuasiArray(L.B))
-@inline function copy(L::Ldiv{<:SubBasisLayouts,<:AbstractBasisLayout}) 
+@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
 end
 
+# default to transform for expanding weights
+copy(L::Ldiv{<:AbstractBasisLayout,WeightLayout}) = transform_ldiv(L.A, L.B)
+
 
 for Bas1 in (:Basis, :WeightedBasis), Bas2 in (:Basis, :WeightedBasis)
     @eval ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2}) =
@@ -185,7 +188,7 @@ associates a planned inverse transform with a grid. That is, if `F` is a `Transf
 """
 TransformFactorization(grid, ::Nothing, iplan) = TransformFactorization{promote_type(eltype(eltype(grid)),eltype(iplan))}(grid, nothing, iplan)
 
-grid(T::TransformFactorization) = T.grid    
+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])
@@ -259,13 +262,20 @@ 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)
+# A BroadcastLayout of unknown function is only knowable pointwise
+transform_ldiv_if_columns(A, B, _) = ApplyQuasiArray(\, A, B)
+transform_ldiv_if_columns(A, B, ::Base.OneTo) = transform_ldiv(A,B)
+copy(L::Ldiv{<:AbstractBasisLayout,<:BroadcastLayout}) = transform_ldiv_if_columns(L.A, L.B, axes(L.B,2))
+# Inclusion are QuasiArrayLayout
+copy(L::Ldiv{<:AbstractBasisLayout,QuasiArrayLayout}) = transform_ldiv(L.A, L.B)
+# Otherwise keep lazy to support, e.g., U\D*T
+copy(L::Ldiv{<:AbstractBasisLayout,<:AbstractLazyLayout}) = ApplyQuasiArray(\, L.A, L.B)
 copy(L::Ldiv{<:AbstractBasisLayout,ZerosLayout}) = Zeros{eltype(L)}(axes(L)...)
 
 """
     WeightedFactorization(w, F)
 
-weights a factorization `F` by `w`. 
+weights a factorization `F` by `w`.
 """
 struct WeightedFactorization{T, WW, FAC<:Factorization{T}} <: Factorization{T}
     w::WW
@@ -416,7 +426,7 @@ demap(x) = x
 demap(x::BroadcastQuasiArray) = BroadcastQuasiArray(x.f, map(demap, arguments(x))...)
 demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Any,Slice}}) = parent(V)
 demap(V::SubQuasiArray{<:Any,1}) = parent(V)
-function demap(V::SubQuasiArray{<:Any,2}) 
+function demap(V::SubQuasiArray{<:Any,2})
     kr, jr = parentindices(V)
     demap(parent(V)[kr,:])[:,jr]
 end
@@ -466,7 +476,7 @@ end
 ####
 
 
-function __sum(::SubBasisLayout, Vm, dims) 
+function __sum(::SubBasisLayout, Vm, dims)
     @assert dims == 1
     sum(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
 end

src/bases/splines.jl

@@ -7,6 +7,7 @@ const LinearSpline = Spline{1}
 const HeavisideSpline = Spline{0}
 
 Spline{o}(pts::AbstractVector{T}) where {o,T} = Spline{o,float(T)}(pts)
+Spline{o}(S::Spline) where {o} = Spline{o}(S.points)
 
 axes(B::Spline{o}) where o =
     (Inclusion(first(B.points)..last(B.points)), OneTo(length(B.points)+o-1))