Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lowering connection cfs from raising connection cfs #112

Merged
merged 10 commits into from
Jun 24, 2024

Conversation

TSGut
Copy link
Member

@TSGut TSGut commented Jun 21, 2024

Not sure whether we will end up merging this or your PR @DanielVandH but the point of this is to show how a lowering implementation would look. It actually works already but it's not efficient.

julia> @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
Test Summary:  | Pass  Total  Time
basic lowering |    3      3  0.0s

The reason for the inefficiency lies entirely with three lines where I 'cheat'. The iterative connection coefficients above are computed correctly, that's not the problem. The problem is in the computation we need Q.X which the user could either supply always (this would correspond to @dlfivefifty 's suggestion of changing all conventions to be an semiclassical_ldiv(Q,X) format or alternatively we can implement the reverse Cholesky Jacobi matrix which would go in the big marked by TODO section, i.e. this bit

 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)

If the above was reverse Cholesky jacobimatrix then this PR would do everything you want I think. Without that, we have a design decision to make.

@TSGut
Copy link
Member Author

TSGut commented Jun 21, 2024

Guess something else broke as a result, I will fix that later today.

@DanielVandH
Copy link
Member

Thanks a lot @TSGut, this seems to be a better approach to it.

The problem is in the computation we need Q.X which the user could either supply always (this would correspond to @dlfivefifty 's suggestion of changing all conventions to be an semiclassical_ldiv(Q,X) format or alternatively we can implement the reverse Cholesky Jacobi matrix

I think it would be easier to work with just being able to keep Q around. Do you have a preference? Maybe makes this easier in the future too if ever we need to come back to these functions and again have a need for Q.

@TSGut
Copy link
Member Author

TSGut commented Jun 21, 2024

I think either is fine, there are some cost savings here and there depending on the use case with either approach but probably all negligible.

@TSGut TSGut changed the title WIP: How to make lowering work WIP: Lowering connection cfs from raising connection cfs Jun 21, 2024
@TSGut
Copy link
Member Author

TSGut commented Jun 22, 2024

Alright, tests passing (though I may add more).

In this package I am usually happy to just merge things myself but since you guys are actively developing with this, I would prefer not to mess unless this is wanted, so please take a look @dlfivefifty and @DanielVandH.


Changes and new features:

  • Lowering implemented by inverting raising operator assuming Q and P are in hand (which they ought to be, else how are you even calling the \ function...).

  • There is still a fallback option that reconstructs marked with TODO, I think to remove that essentially requires reverse Cholesky / QL based jacobimatrices (which would live in ClassicalOPs.jl I guess). I am happy to eventually get around to that but it may take me some time since doing that right will be a bit of work. Once that is done it just plugs into the TODO segments here. For now this should be fine in 99% of applications but just want to point out there is a lot to be gained there.

  • Changed convention to follow semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi), old version deprecated for now (mainly to give John time to test whether this breaks something for his Annuli package)

  • Some additional docstrings

  • Some clean-up based on discussions in PR Support computing connection matrices when decrementing parameters #109, attempt to make things more readable


src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
@TSGut
Copy link
Member Author

TSGut commented Jun 23, 2024

Any idea why codecov went and died on this repo? It would be nice to add more tests to this but it isn't being displayed and manually checking https://app.codecov.io/gh/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/pull/112 shows some error in the config that started a few PRs ago

@dlfivefifty
Copy link
Member

At some point I had to add keys for codexov, probably I haven’t done this yet. I’ll check when I’m home

src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
# 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],∞)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
return ApplyArray(\, M, Diagonal(Fill(M[1],∞)))
return inv(M/M[1])

Copy link
Member Author

@TSGut TSGut Jun 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

inv doesn't work when applied to this type. ApplyArray(inv,M/M[1]) works similarly to ApplyArray(\, M, Diagonal(Fill(M[1],∞))) though Julia doesn't seem to know it's Triangular so it causes some odd behavior here and there. This feels like something to fix elsewhere and then we can make it just read inv(M/M[1]) here? Altrenatively I can always just return UpperTriangular(ApplyArray(inv,M/M[1]) for now...

