Merge pull request #9 from JuliaApproximation/dl/collocation

Support LazyArrays v0.6

Author: dlfivefifty

examples/collocation.jl

@@ -0,0 +1,13 @@
+using ContinuumArrays, FillArrays, InfiniteArrays
+import ContinuumArrays.QuasiArrays: Inclusion, QuasiDiagonal
+
+P = Legendre()
+X = QuasiDiagonal(Inclusion(-1..1))
+
+@test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
+
+axes(X)
+J = pinv(P)*X*P
+
+J - I
+Vcat(Hcat(1, Zeros(1,∞)), J)

src/ContinuumArrays.jl

@@ -3,15 +3,15 @@ using IntervalSets, LinearAlgebra, LazyArrays, BandedMatrices, InfiniteArrays, D
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear, ==, OneTo, tail
 import Base.Broadcast: materialize
-import LazyArrays: Mul2, MemoryLayout
+import LazyArrays: Mul2, MemoryLayout, Applied
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 
 include("QuasiArrays/QuasiArrays.jl")
 using .QuasiArrays
-import .QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, slice, QSlice, SubQuasiArray,
+import .QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, slice, Inclusion, SubQuasiArray,
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
-                    _PInvQuasiMatrix
+                    ApplyQuasiArray, ApplyQuasiMatrix
 
 export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, JacobiWeight, Jacobi, Legendre,
             fullmaterialize
@@ -24,12 +24,14 @@ struct AlephInfinity{N} <: Integer end
 const ℵ₁ = AlephInfinity{1}()
 
 
+const QMul2{A,B} = Mul{<:Any, <:Tuple{A,B}}
+
 cardinality(::AbstractInterval) = ℵ₁
 *(ℵ::AlephInfinity) = ℵ
 
 
 checkindex(::Type{Bool}, inds::AbstractInterval, i::Real) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
-checkindex(::Type{Bool}, inds::AbstractInterval, i::QSlice) = i.axis ⊆ inds
+checkindex(::Type{Bool}, inds::AbstractInterval, i::Inclusion) = i.axis ⊆ inds
 function checkindex(::Type{Bool}, inds::AbstractInterval, I::AbstractArray)
     @_inline_meta
     b = true
@@ -41,7 +43,7 @@ end
 
 
 # we represent as a Mul with a banded matrix
