Add Quasildiv

Project.toml

@@ -11,16 +11,15 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 
 [compat]
-BandedMatrices = "0.14"
-FillArrays = "0.8"
-IntervalSets = "0.3.1"
-LazyArrays = "0.14"
-QuasiArrays = "0.0.4, 0.0.5"
+BandedMatrices = "0.14.1"
+FillArrays = "0.8.2"
+IntervalSets = "0.3.2"
+LazyArrays = "0.14.7"
+QuasiArrays = "0.0.6"
 julia = "1.1"
 
 [extras]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ForwardDiff"]
+test = ["Test"]

src/ContinuumArrays.jl

@@ -1,20 +1,20 @@
 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, union, minimum, maximum
+                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, diff,
+                first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport,
                         adjointlayout, LdivApplyStyle, arguments, _arguments, call, broadcastlayout, lazy_getindex,
-                        sublayout, ApplyLayout, BroadcastLayout, combine_mul_styles
+                        sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
-import FillArrays: AbstractFill, getindex_value
+import FillArrays: AbstractFill, getindex_value, SquareEye
 
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, quasimulapplystyle,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
-                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle
+                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, quasildivapplystyle
 
 export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform
 
@@ -86,6 +86,15 @@ AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*
 
 axes(A::AffineQuasiVector) = axes(A.x)
 getindex(A::AffineQuasiVector, k::Number) = A.A*A.x[k] .+ A.b
+function getindex(A::AffineQuasiVector, k::Inclusion) 
+    @boundscheck A.x[k] # throws bounds error if k ≠ x
+    A
+end
+
+getindex(A::AffineQuasiVector, ::Colon) = copy(A)
+
+copy(A::AffineQuasiVector) = A
+
 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
@@ -107,21 +116,28 @@ broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AffineQuasiVector, b::Numb
 broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
 broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)
 
-function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AffineQuasiVector{<:Real,<:Real,<:Inclusion{<:Real,<:AbstractInterval}})
+function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AffineQuasiVector{<:Any,<:Any,<:Inclusion{<:Any,<:AbstractInterval}})
     @_propagate_inbounds_meta
     isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
 end
 
+igetindex(d::AffineQuasiVector, x) = d.A \ (x .- d.b)
+
 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)
+    @eval $find(f::Base.Fix2{typeof(isequal)}, d::AffineQuasiVector) = $find(isequal(igetindex(d, f.x)), 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))
 
+const QInfAxes = Union{Inclusion,AffineQuasiVector}
 
+sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes}) = V
+sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,QInfAxes}) = V
+sub_materialize(_, V::AbstractQuasiArray, ::Tuple{<:Any,QInfAxes}) = V
+sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V
 
 
 include("operators.jl")

src/bases/bases.jl

@@ -30,56 +30,53 @@ combine_mul_styles(::AbstractAdjointBasisLayout) = LazyQuasiArrayApplyStyle()
 ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
 pinv(J::Basis) = apply(pinv,J)
 
-_multup(a::Tuple) = Mul(a...)
-_multup(a) = a
-
 
 function ==(A::Basis, B::Basis)
     axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
     false
 end
 
+@inline quasildivapplystyle(::AbstractBasisLayout, ::AbstractBasisLayout) = LdivApplyStyle()
+@inline quasildivapplystyle(::AbstractBasisLayout, _) = LdivApplyStyle()
+@inline quasildivapplystyle(_, ::AbstractBasisLayout) = LdivApplyStyle()
 
-ApplyStyle(::typeof(\), ::Type{<:Basis}, ::Type{<:AbstractQuasiMatrix}) = LdivApplyStyle()
-ApplyStyle(::typeof(\), ::Type{<:Basis}, ::Type{<:AbstractQuasiVector}) = LdivApplyStyle()
-ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiMatrix}) = LdivApplyStyle()
-ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiVector}) = LdivApplyStyle()
 
-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))...)
+@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
+@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) = 
+    transform_ldiv(L.A, L.B)
 
-function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)}})
+@inline 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
+@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)},<:Any,<:AbstractQuasiVector}) =
+    transform_ldiv(L.A, L.B)
 
 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))
+    A == B || throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
+    SquareEye{eltype(P)}(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))
+        throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
+    SquareEye{eltype(P)}(size(A,2))
 end
 
