Support Mapped-Weighted OPs (#69)

* Support Mapped-Weighted OPs

* Improvements to massmatrix and interlace

* Update test_lanczos.jl

* Much faster OffHilbert

* two-interval Hilbert inverse works

* much faster solves

* Update Project.toml

* v0.4.7

* sum for PiecewiseInterlace

* mapped log kernel

* Support both Weighted(::SubArray) and SubArray(::Weighted,...)

* remove old Hilber

* Update stieltjes.jl

* Update test_stieltjes.jl

* bug fixes

* Update Project.toml

* fix plotting for PiecewiseInterlace

* Update Project.toml

* Update Project.toml

* fixes for DomainSets v0.5.5

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.4.6"
+version = "0.4.7"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -25,22 +25,22 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "0.7"
-BandedMatrices = "0.16.5"
-BlockArrays = "0.16"
-BlockBandedMatrices = "0.10.9"
+ArrayLayouts = "0.7.5"
+BandedMatrices = "0.16.11"
+BlockArrays = "0.16.6"
+BlockBandedMatrices = "0.11"
 ContinuumArrays = "0.9.1"
-DomainSets = "0.5"
+DomainSets = "0.5.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
 FastTransforms = "0.12"
 FillArrays = "0.11, 0.12"
 HypergeometricFunctions = "0.3.4"
-InfiniteArrays = "0.11.1, 0.12"
-InfiniteLinearAlgebra = "0.5.8"
+InfiniteArrays = "0.12"
+InfiniteLinearAlgebra = "0.6"
 IntervalSets = "0.5"
-LazyArrays = "0.21.16"
-LazyBandedMatrices = "0.6.6"
+LazyArrays = "0.21.20"
+LazyBandedMatrices = "0.7"
 QuasiArrays = "0.7"
 SpecialFunctions = "1.0"
 julia = "1.6"

src/ClassicalOrthogonalPolynomials.jl

@@ -12,24 +12,26 @@ import Base: @_inline_meta, axes, getindex, unsafe_getindex, convert, prod, *, /
                 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, colsupport, adjointlayout,
+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!
+                AbstractCachedVector, AbstractCachedMatrix, paddeddata, cache_filldata!,
+                simplifiable, PaddedArray
 import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
-                        subdiagonaldata, diagonaldata, supdiagonaldata
+                        subdiagonaldata, diagonaldata, supdiagonaldata, mul, rowsupport, colsupport
 import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, unitblocks, BlockRange1, AbstractLazyBandedLayout
 import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot
 import BandedMatrices: AbstractBandedLayout, AbstractBandedMatrix, _BandedMatrix, bandeddata
 import FillArrays: AbstractFill, getindex_value, SquareEye
-
+import DomainSets: components
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
                     _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector
 
 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,
@@ -60,22 +62,12 @@ include("interlace.jl")
 
 cardinality(::FullSpace{<:AbstractFloat}) = ℵ₁
 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)
 
-function chop!(c::AbstractVector, tol::Real)
-    @assert tol >= 0
-
-    for k=length(c):-1:1
-        if abs(c[k]) > tol
-            resize!(c,k)
-            return c
-        end
-    end
-    resize!(c,0)
-    c
-end
 
 setaxis(c, ::OneToInf) = c
 setaxis(c, ax::BlockedUnitRange) = PseudoBlockVector(c, (ax,))
@@ -185,6 +177,12 @@ singularities(r::Base.RefValue) = r[] # pass through
 orthogonalityweight(P::SubQuasiArray{<:Any,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) =
     orthogonalityweight(parent(P))[parentindices(P)[1]]
 
+function massmatrix(P::SubQuasiArray{<:Any,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}})
+    Q = parent(P)
+    kr,jr = parentindices(P)
+    massmatrix(Q)/kr.A
+end
+
 _weighted(w, P) = w .* P
 weighted(P::AbstractQuasiMatrix) = _weighted(orthogonalityweight(P), P)
 

src/classical/chebyshev.jl