-function materialize(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:QSlice,<:AbstractUnitRange}})
+function materialize(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
     A = parent(V)
     _,jr = parentindices(V)
     first(jr) ≥ 1 || throw(BoundsError())
@@ -51,42 +53,6 @@ end
 
 
 
-
-most(a) = reverse(tail(reverse(a)))
-
-MulQuasiOrArray = Union{MulArray,MulQuasiArray}
-
-_factors(M::MulQuasiOrArray) = M.mul.factors
-_factors(M) = (M,)
-
-function fullmaterialize(M::MulQuasiOrArray)
-    M_mat = materialize(M.mul)
-    typeof(M_mat) <: MulQuasiOrArray || return M_mat
-    typeof(M_mat) == typeof(M) || return(fullmaterialize(M_mat))
-
-    ABC = M_mat.mul.factors
-    length(ABC) ≤ 2 && return M_mat
-
-    AB = most(ABC)
-    Mhead = fullmaterialize(Mul(AB...))
-
-    typeof(_factors(Mhead)) == typeof(AB) ||
-        return fullmaterialize(Mul(_factors(Mhead)..., last(ABC)))
-
-    BC = tail(ABC)
-    Mtail =  fullmaterialize(Mul(BC...))
-    typeof(_factors(Mtail)) == typeof(BC) ||
-        return fullmaterialize(Mul(first(ABC), _factors(Mtail)...))
-
-    first(ABC) * Mtail
-end
-
-fullmaterialize(M::Mul) = fullmaterialize(materialize(M))
-fullmaterialize(M) = M
-
-materialize(M::Mul{<:Tuple,<:Tuple{Vararg{<:Union{Adjoint,QuasiAdjoint,QuasiDiagonal}}}}) =
-    materialize(Mul(reverse(adjoint.(M.factors))))'
-
 include("operators.jl")
 include("bases/bases.jl")
 

src/QuasiArrays/QuasiArrays.jl

@@ -2,9 +2,9 @@ module QuasiArrays
 using Base, LinearAlgebra, LazyArrays
 import Base: getindex, size, axes, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
                 Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, zero,
-                map, eachindex, eltype, first, last, firstindex, lastindex
+                map, eachindex, eltype, first, last, firstindex, lastindex, in
 import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinline_meta,
-                DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type, _default_type,
+                DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type,
                 _maybetail, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,
                 index_shape, to_shape, unsafe_length, @nloops, @ncall, Slice, unalias
 import Base: ViewIndex, Slice, ScalarIndex, RangeIndex, view, viewindexing, ensure_indexable, index_dimsum,
@@ -21,10 +21,11 @@ import Base: Array, Matrix, Vector
 import Base.Broadcast: materialize
 
 import LinearAlgebra: transpose, adjoint, checkeltype_adjoint, checkeltype_transpose, Diagonal,
-                        AbstractTriangular
+                        AbstractTriangular, pinv, inv
 
-import LazyArrays: MemoryLayout, UnknownLayout, Mul2, _materialize, MulLayout, ⋆, rmaterialize,
-                    _rmaterialize, _lmaterialize, flatten, _flatten, AbstractPInv
+import LazyArrays: MemoryLayout, UnknownLayout, Mul2, _materialize, MulLayout, ⋆,
+                    _lmaterialize, InvOrPInv, ApplyStyle,
+                    LayoutApplyStyle, Applied
 
 export AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector, materialize
 
@@ -49,4 +50,8 @@ include("abstractquasiarraymath.jl")
 include("quasiadjtrans.jl")
 include("quasidiagonal.jl")
 
+
+materialize(M::Applied{<:Any,typeof(*),<:Tuple{Vararg{<:Union{Adjoint,QuasiAdjoint,QuasiDiagonal}}}}) =
+    materialize(Mul(reverse(adjoint.(M.args))...))'
+
 end

src/QuasiArrays/indices.jl

@@ -148,32 +148,51 @@ LinearIndices(A::AbstractQuasiArray) = LinearIndices(axes(A))
 
 
 """
-   QSlice(indices)
+   Inclusion(domain)
 
-Represent an axis as a quasi-vector that returns itself.
+Represents the inclusion operator of a domain (that is, a type that overrides in)
+as an AbstractQuasiVector. That is, if `v = Inclusion(domain)`, then
+`v[x] == x` if `x in domain`, otherwise it throws a `DomainError`.
+
+Inclusions are useful for turning domains into axes. They also serve the same
+role as `Slice` does for offset arrays.
 """
-struct QSlice{T,AX} <: AbstractQuasiVector{T}
-    axis::AX
+struct Inclusion{T,AX} <: AbstractQuasiVector{T}
+    domain::AX
+end
+Inclusion(domain) = Inclusion{eltype(domain),typeof(domain)}(domain)
+Inclusion(S::Inclusion) = S
+==(A::Inclusion, B::Inclusion) = A.domain == B.domain
+axes(S::Inclusion) = (S,)
+unsafe_indices(S::Inclusion) = (S,)
+axes1(S::Inclusion) = S
+axes(S::Inclusion{<:Any,<:OneTo}) = (S.domain,)
+unsafe_indices(S::Inclusion{<:Any,<:OneTo}) = (S.domain,)
+axes1(S::Inclusion{<:Any,<:OneTo}) = S.domain
+
+first(S::Inclusion) = first(S.domain)
+last(S::Inclusion) = last(S.domain)
+size(S::Inclusion) = (length(S.domain),)
+length(S::Inclusion) = length(S.domain)
+unsafe_length(S::Inclusion) = unsafe_length(S.domain)
+cardinality(S::Inclusion) = cardinality(S.domain)
+getindex(S::Inclusion{T}, i::Real) where T =
+    (@_inline_meta; @boundscheck checkbounds(S, i); convert(T,i))
+getindex(S::Inclusion{T}, i::AbstractVector{<:Real}) where T =
+    (@_inline_meta; @boundscheck checkbounds(S, i); convert(AbstractVector{T},i))
+show(io::IO, r::Inclusion) = print(io, "Inclusion(", r.domain, ")")
+iterate(S::Inclusion, s...) = iterate(S.domain, s...)
+
+in(x, S::Inclusion) = x in S.domain
+
+checkindex(::Type{Bool}, inds::Inclusion, i::Real) = i ∈ inds.domain
+checkindex(::Type{Bool}, inds::Inclusion, ::Colon) = true
+checkindex(::Type{Bool}, inds::Inclusion, ::Inclusion) = true
+function checkindex(::Type{Bool}, inds::Inclusion, I::AbstractArray)
+    @_inline_meta
+    b = true
+    for i in I
+        b &= checkindex(Bool, inds, i)
+    end
+    b
 end
-QSlice(axis) = QSlice{eltype(axis),typeof(axis)}(axis)
-QSlice(S::QSlice) = S
-==(A::QSlice, B::QSlice) = A.axis == B.axis
-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)
-cardinality(S::QSlice) = cardinality(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...)
-
-checkindex(::Type{Bool}, inds::QSlice, i::Real) = i ∈ inds.axis