Extra Legendre constructors, sum and dot for generic quasi-vectors (#137)

* Extra Legendre constructors

* v0.8.2

* Add routines for non-expanded functions

* singularities on onclusion

* dot tests

* tests pass

* Update show

* sum_layout

* v0.10

* Update documentation.yml

Author: dlfivefifty

.github/workflows/documentation.yml

@@ -16,7 +16,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.7'
+          version: '1.9'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.9.0"
+version = "0.10"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,19 +29,19 @@ ArrayLayouts = "1.0.1"
 BandedMatrices = "0.17.17"
 BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.12"
-ContinuumArrays = "0.12.4"
-DomainSets = "0.5.6, 0.6"
+ContinuumArrays = "0.13"
+DomainSets = "0.6"
 FFTW = "1.1"
-FastGaussQuadrature = "0.4.3, 0.5"
+FastGaussQuadrature = "0.5"
 FastTransforms = "0.15.2"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.11"
 InfiniteLinearAlgebra = "0.6.16"
-IntervalSets = "0.5, 0.6, 0.7"
+IntervalSets = "0.7"
 LazyArrays = "1.0.1"
 LazyBandedMatrices = "0.8.5"
-QuasiArrays = "0.9.6"
+QuasiArrays = "0.10"
 SpecialFunctions = "1.0, 2"
 julia = "1.9"
 

src/ClassicalOrthogonalPolynomials.jl

@@ -34,10 +34,10 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!
-import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
+import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     grid, plotgrid, _plotgrid, _grid, 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,
+    checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, _broadcastbasis,
     plan_transform, plan_grid_transform, MAX_PLOT_POINTS, MulPlan
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg
@@ -94,7 +94,7 @@ _equals(::MappedBasisLayouts, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) &&
 _broadcastbasis(::typeof(+), ::MappedOPLayout, ::MappedOPLayout, P, Q) = _broadcastbasis(+, MappedBasisLayout(), MappedBasisLayout(), P, Q)
 _broadcastbasis(::typeof(+), ::MappedOPLayout, M::MappedBasisLayout, P, Q) = _broadcastbasis(+, MappedBasisLayout(), M, P, Q)
 _broadcastbasis(::typeof(+), L::MappedBasisLayout, ::MappedOPLayout, P, Q) = _broadcastbasis(+, L, MappedBasisLayout(), P, Q)
-__sum(::MappedOPLayout, A, dims) = __sum(MappedBasisLayout(), A, dims)
+sum_layout(::MappedOPLayout, A, dims) = sum_layout(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)
@@ -161,13 +161,16 @@ gives the singularity structure of an expansion, e.g.,
 `JacobiWeight`.
 """
 singularities(::WeightLayout, w) = w
-singularities(lay::BroadcastLayout, a) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, a))...)
+singularities(lay::BroadcastLayout, a::AbstractQuasiVector) = 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(::ExpansionLayout, f) = singularities(basis(f))
 
-singularities(S::SubQuasiArray) = singularities(parent(S))[parentindices(S)[1]]
+singularitiesview(w, ::Inclusion) = w # for now just assume it doesn't change
+singularities(S::SubQuasiArray) = singularitiesview(singularities(parent(S)), parentindices(S)[1])
+
+basis_axes(::Inclusion{<:Any,<:AbstractInterval}, v) = basis_singularities(singularities(v))
 
 struct NoSingularities end
 

src/classical/jacobi.jl

@@ -83,6 +83,8 @@ jacobi(a,b, d::AbstractInterval{T}) where T = Jacobi{float(promote_type(eltype(a
 
 Jacobi(P::Legendre{T}) where T = Jacobi(zero(T), zero(T))
 
+basis_singularities(w::JacobiWeight) = Weighted(Jacobi(w.a, w.b))
+
 """
      jacobip(n, a, b, z)
 
@@ -154,7 +156,8 @@ axes(::AbstractJacobi{T}) where T = (Inclusion{T}(ChebyshevInterval{real(T)}()),
 ==(A::Legendre, B::Weighted{<:Any,<:AbstractJacobi}) = A == B.P
 ==(A::Weighted{<:Any,<:AbstractJacobi}, B::Legendre) = A.P == B
 
-
+show(io::IO, w::AbstractJacobiWeight) = summary(io, w)
+show(io::IO, P::AbstractJacobi) = summary(io, P)
 summary(io::IO, P::Jacobi) = print(io, "Jacobi($(P.a), $(P.b))")
 
 ###

src/classical/legendre.jl

@@ -35,10 +35,13 @@ singularitiesbroadcast(::typeof(^), L::LegendreWeight, ::NoSingularities) = L
 singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L # can't find roots
 
 singularities(::AbstractJacobi{T}) where T = LegendreWeight{T}()
-singularities(::Inclusion{T,<:AbstractInterval}) where T = LegendreWeight{T}()
-singularities(d::Inclusion{T,<:Interval}) where T = LegendreWeight{T}()[affine(d,ChebyshevInterval{T}())]
+singularities(::Inclusion{T,<:ChebyshevInterval}) where T = LegendreWeight{T}()
+singularities(d::Inclusion{T,<:AbstractInterval}) where T = LegendreWeight{T}()[affine(d,ChebyshevInterval{T}())]
 singularities(::AbstractFillLayout, P) = LegendreWeight{eltype(P)}()
 
+basis_singularities(::LegendreWeight{T}) where T = Legendre{T}()
+basis_singularities(v::SubQuasiArray) = view(basis_singularities(parent(v)), parentindices(v)[1], :)
+
 struct Legendre{T} <: AbstractJacobi{T} end
 Legendre() = Legendre{Float64}()
 
@@ -49,6 +52,8 @@ weighted(P::SubQuasiArray{<:Any,2,<:Normalized{<:Any,<:Legendre}}) = P
 
 legendre() = Legendre()
 legendre(d::AbstractInterval{T}) where T = Legendre{float(T)}()[affine(d,ChebyshevInterval{T}()), :]
+legendre(d::ChebyshevInterval{T}) where T = Legendre{float(T)}()
+legendre(d::Inclusion) = legendre(d.domain)
 
 """
      legendrep(n, z)
@@ -161,14 +166,3 @@ function _sum(P::Legendre{T}, dims) where T
 end
 
 _sum(p::SubQuasiArray{T,1,Legendre{T},<:Tuple{Inclusion,Int}}, ::Colon) where T = parentindices(p)[2] == 1 ? convert(T, 2) : zero(T)
-
-
-###
-# dot
-###
-
-_dot(::Inclusion{<:Any,<:ChebyshevInterval}, a, b) = _dot(singularities(a), singularities(b), a, b)
-function _dot(::LegendreWeight, ::LegendreWeight, a, b)
-    P = Legendre{promote_type(eltype(a),eltype(b))}()
-    dot(P\a, (massmatrix(P) * (P\b)))
-end

src/classical/ultraspherical.jl

@@ -13,6 +13,8 @@ UltrasphericalWeight{T}(λ) where T = UltrasphericalWeight{T,typeof(λ)}(λ)
 UltrasphericalWeight(λ) = UltrasphericalWeight{float(typeof(λ)),typeof(λ)}(λ)
 UltrasphericalWeight(::LegendreWeight{T}) where T = UltrasphericalWeight(one(T)/2)
 
+summary(io::IO, w::UltrasphericalWeight) = print(io, "UltrasphericalWeight($(w.λ))")
+
 ==(a::UltrasphericalWeight, b::UltrasphericalWeight) = a.λ == b.λ
 
 function getindex(w::UltrasphericalWeight, x::Number)
@@ -36,6 +38,8 @@ end
 
 Ultraspherical(::ChebyshevU{T}) where T = Ultraspherical{T}(1)
 
+summary(io::IO, w::Ultraspherical) = print(io, "Ultraspherical($(w.λ))")
+
 const WeightedUltraspherical{T} = WeightedBasis{T,<:UltrasphericalWeight,<:Ultraspherical}
 
 orthogonalityweight(C::Ultraspherical) = UltrasphericalWeight(C.λ)

src/normalized.jl

@@ -298,10 +298,10 @@ function mul(A::Derivative, B::Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:Abstract
 end
 mul(Ac::QuasiAdjoint{<:Any, Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = mul(Bc', Ac')'    
 
-summary(io::IO, Q::Weighted) = print(io, "Weighted($(Q.P))")
+show(io::IO, Q::Weighted) = print(io, "Weighted($(Q.P))")
 
-__sum(::AbstractNormalizedOPLayout, A, dims) = __sum(ApplyLayout{typeof(*)}(), A, dims)
-function __sum(::WeightedOPLayout, A, dims)
+sum_layout(::AbstractNormalizedOPLayout, A, dims) = sum_layout(ApplyLayout{typeof(*)}(), A, dims)
+function sum_layout(::WeightedOPLayout, A, dims)
     @assert dims == 1
     Hcat(sum(weight(A)), Zeros{eltype(A)}(1,∞))
 end

test/runtests.jl

@@ -56,4 +56,20 @@ include("test_roots.jl")
 
     f = x -> ChebyshevT{eltype(x)}()[x,1:5]
     @test ForwardDiff.derivative(f,0.1) ≈ [0;(1:4).*U[0.1,1:4]]
+end
+
+@testset "basis" begin
+    for x in (Inclusion(ChebyshevInterval()), Inclusion(1 .. 2))
+        a,b = first(x),last(x)
+        @test sum(x) == b-a
+        @test sum(x .^ 2) ≈ (b^3 - a^3)/3
+        @test sum(exp.(x)) ≈ exp(b) - exp(a)
+        @test dot(x, x) ≈ sum(expand(x .^2))
+        @test dot(x.^2, x.^2) ≈ sum(expand(x .^4))
+        @test dot(exp.(x), x.^2) ≈ sum(expand(x .^2 .* exp.(x)))
+        @test dot(x, exp.(x)) ≈ dot(exp.(x), x)
+    end
+
+    # A = x .^ (0:2)'
+    # sum(A; dims=1)
 end

test/test_legendre.jl

@@ -86,10 +86,17 @@ import QuasiArrays: MulQuasiArray
         @test W isa Clenshaw
         @test W * [1; 2; zeros(∞)] ≈ P \ (w .* (P[:,1:2] * [1,2]))
 
+        f = expand(P, exp)
+        @test (w .* f)[0.1] ≈ w[0.1]exp(0.1)
+
         M = P'P
         @test M isa Diagonal
         @test P'x ≈ [0; 2/3; zeros(∞)]
-        @test P'exp.(x) ≈ M * (P\exp.(x))
+        @test P'f ≈ M * (P\f)
+
+        v = f.args[2]
+        @test v'*M isa Adjoint
+        @test v'*M*v ≈ dot(f,f) ≈ f'f
     end
 
     @testset "test on functions" begin

test/test_ultraspherical.jl

@@ -145,4 +145,9 @@ import LazyArrays: rowsupport, colsupport
         c = [1; 2; 3; zeros(BigFloat,∞)]
         @test C³[big(1)/10,:]'*(C³ \ U) * c ≈ U[big(1)/10,:]'c
     end
+
+    @testset "show" begin
+        @test stringmime("text/plain",UltrasphericalWeight(1)) == "UltrasphericalWeight(1)"
+        @test stringmime("text/plain",Ultraspherical(1)) == "Ultraspherical(1)"
+    end
 end