@@ -73,6 +73,9 @@ chebysevuweight(d::AbstractInterval{T}) where T = ChebyshevUWeight{float(T)}[aff
 ==(::Chebyshev, ::Legendre) = false
 ==(::Legendre, ::Chebyshev) = false
 
+summary(io::IO, w::ChebyshevT{Float64}) = print(io, "ChebyshevT()")
+summary(io::IO, w::ChebyshevU{Float64}) = print(io, "ChebyshevU()")
+
 OrthogonalPolynomial(w::ChebyshevWeight{kind,T}) where {kind,T} = Chebyshev{kind,T}()
 orthogonalityweight(P::Chebyshev{kind,T}) where {kind,T} = ChebyshevWeight{kind,T}()
 
@@ -158,6 +161,12 @@ end
     ApplyQuasiMatrix(*, ChebyshevU{T}(), A)
 end
 
+@simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, W::Weighted{<:Any,<:ChebyshevU})
+    T = promote_type(eltype(D),eltype(W))
+    Weighted(ChebyshevT{T}()) * _BandedMatrix((-one(T):-one(T):(-∞))', ℵ₀, 1,-1)
+end
+
+
 #####
 # Conversion
 #####

src/classical/fourier.jl

@@ -61,7 +61,7 @@ end
 end
 
 
-function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::BroadcastQuasiVector{<:Any,typeof(cos),<:Tuple{<:Inclusion{<:Any,<:FullSpace}}}, F::Fourier)
+function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::BroadcastQuasiVector{<:Any,typeof(cos),<:Tuple{<:Inclusion{<:Any,RealNumbers}}}, F::Fourier)
     axes(c,1) == axes(F,1) || throw(DimensionMismatch())
     T = promote_type(eltype(c), eltype(F))
     # Use LinearAlgebra.Tridiagonal for now since broadcasting support not complete for LazyBandedMatrices.Tridiagonal
@@ -71,7 +71,7 @@ function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::BroadcastQuasiVec
 end
 
 
-function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), s::BroadcastQuasiVector{<:Any,typeof(sin),<:Tuple{<:Inclusion{<:Any,<:FullSpace}}}, F::Fourier)
+function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), s::BroadcastQuasiVector{<:Any,typeof(sin),<:Tuple{<:Inclusion{<:Any,RealNumbers}}}, F::Fourier)
     axes(s,1) == axes(F,1) || throw(DimensionMismatch())
     T = promote_type(eltype(s), eltype(F))
     # Use LinearAlgebra.Tridiagonal for now since broadcasting support not complete for LazyBandedMatrices.Tridiagonal

src/classical/jacobi.jl

@@ -10,7 +10,7 @@ broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), w::AbstractJacobiWeight, v::A
 broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(sqrt), w::AbstractJacobiWeight) =
     JacobiWeight(w.a/2, w.b/2)
 
-broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, w::AbstractJacobiWeight, ::Base.RefValue{Val{k}}) where k = 
+broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, w::AbstractJacobiWeight, ::Base.RefValue{Val{k}}) where k =
     JacobiWeight(k * w.a, k * w.b)
 
 struct JacobiWeight{T} <: AbstractJacobiWeight{T}
@@ -136,7 +136,7 @@ Jacobi(P::Legendre{T}) where T = Jacobi(zero(T), zero(T))
 """
      jacobip(n, a, b, z)
 
-computes the `n`-th Jacobi polynomial, orthogonal with 
+computes the `n`-th Jacobi polynomial, orthogonal with
 respec to `(1-x)^a*(1+x)^b`, at `z`.
 """
 jacobip(n::Integer, a, b, z::Number) = Base.unsafe_getindex(Jacobi{promote_type(typeof(a), typeof(b), typeof(z))}(a,b), z, n+1)
@@ -237,25 +237,40 @@ end
 # Mass Matrix
 #########
 
-legendre_massmatrix(::Type{T}) where T = Diagonal(convert(T,2) ./ (2(0:∞) .+ 1))
+# massmatrix(P) = Weighted(P)'P
+massmatrix(P::Legendre{T}) where T = Diagonal(convert(T,2) ./ (2(0:∞) .+ 1))
 
+
+
+"""
+    legendre_massmatrix
+
+computes the massmatrix by first re-expanding in Legendre
+"""
 function legendre_massmatrix(Ac, B)
     A = parent(Ac)
     P = Legendre{eltype(B)}()
-    (P\A)'*legendre_massmatrix(eltype(P))*(P\B)
+    (P\A)'*massmatrix(P)*(P\B)
 end
 
-@simplify *(Ac::QuasiAdjoint{<:Any,<:Legendre}, B::Legendre) = legendre_massmatrix(Ac, B)
+@simplify *(Ac::QuasiAdjoint{<:Any,<:Legendre}, B::Legendre) = massmatrix(Legendre{promote_type(eltype(Ac), eltype(B))}())
 @simplify *(Ac::QuasiAdjoint{<:Any,<:AbstractJacobi}, B::AbstractJacobi) = legendre_massmatrix(Ac,B)
+@simplify *(Ac::QuasiAdjoint{<:Any,<:AbstractJacobi}, B::Weighted{<:Any,<:AbstractJacobi}) = legendre_massmatrix(Ac,B)
 