-function copy(P::Ldiv{MappedBasisLayout,MappedBasisLayout})
+@inline 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
+
+@inline copy(L::Ldiv{BasisLayout,SubBasisLayout}) = apply(\, L.A, ApplyQuasiArray(L.B))
+@inline function copy(L::Ldiv{SubBasisLayout,BasisLayout}) 
+    P = parent(L.A)
+    kr, jr = parentindices(L.A)
+    lazy_getindex(apply(\, P, L.B), jr, :) # avoid sparse arrays
+end
+
 
 for Bas1 in (:Basis, :WeightedBasis), Bas2 in (:Basis, :WeightedBasis)
     @eval ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2}) =
@@ -88,23 +85,39 @@ end
 
 
 # expansion
-grid(P) = error("Overload Grid")
+_grid(_, P) = error("Overload Grid")
+_grid(::MappedBasisLayout, P) = igetindex.(Ref(parentindices(P)[1]), grid(demap(P)))
+_grid(::SubBasisLayout, P) = grid(parent(P))
+grid(P) = _grid(MemoryLayout(typeof(P)), P)
+
 function transform(L)
     p = grid(L)
     p,L[p,:]
 end
 
-function copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector})
-    p,T = transform(L.A)
-    T \ L.B[p]
+function transform_ldiv(A, B, _)
+    p,T = transform(A)
+    T \ convert(Array, B[p])
 end
 
+transform_ldiv(A, B) = transform_ldiv(A, B, axes(A))
+
+copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector}) =
+    transform_ldiv(L.A, L.B)
+
 copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
-    copy(Ldiv{LazyLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
+    transform_ldiv(L.A, L.B)
+
+function copy(L::Ldiv{ApplyLayout{typeof(*)},<:AbstractBasisLayout})
+    args = arguments(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 = transform(L.A)
-     T \ L.B[p]
+    T \ L.B[p]
 end
 
 ## materialize views
@@ -144,6 +157,9 @@ end
 sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
 sublayout(::BasisLayout, ::Type{<:Tuple{<:AffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()
 
+@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)
 function demap(V::SubQuasiArray{<:Any,2}) 
@@ -172,8 +188,28 @@ function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:Abstract
     A,P
 end
 
-copy(L::Ldiv{BasisLayout,SubBasisLayout}) = apply(\, L.A, ApplyQuasiArray(L.B))
+####
+# sum
+####
 
+_sum(V::AbstractQuasiArray, dims) = __sum(MemoryLayout(typeof(V)), V, dims)
+_sum(V::AbstractQuasiArray, ::Colon) = __sum(MemoryLayout(typeof(V)), V, :)
+sum(V::AbstractQuasiArray; dims=:) = _sum(V, dims)
 
+__sum(L, Vm, _) = error("Override for $L")
+function __sum(::SubBasisLayout, Vm, dims) 
+    @assert dims == 1
+    sum(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
+end
+function __sum(::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
+    a = arguments(V)
+    first(apply(*, sum(a[1]; dims=1), tail(a)...))
+end
+
+function __sum(::MappedBasisLayout, V::AbstractQuasiArray, dims)
+    kr, jr = parentindices(V)
+    @assert kr isa AffineQuasiVector
+    sum(demap(V); dims=dims)/kr.A
+end
 
 include("splines.jl")

src/bases/splines.jl

@@ -39,6 +39,10 @@ end
 
 grid(L::LinearSpline) = L.points
 transform(L::LinearSpline{T}) where T = grid(L),Eye{T}(size(L,2))
+function transform(V::SubQuasiArray{<:Any,2,<:LinearSpline,<:Tuple{<:Inclusion,<:Any}})
+    g, T = transform(parent(V))
+    g, qr(lazy_getindex(T,:,parentindices(V)[2])) # avoid sparse matrices of sub-diagonal
+end
 
 ## Sub-bases
 
@@ -124,4 +128,14 @@ ApplyStyle(::typeof(*), ::Type{<:QuasiAdjoint{<:Any,<:LinearSpline}}, ::Type{<:Q
 function copy(M::QMul2{<:QuasiAdjoint{<:Any,<:LinearSpline},<:QuasiAdjoint{<:Any,<:Derivative}})
     Bc,Ac = M.args
     apply(*,Ac',Bc')'
+end
+
+
+##
+# sum
+##
+
+function _sum(A::HeavisideSpline, dims)
+    @assert dims == 1
+    permutedims(diff(A.points))
 end

src/operators.jl

@@ -153,6 +153,10 @@ axes(D::Derivative) = (D.axis, D.axis)
 ==(a::Derivative, b::Derivative) = a.axis == b.axis
 copy(D::Derivative) = Derivative(copy(D.axis))
 
+function diff(d::AbstractQuasiVector)
+    x = axes(d,1)
+    Derivative(x)*d
+end
 
 
 # struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}