Add AffineMap

Project.toml

@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.0.1"
+version = "0.0.2"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
@@ -11,11 +11,11 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 
 [compat]
-BandedMatrices = "0.12, 0.13"
-FillArrays = "0.7"
+BandedMatrices = "0.14"
+FillArrays = "0.8"
 IntervalSets = "0.3.1"
-LazyArrays = "0.12, 0.14"
-QuasiArrays = "0.0.2, 0.0.3, 0.0.4"
+LazyArrays = "0.14"
+QuasiArrays = "0.0.4"
 julia = "1.1"
 
 [extras]

examples/heat.jl

@@ -0,0 +1,37 @@
+using ContinuumArrays, DifferentialEquations, Plots
+
+##
+# natural boundary conditions
+##
+
+L = LinearSpline(range(-1,1; length=20))
+x = axes(L,1)
+D = Derivative(x)
+M = L'L
+Δ = -((D*L)'D*L)
+u0 = copy(L \ exp.(x))
+
+heat(u,(M,Δ),t) = M\(Δ*u)
+prob = ODEProblem(heat,u0,(0.0,1.0),(M,Δ))
+sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
+
+g = range(-1,1;length=1000)
+@gif for t in 0.0:0.01:1.0
+    plot(g, L[g,:]*sol(t); ylims=(0,3))
+end
+
+##
+# Dirichlet boundary conditions
+##
+
+S = L[:,2:end-1]
+M = S'S
+Δ = -((D*S)'D*S)
+u0 = (L \ broadcast(x -> (1-x^2)*exp(x), x))[2:end-1]
+prob = ODEProblem(heat,u0,(0.0,1.0),(M,Δ))
+sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
+
+g = range(-1,1;length=1000)
+@gif for t in 0.0:0.01:1.0
+    plot(g, S[g,:]*sol(t); ylims=(0,3))
+end

src/ContinuumArrays.jl

@@ -1,20 +1,22 @@
 module ContinuumArrays
 using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays
-import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
+import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
-                first, last, show
-import Base.Broadcast: materialize, BroadcastStyle
-import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout, LdivApplyStyle
+                first, last, show, isempty, findfirst, findlast, findall, Slice
+import Base.Broadcast: materialize, BroadcastStyle, broadcasted
+import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, 
+                        adjointlayout, LdivApplyStyle, arguments, broadcastlayout, lazy_getindex,
+                        sublayout, ApplyLayout, BroadcastLayout, combine_mul_styles
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 import FillArrays: AbstractFill, getindex_value
 
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
-                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
+                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, quasimulapplystyle,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid
 
 ####
 # Interval indexing support
@@ -38,35 +40,85 @@ const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}
 cardinality(::AbstractInterval) = ℵ₁
 *(ℵ::AlephInfinity) = ℵ
 
+Inclusion(d::AbstractInterval{T}) where T = Inclusion{float(T)}(d)
 first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
 last(S::Inclusion{<:Any,<:AbstractInterval}) = rightendpoint(S.domain)
 
+for find in (:findfirst, :findlast)
+    @eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
+end
 
-checkindex(::Type{Bool}, inds::AbstractInterval, i::Number) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
-checkindex(::Type{Bool}, inds::AbstractInterval, i::Inclusion) = i.domain ⊆ inds
-function checkindex(::Type{Bool}, inds::AbstractInterval, I::AbstractArray)
-    @_inline_meta
-    b = true
-    for i in I
-        b &= checkindex(Bool, inds, i)
-    end
-    b
+function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion) 
+    r = findfirst(f,d)
+    r === nothing ? eltype(d)[] : [r]
 end
 
 
-# we represent as a Mul with a banded matrix
-function materialize(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
-    A = parent(V)
-    _,jr = parentindices(V)
-    first(jr) ≥ 1 || throw(BoundsError())
-    P = _BandedMatrix(Ones{Int}(1,length(jr)), axes(A,2), first(jr)-1,1-first(jr))
-    A*P
+function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval})
+    @_propagate_inbounds_meta
+    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
 end
 
+
 BroadcastStyle(::Type{<:Inclusion{<:Any,<:AbstractInterval}}) = LazyQuasiArrayStyle{1}()
 BroadcastStyle(::Type{<:QuasiAdjoint{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()
 BroadcastStyle(::Type{<:QuasiTranspose{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()
 
+
+###
+# Maps
+###
+
+# Affine map represents A*x .+ b
+struct AffineQuasiVector{T,AA,X,B} <: AbstractQuasiVector{T}
+    A::AA
+    x::X
+    b::B
+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(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    
+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
+
+BroadcastStyle(::Type{<:AffineQuasiVector}) = LazyQuasiArrayStyle{1}()
+
+for op in(:*, :\, :+, :-)
+    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), a::Number, x::Inclusion) = broadcast($op, a, AffineQuasiVector(x))
+end
+for op in(:/, :+, :-)
+    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Inclusion, a::Number) = broadcast($op, AffineQuasiVector(x), a)
+end
+
+broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Number, x::AffineQuasiVector) = AffineQuasiVector(a, x)
+broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(\), a::Number, x::AffineQuasiVector) = AffineQuasiVector(inv(a), x)
+broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(/), x::AffineQuasiVector, a::Number) = AffineQuasiVector(inv(a), x)
+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)
+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}})
+    @_propagate_inbounds_meta
+    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
+end
+
+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
+
+
+
+
 include("operators.jl")
 include("bases/bases.jl")
 