-# 2^{a + b + 1} {\Gamma(n+a+1) \Gamma(n+b+1) \over (2n+a+b+1) \Gamma(n+a+b+1) n!}.
+@simplify function *(Ac::QuasiAdjoint{<:Any,<:AbstractJacobi}, B::AbstractQuasiVector)
+    P = Legendre{eltype(Ac)}()
+    (Ac * P) * (P \ B)
+end
 
-function jacobi_massmatrix(a, b)
+# 2^{a + b + 1} {\Gamma(n+a+1) \Gamma(n+b+1) \over (2n+a+b+1) \Gamma(n+a+b+1) n!}.
+function massmatrix(P::Jacobi)
+    a,b = P.a,P.b
     n = 0:∞
     Diagonal(2^(a+b+1) .* (exp.(loggamma.(n .+ (a+1)) .+ loggamma.(n .+ (b+1)) .- loggamma.(n .+ (a+b+1)) .- loggamma.(n .+ 1)) ./ (2n .+ (a+b+1))))
 end
 
 
+
 @simplify function *(wAc::QuasiAdjoint{<:Any,<:WeightedBasis{<:Any,<:AbstractJacobiWeight}}, wB::WeightedBasis{<:Any,<:AbstractJacobiWeight})
     w,A = arguments(parent(wAc))
     v,B = arguments(wB)
@@ -265,7 +280,7 @@ end
     A = parent(Ac)
     w,B = arguments(wB)
     P = Jacobi(w.a, w.b)
-    (P\A)' * jacobi_massmatrix(w.a, w.b) * (P \ B)
+    (P\A)' * massmatrix(P) * (P \ B)
 end
 
 ########

src/clenshaw.jl

@@ -133,7 +133,9 @@ Base.@propagate_inbounds function _clenshaw_next!(n, A::AbstractVector, ::Zeros,
 end
 
 Base.@propagate_inbounds function _clenshaw_next!(n, A::AbstractVector, B::AbstractVector, C::AbstractVector, x::AbstractMatrix, c, bn1::AbstractMatrix{T}, bn2::AbstractMatrix{T}) where T
-    bn2 .= B[n] .* bn1 .- C[n+1] .* bn2
+    # bn2 .= B[n] .* bn1 .- C[n+1] .* bn2
+    lmul!(-C[n+1], bn2)
+    BLAS.axpy!(B[n], bn1, bn2)
     muladd!(A[n], x, bn1, one(T), bn2)
     view(bn2,band(0)) .+= c[n]
     bn2
@@ -166,7 +168,8 @@ Base.@propagate_inbounds function _clenshaw_first!(A, ::Zeros, C, X, c, bn1, bn2
 end
 
 Base.@propagate_inbounds function _clenshaw_first!(A, B, C, X, c, bn1, bn2) 
-    bn2 .= B[1] .* bn1 .- C[2] .* bn2
+    lmul!(-C[2], bn2)
+    BLAS.axpy!(B[1], bn1, bn2)
     muladd!(A[1], X, bn1, one(eltype(bn2)), bn2)
     view(bn2,band(0)) .+= c[1]
     bn2
@@ -283,17 +286,17 @@ function _BandedMatrix(::ClenshawLayout, V::SubArray{<:Any,2})
     M = parent(V)
     kr,jr = parentindices(V)
     b = bandwidth(M,1)
-    jkr=max(1,min(jr[1],kr[1])-b÷2):max(jr[end],kr[end])+b÷2
+    jkr=max(1,min(first(jr),first(kr))-b÷2):max(last(jr),last(kr))+b÷2
     # relationship between jkr and kr, jr
-    kr2,jr2 = kr.-jkr[1].+1,jr.-jkr[1].+1
+    kr2,jr2 = kr.-first(jkr).+1,jr.-first(jkr).+1
     lmul!(M.p0, clenshaw(M.c, M.A, M.B, M.C, M.X[jkr, jkr])[kr2,jr2])
 end
 
 function getindex(M::Clenshaw{T}, kr::AbstractUnitRange, j::Integer) where T
     b = bandwidth(M,1)
-    jkr=max(1,min(j,kr[1])-b÷2):max(j,kr[end])+b÷2
+    jkr=max(1,min(j,first(kr))-b÷2):max(j,last(kr))+b÷2
     # relationship between jkr and kr, jr
-    kr2,j2 = kr.-jkr[1].+1,j-jkr[1]+1
+    kr2,j2 = kr.-first(jkr).+1,j-first(jkr)+1
     f = [Zeros{T}(j2-1); one(T); Zeros{T}(length(jkr)-j2)]
     lmul!(M.p0, clenshaw(M.c, M.A, M.B, M.C, M.X[jkr, jkr], f)[kr2])
 end
@@ -322,11 +325,19 @@ function materialize!(M::MatMulVecAdd{<:ClenshawLayout,<:PaddedLayout,<:PaddedLa
     y
 end
 
+# TODO: generalise this to be trait based
 function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:OrthogonalPolynomial}, P::OrthogonalPolynomial)
     axes(a,1) == axes(P,1) || throw(DimensionMismatch())
     P * Clenshaw(a, P)
 end
 
+function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{AbstractAffineQuasiVector,Slice}}}, V::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{AbstractAffineQuasiVector,Any}})
+    axes(a,1) == axes(V,1) || throw(DimensionMismatch())
+    kr,jr = parentindices(V)
+    P = view(parent(V),kr,:)
+    P * Clenshaw(a, P)[:,jr]
+end
+
 function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:WeightedOrthogonalPolynomial}, P::OrthogonalPolynomial)
     axes(a,1) == axes(P,1) || throw(DimensionMismatch())
     wQ,c = arguments(a)

