From da33c81d2a9916bf0e2879323be6fe243de03161 Mon Sep 17 00:00:00 2001 From: "Timon S. Gutleb" Date: Fri, 21 Jun 2024 04:28:29 -0700 Subject: [PATCH 01/10] lowering works --- src/SemiclassicalOrthogonalPolynomials.jl | 45 +++++++++++++++++------ test/runtests.jl | 21 ++++++++++- 2 files changed, 54 insertions(+), 12 deletions(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index b562b8a..1c92b59 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -180,19 +180,27 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c) cholesky_jacobimatrix(I-Q.X,Q) elseif iszero(Δa) && iszero(Δb) && isone(Δc) cholesky_jacobimatrix(Q.t*I-Q.X,Q) + elseif isone(-Δa) && iszero(Δb) && iszero(Δc) # raising by 1 + # TODO: This is re-constructing. It should instead use reverse Cholesky (or an alternative)! + semiclassical_jacobimatrix(Q.t,a,b,c) + elseif iszero(Δa) && isone(-Δb) && iszero(Δc) + # TODO: This is re-constructing. It should instead use reverse Cholesky (or an alternative)! + semiclassical_jacobimatrix(Q.t,a,b,c) + elseif iszero(Δa) && iszero(Δb) && isone(-Δc) + # TODO: This is re-constructing. It should instead use reverse Cholesky (or an alternative)! + semiclassical_jacobimatrix(Q.t,a,b,c) elseif a > Q.a # iterative raising by 1 semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a+1, Q.b, Q.c, Q), a, b, c) elseif b > Q.b semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b+1, Q.c, Q), a, b, c) elseif c > Q.c semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1, Q), a, b, c) - # TODO: Implement lowering via QL/reverse Cholesky or via inverting R? - # elseif b < Q.b # iterative lowering by 1 - # semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c) - # elseif c < Q.c - # semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1, Q), a, b, c) - # elseif a < Q.a - # semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a-1, Q.b, Q.c, Q), a, b, c) + elseif a < Q.a # iterative lowering by 1 + semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a-1, Q.b, Q.c, Q), a, b, c) + elseif b < Q.b + semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c) + elseif c < Q.c + semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1, Q), a, b, c) else error("Not Implemented") end @@ -268,7 +276,7 @@ function semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) Δb = Qb-P.b Δc = Qc-P.c M = Diagonal(Ones(∞)) - if iseven(Δa) && iseven(Δb) && iseven(Δc) + if iseven(Δa) && iseven(Δb) && iseven(Δc) && (Δa > 0) && (Δb > 0) && (Δc > 0) M = qr(P.X^(Δa÷2)*(I-P.X)^(Δb÷2)*(Qt*I-P.X)^(Δc÷2)).R return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) elseif isone(Δa) && iszero(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (1,0,0) @@ -280,7 +288,16 @@ function semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) elseif iszero(Δa) && iszero(Δb) && isone(Δc) # special case (Δa,Δb,Δc) = (0,0,1) M = cholesky(Qt*I-P.X).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) - elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) + elseif isone(-Δa) && iszero(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (-1,0,0) + M = cholesky(semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U + return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) + elseif iszero(Δa) && isone(-Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (0,-1,0) + M = cholesky(I-semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U + return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) + elseif iszero(Δa) && iszero(Δb) && isone(-Δc) # special case (Δa,Δb,Δc) = (0,0,-1) + M = cholesky(Qt*I-semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U + return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) + elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa > 0) && (Δb > 0) && (Δc > 0) M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Qt*I-P.X)^(Δc))).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1 else @@ -297,14 +314,20 @@ function semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) Δb = Qb-P.b Δc = Qc-P.c if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively - if ((isone(Δa)||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(Δb)||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(Δc)||isone(Δc/2))) + if ((isone(abs(Δa))||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(abs(Δb))||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(abs(Δc))||isone(Δc/2))) M = semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P) - elseif Δa > 0 # iterative modification by 1 + elseif Δa > 0 # iterative raising by 1 M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa-1-iseven(Δa), Qb, Qc, P)),semijacobi_ldiv(Qt, Qa-1-iseven(Δa), Qb, Qc, P)) elseif Δb > 0 M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb-1-iseven(Δb), Qc, P)),semijacobi_ldiv(Qt, Qa, Qb-1-iseven(Δb), Qc, P)) elseif Δc > 0 M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb, Qc-1-iseven(Δc), P)),semijacobi_ldiv(Qt, Qa, Qb, Qc-1-iseven(Δc), P)) + elseif Δa < 0 # iterative lowering by 1 + M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa+1, Qb, Qc, P)),semijacobi_ldiv(Qt, Qa+1, Qb, Qc, P)) + elseif Δb < 0 + M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb+1, Qc, P)),semijacobi_ldiv(Qt, Qa, Qb+1, Qc, P)) + elseif Δc < 0 + M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb, Qc+1, P)),semijacobi_ldiv(Qt, Qa, Qb, Qc+1, P)) end else # fallback to Lancos R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1)) diff --git a/test/runtests.jl b/test/runtests.jl index 2791b13..7f1073d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -619,5 +619,24 @@ end @test sprint(show, W) == "Weighted(SemiclassicalJacobi with weight x^2.3 * (1-x)^5.3 * (2.0-x)^0.4 on 0..1)" end +@testset "basic lowering" begin + # lowering + t = 2.2 + P = SemiclassicalJacobi(t, 2, 2, 2) + R_0 = SemiclassicalJacobi(t, 1, 2, 2) \ P + R_1 = SemiclassicalJacobi(t, 2, 1, 2) \ P + R_t = SemiclassicalJacobi(t, 2, 2, 1) \ P + + R_0inv = P \ SemiclassicalJacobi(t, 1, 2, 2) + R_1inv = P \ SemiclassicalJacobi(t, 2, 1, 2) + R_tinv = P \ SemiclassicalJacobi(t, 2, 2, 1) + + @test ApplyArray(inv,R_0inv)[1:20,1:20] ≈ R_0[1:20,1:20] + @test ApplyArray(inv,R_1inv)[1:20,1:20] ≈ R_1[1:20,1:20] + @test ApplyArray(inv,R_tinv)[1:20,1:20] ≈ R_t[1:20,1:20] +end + include("test_derivative.jl") -include("test_neg1c.jl") \ No newline at end of file +include("test_neg1c.jl") + + \ No newline at end of file From 04109ccbddd9e520cf8bdb68d26a10e16d95bba9 Mon Sep 17 00:00:00 2001 From: "Timon S. Gutleb" Date: Fri, 21 Jun 2024 12:17:51 -0700 Subject: [PATCH 02/10] fix bugs, add QR lowering too --- src/SemiclassicalOrthogonalPolynomials.jl | 46 +++++++++++++++----- test/runtests.jl | 52 +++++++++++++++-------- 2 files changed, 70 insertions(+), 28 deletions(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 1c92b59..18bedc8 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -275,31 +275,57 @@ function semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) Δa = Qa-P.a Δb = Qb-P.b Δc = Qc-P.c - M = Diagonal(Ones(∞)) - if iseven(Δa) && iseven(Δb) && iseven(Δc) && (Δa > 0) && (Δb > 0) && (Δc > 0) - M = qr(P.X^(Δa÷2)*(I-P.X)^(Δb÷2)*(Qt*I-P.X)^(Δc÷2)).R + # special case (Δa,Δb,Δc) = (2,0,0) + if isone(Δa÷2) && iszero(Δb) && iszero(Δc) + M = qr(P.X^(Δa÷2)).R + return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # match normalization choice P_0(x) = 1 + # special case (Δa,Δb,Δc) = (0,2,0) + elseif iszero(Δa) && isone(Δb÷2) && iszero(Δc) + M = qr((I-P.X)^(Δb÷2)).R return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) - elseif isone(Δa) && iszero(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (1,0,0) + # special case (Δa,Δb,Δc) = (0,0,2) + elseif iszero(Δa) && iszero(Δb) && isone(Δc÷2) + M = qr((Qt*I-P.X)^(Δc÷2)).R + return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) + # special case (Δa,Δb,Δc) = (-2,0,0) + elseif isone(Δa÷2) && iszero(Δb) && iszero(Δc) + M = qr(P.X^(Δa÷2)).R + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) + # special case (Δa,Δb,Δc) = (0,-2,0) + elseif iszero(Δa) && isone(Δb÷2) && iszero(Δc) + M = qr((I-P.X)^(Δb÷2)).R + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) + # special case (Δa,Δb,Δc) = (0,0,-2) + elseif iszero(Δa) && iszero(Δb) && isone(Δc÷2) + M = qr((Qt*I-P.X)^(Δc÷2)).R + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) + # special case (Δa,Δb,Δc) = (1,0,0) + elseif isone(Δa) && iszero(Δb) && iszero(Δc) M = cholesky(P.X).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) - elseif iszero(Δa) && isone(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (0,1,0) + # special case (Δa,Δb,Δc) = (0,1,0) + elseif iszero(Δa) && isone(Δb) && iszero(Δc) M = cholesky(I-P.X).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) - elseif iszero(Δa) && iszero(Δb) && isone(Δc) # special case (Δa,Δb,Δc) = (0,0,1) + # special case (Δa,Δb,Δc) = (0,0,1) + elseif iszero(Δa) && iszero(Δb) && isone(Δc) M = cholesky(Qt*I-P.X).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) - elseif isone(-Δa) && iszero(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (-1,0,0) + # special case (Δa,Δb,Δc) = (-1,0,0) + elseif isone(-Δa) && iszero(Δb) && iszero(Δc) M = cholesky(semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) - elseif iszero(Δa) && isone(-Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (0,-1,0) + # special case (Δa,Δb,Δc) = (0,-1,0) + elseif iszero(Δa) && isone(-Δb) && iszero(Δc) M = cholesky(I-semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) - elseif iszero(Δa) && iszero(Δb) && isone(-Δc) # special case (Δa,Δb,Δc) = (0,0,-1) + # special case (Δa,Δb,Δc) = (0,0,-1) + elseif iszero(Δa) && iszero(Δb) && isone(-Δc) M = cholesky(Qt*I-semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa > 0) && (Δb > 0) && (Δc > 0) M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Qt*I-P.X)^(Δc))).U - return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1 + return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) else error("Implement") end diff --git a/test/runtests.jl b/test/runtests.jl index 7f1073d..01db9f1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -619,24 +619,40 @@ end @test sprint(show, W) == "Weighted(SemiclassicalJacobi with weight x^2.3 * (1-x)^5.3 * (2.0-x)^0.4 on 0..1)" end -@testset "basic lowering" begin - # lowering - t = 2.2 - P = SemiclassicalJacobi(t, 2, 2, 2) - R_0 = SemiclassicalJacobi(t, 1, 2, 2) \ P - R_1 = SemiclassicalJacobi(t, 2, 1, 2) \ P - R_t = SemiclassicalJacobi(t, 2, 2, 1) \ P - - R_0inv = P \ SemiclassicalJacobi(t, 1, 2, 2) - R_1inv = P \ SemiclassicalJacobi(t, 2, 1, 2) - R_tinv = P \ SemiclassicalJacobi(t, 2, 2, 1) - - @test ApplyArray(inv,R_0inv)[1:20,1:20] ≈ R_0[1:20,1:20] - @test ApplyArray(inv,R_1inv)[1:20,1:20] ≈ R_1[1:20,1:20] - @test ApplyArray(inv,R_tinv)[1:20,1:20] ≈ R_t[1:20,1:20] +@testset "Lowering" begin + @testset "Decrements of 1 via Cholesky" begin + # Cholesky lowering + t = 2.2 + P = SemiclassicalJacobi(t, 2, 2, 2) + R_0 = SemiclassicalJacobi(t, 1, 2, 2) \ P + R_1 = SemiclassicalJacobi(t, 2, 1, 2) \ P + R_t = SemiclassicalJacobi(t, 2, 2, 1) \ P + + R_0inv = P \ SemiclassicalJacobi(t, 1, 2, 2) + R_1inv = P \ SemiclassicalJacobi(t, 2, 1, 2) + R_tinv = P \ SemiclassicalJacobi(t, 2, 2, 1) + + @test ApplyArray(inv,R_0inv)[1:20,1:20] ≈ R_0[1:20,1:20] + @test ApplyArray(inv,R_1inv)[1:20,1:20] ≈ R_1[1:20,1:20] + @test ApplyArray(inv,R_tinv)[1:20,1:20] ≈ R_t[1:20,1:20] + end + @testset "Decrements of 2 via QR" begin + # Cholesky lowering + t = 2.2 + P = SemiclassicalJacobi(t, 3, 3, 3) + R_0 = SemiclassicalJacobi(t, 1, 3, 3) \ P + R_1 = SemiclassicalJacobi(t, 3, 1, 3) \ P + R_t = SemiclassicalJacobi(t, 3, 3, 1) \ P + + R_0inv = P \ SemiclassicalJacobi(t, 1, 3, 3) + R_1inv = P \ SemiclassicalJacobi(t, 3, 1, 3) + R_tinv = P \ SemiclassicalJacobi(t, 3, 3, 1) + + @test ApplyArray(inv,R_0inv)[1:20,1:20] ≈ R_0[1:20,1:20] + @test ApplyArray(inv,R_1inv)[1:20,1:20] ≈ R_1[1:20,1:20] + @test ApplyArray(inv,R_tinv)[1:20,1:20] ≈ R_t[1:20,1:20] + end end include("test_derivative.jl") -include("test_neg1c.jl") - - \ No newline at end of file +include("test_neg1c.jl") \ No newline at end of file From d77c8981e3c75bfbefae542277dc0f533b78c94c Mon Sep 17 00:00:00 2001 From: "Timon S. Gutleb" Date: Fri, 21 Jun 2024 16:25:52 -0700 Subject: [PATCH 03/10] change convention, implement lowering --- src/SemiclassicalOrthogonalPolynomials.jl | 98 +++++++++++++---------- src/deprecated.jl | 3 + 2 files changed, 59 insertions(+), 42 deletions(-) create mode 100644 src/deprecated.jl diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 18bedc8..5ae882c 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -180,7 +180,7 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c) cholesky_jacobimatrix(I-Q.X,Q) elseif iszero(Δa) && iszero(Δb) && isone(Δc) cholesky_jacobimatrix(Q.t*I-Q.X,Q) - elseif isone(-Δa) && iszero(Δb) && iszero(Δc) # raising by 1 + elseif isone(-Δa) && iszero(Δb) && iszero(Δc) # in these cases we currently have to reconstruct # TODO: This is re-constructing. It should instead use reverse Cholesky (or an alternative)! semiclassical_jacobimatrix(Q.t,a,b,c) elseif iszero(Δa) && isone(-Δb) && iszero(Δc) @@ -214,7 +214,6 @@ LanczosPolynomial(P::SemiclassicalJacobi{T}) where T = gives either a mapped `Jacobi` or `LanczosPolynomial` version of `P`. """ -# TODO: Use ConvertedOPs for integer special cases? toclassical(P::SemiclassicalJacobi{T}) where T = iszero(P.c) ? Normalized(jacobi(P.b, P.a, UnitInterval{T}())) : LanczosPolynomial(P) copy(P::SemiclassicalJacobi) = P @@ -239,6 +238,7 @@ jacobimatrix(P::SemiclassicalJacobi) = P.X """ op_lowering(Q, y) + Gives the Lowering operator from OPs w.r.t. (x-y)*w(x) to Q as constructed from Chistoffel–Darboux """ @@ -268,36 +268,40 @@ function semijacobi_ldiv(P::SemiclassicalJacobi, Q) (P \ R) * _p0(R̃) * (R̃ \ Q) end -# returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q in a single step via decomposition. -semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) = semijacobi_ldiv_direct(Q.t, Q.a, Q.b, Q.c, P) -function semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) - (Qt ≈ P.t) && (Qa ≈ P.a) && (Qb ≈ P.b) && (Qc ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞)) - Δa = Qa-P.a - Δb = Qb-P.b - Δc = Qc-P.c +""" + semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) + +Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q` in a single step via decomposition. +Numerically unstable if the parameter modification is large. Typically one should instead use `P \\ Q` which is equivalent to `semijacobi_ldiv(P,Q)` and proceeds step by step. +""" +function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) + (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞)) + Δa = Q.a-P.a + Δb = Q.b-P.b + Δc = Q.c-P.c # special case (Δa,Δb,Δc) = (2,0,0) - if isone(Δa÷2) && iszero(Δb) && iszero(Δc) + if isone(Δa/2) && iszero(Δb) && iszero(Δc) M = qr(P.X^(Δa÷2)).R return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # match normalization choice P_0(x) = 1 # special case (Δa,Δb,Δc) = (0,2,0) - elseif iszero(Δa) && isone(Δb÷2) && iszero(Δc) + elseif iszero(Δa) && isone(Δb/2) && iszero(Δc) M = qr((I-P.X)^(Δb÷2)).R return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # special case (Δa,Δb,Δc) = (0,0,2) - elseif iszero(Δa) && iszero(Δb) && isone(Δc÷2) - M = qr((Qt*I-P.X)^(Δc÷2)).R + elseif iszero(Δa) && iszero(Δb) && isone(Δc/2) + M = qr((Q.t*I-P.X)^(Δc÷2)).R return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # special case (Δa,Δb,Δc) = (-2,0,0) - elseif isone(Δa÷2) && iszero(Δb) && iszero(Δc) - M = qr(P.X^(Δa÷2)).R + elseif isone(-Δa/2) && iszero(Δb) && iszero(Δc) + M = qr(Q.X^(-Δa÷2)).R return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) # special case (Δa,Δb,Δc) = (0,-2,0) - elseif iszero(Δa) && isone(Δb÷2) && iszero(Δc) - M = qr((I-P.X)^(Δb÷2)).R + elseif iszero(Δa) && isone(-Δb/2) && iszero(Δc) + M = qr((I-Q.X)^(-Δb÷2)).R return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) # special case (Δa,Δb,Δc) = (0,0,-2) - elseif iszero(Δa) && iszero(Δb) && isone(Δc÷2) - M = qr((Qt*I-P.X)^(Δc÷2)).R + elseif iszero(Δa) && iszero(Δb) && isone(-Δc/2) + M = qr((Q.t*I-Q.X)^(-Δc÷2)).R return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) # special case (Δa,Δb,Δc) = (1,0,0) elseif isone(Δa) && iszero(Δb) && iszero(Δc) @@ -309,56 +313,65 @@ function semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # special case (Δa,Δb,Δc) = (0,0,1) elseif iszero(Δa) && iszero(Δb) && isone(Δc) - M = cholesky(Qt*I-P.X).U + M = cholesky(Q.t*I-P.X).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # special case (Δa,Δb,Δc) = (-1,0,0) elseif isone(-Δa) && iszero(Δb) && iszero(Δc) - M = cholesky(semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U + M = cholesky(Q.X).U return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) # special case (Δa,Δb,Δc) = (0,-1,0) elseif iszero(Δa) && isone(-Δb) && iszero(Δc) - M = cholesky(I-semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U + M = cholesky(I-Q.X).U return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) # special case (Δa,Δb,Δc) = (0,0,-1) elseif iszero(Δa) && iszero(Δb) && isone(-Δc) - M = cholesky(Qt*I-semiclassical_jacobimatrix(P, Qa, Qb, Qc)).U + M = cholesky(Q.t*I-Q.X).U return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) - elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa > 0) && (Δb > 0) && (Δc > 0) - M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Qt*I-P.X)^(Δc))).U + elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa >= 0) && (Δb >= 0) && (Δc >= 0) + M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Q.t*I-P.X)^(Δc))).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) else - error("Implement") + error("Implement modification by ($Δa,$Δb,$Δc)") end end -# returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q. -semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) = semijacobi_ldiv(Q.t, Q.a, Q.b, Q.c, P) -function semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) - @assert Qt ≈ P.t - (Qt ≈ P.t) && (Qa ≈ P.a) && (Qb ≈ P.b) && (Qc ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞)) - Δa = Qa-P.a - Δb = Qb-P.b - Δc = Qc-P.c +""" + semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) + +Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q`. Integer distances are covered by decomposition methods, for non-integer cases a Lanczos fallback is attempted. +""" +function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) + @assert Q.t ≈ P.t + (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞)) + Δa = Q.a-P.a + Δb = Q.b-P.b + Δc = Q.c-P.c if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively if ((isone(abs(Δa))||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(abs(Δb))||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(abs(Δc))||isone(Δc/2))) - M = semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P) + return semijacobi_ldiv_direct(Q, P) elseif Δa > 0 # iterative raising by 1 - M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa-1-iseven(Δa), Qb, Qc, P)),semijacobi_ldiv(Qt, Qa-1-iseven(Δa), Qb, Qc, P)) + QQ = SemiclassicalJacobi(Q.t, Q.a-1-iseven(Δa), Q.b, Q.c, P) + return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P)) elseif Δb > 0 - M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb-1-iseven(Δb), Qc, P)),semijacobi_ldiv(Qt, Qa, Qb-1-iseven(Δb), Qc, P)) + QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b-1-iseven(Δb), Q.c, P) + return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P)) elseif Δc > 0 - M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb, Qc-1-iseven(Δc), P)),semijacobi_ldiv(Qt, Qa, Qb, Qc-1-iseven(Δc), P)) + QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1-iseven(Δc), P) + return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P)) elseif Δa < 0 # iterative lowering by 1 - M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa+1, Qb, Qc, P)),semijacobi_ldiv(Qt, Qa+1, Qb, Qc, P)) + QQ = SemiclassicalJacobi(Q.t, Q.a+1+iseven(Δa), Q.b, Q.c, P) + return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P)) elseif Δb < 0 - M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb+1, Qc, P)),semijacobi_ldiv(Qt, Qa, Qb+1, Qc, P)) + QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b+1+iseven(Δb), Q.c, P) + return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P)) elseif Δc < 0 - M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb, Qc+1, P)),semijacobi_ldiv(Qt, Qa, Qb, Qc+1, P)) + QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1+iseven(Δc), P) + return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P)) end else # fallback to Lancos R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1)) R̃ = toclassical(R) - return (P \ R) * _p0(R̃) * (R̃ \ SemiclassicalJacobi(Qt, Qa, Qb, Qc)) + return (P \ R) * _p0(R̃) * (R̃ \ Q) end end @@ -645,5 +658,6 @@ function sumquotient(wP::SemiclassicalJacobiWeight{T},wQ::SemiclassicalJacobiWei end include("neg1c.jl") +include("deprecated.jl") end diff --git a/src/deprecated.jl b/src/deprecated.jl new file mode 100644 index 0000000..a87a92f --- /dev/null +++ b/src/deprecated.jl @@ -0,0 +1,3 @@ +# The old ldiv variants are still supported for now but deprecated and implemented using non-efficient constructors +Base.@deprecate semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) semijacobi_ldiv_direct(SemiclassicalJacobi(Qt, Qa, Qb, Qc), P) +Base.@deprecate semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) semijacobi_ldiv(SemiclassicalJacobi(Qt, Qa, Qb, Qc), P) From dc70b6b63a54b00c9522426781d3b05f4d2675c2 Mon Sep 17 00:00:00 2001 From: "Timon S. Gutleb" Date: Sun, 23 Jun 2024 03:15:08 -0700 Subject: [PATCH 04/10] X^1 = X --- src/SemiclassicalOrthogonalPolynomials.jl | 42 +++++++++++------------ 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 5ae882c..025de65 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -144,7 +144,7 @@ function semiclassical_jacobimatrix(t, a, b, c) iszero(c) && return jacobimatrix(P) if isone(c) return cholesky_jacobimatrix(Symmetric(P \ ((t.-axes(P,1)).*P)), P) - elseif isone(c/2) + elseif c == 2 return qr_jacobimatrix(Symmetric(P \ ((t.-axes(P,1)).*P)), P) elseif isinteger(c) && c ≥ 0 # reduce other integer c cases to hierarchy return SemiclassicalJacobi.(t, a, b, 0:Int(c))[end].X @@ -280,29 +280,29 @@ function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) Δb = Q.b-P.b Δc = Q.c-P.c # special case (Δa,Δb,Δc) = (2,0,0) - if isone(Δa/2) && iszero(Δb) && iszero(Δc) - M = qr(P.X^(Δa÷2)).R - return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # match normalization choice P_0(x) = 1 + if (Δa == 2) && iszero(Δb) && iszero(Δc) + M = qr(P.X).R + return ApplyArray(*, Diagonal(sign.(view(M,band(0)))/abs(M[1])), M) # special case (Δa,Δb,Δc) = (0,2,0) - elseif iszero(Δa) && isone(Δb/2) && iszero(Δc) - M = qr((I-P.X)^(Δb÷2)).R - return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) + elseif iszero(Δa) && (Δb == 2) && iszero(Δc) + M = qr(I-P.X).R + return ApplyArray(*, Diagonal(sign.(view(M,band(0)))/abs(M[1])), M) # special case (Δa,Δb,Δc) = (0,0,2) - elseif iszero(Δa) && iszero(Δb) && isone(Δc/2) - M = qr((Q.t*I-P.X)^(Δc÷2)).R - return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) + elseif iszero(Δa) && iszero(Δb) && (Δc == 2) + M = qr(Q.t*I-P.X).R + return ApplyArray(*, Diagonal(sign.(view(M,band(0)))/abs(M[1])), M) # special case (Δa,Δb,Δc) = (-2,0,0) - elseif isone(-Δa/2) && iszero(Δb) && iszero(Δc) - M = qr(Q.X^(-Δa÷2)).R - return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) + elseif (Δa == -2) && iszero(Δb) && iszero(Δc) + M = qr(Q.X).R + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs.(M[1]))) # special case (Δa,Δb,Δc) = (0,-2,0) - elseif iszero(Δa) && isone(-Δb/2) && iszero(Δc) - M = qr((I-Q.X)^(-Δb÷2)).R - return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) + elseif iszero(Δa) && (Δb == -2) && iszero(Δc) + M = qr(I-Q.X).R + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs.(M[1]))) # special case (Δa,Δb,Δc) = (0,0,-2) - elseif iszero(Δa) && iszero(Δb) && isone(-Δc/2) - M = qr((Q.t*I-Q.X)^(-Δc÷2)).R - return ApplyArray(\, M, Diagonal(sign.(view(M,band(0))).*Fill(abs.(M[1]),∞))) + elseif iszero(Δa) && iszero(Δb) && (Δc == -2) + M = qr(Q.t*I-Q.X).R + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs.(M[1]))) # special case (Δa,Δb,Δc) = (1,0,0) elseif isone(Δa) && iszero(Δb) && iszero(Δc) M = cholesky(P.X).U @@ -327,7 +327,7 @@ function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) elseif iszero(Δa) && iszero(Δb) && isone(-Δc) M = cholesky(Q.t*I-Q.X).U return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) - elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa >= 0) && (Δb >= 0) && (Δc >= 0) + elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa ≥ 0) && (Δb ≥ 0) && (Δc ≥ 0) M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Q.t*I-P.X)^(Δc))).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) else @@ -347,7 +347,7 @@ function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) Δb = Q.b-P.b Δc = Q.c-P.c if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively - if ((isone(abs(Δa))||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(abs(Δb))||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(abs(Δc))||isone(Δc/2))) + if ((isone(abs(Δa))||(Δa == 2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(abs(Δb))||(Δb == 2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(abs(Δc))||(Δc == 2))) return semijacobi_ldiv_direct(Q, P) elseif Δa > 0 # iterative raising by 1 QQ = SemiclassicalJacobi(Q.t, Q.a-1-iseven(Δa), Q.b, Q.c, P) From 5279d290628078243bb02232a813aa33f0d5f8fe Mon Sep 17 00:00:00 2001 From: Timon Salar Gutleb Date: Mon, 24 Jun 2024 07:38:04 -0700 Subject: [PATCH 05/10] Update src/SemiclassicalOrthogonalPolynomials.jl Co-authored-by: Sheehan Olver --- src/SemiclassicalOrthogonalPolynomials.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 025de65..29d7f9c 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -342,7 +342,8 @@ Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi """ function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) @assert Q.t ≈ P.t - (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞)) + T = promote_type(eltype(Q), eltype(P)) + (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SquareEye{T}(∞) Δa = Q.a-P.a Δb = Q.b-P.b Δc = Q.c-P.c From d0b95ddc7e3bfc0c70192b9bee3525ef58e9957e Mon Sep 17 00:00:00 2001 From: Timon Salar Gutleb Date: Mon, 24 Jun 2024 07:38:15 -0700 Subject: [PATCH 06/10] Update src/SemiclassicalOrthogonalPolynomials.jl Co-authored-by: Sheehan Olver --- src/SemiclassicalOrthogonalPolynomials.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 29d7f9c..1796c6d 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -314,7 +314,7 @@ function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) # special case (Δa,Δb,Δc) = (0,0,1) elseif iszero(Δa) && iszero(Δb) && isone(Δc) M = cholesky(Q.t*I-P.X).U - return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) + return M/M[1] # special case (Δa,Δb,Δc) = (-1,0,0) elseif isone(-Δa) && iszero(Δb) && iszero(Δc) M = cholesky(Q.X).U From 24d7dd12b8c05f72708eb6f995b3dba6bf5b1d64 Mon Sep 17 00:00:00 2001 From: "Timon S. Gutleb" Date: Mon, 24 Jun 2024 07:39:47 -0700 Subject: [PATCH 07/10] consistency change --- src/SemiclassicalOrthogonalPolynomials.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 1796c6d..11c5d22 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -294,23 +294,23 @@ function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) # special case (Δa,Δb,Δc) = (-2,0,0) elseif (Δa == -2) && iszero(Δb) && iszero(Δc) M = qr(Q.X).R - return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs.(M[1]))) + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs(M[1]))) # special case (Δa,Δb,Δc) = (0,-2,0) elseif iszero(Δa) && (Δb == -2) && iszero(Δc) M = qr(I-Q.X).R - return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs.(M[1]))) + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs(M[1]))) # special case (Δa,Δb,Δc) = (0,0,-2) elseif iszero(Δa) && iszero(Δb) && (Δc == -2) M = qr(Q.t*I-Q.X).R - return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs.(M[1]))) + return ApplyArray(\, M, Diagonal(sign.(view(M,band(0)))*abs(M[1]))) # special case (Δa,Δb,Δc) = (1,0,0) elseif isone(Δa) && iszero(Δb) && iszero(Δc) M = cholesky(P.X).U - return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) + return M/M[1] # special case (Δa,Δb,Δc) = (0,1,0) elseif iszero(Δa) && isone(Δb) && iszero(Δc) M = cholesky(I-P.X).U - return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) + return M/M[1] # special case (Δa,Δb,Δc) = (0,0,1) elseif iszero(Δa) && iszero(Δb) && isone(Δc) M = cholesky(Q.t*I-P.X).U From 6bafbd9113b8a0542bf88ed3404363ef5724b6d2 Mon Sep 17 00:00:00 2001 From: "Timon S. Gutleb" Date: Mon, 24 Jun 2024 07:49:04 -0700 Subject: [PATCH 08/10] bugfix --- src/SemiclassicalOrthogonalPolynomials.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 11c5d22..6541f2d 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -275,7 +275,7 @@ Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi Numerically unstable if the parameter modification is large. Typically one should instead use `P \\ Q` which is equivalent to `semijacobi_ldiv(P,Q)` and proceeds step by step. """ function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) - (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞)) + (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SquareEye{eltype(Q.t)}((axes(P,2),)) Δa = Q.a-P.a Δb = Q.b-P.b Δc = Q.c-P.c @@ -318,15 +318,15 @@ function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) # special case (Δa,Δb,Δc) = (-1,0,0) elseif isone(-Δa) && iszero(Δb) && iszero(Δc) M = cholesky(Q.X).U - return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) + return UpperTriangular(ApplyArray(inv,M/M[1])) # special case (Δa,Δb,Δc) = (0,-1,0) elseif iszero(Δa) && isone(-Δb) && iszero(Δc) M = cholesky(I-Q.X).U - return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) + return UpperTriangular(ApplyArray(inv,M/M[1])) # special case (Δa,Δb,Δc) = (0,0,-1) elseif iszero(Δa) && iszero(Δb) && isone(-Δc) M = cholesky(Q.t*I-Q.X).U - return ApplyArray(\, M, Diagonal(Fill(M[1],∞))) + return UpperTriangular(ApplyArray(inv,M/M[1])) elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) && (Δa ≥ 0) && (Δb ≥ 0) && (Δc ≥ 0) M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Q.t*I-P.X)^(Δc))).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) @@ -343,7 +343,7 @@ Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) @assert Q.t ≈ P.t T = promote_type(eltype(Q), eltype(P)) - (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return SquareEye{T}(∞) + (Q.t ≈ P.t) && (Q.a ≈ P.a) && (Q.b ≈ P.b) && (Q.c ≈ P.c) && return return SquareEye{eltype(Q.t)}((axes(P,2),)) Δa = Q.a-P.a Δb = Q.b-P.b Δc = Q.c-P.c From e3be3fbffbd733878b25c75919b9435545fbeef2 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 24 Jun 2024 16:34:53 +0100 Subject: [PATCH 09/10] Update ci.yml --- .github/workflows/ci.yml | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 588a137..6108091 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -20,7 +20,7 @@ jobs: fail-fast: false matrix: version: - - '1' + - '1.10' os: - ubuntu-latest - macOS-latest @@ -28,12 +28,12 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v1 + - uses: actions/cache@v4 env: cache-name: cache-artifacts with: @@ -46,6 +46,7 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v4 with: + token: ${{ secrets.CODECOV_TOKEN }} file: lcov.info From f4e9683cc1aa5402f4592f69fdf1f3ff32f04384 Mon Sep 17 00:00:00 2001 From: "Timon S. Gutleb" Date: Mon, 24 Jun 2024 09:00:53 -0700 Subject: [PATCH 10/10] increase test coverage --- test/runtests.jl | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 01db9f1..5a141b8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -620,6 +620,25 @@ end end @testset "Lowering" begin + @testset "Constructing Lowered Polynomials" begin + t = 1.3184 + P = SemiclassicalJacobi(t, 3, 3, 3) + Q_1 = SemiclassicalJacobi(t, 2, 3, 3, P) + Q_2 = SemiclassicalJacobi(t, 3, 3, 2, P) + Q_3 = SemiclassicalJacobi(t, 3, 2, 3, P) + Q_4 = SemiclassicalJacobi(t, 1, 2, 3, P) + Q_5 = SemiclassicalJacobi(t, 3, 2, 1, P) + Q_6 = SemiclassicalJacobi(t, 2, 2, 2, P) + Q_7 = SemiclassicalJacobi(t, 1, 1, 1, P) + + @test Q_1.X[1:20,1:20] ≈ SemiclassicalJacobi(t, 2, 3, 3).X[1:20,1:20] + @test Q_2.X[1:20,1:20] ≈ SemiclassicalJacobi(t, 3, 3, 2).X[1:20,1:20] + @test Q_3.X[1:20,1:20] ≈ SemiclassicalJacobi(t, 3, 2, 3).X[1:20,1:20] + @test Q_4.X[1:20,1:20] ≈ SemiclassicalJacobi(t, 1, 2, 3).X[1:20,1:20] + @test Q_5.X[1:20,1:20] ≈ SemiclassicalJacobi(t, 3, 2, 1).X[1:20,1:20] + @test Q_6.X[1:20,1:20] ≈ SemiclassicalJacobi(t, 2, 2, 2).X[1:20,1:20] + @test Q_7.X[1:20,1:20] ≈ SemiclassicalJacobi(t, 1, 1, 1).X[1:20,1:20] + end @testset "Decrements of 1 via Cholesky" begin # Cholesky lowering t = 2.2 @@ -655,4 +674,5 @@ end end include("test_derivative.jl") -include("test_neg1c.jl") \ No newline at end of file +include("test_neg1c.jl") +