src/bases/bases.jl

@@ -1,27 +1,146 @@
 abstract type Basis{T} <: LazyQuasiMatrix{T} end
+abstract type Weight{T} <: LazyQuasiVector{T} end
 
 
 const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}
 
-MemoryLayout(::Type{<:Basis}) = LazyLayout()
+struct WeightLayout <: MemoryLayout end
+struct BasisLayout <: MemoryLayout end
+struct AdjointBasisLayout <: MemoryLayout end
+
+MemoryLayout(::Type{<:Basis}) = BasisLayout()
+MemoryLayout(::Type{<:Weight}) = WeightLayout()
+
+adjointlayout(::Type, ::BasisLayout) = AdjointBasisLayout()
+transposelayout(::Type{<:Real}, ::BasisLayout) = AdjointBasisLayout()
+broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = BasisLayout()
+
+combine_mul_styles(::BasisLayout) = LazyQuasiArrayApplyStyle()
+combine_mul_styles(::AdjointBasisLayout) = LazyQuasiArrayApplyStyle()
+
 ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
 pinv(J::Basis) = apply(pinv,J)
 
 _multup(a::Tuple) = Mul(a...)
 _multup(a) = a
 
 
-==(A::Basis, B::Basis) = axes(A) ≠ axes(B) ||
-    throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
+function ==(A::Basis, B::Basis) 
+    axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
+    false
+end
+
 
 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()
 
-function copy(P::Ldiv{<:Any,<:Any,<:Basis,<:Basis})
-    A, B = P.A, P.B
-    A == B || throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
-    Eye(size(A,2))
+copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
+function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(-)}}) 
+    a,b = arguments(L.B)
+    (L.A\a)-(L.A\b)
+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
+end
+
+
+# expansion
+grid(P) = error("Overload Grid")
+function transform(L)
+    p = grid(L)
+    p,L[p,:]
+end
+
+function copy(L::Ldiv{BasisLayout,<:Any,<:Any,<:AbstractQuasiVector})
+    p,T = transform(L.A)
+    T \ L.B[p]
+end
+
+copy(L::Ldiv{BasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
+    copy(Ldiv{LazyLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
+
+function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
+    p,T = transform(L.A)
+     T \ L.B[p]
+end
+
+## materialize views
+
+# materialize(S::SubQuasiArray{<:Any,2,<:ApplyQuasiArray{<:Any,2,typeof(*),<:Tuple{<:Basis,<:Any}}}) =
+#     *(arguments(S)...)
+
+
+# Differentiation of sub-arrays 
+function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}})
+    A, B = M.args
+    P = parent(B)
+    (Derivative(axes(P,1))*P)[parentindices(B)...]
+end
+
+function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AffineQuasiVector,<:Any}}})
+    A, B = M.args
+    P = parent(B)
+    kr,jr = parentindices(B)
+    (Derivative(axes(P,1))*P*kr.A)[kr,jr]
+end
+
+function copy(L::Ldiv{BasisLayout,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(::BasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = ApplyLayout{typeof(*)}()
+sublayout(::BasisLayout, ::Type{<:Tuple{<:AffineQuasiVector,<:AbstractUnitRange}}) = BasisLayout()
+function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
+    A = parent(V)
+    _,jr = parentindices(V)
+    first(jr) ≥ 1 || throw(BoundsError())
+    P = _BandedMatrix(Ones{Int}(1,length(jr)), axes(A,2), first(jr)-1,1-first(jr))
+    A,P
 end
 
 
 include("splines.jl")
+
+
+

src/bases/splines.jl

@@ -1,9 +1,10 @@
-struct Spline{order,T} <: Basis{T}
-    points::Vector{T}
+struct Spline{order,T,P<:AbstractVector} <: Basis{T}
+    points::P
 end
+Spline{o,T}(pts::P) where {o,T,P} = Spline{o,T,P}(pts)
 
-const LinearSpline{T} = Spline{1,T}
-const HeavisideSpline{T} = Spline{0,T}
+const LinearSpline = Spline{1}
+const HeavisideSpline = Spline{0}
 
 Spline{o}(pts::AbstractVector{T}) where {o,T} = Spline{o,float(T)}(pts)
 
@@ -36,6 +37,9 @@ function getindex(B::HeavisideSpline{T}, x::Number, k::Int) where T
     return zero(T)
 end
 
+grid(L::LinearSpline) = L.points
+transform(L::LinearSpline{T}) where T = grid(L),Eye{T}(size(L,2))
+
 ## Sub-bases
 
 

src/operators.jl

@@ -154,12 +154,6 @@ axes(D::Derivative) = (D.axis, D.axis)
 copy(D::Derivative) = Derivative(copy(D.axis))
 
 
-function copy(M::QMul2{<:Derivative,<:SubQuasiArray})
-    A, B = M.args
-    P = parent(B)
-    (Derivative(axes(P,1))*P)[parentindices(B)...]
-end
-
 
 # struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}
 #     f::F