Add fullmaterialize

Author: dlfivefifty

src/ContinuumArrays.jl

@@ -1,18 +1,20 @@
 module ContinuumArrays
 using IntervalSets, LinearAlgebra, LazyArrays, BandedMatrices, InfiniteArrays, DomainSets
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
-                IndexStyle, IndexLinear, ==, OneTo
+                IndexStyle, IndexLinear, ==, OneTo, tail
 import Base.Broadcast: materialize
-import LazyArrays: Mul2
+import LazyArrays: Mul2, MemoryLayout
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 
 include("QuasiArrays/QuasiArrays.jl")
 using .QuasiArrays
 import .QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, slice, QSlice, SubQuasiArray,
-                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat
+                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
+                    _PInvQuasiMatrix
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, JacobiWeight, Jacobi, Legendre
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, JacobiWeight, Jacobi, Legendre,
+            fullmaterialize
 
 ####
 # Interval indexing support
@@ -48,6 +50,43 @@ function materialize(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:QSlice,<:AbstractU
 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

@@ -24,7 +24,7 @@ import LinearAlgebra: transpose, adjoint, checkeltype_adjoint, checkeltype_trans
                         AbstractTriangular
 
 import LazyArrays: MemoryLayout, UnknownLayout, Mul2, _materialize, MulLayout, ⋆, rmaterialize,
-                    _rmaterialize, _lmaterialize, flatten, _flatten
+                    _rmaterialize, _lmaterialize, flatten, _flatten, AbstractPInv
 
 export AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector, materialize
 

src/QuasiArrays/matmul.jl

@@ -35,15 +35,14 @@ end
 QuasiMatMulMat{styleA, styleB, T, V} = QuasiArrayMulArray{styleA, styleB, 2, 2, T, V}
 QuasiMatMulQuasiMat{styleA, styleB, T, V} = QuasiArrayMulQuasiArray{styleA, styleB, 2, 2, T, V}
 
-*(A::AbstractQuasiArray, B::AbstractQuasiArray, C::AbstractQuasiArray, D::AbstractQuasiArray...) = materialize(Mul(A,B,C,D...))
-*(A::AbstractQuasiArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
-*(A::AbstractQuasiArray, B::AbstractArray) = materialize(Mul(A,B))
-*(A::AbstractArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
+*(A::AbstractQuasiArray, B...) = materialize(Mul(A,B...))
+*(A::AbstractQuasiArray, B::AbstractQuasiArray, C...) = materialize(Mul(A,B,C...))
+*(A::AbstractArray, B::AbstractQuasiArray, C...) = materialize(Mul(A,B,C...))
 pinv(A::AbstractQuasiArray) = materialize(PInv(A))
 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, C...) = materialize(Mul(A, B.factors..., C...))
+*(A::Mul, B::AbstractQuasiArray, C...) = materialize(Mul(A.factors..., B, C...))
 
 
 ####
@@ -99,6 +98,31 @@ transpose(A::MulQuasiArray) = MulQuasiArray(reverse(transpose.(A.mul.factors))..
 MemoryLayout(M::MulQuasiArray) = MulLayout(MemoryLayout.(M.mul.factors))
 
 
+## PInvQuasiMatrix
+
+function _PInvQuasiMatrix end
+
+struct PInvQuasiMatrix{T, PINV<:AbstractPInv} <: AbstractQuasiMatrix{T}
+    pinv::PINV
+    global _PInvQuasiMatrix(pinv::PINV) where PINV = new{eltype(pinv), PINV}(pinv)
+end
+
+const InvQuasiMatrix{T, INV<:Inv} = PInvQuasiMatrix{T,INV}
+
+PInvQuasiMatrix(M) = _PInvQuasiMatrix(PInv(M))
+InvQuasiMatrix(M) = _PInvQuasiMatrix(Inv(M))
+
+axes(A::PInvQuasiMatrix) = axes(A.pinv)
+size(A::PInvQuasiMatrix) = map(length, axes(A))
+
+@propagate_inbounds getindex(A::PInvQuasiMatrix{T}, k::Int, j::Int) where T =
+    (A.pinv*[Zeros(j-1); one(T); Zeros(size(A,2) - j)])[k]
+
+MemoryLayout(M::PInvQuasiMatrix) = MemoryLayout(M.pinv)
+
+*(A::PInvQuasiMatrix, B::AbstractQuasiMatrix, C...) = *(A.pinv, B, C...)
+*(A::PInvQuasiMatrix, B::MulQuasiArray, C...) = *(A.pinv, B.mul, C...)
+
 
 ####
 # Matrix * Array
@@ -130,7 +154,7 @@ end
 
 function _materialize(M::Mul2{<:Any,<:Any,<:MulQuasiArray,<:AbstractQuasiArray}, _)
     As, B = M.factors
-    ⋆(As.mul.factors..., B)
+    rmaterialize(Mul(As.mul.factors..., B))
 end
 
 function _materialize(M::Mul2{<:Any,<:Any,<:AbstractQuasiArray,<:MulQuasiArray}, _)
@@ -141,7 +165,7 @@ end
 # A MulQuasiArray can't be materialized further left-to-right, so we do right-to-left
 function _materialize(M::Mul2{<:Any,<:Any,<:MulQuasiArray,<:AbstractArray}, _)
     As, B = M.factors
-    ⋆(As.mul.factors..., B)
+    rmaterialize(Mul(As.mul.factors..., B))
 end
 
 function _lmaterialize(A::MulQuasiArray, B, C...)

src/bases/bases.jl

@@ -1,2 +1,21 @@
+abstract type Basis{T} <: AbstractQuasiMatrix{T} end
+struct BasisLayout <: MemoryLayout end
+
+MemoryLayout(::Basis) = BasisLayout()
+
+pinv(J::Basis) = materialize(PInv(J))
+materialize(P::PInv{BasisLayout}) = _PInvQuasiMatrix(P)
+
+==(A::Basis, B::Basis) = axes(A) ≠ axes(B) ||
+    throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
+
+function materialize(P::Ldiv{BasisLayout,BasisLayout})
+    Ai, B = P.factors
+    A = parent(Ai)
+    axes(A) == axes(B) || throw(DimensionMismatch("axes of bases must match"))
+    A == B || throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
+    Eye(size(A,2))
+end
+
 include("splines.jl")
 include("jacobi.jl")

src/bases/jacobi.jl

@@ -9,13 +9,13 @@ function getindex(w::JacobiWeight, x::Real)
     (1-x)^w.a * (1+x)^w.b
 end
 
-abstract type AbstractJacobi{T} <: AbstractQuasiMatrix{T} end
-
-pinv(J::AbstractJacobi) = PInv(J)
+abstract type AbstractJacobi{T} <: Basis{T} end
 
 struct Legendre{T} <: AbstractJacobi{T} end
 Legendre() = Legendre{Float64}()
 
+==(::Legendre, ::Legendre) = true
+
 struct Jacobi{T} <: AbstractJacobi{T}
     b::T
     a::T

src/bases/splines.jl

@@ -1,5 +1,5 @@
 
-struct Spline{order,T} <: AbstractQuasiMatrix{T}
+struct Spline{order,T} <: Basis{T}
     points::Vector{T}
 end
 

src/operators.jl

@@ -48,3 +48,7 @@ function materialize(M::Mul2{<:Any,<:Any,<:Derivative,<:SubQuasiArray})
     P = parent(B)
     (Derivative(axes(P,1))*P)[parentindices(B)...]
 end
+
+
+materialize(M::Mul2{<:Any,<:Any,<:Derivative}) = MulQuasiArray(M)
+materialize(M::Mul2{<:Any,<:Any,<:Any,<:Derivative}) = MulQuasiArray(M)