cumsum for Legendre (#248)

* cumsum for Legendre

* Update test_roots.jl

* Create detpointprocess.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.15.5"
+version = "0.15.6"
 
 
 [deps]

examples/detpointprocess.jl

@@ -0,0 +1,35 @@
+using ClassicalOrthogonalPolynomials, Plots
+using ClassicalOrthogonalPolynomials: sample
+
+x = Inclusion(ChebyshevInterval())
+
+function christoffel(A)
+    Q,R = qr(A)
+    n = size(A,2)
+    sum(expand(Q[:,k] .^2) for k=1:n)/n
+end
+
+function dpp(A)
+    m = size(A,2)
+    Q,R = qr(A)
+    r = Float64[]
+    for n = m:-1:1
+        Kₙ = sum(expand(Q[:,k] .^2) for k=1:n)/n
+        r₁ = sample(Kₙ)
+        push!(r, r₁)
+        Q = Q*nullspace(Q[r₁, :]')
+    end
+    r
+end
+
+m = 10
+A = cos.((0:m)' .* x)
+r = union([dpp(A) for _=1:1000]...)
+histogram(r; nbins=50, normalized=true)
+plot!(christoffel(A); ylims=(0,1))
+
+## DPPs are much better condtioned
+Q,R = qr(A)
+@test cond(Q[dpp(A),:]) ≤ 100
+@test cond(Q[sample(christoffel(A),11),:]) ≥ 1000
+@test cond(Q[range(-1,1,11),:]) > 1E13

src/classical/ultraspherical.jl

@@ -185,6 +185,12 @@ function diff(WS::WeightedUltraspherical{T}; dims=1) where T
     end
 end
 
+function _cumsum(P::Legendre{V}, dims) where V
+    @assert dims == 1
+    Σ = Bidiagonal(Vcat(1, Zeros{V}(∞)), Fill(-one(V), ∞), :L)
+    ApplyQuasiArray(*, Ultraspherical(-one(V)/2), Σ)
+end
+
 
 ##########
 # Conversion

test/test_legendre.jl

@@ -229,4 +229,12 @@ import QuasiArrays: MulQuasiArray
     @testset "ChebyshevInterval constructor" begin
         @test legendre(ChebyshevInterval()) ≡ Legendre()
     end
+
+    @testset "cumsum" begin
+        P = Legendre()
+        x = axes(P,1)
+        @test (P \ cumsum(P)) * (P \ exp.(x)) ≈ P \ (exp.(x) .- exp(-1))
+        f = P * (P \ exp.(x))
+        @test P \ cumsum(f) ≈ P \ (exp.(x) .- exp(-1))
+    end
 end

test/test_roots.jl

@@ -30,7 +30,7 @@ end
     Q,R = qr(A)
     @test Q[0.1,:]'R ≈ A[0.1,:]'
 
-    # sum(expand(Q[:,k] .^2) for k=axes(Q,2))
+    @test abs(sum(sample(sum(expand(Q[:,k] .^2) for k=axes(Q,2)), 1000))) ≤ 100 # mean is (numerically) zero
 end
 
 @testset "minimum/maximum/extrema (#242)" begin