src/interlace.jl

@@ -96,7 +96,7 @@ end
 
 
 abstract type AbstractInterlaceBasis{T} <: Basis{T} end
-copy(A::AbstractInterlaceBasis) = interlacebasis(A, map(copy, A.args))
+copy(A::AbstractInterlaceBasis) = interlacebasis(A, map(copy, A.args)...)
 
 """
     PiecewiseInterlace(args...)
@@ -145,6 +145,18 @@ end
     interlacebasis(S, map(first, args)...) * BlockBroadcastArray{promote_type(eltype(D),eltype(S))}(Diagonal, unitblocks.(last.(args))...)
 end
 
+@simplify function *(Ac::QuasiAdjoint{<:Any,<:AbstractInterlaceBasis}, B::AbstractInterlaceBasis)
+    axes(Ac,2) == axes(B,1) || throw(DimensionMismatch())
+    BlockBroadcastArray{promote_type(eltype(Ac),eltype(B))}(Diagonal, unitblocks.(adjoint.(parent(Ac).args) .* B.args)...)
+end
+
+@simplify function *(Ac::QuasiAdjoint{<:Any,<:AbstractInterlaceBasis}, B::AbstractQuasiVector)
+    axes(Ac,2) == axes(B,1) || throw(DimensionMismatch())
+    args = (Ac').args
+    cs = (adjoint.(args) .* getindex.(Ref(B), axes.(args,1)))
+    BlockBroadcastArray(vcat, unitblocks.(cs)...)
+end
+
 struct PiecewiseFactorization{T,FF,Ax} <: Factorization{T}
     factorizations::FF
     axes::Ax
@@ -168,4 +180,45 @@ function factorize(V::SubQuasiArray{<:Any,2,<:AbstractInterlaceBasis,<:Tuple{Inc
     _,jr = parentindices(V)
     J = findblock(axes(P,2),maximum(jr))
     ProjectionFactorization(factorize(P[:,Block.(OneTo(Int(J)))]), jr)
-end
+end
+
+###
+# sum
+###
+
+_sum(P::PiecewiseInterlace, dims) = BlockBroadcastArray(hcat, unitblocks.(_sum.(P.args, dims))...)
+
+# blockvector2vectortuple
+
+function components(f::Expansion{<:Any,<:PiecewiseInterlace})
+    P,c = arguments(*, f)
+    P.args .* blockvector2vectortuple(c)
+end
+
+##
+# summary
+###
+
+function summary(io::IO, w::PiecewiseInterlace)
+    print(io, "PiecewiseInterlace(")
+    args = w.args
+    summary(io, first(args))
+    for a in tail(args)
+        print(io, ", ")
+        summary(io, a)
+    end
+    print(io, ")")
+end
+
+###
+# plot
+###
+
+
+function plotgrid(Pn::SubQuasiArray{T,2,<:PiecewiseInterlace,<:Tuple{Inclusion,Any}}) where T
+    kr,jr = parentindices(Pn)
+    P = parent(Pn)
+    N = findblock(axes(P,2),maximum(jr))
+    sort!(vcat(plotgrid.(getindex.(P.args, :, Ref(oneto(Int(N)))))...))
+end
+

src/lanczos.jl

@@ -112,7 +112,7 @@ triangulardata(R::LanczosConversion) = R
 sublayout(::LanczosConversionLayout, ::Type{<:Tuple{KR,Integer}}) where KR = 
     sublayout(PaddedLayout{UnknownLayout}(), Tuple{KR})
 
-function sub_paddeddata(::LanczosConversionLayout, S::SubArray{<:Any,1,<:AbstractMatrix})
+function sub_paddeddata(::LanczosConversionLayout, S::SubArray{<:Any,1,<:AbstractMatrix,<:Tuple{Any,Integer}})
     P = parent(S)
     (kr,j) = parentindices(S)
     resizedata!(P.data, j)