For the record, with ApplyArray(inv,M/M[1]) we get

MemoryLayout(R_0)
LazyArrays.InvLayout{TriangularLayout{'U', 'N', LazyArrays.LazyLayout}}()

with inv(M/M[1]) it errors:

julia>         R_0 = SemiclassicalJacobi(t, 1, 2, 2) \ P
ERROR: MethodError: no method matching Matrix{Float64}(::UniformScaling{Bool}, ::Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}})

Closest candidates are:
  Matrix{T}(::UndefInitializer, ::Tuple{Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}, Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}}) where T
   @ InfiniteArrays C:\Users\timon\.julia\packages\InfiniteArrays\heAfT\src\infarrays.jl:7
  Matrix{T}(::UndefInitializer, ::Tuple{Integer, Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}}) where T
   @ InfiniteArrays C:\Users\timon\.julia\packages\InfiniteArrays\heAfT\src\infarrays.jl:4
  Matrix{T}(::UndefInitializer, ::Tuple{Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}, Integer}) where T
   @ InfiniteArrays C:\Users\timon\.julia\packages\InfiniteArrays\heAfT\src\infarrays.jl:5
  ...

Stacktrace:
 [1] inv(A::UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{…}}})
   @ LinearAlgebra C:\Users\timon\AppData\Local\Programs\Julia-1.10.0\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:797    
 [2] semijacobi_ldiv_direct(Q::SemiclassicalJacobi{Float64}, P::SemiclassicalJacobi{Float64})
   @ SemiclassicalOrthogonalPolynomials c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:322
 [3] semijacobi_ldiv(Q::SemiclassicalJacobi{Float64}, P::SemiclassicalJacobi{Float64})
   @ SemiclassicalOrthogonalPolynomials c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:352
 [4] \(A::SemiclassicalJacobi{Float64}, B::SemiclassicalJacobi{Float64})
   @ SemiclassicalOrthogonalPolynomials c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:411
 [5] top-level scope
   @ c:\Users\timon\OneDrive\Documents\Work projects\Archive\SemiclassicalOrthogonalPolynomials.jl\test\runtests.jl:663
Some type information was truncated. Use `show(err)` to see complete types.

src/SemiclassicalOrthogonalPolynomials.jl Outdated Show resolved Hide resolved
@dlfivefifty
Copy link
Member

codecov shouldb e fixded

@dlfivefifty
Copy link
Member

can you remove the WIP when this is ready?

Copy link

codecov bot commented Jun 24, 2024

Codecov Report

Attention: Patch coverage is 98.76543% with 1 line in your changes missing coverage. Please review.

Project coverage is 92.52%. Comparing base (5d48cc5) to head (f4e9683).
Report is 3 commits behind head on master.

Files Patch % Lines
src/SemiclassicalOrthogonalPolynomials.jl 98.76% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master     #112      +/-   ##
==========================================
+ Coverage   92.40%   92.52%   +0.11%     
==========================================
  Files           3        3              
  Lines         527      575      +48     
==========================================
+ Hits          487      532      +45     
- Misses         40       43       +3     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@dlfivefifty
Copy link
Member

also get the codecov up

@TSGut TSGut changed the title WIP: Lowering connection cfs from raising connection cfs Lowering connection cfs from raising connection cfs Jun 24, 2024
@TSGut
Copy link
Member Author

TSGut commented Jun 24, 2024

Ok I think this is ready if you're happy with it @dlfivefifty. Codecov is good and WIP tag removed. Version was already previously bumped.

Once reverse Cholesky jacobimatrix is implemented, which if nobody beats me to it I will get to this summer, I will also update this package to remove the fallback option.

I am excited to see the negative parameter stuff @DanielVandH! Let me know if I can help any more regarding this package.

@dlfivefifty dlfivefifty merged commit 52bbce9 into JuliaApproximation:master Jun 24, 2024
5 checks passed
@DanielVandH
Copy link
Member

Great! Thanks so much @TSGut for the help here.

@DanielVandH DanielVandH mentioned this pull request Jul 1, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants