Support vector/matrix bases (#77)

* Support vector/matrix bases

* vector/matrix expansions

* Support operators

* Update Project.toml

* Fourier tests

* start Fill

* Support matrix transforms

* start matrix Fourier

* Support Fourier matrix transform

* Update test_chebyshev.jl

* Expansion -> ExpansionLayout

* Update lanczos.jl

* Remove WeightedOrthogonalPolynomial, no longer needed broadcasted code

* improve Weighted(::Jacobi)

* WeightedBassis takes Layout

* Work on SetindexInterlace evaluation

* improve adaptive block constructor

* Fix mapped broadcastbasis

* tests pass

* Increase coverage

* Update test_chebyshev.jl

* Update test_legendre.jl

* increase coverage

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.4.10"
+version = "0.5.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -25,31 +25,32 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "0.7.5"
+ArrayLayouts = "0.7.7"
 BandedMatrices = "0.16.11"
-BlockArrays = "0.16.6"
+BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.11"
-ContinuumArrays = "0.9.4"
+ContinuumArrays = "0.10"
 DomainSets = "0.5.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
-FastTransforms = "0.12"
+FastTransforms = "0.13"
 FillArrays = "0.12"
 HypergeometricFunctions = "0.3.4"
-InfiniteArrays = "0.12"
-InfiniteLinearAlgebra = "0.6"
+InfiniteArrays = "0.12.3"
+InfiniteLinearAlgebra = "0.6.5"
 IntervalSets = "0.5"
 LazyArrays = "0.22"
 LazyBandedMatrices = "0.7"
-QuasiArrays = "0.8"
+QuasiArrays = "0.9"
 SpecialFunctions = "1.0"
 julia = "1.6"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices"]
+test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays"]

src/ClassicalOrthogonalPolynomials.jl

@@ -8,19 +8,19 @@ using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, Bloc
     LazyBandedMatrices, HypergeometricFunctions
 
 import Base: @_inline_meta, axes, getindex, unsafe_getindex, convert, prod, *, /, \, +, -,
-                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
+                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, setindex,
                 first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum, cumsum,
                 to_indices, _maybetail, tail, getproperty, inv, show, isapprox, summary
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, adjointlayout,
                 sub_materialize, arguments, sub_paddeddata, paddeddata, PaddedLayout, resizedata!, LazyVector, ApplyLayout, call,
                 _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout,
                 AbstractCachedVector, AbstractCachedMatrix, paddeddata, cache_filldata!,
-                simplifiable, PaddedArray
+                simplifiable, PaddedArray, converteltype
 import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
                         subdiagonaldata, diagonaldata, supdiagonaldata, mul, rowsupport, colsupport
 import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, unitblocks, BlockRange1, AbstractLazyBandedLayout
-import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot
+import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot, mul!
 import BandedMatrices: AbstractBandedLayout, AbstractBandedMatrix, _BandedMatrix, bandeddata
 import FillArrays: AbstractFill, getindex_value, SquareEye
 import DomainSets: components
@@ -29,20 +29,20 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
                     _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector,
-                    AbstractQuasiFill, _dot
+                    AbstractQuasiFill, _dot, _equals, QuasiArrayLayout, PolynomialLayout
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
-    inbounds_getindex, grid, plotgrid, transform, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, Expansion,
-    AffineQuasiVector, AffineMap, WeightLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
-    checkpoints, weight, unweightedbasis, MappedBasisLayouts, __sum, invmap
+    inbounds_getindex, grid, plotgrid, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
+    AffineQuasiVector, AffineMap, WeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
+    checkpoints, weight, unweighted, MappedBasisLayouts, __sum, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, _broadcastbasis
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
-                        _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints
+                        _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan
 
 import FastGaussQuadrature: jacobimoment
 
-import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray, block, blockindex, BlockSlice
+import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray, block, blockindex, BlockSlice, blockvec
 import BandedMatrices: bandwidths
 
 export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial,
@@ -66,42 +66,53 @@ cardinality(::EuclideanDomain) = ℵ₁
 cardinality(::Union{DomainSets.RealNumbers,DomainSets.ComplexNumbers}) = ℵ₁
 cardinality(::Union{DomainSets.Integers,DomainSets.Rationals,DomainSets.NaturalNumbers}) = ℵ₀
 
-transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V}  =
-    adaptivetransform_ldiv(convert(AbstractQuasiArray{promote_type(T,V)}, A), f)
+include("adaptivetransform.jl")
 
