cumsum for Weighted(ChebyshevT()) (#170)

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -38,7 +38,7 @@ import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, Abstract
     grid, plotgrid, plotgrid_layout, plotvalues_layout, grid_layout, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
     AffineQuasiVector, AffineMap, AbstractWeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
     checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, broadcastbasis_layout,
-    plan_grid_transform, plan_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout, plan_transform_layout
+    plan_grid_transform, plan_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout, plan_transform_layout, _cumsum
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg
 

src/classical/chebyshev.jl

@@ -96,8 +96,14 @@ function getindex(w::ChebyshevUWeight, x::Number)
     sqrt(1-x^2)
 end
 
-sum(::ChebyshevWeight{1,T}) where T = convert(T,π)
-sum(::ChebyshevWeight{2,T}) where T = convert(T,π)/2
+sum(::ChebyshevTWeight{T}) where T = convert(T,π)
+sum(::ChebyshevUWeight{T}) where T = convert(T,π)/2
+
+function _cumsum(w::ChebyshevTWeight{V}, dims) where V
+    @assert dims == 1
+    x = axes(w,1)
+    (V(π)/2 .+ asin.(x))
+end
 
 normalizationconstant(::ChebyshevT{T}) where T = Vcat(sqrt(inv(convert(T,π))), Fill(sqrt(2/convert(T,π)),∞))
 
@@ -298,12 +304,17 @@ function _sum(::Weighted{T,<:Chebyshev}, dims) where T
     Hcat(convert(T, π), Zeros{T}(1,∞))
 end
 
-function cumsum(T::ChebyshevT{V}; dims::Integer) where V
+function _cumsum(T::ChebyshevT{V}, dims) where V
     @assert dims == 1
     Σ = _BandedMatrix(Vcat(-one(V) ./ (-2:2:∞)', Zeros{V}(1,∞), Hcat(one(V), one(V) ./ (4:2:∞)')), ℵ₀, 0, 2)
     ApplyQuasiArray(*, T, Vcat((-1).^(0:∞)'* Σ, Σ))
 end
 
+function _cumsum(W::Weighted{V, ChebyshevT{V}}, dims) where V
+    @assert dims == 1
+    [cumsum(ChebyshevTWeight{V}()) Weighted(ChebyshevU{V}())] * Diagonal(Vcat(one(V), -inv.(one(V):∞)))
+end
+
 ####
 # algebra
 ####

test/test_chebyshev.jl

@@ -390,6 +390,9 @@ import BandedMatrices: isbanded
         @test (T \ cumsum(T; dims=1)) * (T \ exp.(x)) ≈ T \ (exp.(x) .- exp(-1))
         f = T * (T \ exp.(x))
         @test T \ cumsum(f) ≈ T \ (exp.(x) .- exp(-1))
+
+        f = expand(Weighted(T), x -> exp(x)/sqrt(1-x^2))
+        @test cumsum(f)[0.1] ≈ 0.9784356617497387
     end
 
     @testset "algebra" begin