Support Infinities v0.1 (#83)

* Support Infinities v0.1

* DiracDelta fixes

* avoid stack overflow in derivative * basis view materialize

* Update bases.jl

* Fix B'D' simplify

* Check axes match

* better support for vector sub-bases

* Update runtests.jl

* v0.6.4

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.6.3"
+version = "0.6.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,10 +18,10 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 [compat]
 ArrayLayouts = "0.5, 0.6"
 BandedMatrices = "0.16"
-BlockArrays = "0.14, 0.15"
+BlockArrays = "0.15.1"
 FillArrays = "0.11"
 InfiniteArrays = "0.9, 0.10"
-Infinities = "0.0.1, 0.0.2"
+Infinities = "0.0.1, 0.0.2, 0.1"
 IntervalSets = "0.5"
 LazyArrays = "0.20, 0.21"
 QuasiArrays = "0.4.1"

src/ContinuumArrays.jl

@@ -3,7 +3,7 @@ using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, Quasi
 import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==, ^,
                 IndexStyle, IndexLinear, ==, OneTo, _maybetail, tail, similar, copyto!, copy, diff,
                 first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum,
-                getproperty, isone, iszero, zero, abs, <, ≤, >, ≥, string, summary, to_indices
+                getproperty, isone, iszero, zero, abs, <, ≤, >, ≥, string, summary, to_indices, view
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, most, combine_mul_styles, AbstractArrayApplyStyle,
                         adjointlayout, arguments, _mul_arguments, call, broadcastlayout, layout_getindex, UnknownLayout,

src/bases/bases.jl

@@ -34,7 +34,7 @@ broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::MappedBasisLayouts) = Mappe
 
 # A sub of a weight is still a weight
 sublayout(::WeightLayout, _) = WeightLayout()
-sublayout(::AbstractBasisLayout, ::Type{<:Tuple{Map,AbstractUnitRange}}) = MappedBasisLayout()
+sublayout(::AbstractBasisLayout, ::Type{<:Tuple{Map,AbstractVector}}) = MappedBasisLayout()
 
 
 ## Weighted basis interface
@@ -290,23 +290,37 @@ end
 
 
 # Differentiation of sub-arrays
-@simplify function *(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}})
+
+# avoid stack overflow from unmaterialize Derivative() * parent()
+_der_sub(DP, inds...) = DP[inds...]
+_der_sub(DP::ApplyQuasiMatrix{T,typeof(*),<:Tuple{Derivative,Any}}, kr, jr) where T = ApplyQuasiMatrix{T}(*, DP.args[1], view(DP.args[2], kr, jr))
+
+# need to customise simplifiable so can't use @simplify
+simplifiable(::typeof(*), A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}})= simplifiable(*, A, parent(B))
+simplifiable(::typeof(*), Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = simplifiable(*, Bc', Ac')
+function mul(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}})
+    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
     P = parent(B)
-    (Derivative(axes(P,1))*P)[parentindices(B)...]
+    _der_sub(Derivative(axes(P,1))*P, parentindices(B)...)
 end
+mul(Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = mul(Bc', Ac')'
 
-@simplify function *(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
+simplifiable(::typeof(*), A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) = simplifiable(*, A, parent(B))
+simplifiable(::typeof(*), Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = simplifiable(*, Bc', Ac')
+function mul(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
+    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
     P = parent(B)
     kr,jr = parentindices(B)
     (Derivative(axes(P,1))*P*kr.A)[kr,jr]
 end
+mul(Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = mul(Bc', Ac')'
 
 # we represent as a Mul with a banded matrix
-sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
-sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()
-sublayout(::WeightedBasisLayouts, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedWeightedBasisLayout()
-sublayout(::WeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubWeightedBasisLayout()
-sublayout(::MappedWeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = MappedWeightedBasisLayout()
+sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) = SubBasisLayout()
+sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractVector}}) = MappedBasisLayout()
+sublayout(::WeightedBasisLayouts, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractVector}}) = MappedWeightedBasisLayout()
+sublayout(::WeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) = SubWeightedBasisLayout()
+sublayout(::MappedWeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) = MappedWeightedBasisLayout()
 
 @inline sub_materialize(::AbstractBasisLayout, V::AbstractQuasiArray) = V
 @inline sub_materialize(::AbstractBasisLayout, V::AbstractArray) = V
@@ -348,6 +362,14 @@ function arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2,<:Any,<:Tu
     A,P
 end
 
+function arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractVector}})
+    A = parent(V)
+    _,jr = parentindices(V)
+    first(jr) ≥ 1 || throw(BoundsError())
+    P = Eye{Int}((axes(A,2),))[:,jr]
+    A,P
+end
+
 ####
 # sum
 ####

