Fixe constant for ConvertedOrthogonalPolynomial (#226)

* Introduce GramMatrixPolynomial

* Fixes to ConvertedOP

* Update choleskyQR.jl

* make ConvertedOP normalized again

* Update Project.toml

* Update choleskyQR.jl

* Update test_choleskyQR.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.14.2"
+version = "0.14.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/choleskyQR.jl

@@ -1,5 +1,5 @@
 """
-Represent an Orthogonal polynomial which has a conversion operator from P, that is, Q = P * inv(U).
+Represents orthonormal polynomials defined via a conversion operator from P, that is, Q = P * inv(U).
 """
 struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP} <: OrthonormalPolynomial{T}
     weight::WW
@@ -8,15 +8,22 @@ struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP}
     P::PP
 end
 
-_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P)
+_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P)/Q.U[1,1]
 
 axes(Q::ConvertedOrthogonalPolynomial) = axes(Q.P)
+
+
+struct ConvertedOPLayout <: AbstractNormalizedOPLayout end
 MemoryLayout(::Type{<:ConvertedOrthogonalPolynomial}) = ConvertedOPLayout()
+
+
+
+
 jacobimatrix(Q::ConvertedOrthogonalPolynomial) = Q.X
 orthogonalityweight(Q::ConvertedOrthogonalPolynomial) = Q.weight
 
 
-# transform to P * U if needed for differentiation, etc.
+# transform to P * inv(U) if needed for differentiation, etc.
 arguments(::ApplyLayout{typeof(*)}, Q::ConvertedOrthogonalPolynomial) = Q.P, ApplyArray(inv, Q.U)
 
 OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthogonalpolynomial(singularities(w)))

src/normalized.jl

@@ -35,7 +35,6 @@ isnormalized(_) = false
 represents OPs that are of the form P * R where P is another family of OPs and R is upper-triangular.
 """
 abstract type AbstractNormalizedOPLayout <: AbstractOPLayout end
-struct ConvertedOPLayout <: AbstractNormalizedOPLayout end
 struct NormalizedOPLayout{LAY<:AbstractBasisLayout} <: AbstractNormalizedOPLayout end
 
 MemoryLayout(::Type{<:Normalized{<:Any, OPs}}) where OPs = NormalizedOPLayout{typeof(MemoryLayout(OPs))}()

test/test_choleskyQR.jl

@@ -1,5 +1,5 @@
 using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices, InfiniteLinearAlgebra
-import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial
+import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial, _p0
 import LazyArrays: AbstractCachedMatrix, resizedata!
 
 @testset "CholeskyQR" begin
@@ -228,9 +228,8 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
         @test Q == Q̃
         @test Q̃ == Q
 
-        @test Q[0.1,1] ≈ 1/sqrt(2)
+        @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(2)
         @test Q[0.1,1:10] ≈ Q̃[0.1,1:10]
-        # AWESOME, thanks TSGUT!!
         @test Q[0.1,10_000] ≈ Q̃[0.1,10_000]
 
         R = P \ Q
@@ -247,6 +246,9 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
     @testset "Chebyshev" begin
          U = ChebyshevU()
          Q = orthogonalpolynomial(x -> (1+x^2)*sqrt(1-x^2), U)
+         x = axes(U,1)
+         
+         @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(sum(expand(Weighted(U),x -> (1+x^2)*sqrt(1-x^2))))
          @test bandwidths(Q\U) == (0,2)
 
          Q̃ = OrthogonalPolynomial(x -> (1+x^2)*sqrt(1-x^2), U)

test/test_gram.jl

@@ -0,0 +1,10 @@
+using ClassicalOrthogonalPolynomials, FastTransforms
+
+P = Legendre()
+x = axes(P,1)
+w = @.(1-x^2)
+
+μ = P'w
+X = jacobimatrix(P)
+n = 20
+@test GramMatrix(μ[1:2n], X[1:2n,1:2n]) ≈ (P' * (w .* P))[1:n,1:n]