+const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}
+abstract type OrthogonalPolynomial{T} <: Basis{T} end
+abstract type AbstractOPLayout <: AbstractBasisLayout end
+struct OPLayout <: AbstractOPLayout end
+MemoryLayout(::Type{<:OrthogonalPolynomial}) = OPLayout()
 
-setaxis(c, ::OneToInf) = c
-setaxis(c, ax::BlockedUnitRange) = PseudoBlockVector(c, (ax,))
-
-function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiArray{V}) where {U,V}
-    T = promote_type(U,V)
 
-    r = checkpoints(A)
-    fr = f[r]
-    maxabsfr = norm(fr,Inf)
 
-    tol = 20eps(real(T))
+sublayout(::AbstractOPLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) = MappedOPLayout()
 
-    for n = 2 .^ (4:∞)
-        An = A[:,oneto(n)]
-        cfs = An \ f
-        maxabsc = maximum(abs, cfs)
-        if maxabsc == 0 && maxabsfr == 0
-            return zeros(T,∞)
-        end
+struct MappedOPLayout <: AbstractOPLayout end
+struct WeightedOPLayout{Lay<:AbstractOPLayout} <: AbstractWeightedBasisLayout end
 
-        un = A * [cfs; Zeros{T}(∞)]
-        # we allow for transformed coefficients being a different size
-        ##TODO: how to do scaling for unnormalized bases like Jacobi?
-        if maximum(abs,@views(cfs[n-2:end])) < 10tol*maxabsc &&
-                all(norm.(un[r] - fr, 1) .< tol * n * maxabsfr*1000)
-            return setaxis([chop!(cfs, tol); zeros(T,∞)], axes(A,2))
-        end
-    end
-    error("Have not converged")
+isorthogonalityweighted(::WeightedOPLayout, _) = true
+function isorthogonalityweighted(::AbstractWeightedBasisLayout, wS)
+    w,S = arguments(wS)
+    w == orthogonalityweight(S)
 end
 