src/basisconcat.jl

@@ -15,9 +15,10 @@ end
 
 
 """
-    PiecewiseBasis
+    PiecewiseBasis(args...)
 
-is an analogue of `Basis` that takes the union of the first axis.
+is an analogue of `Basis` that takes the union of the first axis,
+and the second axis is a blocked concatenatation of args.
 If there is overlap, it uses the first in order.
 """
 struct PiecewiseBasis{T, Args} <: AbstractConcatBasis{T}
@@ -31,13 +32,12 @@ PiecewiseBasis(args::AbstractVector) = PiecewiseBasis{eltype(eltype(args))}(args
 
 concatbasis(::PiecewiseBasis, args...) = PiecewiseBasis(args...)
 
-_blockedrange(a::Tuple) = blockedrange(SVector(a...))
-_blockedrange(a::AbstractVector) = blockedrange(a)
 
-axes(A::PiecewiseBasis) = (union(axes.(A.args,1)...), _blockedrange(size.(A.args,2)))
+axes(A::PiecewiseBasis) = (union(axes.(A.args,1)...), blockedrange(size.(A.args,2)))
 
 ==(A::PiecewiseBasis, B::PiecewiseBasis) = all(A.args .== B.args)
 
+
 function QuasiArrays._getindex(::Type{IND}, A::PiecewiseBasis{T}, (x,j)::IND) where {IND,T}
     Jj = findblockindex(axes(A,2), j)
     J = Int(block(Jj))
@@ -67,7 +67,7 @@ VcatBasis(args::Vararg{Any,N}) where N = VcatBasis{SVector{N,mapreduce(eltype,pr
 
 concatbasis(::VcatBasis, args...) = VcatBasis(args...)
 
-axes(A::VcatBasis{<:Any,<:Tuple}) = (axes(A.args[1],1), _blockedrange(size.(A.args,2)))
+axes(A::VcatBasis{<:Any,<:Tuple}) = (axes(A.args[1],1), blockedrange(size.(A.args,2)))
 
 ==(A::VcatBasis, B::VcatBasis) = all(A.args .== B.args)
 
@@ -99,12 +99,12 @@ HvcatBasis(n::Int, args::Vararg{Any,N}) where N = HvcatBasis{Matrix{mapreduce(el
 
 concatbasis(S::HvcatBasis, args...) = HvcatBasis(S.n, args...)
 
-axes(A::HvcatBasis{<:Any,<:Tuple}) = (axes(A.args[1],1), _blockedrange(size.(A.args,2)))
+axes(A::HvcatBasis{<:Any,<:Tuple}) = (axes(A.args[1],1), blockedrange(size.(A.args,2)))
 
 ==(A::HvcatBasis, B::HvcatBasis) = all(A.args .== B.args)
 
 function QuasiArrays._getindex(::Type{IND}, A::HvcatBasis{T}, (x,j)::IND) where {T,IND}
     Jj = findblockindex(axes(A,2), j)
     J = Int(block(Jj))
     hvcat(A.n, zeros(J-1)..., A.args[J][x, blockindex(Jj)], zeros(length(A.args)-J)...)::T
-end
+end

src/operators.jl

@@ -67,9 +67,12 @@ end
 
 DiracDelta{T}(x, axis::A) where {T,A<:AbstractQuasiVector} = DiracDelta{T,A}(x, axis)
 DiracDelta{T}(x, domain) where T = DiracDelta{T}(x, Inclusion(domain))
-DiracDelta(x, domain) = DiracDelta{eltype(x)}(x, Inclusion(domain))
+DiracDelta(x, domain) = DiracDelta{float(eltype(x))}(x, Inclusion(domain))
 DiracDelta(axis::AbstractQuasiVector) = DiracDelta(zero(float(eltype(axis))), axis)
-DiracDelta(domain) = DiracDelta(Inclusion(domain))
+DiracDelta(x) = DiracDelta(x, Inclusion(x))
+DiracDelta{T}() where T = DiracDelta(zero(T))
+DiracDelta() = DiracDelta{Float64}()
+
 
 axes(δ::DiracDelta) = (δ.axis,)
 IndexStyle(::Type{<:DiracDelta}) = IndexLinear()
@@ -78,7 +81,7 @@ IndexStyle(::Type{<:DiracDelta}) = IndexLinear()
 
 function getindex(δ::DiracDelta{T}, x::Number) where T
     x ∈ δ.axis || throw(BoundsError())
-    x == δ.x ? inv(zero(T)) : zero(T)
+    convert(T, x == δ.x ? inv(zero(T)) : zero(T))::T
 end
 
 
@@ -111,6 +114,12 @@ end
 
 ^(D::Derivative, k::Integer) = ApplyQuasiArray(^, D, k)
 
+
+function view(D::Derivative, kr::Inclusion, jr::Inclusion)
+    @boundscheck axes(D,1) == kr == jr || throw(BoundsError(D,(kr,jr)))
+    D
+end
+
 # struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}
 #     f::F
 #     axis::A

test/runtests.jl

@@ -31,7 +31,16 @@ end
 include("test_maps.jl")
 
 @testset "DiracDelta" begin
-    δ = DiracDelta(-1..3)
+    δ = DiracDelta()
+    @test δ == DiracDelta(0, Inclusion(0))
+    @test axes(δ) ≡ (axes(δ,1),) ≡ (Inclusion(0.0),)
+    @test size(δ) ≡ (length(δ),) ≡ (1,)
+    @test_throws BoundsError δ[1.1]
+    @test δ[0.0] ≡ Inf
+    @test Base.IndexStyle(δ) ≡ Base.IndexLinear()
+
+    δ = DiracDelta(0, -1..3)
+    @test δ == DiracDelta{Float64}(0, -1..3)
     @test axes(δ) ≡ (axes(δ,1),) ≡ (Inclusion(-1..3),)
     @test size(δ) ≡ (length(δ),) ≡ (ℵ₁,)
     @test δ[1.1] ≡ 0.0
@@ -203,6 +212,33 @@ end
         @test length(fp.args) == 2
         @test fp[1.1] ≈ 1
         @test fp[2.2] ≈ 2
+
+
+        @testset "View derivatives" begin
+            L = LinearSpline(1:5)
+            H = HeavisideSpline(1:5)
+            x = axes(L,1)
+            D = Derivative(x)
+
+            @test view(D, x, x) ≡ D
+            @test_throws BoundsError view(D, Inclusion(0..1), x)
+            jr = [1,3,4]
+            @test H \ (D*L[:,jr]) == H\(L[:,jr]'D')' == (H \ (D*L))[:,jr]
+
+            @test D*L[:,jr]*[1,2,3] == D * L * [1,0,2,3,0]
+            @test L[:,jr]'D'D*L[:,jr] == (L'D'D*L)[jr,jr]
+
+            a = affine(0..1, 1..5)
+            D̃ = Derivative(axes(a,1))
+            @test H[a,:] \ (D̃ * L[a,:]) == H[a,:] \ (L[a,:]'D̃')' ==  4*(H\(D*L))
+            @test_throws DimensionMismatch D * L[a,:]
+            @test_throws DimensionMismatch D̃ * L[:,jr]
+
+            @test ContinuumArrays.simplifiable(*, D, L[:,jr]) isa Val{true}
+            @test ContinuumArrays.simplifiable(*, L[:,jr]', D') isa Val{true}
+            @test ContinuumArrays.simplifiable(*, D̃, L[a,jr]) isa Val{true}
+            @test ContinuumArrays.simplifiable(*, L[a,jr]', D̃') isa Val{true}
+        end
     end
 
     @testset "Weak Laplacian" begin