views almost working

Author: dlfivefifty

src/ContinuumArrays.jl

@@ -1,14 +1,14 @@
 module ContinuumArrays
 using IntervalSets, LinearAlgebra, LazyArrays, BandedMatrices
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
-                IndexStyle, IndexLinear
+                IndexStyle, IndexLinear, ==
 import Base.Broadcast: materialize
 import LazyArrays: Mul2
 import BandedMatrices: AbstractBandedLayout
 
 include("QuasiArrays/QuasiArrays.jl")
 using .QuasiArrays
-import .QuasiArrays: _length, checkindex, Adjoint, Transpose
+import .QuasiArrays: _length, checkindex, Adjoint, Transpose, slice, QSlice
 
 export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative
 
@@ -25,6 +25,7 @@ _length(::AbstractInterval) = ℵ₁
 
 
 checkindex(::Type{Bool}, inds::AbstractInterval, i::Real) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
+checkindex(::Type{Bool}, inds::AbstractInterval, i::QSlice) = i.axis ⊆ inds
 function checkindex(::Type{Bool}, inds::AbstractInterval, I::AbstractArray)
     @_inline_meta
     b = true

src/QuasiArrays/QuasiArrays.jl

@@ -6,8 +6,9 @@ import Base: getindex, size, axes, length, ==, isequal, iterate, CartesianIndice
 import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinline_meta,
                 DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type, _default_type,
                 _maybetail, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,
-                index_shape, to_shape, unsafe_length, @nloops, @ncall, Slice
-import Base: ViewIndex, Slice, ScalarIndex, RangeIndex
+                index_shape, to_shape, unsafe_length, @nloops, @ncall, Slice, unalias
+import Base: ViewIndex, Slice, ScalarIndex, RangeIndex, view, viewindexing, ensure_indexable, index_dimsum,
+                check_parent_index_match, reindex, _isdisjoint
 import Base: *, /, \, +, -, inv
 
 import Base.Broadcast: materialize

src/QuasiArrays/abstractquasiarray.jl