-abstract type OrthogonalPolynomial{T} <: Basis{T} end
+isorthogonalityweighted(wS) = isorthogonalityweighted(MemoryLayout(wS), wS)
+
+
+_equals(::MappedOPLayout, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
+_equals(::MappedOPLayout, ::MappedBasisLayouts, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
+_equals(::MappedBasisLayouts, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
+
+_broadcastbasis(::typeof(+), ::MappedOPLayout, ::MappedOPLayout, P, Q) where {L,M} = _broadcastbasis(+, MappedBasisLayout(), MappedBasisLayout(), P, Q)
+_broadcastbasis(::typeof(+), ::MappedOPLayout, M::MappedBasisLayout, P, Q) where L = _broadcastbasis(+, MappedBasisLayout(), M, P, Q)
+_broadcastbasis(::typeof(+), L::MappedBasisLayout, ::MappedOPLayout, P, Q) where M = _broadcastbasis(+, L, MappedBasisLayout(), P, Q)
+__sum(::MappedOPLayout, A, dims) = __sum(MappedBasisLayout(), A, dims)
+
+# demap to avoid Golub-Welsch fallback
+ContinuumArrays.transform_ldiv_if_columns(L::Ldiv{MappedOPLayout,Lay}, ax::OneTo) where Lay = ContinuumArrays.transform_ldiv_if_columns(Ldiv{MappedBasisLayout,Lay}(L.A,L.B), ax)
+ContinuumArrays.transform_ldiv_if_columns(L::Ldiv{MappedOPLayout,ApplyLayout{typeof(hcat)}}, ax::OneTo) = ContinuumArrays.transform_ldiv_if_columns(Ldiv{MappedBasisLayout,UnknownLayout}(L.A,L.B), ax)
+
+_equals(::AbstractOPLayout, ::AbstractWeightedBasisLayout, _, _) = false # Weighedt-Legendre doesn't exist
+_equals(::AbstractWeightedBasisLayout, ::AbstractOPLayout, _, _) = false # Weighedt-Legendre doesn't exist
+
+_equals(::WeightedOPLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ)
+_equals(::WeightedOPLayout, ::WeightedBasisLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
+_equals(::WeightedBasisLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
+_equals(::WeightedBasisLayout{<:AbstractOPLayout}, ::WeightedBasisLayout{<:AbstractOPLayout}, A, B) = A.f == B.f && all(A.args .== B.args)
+    
+
+copy(L::Ldiv{MappedOPLayout,Lay}) where Lay<:MappedBasisLayouts = copy(Ldiv{MappedBasisLayout,Lay}(L.A,L.B))
 
 # OPs are immutable
 copy(a::OrthogonalPolynomial) = a
@@ -146,12 +157,6 @@ function recurrencecoefficients(Q::AbstractQuasiMatrix{T}) where T
 end
 
 
-const WeightedOrthogonalPolynomial{T, A<:AbstractQuasiVector, B<:OrthogonalPolynomial} = WeightedBasis{T, A, B}
-
-function isorthogonalityweighted(wS::WeightedOrthogonalPolynomial)
-    w,S = wS.args
-    w == orthogonalityweight(S)
-end
 
 """
     singularities(f)
@@ -162,9 +167,9 @@ gives the singularity structure of an expansion, e.g.,
 singularities(::WeightLayout, w) = w
 singularities(lay::BroadcastLayout, a) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, a))...)
 singularities(::WeightedBasisLayouts, a) = singularities(BroadcastLayout{typeof(*)}(), a)
+singularities(::WeightedOPLayout, a) = singularities(weight(a))
 singularities(w) = singularities(MemoryLayout(w), w)
-singularities(f::Expansion) = singularities(basis(f))
-singularities(S::WeightedOrthogonalPolynomial) = singularities(S.args[1])
+singularities(::ExpansionLayout, f) = singularities(basis(f))
 
 singularities(S::SubQuasiArray) = singularities(parent(S))[parentindices(S)[1]]
 
@@ -184,83 +189,32 @@ function massmatrix(P::SubQuasiArray{<:Any,2,<:Any,<:Tuple{AbstractAffineQuasiVe
     massmatrix(Q)/kr.A
 end
 
-_weighted(w, P) = w .* P
-weighted(P::AbstractQuasiMatrix) = _weighted(orthogonalityweight(P), P)
+weighted(P::AbstractQuasiMatrix) = Weighted(P)
 
 OrthogonalPolynomial(w::Weight) =error("Override for $(typeof(w))")
 
 @simplify *(B::Identity, C::OrthogonalPolynomial) = C*jacobimatrix(C)
 
-function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::OrthogonalPolynomial)
+function layout_broadcasted(::Tuple{PolynomialLayout,AbstractOPLayout}, ::typeof(*), x::Inclusion, C)
     x == axes(C,1) || throw(DimensionMismatch())
     C*jacobimatrix(C)
 end
 
+
+
 # function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::BroadcastQuasiVector, C::OrthogonalPolynomial)
 #     axes(a,1) == axes(C,1) || throw(DimensionMismatch())
 #     # re-expand in OP basis
 #     broadcast(*, C * (C \ a), C)
 # end
 
-function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::AbstractAffineQuasiVector, C::OrthogonalPolynomial)
-    x = axes(C,1)
-    axes(a,1) == x || throw(DimensionMismatch())
-    broadcast(*, C * (C \ a), C)
-end
-
-function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::WeightedOrthogonalPolynomial)
-    x == axes(C,1) || throw(DimensionMismatch())
-    w,P = C.args
-    P2, J = (x .* P).args
-    @assert P == P2
-    (w .* P) * J
-end
-
-##
-# Multiplication for mapped and subviews x .* view(P,...)
-##
-
-function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
-    T = promote_type(eltype(x), eltype(C))
-    x == axes(C,1) || throw(DimensionMismatch())
-    P = parent(C)
-    kr,jr = parentindices(C)
-    y = axes(P,1)
-    Y = P \ (y .* P)
-    X = kr.A \ (Y     - kr.b * Eye{T}(∞))
-    P[kr, :] * view(X,:,jr)
-end
-
-function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
-    T = promote_type(eltype(x), eltype(C))
-    x == axes(C,1) || throw(DimensionMismatch())
-    P = parent(C)
-    kr,_ = parentindices(C)
-    y = axes(P,1)
-    Y = P \ (y .* P)
-    X = kr.A \ (Y     - kr.b * Eye{T}(∞))
-    P[kr, :] * X
-end
-
-
-# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
-#     T = promote_type(eltype(f), eltype(C))
-#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
-#     P = parent(C)
-#     kr,jr = parentindices(C)
-#     (f[invmap(kr)] .* P)[kr,jr]
+# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::AbstractAffineQuasiVector, C::OrthogonalPolynomial)
+#     x = axes(C,1)
+#     axes(a,1) == x || throw(DimensionMismatch())
+#     broadcast(*, C * (C \ a), C)
 # end
 
-# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
-#     T = promote_type(eltype(f), eltype(C))
-#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
-#     P = parent(C)
-#     kr,jr = parentindices(C)
-#     (f[invmap(kr)] .* P)[kr,jr]
-# end
 
-broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::Broadcasted, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) =
-    broadcast(*, materialize(f), C)
 
 function jacobimatrix(C::SubQuasiArray{T,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) where T
     P = parent(C)
@@ -315,22 +269,22 @@ function golubwelsch(V::SubQuasiArray)
     x,w
 end
 
-function factorize(L::SubQuasiArray{T,2,<:Normalized,<:Tuple{Inclusion,OneTo}}) where T
+function factorize(L::SubQuasiArray{T,2,<:Normalized,<:Tuple{Inclusion,OneTo}}, dims...; kws...) where T
     x,w = golubwelsch(L)
     TransformFactorization(x, L[x,:]'*Diagonal(w))
 end
 
 
-function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}) where T
+function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}, dims...; kws...) where T
     Q = Normalized(parent(L))[parentindices(L)...]
     D = L \ Q
-    F = factorize(Q)
+    F = factorize(Q, dims...; kws...)
     TransformFactorization(F.grid, D*F.plan)
 end
 
-function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
+function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:AbstractUnitRange}}, dims...; kws...) where T
     _,jr = parentindices(L)
-    ProjectionFactorization(factorize(parent(L)[:,oneto(maximum(jr))]), jr)
+    ProjectionFactorization(factorize(parent(L)[:,oneto(maximum(jr))], dims...; kws...), jr)
 end
 
 function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial})
@@ -346,11 +300,11 @@ function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}},
 end
 
 # assume we can expand w_B in wA to reduce to polynomial multiplication
-function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
-    _,A = arguments(wA)
-    w_B,B = arguments(wB)
-    A \ ((A * (wA \ w_B)) .* B)
-end
+# function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
+#     _,A = arguments(wA)
+#     w_B,B = arguments(wB)
+#     A \ ((A * (wA \ w_B)) .* B)
+# end
 
 ## special expansion routines for constants and x
 function _op_ldiv(P::AbstractQuasiMatrix{V}, f::AbstractQuasiFill{T,1}) where {T,V}