@@ -560,8 +560,9 @@ error_if_canonical_getindex(::IndexCartesian, A::AbstractQuasiArray{T,N}, ::Vara
 error_if_canonical_getindex(::IndexStyle, ::AbstractQuasiArray, ::Any...) = nothing
 
 ## Internal definitions
-_getindex(::IndexStyle, A::AbstractQuasiArray, I...) =
-    error("getindex for $(typeof(A)) with types $(typeof(I)) is not supported")
+function _getindex(::IndexStyle, A::AbstractQuasiArray, I...)
+    materialize(view(A, I...))
+end
 
 ## IndexLinear Scalar indexing: canonical method is one Int
 _getindex(::IndexLinear, A::AbstractQuasiArray, i::Real) = (@_propagate_inbounds_meta; getindex(A, i))

src/QuasiArrays/adjtrans.jl

@@ -218,3 +218,5 @@ end
 function adjoint(M::Mul)
     Mul(reverse(adjoint.(M.factors))...)
 end
+
+==(A::Adjoint, B::Adjoint) = parent(A) == parent(B)

src/QuasiArrays/indices.jl

@@ -144,3 +144,32 @@ to_index(I::AbstractQuasiArray{<:Union{AbstractArray, Colon}}) =
     throw(ArgumentError("invalid index: $I of type $(typeof(I))"))
 
 LinearIndices(A::AbstractQuasiArray) = LinearIndices(axes(A))
+
+
+
+"""
+   QSlice(indices)
+
+Represent an axis as a quasi-vector that returns itself.
+"""
+struct QSlice{T,AX} <: AbstractQuasiVector{T}
+    axis::AX
+end
+QSlice(axis) = QSlice{eltype(axis),typeof(axis)}(axis)
+QSlice(S::QSlice) = S
+axes(S::QSlice) = (S,)
+unsafe_indices(S::QSlice) = (S,)
+axes1(S::QSlice) = S
+axes(S::QSlice{<:OneTo}) = (S.axis,)
+unsafe_indices(S::QSlice{<:OneTo}) = (S.axis,)
+axes1(S::QSlice{<:OneTo}) = S.axis
+
+first(S::QSlice) = first(S.axis)
+last(S::QSlice) = last(S.axis)
+size(S::QSlice) = (length(S.axis),)
+length(S::QSlice) = length(S.axis)
+unsafe_length(S::QSlice) = unsafe_length(S.axis)
+getindex(S::QSlice, i::Real) = (@_inline_meta; @boundscheck checkbounds(S, i); i)
+getindex(S::QSlice, i::AbstractVector{<:Real}) = (@_inline_meta; @boundscheck checkbounds(S, i); i)
+show(io::IO, r::QSlice) = print(io, "QSlice(", r.axis, ")")
+iterate(S::QSlice, s...) = iterate(S.axis, s...)

src/QuasiArrays/matmul.jl

@@ -28,5 +28,5 @@ QuasiMatMulQuasiMat{styleA, styleB, T, V} = QuasiArrayMulQuasiArray{styleA, styl
 *(A::AbstractArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
 inv(A::AbstractQuasiArray) = materialize(Inv(A))
 
-*(A::AbstractQuasiArray, B::Mul) = materialize(Mul(A, B.factors...))
-*(A::Mul, B::AbstractQuasiArray) = materialize(Mul(A.factors..., B))
+*(A::AbstractQuasiArray, B::Mul) = Mul(A, B.factors...)
+*(A::Mul, B::AbstractQuasiArray) = Mul(A.factors..., B)

src/QuasiArrays/multidimensional.jl

@@ -7,8 +7,11 @@ end
 @inline to_indices(A, inds, I::Tuple{Colon, Vararg{Any}}) =
     (uncolon(inds, I), to_indices(A, _maybetail(inds), tail(I))...)
 
-uncolon(inds::Tuple{},    I::Tuple{Colon, Vararg{Any}}) = Slice(OneTo(1))
-uncolon(inds::Tuple,      I::Tuple{Colon, Vararg{Any}}) = Slice(inds[1])
+slice(d::AbstractVector) = Slice(d)
+slice(d) = QSlice(d)
+
+uncolon(inds::Tuple{},    I::Tuple{Colon, Vararg{Any}}) = slice(OneTo(1))
+uncolon(inds::Tuple,      I::Tuple{Colon, Vararg{Any}}) = slice(inds[1])
 
 
 _maybe_reshape(::IndexLinear, A::AbstractQuasiArray, I...) = A

src/QuasiArrays/subarray.jl

@@ -69,7 +69,9 @@ function view(A::AbstractQuasiArray, I::Vararg{Any,N}) where {N}
     unsafe_view(_maybe_reshape_parent(A, index_ndims(J...)), J...)
 end
 
-function unsafe_view(A::AbstractQuasiArray, I::Vararg{ViewIndex,N}) where {N}
+const QViewIndex = Union{ViewIndex,AbstractQuasiArray}
+
+function unsafe_view(A::AbstractQuasiArray, I::Vararg{QViewIndex,N}) where {N}
     @_inline_meta
     SubQuasiArray(A, I)
 end

src/bases/splines.jl

@@ -9,6 +9,7 @@ const HeavisideSpline{T} = Spline{0,T}
 Spline{o}(pts::AbstractVector{T}) where {o,T} = Spline{o,float(T)}(pts)
 
 axes(B::Spline{o}) where o = (first(B.points)..last(B.points), Base.OneTo(length(B.points)+o-1))
+==(A::Spline{o}, B::Spline{o}) where o = A.points == B.points
 
 function getindex(B::LinearSpline{T}, x::Real, k::Int) where T
     x ∈ axes(B,1) && 1 ≤ k ≤ size(B,2)|| throw(BoundsError())

src/operators.jl

@@ -11,6 +11,8 @@ DiracDelta(axis) = DiracDelta(zero(float(eltype(axis))), axis)
 axes(δ::DiracDelta) = (δ.axis,)
 IndexStyle(::Type{<:DiracDelta}) = IndexLinear()
 
+==(a::DiracDelta, b::DiracDelta) = a.axis == b.axis && a.x == b.x
+
 function getindex(δ::DiracDelta{T}, x::Real) where T
     x ∈ δ.axis || throw(BoundsError())
     x == δ.x ? inv(zero(T)) : zero(T)
@@ -37,3 +39,4 @@ Derivative{T}(axis::A) where {T,A} = Derivative{T,A}(axis)
 Derivative(axis) = Derivative{Float64}(axis)
 
 axes(D::Derivative) = (D.axis, D.axis)
+==(a::Derivative, b::Derivative) = a.axis == b.axis

test/runtests.jl

@@ -1,5 +1,5 @@
-using ContinuumArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, Test
-    import ContinuumArrays: ℵ₁
+using ContinuumArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test
+    import ContinuumArrays: ℵ₁, materialize
 
 @testset "DiracDelta" begin
     δ = DiracDelta(-1..3)
@@ -80,23 +80,27 @@ end
     @test fp[2.2] ≈ 2
 end
 
+@testset "Weak Laplacian" begin
+    H = HeavisideSpline(0:2)
+    L = LinearSpline(0:2)
 
+    D = Derivative(axes(L,1))
+    M = materialize(Mul(D',D,L))
+    DL = D*L
+    @test M.factors == tuple(D', (D*L).factors...)
 
+    @test materialize(Mul(L', D', D, L)) == materialize(L'D'*D*L) ==
+        [1.0 -1 0; -1.0 2.0 -1.0; 0.0 -1.0 1.0]
 
-L = LinearSpline([1,2,3])
-f = L*[1,2,4]
-D = Derivative(axes(L,1))
-
-H = HeavisideSpline([1,2,3])
-
-D*L
-M = (L'D')*(D*L)
-
-
-*(M.factors...)
-
-
+    @test materialize(Mul(L', D', D, L)) isa BandedMatrix
+    @test materialize(L'D'*D*L) isa BandedMatrix
 
-LazyArrays.MemoryLayout(Diagonal(randn(5)))
+    @test bandwidths(materialize(L'D'*D*L)) == (1,1)
+end
 
-M.factors[2] * M.factors[3]
+@testset "Views" begin
+    L = LinearSpline(0:2)
+    @test view(L,0.1,1)[1] == L[0.1,1]
+end
+L = LinearSpline(0:2)
+L[:,1]