Opened 3 years ago
Closed 3 years ago
#23175 closed enhancement (fixed)
Implement computation of Riemann period matrices etc.
Reported by:  nbruin  Owned by:  

Priority:  major  Milestone:  sage8.0 
Component:  algebraic geometry  Keywords:  sd86.5 
Cc:  aly.deines  Merged in:  
Authors:  Nils Bruin, Alexandre Zotine  Reviewers:  Julian Rüth 
Report Upstream:  N/A  Work issues:  
Branch:  6e689c2 (Commits)  Commit:  6e689c2aa87ec16fa6ce335e171aa4bce1885120 
Dependencies:  #23140  Stopgaps: 
Description
Include a class that supports analytic computation of period matrices and computation of endomorphism matrices.
Change History (32)
comment:1 Changed 3 years ago by
 Branch set to u/nbruin/riemann_surface
comment:2 Changed 3 years ago by
 Commit set to 3ffce6a66ed253e5592f88010cb98beddeca794c
 Dependencies set to #23140
 Keywords sd86.5 added
comment:3 Changed 3 years ago by
 Commit changed from 3ffce6a66ed253e5592f88010cb98beddeca794c to ce9fca4f7042629ad0e733e9bf87044cfa668597
Branch pushed to git repo; I updated commit sha1. New commits:
ce9fca4  some further imports

comment:4 Changed 3 years ago by
 Commit changed from ce9fca4f7042629ad0e733e9bf87044cfa668597 to 0143571215b3720f77a4f62db5d1e9e235758d08
Branch pushed to git repo; I updated commit sha1. New commits:
0143571  renamed file

comment:5 Changed 3 years ago by
 Commit changed from 0143571215b3720f77a4f62db5d1e9e235758d08 to c9edca9a447d5f9fcfa42f939c70cbaa59e126e1
comment:6 Changed 3 years ago by
 Commit changed from c9edca9a447d5f9fcfa42f939c70cbaa59e126e1 to c26f623ce8e5cf57cf8402d84a139220f75ea644
Branch pushed to git repo; I updated commit sha1. New commits:
c26f623  avoid integer division in python code

comment:7 Changed 3 years ago by
 Commit changed from c26f623ce8e5cf57cf8402d84a139220f75ea644 to 8f8dbb1d48834aa301bc8f6bd82b7c122ffcf890
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
8f8dbb1  avoid integer division in python code

comment:8 Changed 3 years ago by
 Commit changed from 8f8dbb1d48834aa301bc8f6bd82b7c122ffcf890 to 54bdd617f2bf7c6f177123bbdc7adffd5673854f
comment:9 Changed 3 years ago by
 Commit changed from 54bdd617f2bf7c6f177123bbdc7adffd5673854f to 5b1c2d61c9f88c175b8edb15ad69ca2307163ae3
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
5b1c2d6  Initial Checkin of RiemannSurface

comment:10 Changed 3 years ago by
 Cc aly.deines added
 Status changed from new to needs_review
Documentation should now be more or less up to standards.
comment:11 followup: ↓ 14 Changed 3 years ago by
 Reviewers set to Julian Rüth
 Status changed from needs_review to needs_work
I am making some style changes. I am still working on this, but here are already some comments.
Please fix/comment the following:
 I guess you want to make this more accessible than
from sage.schemes.riemann_surfaces.riemann_surface import RiemannSurface
?  Why so much precision in the examples? Having 200 bits of precision seems like a lot. I tried a couple of the examples with less precision and found the output much easier to parse.
TODO: have an example where self.genus < 0 to raise error
 Comment on the constant 100 in
_determine_new_w
and_newton_iteration
(nothing fancy necessary, just say that this is a random bound that worked well in practice?)  In
period_matrix
:Computes the period matrix of f given the homology basis and cohomology basis.
What is "given" here?  Why do you perform your own caching in
period_matrix
(self._PM
) andriemann_matrix
(self._RM
),cohomology_basis
(self._differentials
) ?  The
self._differentials[1]
is very confusing. Why do you need it? raise ValueError('Length of differentials list is not equal to genus.')
should this not be aNotImplementedError
? Does
simple_vector_line_integral
require a previous call tocohomology_basis
to work? (The docstring makes it seem so.) If this is the case, then the method should probably not be public. In other words, theNOTE
there is confusing.  The
if/else
inhomology_basis
…could you add a comment on why you can be sure thatcds[0] != cds[1]
? If this could go wrong in low precision, you might want to replace the assert with something else.  Should
RiemannSurface
really be anobject
and not at least aSageObject
? You could probably add a comment about why this is not a scheme in the hierarchy of the namespace there.
There are some confusing comments in the source code, please try to make them more understandable to people who did not write this code:
# Using singular to compute the genus quickly. Later, use rank/2 of # gram matrix?
# TODO: can we get rid of the comments below? # I'm not sure if this check will work, but could be useful? # I think singular's genus computation returns a negative genus # if it's reducible. Either way, if the genus is negative we should # terminate.
#TODO: Why not make a list in the first place with the edges?
#mt = Matrix(ZZ,[[1 if i==j else 0 for j in range(4*g**2)] + # [((S*w.monomial_coefficient(vars[i]).real_part()).round() for w in W] + # [S*w.monomial_coefficient(vars[i]).imag_part()).round() for w in W] for i in range(len(vars))])
comment:12 Changed 3 years ago by
 Branch changed from u/nbruin/riemann_surface to u/saraedum/riemann_surface
comment:13 Changed 3 years ago by
 Commit changed from 5b1c2d61c9f88c175b8edb15ad69ca2307163ae3 to 4b3eeea12060e9d2dc38ab4c1ee70ae8c2cacd45
I hope you do not mind my changes. And sorry for writing so many comments. You don't really have to fix all of them, it's just things I stumbled upon so I thought I could as well mention them.
New commits:
415e21d  Merge branch 'develop' into t/23175/riemann_surface

4b3eeea  Some reviewer style changes

comment:14 in reply to: ↑ 11 Changed 3 years ago by
Replying to saraedum:
I am making some style changes.
Thank you for style formatting!
 I guess you want to make this more accessible than
from sage.schemes.riemann_surfaces.riemann_surface import RiemannSurface
?
Probably affine and projective plane curves will eventually grow a method "riemann_surface", so I don't think accessibility is an issue. RiemannSurface
can mean more than just a compact one, so I'd be careful exposing it in the global namespace for that reason under that name (and other names get unappetizing pretty quickly)
 Why so much precision in the examples? Having 200 bits of precision seems like a lot. I tried a couple of the examples with less precision and found the output much easier to parse.
Turned it back a little
TODO: have an example where self.genus < 0 to raise error
Removed it. Actually, surfaces with multiple components shouldn't be an issue. The singular backend doesn't like them, though.
 Comment on the constant 100 in
_determine_new_w
and_newton_iteration
(nothing fancy necessary, just say that this is a random bound that worked well in practice?)
Done. This will never be triggered in practice.
 In
period_matrix
:Computes the period matrix of f given the homology basis and cohomology basis.
What is "given" here?
Fixed
 Why do you perform your own caching in
period_matrix
(self._PM
) andriemann_matrix
(self._RM
),
We did this to accommodate for the different options in the cohomology basis. However, in practice nobody would change that, and in actuality there's a benefit in providing the user a hook to specify this him/herself. This is implemented now, so the period matrix can now be cached in a vanilla way (and the riemann matrix is easily derived, so doesn't need to be cached)
cohomology_basis
(self._differentials
) ?
That's because the attribute can be specified in multiple ways.
 The
self._differentials[1]
is very confusing. Why do you need it?
Removed now, in favour of a userconfigurable "differentials" value.
raise ValueError('Length of differentials list is not equal to genus.')
should this not be aNotImplementedError
?
No, that's an error. "Genus" here should be H^{0(k,Omega}1_S). Thus the genus of the union of two elliptic curves should be 2.
 Does
simple_vector_line_integral
require a previous call tocohomology_basis
to work? (The docstring makes it seem so.) If this is the case, then the method should probably not be public. In other words, theNOTE
there is confusing.
No, it's a useful routine in general, but it does need stuff preinitialized (you'll get an error otherwise, no incorrect results). It's pretty innerloop otherwise, so it shouldn't be weighed down with an expensive call to homology_basis.
 The
if/else
inhomology_basis
…could you add a comment on why you can be sure thatcds[0] != cds[1]
? If this could go wrong in low precision, you might want to replace the assert with something else.
There are much worse ways in which that code can fail than the assert. I don't think it's worth worrying about.
 Should
RiemannSurface
really be anobject
and not at least aSageObject
?
I don't see a benefit from making it a SageObject
. I don't see how this object would benefit from interacting more closely with sage objects. It's just a container to spit out Riemann matrices.
You could probably add a comment about why this is not a scheme in the hierarchy of the namespace there.
I don't think I have a comment on that.
There are some confusing comments in the source code, please try to make them more understandable to people who did not write this code:
# Using singular to compute the genus quickly. Later, use rank/2 of # gram matrix?
Gone
# TODO: can we get rid of the comments below? # I'm not sure if this check will work, but could be useful? # I think singular's genus computation returns a negative genus # if it's reducible. Either way, if the genus is negative we should # terminate.
fixed
#TODO: Why not make a list in the first place with the edges?
fixed
#mt = Matrix(ZZ,[[1 if i==j else 0 for j in range(4*g**2)] + # [((S*w.monomial_coefficient(vars[i]).real_part()).round() for w in W] + # [S*w.monomial_coefficient(vars[i]).imag_part()).round() for w in W] for i in range(len(vars))])
removed
comment:15 Changed 3 years ago by
 Branch changed from u/saraedum/riemann_surface to u/nbruin/riemann_surface
comment:16 Changed 3 years ago by
 Commit changed from 4b3eeea12060e9d2dc38ab4c1ee70ae8c2cacd45 to 358f995077d983e1efedfe4de3950973910fea42
 Status changed from needs_work to needs_review
New commits:
358f995  Bugfixes and review suggestions.

comment:17 Changed 3 years ago by
 Commit changed from 358f995077d983e1efedfe4de3950973910fea42 to f8b255b2abb1d6091dd783db35232d5b5ec1083a
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
f8b255b  Bugfixes and review suggestions.

comment:18 Changed 3 years ago by
Ok. Feel free to set this to positive review when the patchbot is happy.
comment:19 Changed 3 years ago by
 Commit changed from f8b255b2abb1d6091dd783db35232d5b5ec1083a to 34ac9e1ee7523a9e926fa1582f4d1ace7238e939
Branch pushed to git repo; I updated commit sha1. New commits:
34ac9e1  Several fixes:

comment:20 Changed 3 years ago by
Sorry! I ran into one case where the Voronoi cells as computed did not give rise to a homology basis (because we were discarding edges in a way that led to an nonconnected diagram). I rewrote another routine to be better documented and clearer as well.
Knockon effect of the changed voronoi cells is that a lot of doctests with arbitrary output changed, but those are not very insightful. Nonetheless, probably good form if the reviewer gives a thumbsup to these changes before we set this to positive.
comment:21 Changed 3 years ago by
 Commit changed from 34ac9e1ee7523a9e926fa1582f4d1ace7238e939 to ce5097a4dcb6bf45a11078b3a91fc926c2b73592
Branch pushed to git repo; I updated commit sha1. New commits:
ce5097a  Two more fixes:

comment:22 Changed 3 years ago by
 Commit changed from ce5097a4dcb6bf45a11078b3a91fc926c2b73592 to 32856842d33c8b378d354331ddc12dfd7f5884bc
Branch pushed to git repo; I updated commit sha1. New commits:
3285684  Use sage_mantissa_exponent to get a cheaper log2

comment:23 Changed 3 years ago by
 Commit changed from 32856842d33c8b378d354331ddc12dfd7f5884bc to 4eaaff07e7c64df6f205d9b9293458d1d755b1a9
Branch pushed to git repo; I updated commit sha1. New commits:
4eaaff0  correction to implementation of voronoi_ghost

comment:24 Changed 3 years ago by
Feel free to set this back to positive review if you're confident that tests are going to pass.
comment:25 Changed 3 years ago by
OK. the bots seem busy. I've tested on 8.0beta11 and all tests pass. A previous complaint by a plugin about "EXAMPLE:" rather than "EXAMPLES:" has been fixed too, so following Julian's assessment: positive review. Preferably new issues get their own ticket.
comment:26 Changed 3 years ago by
 Status changed from needs_review to positive_review
comment:27 Changed 3 years ago by
 Status changed from positive_review to needs_work
Fails on 32bit:
sage t long src/sage/schemes/riemann_surfaces/riemann_surface.py ********************************************************************** File "src/sage/schemes/riemann_surfaces/riemann_surface.py", line 806, in sage.schemes.riemann_surfaces.riemann_surface.RiemannSurface.upstairs_edges Failed example: edgeset = S.upstairs_edges(); edgeset Expected: [[(0, 0), (1, 0)], [(0, 1), (1, 1)], [(0, 0), (5, 1)], [(0, 1), (5, 0)], [(1, 0), (4, 1)], [(1, 1), (4, 0)], [(2, 0), (3, 1)], [(2, 1), (3, 0)], [(2, 0), (4, 0)], [(2, 1), (4, 1)], [(3, 0), (5, 0)], [(3, 1), (5, 1)], [(4, 0), (5, 0)], [(4, 1), (5, 1)]] Got: [[(0, 0), (1, 1)], [(0, 1), (1, 0)], [(0, 0), (5, 0)], [(0, 1), (5, 1)], [(1, 0), (4, 0)], [(1, 1), (4, 1)], [(2, 0), (3, 0)], [(2, 1), (3, 1)], [(2, 0), (4, 1)], [(2, 1), (4, 0)], [(3, 0), (5, 1)], [(3, 1), (5, 0)], [(4, 0), (5, 0)], [(4, 1), (5, 1)]] ********************************************************************** File "src/sage/schemes/riemann_surfaces/riemann_surface.py", line 1423, in sage.schemes.riemann_surfaces.riemann_surface.RiemannSurface.period_matrix Failed example: S.period_matrix() # abs tol 0.000001 Expected: [0.77519780590366792 + 0.61819962132820207*I 0.19145804690473785  1.3890819401503317*I 2.1720559850746018 + 0.49575760460334677*I 0.23874427945779023  0.49575760460334677*I 0.43020232636252808  0.89332433554698493*I 0.96665585280840577  2.0072815614785338*I] [ 2.1720559850746018 + 0.49575760460334677*I 3.1387118378830076  0.27512471421878286*I 0.53645352644587770  1.1139572259315488*I 1.3968581791709339 + 1.1139572259315488*I 1.7418536587120737  1.3890819401503317*I 0.96665585280840577  0.77088231882212961*I] [ 0.53645352644587769  1.1139572259315488*I 1.5031093792542835  0.89332433554698493*I 0.77519780590366793 + 0.61819962132820208*I 2.7085095115204795  0.61819962132820208*I 1.2054001322661960  0.27512471421878285*I 0.96665585280840577 + 0.22063289038456391*I] Got: [ 0.43020232636252809  1.3345901163161127*I 0.96665585280840577  0.77088231882212962*I 2.1720559850746018 + 0.49575760460334669*I 1.1926223897340549e18  0.99151520920669346*I 0.96665585280840577  0.22063289038456391*I 2.3852447794681098e18 + 2.2279144518630976*I] [ 1.7418536587120737 + 2.6254811828067359*I 0.96665585280840577 + 0.22063289038456391*I 0.53645352644587770  1.1139572259315489*I 2.1684043449710089e19 + 2.2279144518630977*I 0.96665585280840577 + 2.0072815614785338*I 2.6020852139652106e18  1.2363992426564042*I] [ 1.2054001322661960 + 1.2666399234254763*I 0.96665585280840577  2.0072815614785338*I 0.77519780590366792 + 0.61819962132820217*I 3.2526065174565133e19  1.2363992426564042*I 0.96665585280840577 + 0.77088231882212962*I 4.3368086899420177e18  0.99151520920669343*I] Tolerance exceeded in 33 of 36: 0.77519780590366792 vs 0.43020232636252809, tolerance 1e+00 > 1e06 + 0.61819962132820207 vs  1.3345901163161127, tolerance 2e+00 > 1e06 0.19145804690473785 vs 0.96665585280840577, tolerance 1e+00 > 1e06  1.3890819401503317 vs  0.77088231882212962, tolerance 6e01 > 1e06 2.1720559850746018 vs 2.1720559850746018, tolerance 4e+00 > 1e06 0.23874427945779023 vs 1.1926223897340549e18, tolerance 2e01 > 1e06  0.49575760460334677 vs  0.99151520920669346, tolerance 5e01 > 1e06 0.43020232636252808 vs 0.96665585280840577, tolerance 1e+00 > 1e06  0.89332433554698493 vs  0.22063289038456391, tolerance 7e01 > 1e06 0.96665585280840577 vs 2.3852447794681098e18, tolerance 1e+00 > 1e06  2.0072815614785338 vs + 2.2279144518630976, tolerance 4e+00 > 1e06 2.1720559850746018 vs 1.7418536587120737, tolerance 4e01 > 1e06 + 0.49575760460334677 vs + 2.6254811828067359, tolerance 2e+00 > 1e06 3.1387118378830076 vs 0.96665585280840577, tolerance 4e+00 > 1e06  0.27512471421878286 vs + 0.22063289038456391, tolerance 5e01 > 1e06 0.53645352644587770 vs 0.53645352644587770, tolerance 1e+00 > 1e06 1.3968581791709339 vs 2.1684043449710089e19, tolerance 1e+00 > 1e06 + 1.1139572259315488 vs + 2.2279144518630977, tolerance 1e+00 > 1e06 1.7418536587120737 vs 0.96665585280840577, tolerance 3e+00 > 1e06  1.3890819401503317 vs + 2.0072815614785338, tolerance 3e+00 > 1e06 0.96665585280840577 vs 2.6020852139652106e18, tolerance 1e+00 > 1e06  0.77088231882212961 vs  1.2363992426564042, tolerance 5e01 > 1e06 0.53645352644587769 vs 1.2054001322661960, tolerance 2e+00 > 1e06  1.1139572259315488 vs + 1.2666399234254763, tolerance 2e+00 > 1e06 1.5031093792542835 vs 0.96665585280840577, tolerance 2e+00 > 1e06  0.89332433554698493 vs  2.0072815614785338, tolerance 1e+00 > 1e06 0.77519780590366793 vs 0.77519780590366792, tolerance 2e+00 > 1e06 2.7085095115204795 vs 3.2526065174565133e19, tolerance 3e+00 > 1e06  0.61819962132820208 vs  1.2363992426564042, tolerance 6e01 > 1e06 1.2054001322661960 vs 0.96665585280840577, tolerance 2e01 > 1e06  0.27512471421878285 vs + 0.77088231882212962, tolerance 1e+00 > 1e06 0.96665585280840577 vs 4.3368086899420177e18, tolerance 1e+00 > 1e06 + 0.22063289038456391 vs  0.99151520920669343, tolerance 1e+00 > 1e06 ********************************************************************** File "src/sage/schemes/riemann_surfaces/riemann_surface.py", line 1454, in sage.schemes.riemann_surfaces.riemann_surface.RiemannSurface.riemann_matrix Failed example: S.riemann_matrix() #abs tol 0.0000001 Expected: [ 0.50000000000000000 + 1.3228756555322953*I 0.25000000000000000 + 0.66143782776614765*I 6.7220534694101275e18  2.1684043449710089e18*I] [ 0.25000000000000000 + 0.66143782776614764*I 0.25000000000000000 + 0.66143782776614765*I 0.50000000000000000  1.7347234759768071e18*I] [6.5052130349130266e19 + 3.4694469519536142e18*I 0.50000000000000000 + 3.4694469519536142e18*I 0.25000000000000000 + 0.66143782776614765*I] Got: [ 0.37500000000000000 + 0.33071891388307383*I 0.50000000000000000  6.9388939039072284e18*I 0.50000000000000001  7.9146758591441824e18*I] [ 0.50000000000000001  4.3368086899420177e18*I 0.24999999999999999 + 0.66143782776614764*I 1.1275702593849246e17 + 1.3877787807814457e17*I] [ 0.50000000000000000 + 1.5178830414797062e17*I 2.7105054312137611e18  5.4210108624275222e18*I 0.24999999999999998 + 0.66143782776614762*I] Tolerance exceeded in 10 of 18: 0.50000000000000000 vs 0.37500000000000000, tolerance 1e01 > 1e07 + 1.3228756555322953 vs + 0.33071891388307383, tolerance 1e+00 > 1e07 0.25000000000000000 vs 0.50000000000000000, tolerance 8e01 > 1e07 + 0.66143782776614765 vs  6.9388939039072284e18, tolerance 7e01 > 1e07 6.7220534694101275e18 vs 0.50000000000000001, tolerance 5e01 > 1e07 0.25000000000000000 vs 0.50000000000000001, tolerance 8e01 > 1e07 + 0.66143782776614764 vs  4.3368086899420177e18, tolerance 7e01 > 1e07 0.50000000000000000 vs 1.1275702593849246e17, tolerance 5e01 > 1e07 6.5052130349130266e19 vs 0.50000000000000000, tolerance 5e01 > 1e07 0.50000000000000000 vs 2.7105054312137611e18, tolerance 5e01 > 1e07 ********************************************************************** 3 items had failures: 1 of 6 in sage.schemes.riemann_surfaces.riemann_surface.RiemannSurface.period_matrix 1 of 6 in sage.schemes.riemann_surfaces.riemann_surface.RiemannSurface.riemann_matrix 1 of 6 in sage.schemes.riemann_surfaces.riemann_surface.RiemannSurface.upstairs_edges [206 tests, 3 failures, 28.84 s]
comment:28 Changed 3 years ago by
 Commit changed from 4eaaff07e7c64df6f205d9b9293458d1d755b1a9 to 6e689c2aa87ec16fa6ce335e171aa4bce1885120
Branch pushed to git repo; I updated commit sha1. New commits:
6e689c2  make doctests more robust

comment:29 Changed 3 years ago by
 Status changed from needs_work to needs_review
Sigh ... doctesting doctests with nonuniquely represented objects is a pain. It would be nice if we'd have a way of getting wider architecture exposure without Volker having to intervene. Let's try this.
comment:30 Changed 3 years ago by
 Status changed from needs_review to positive_review
comment:31 Changed 3 years ago by
Looks good. (Yes, I guess a better CI/CD infrastructure would be helpful. As I am just transferring the company I work for to gitlab CI/CD, I discussed this quite a bit with roed during the past few days actually. But I guess it would be quite some effort to swap our patchbot/buildbot out for a more standardized solution…)
comment:32 Changed 3 years ago by
 Branch changed from u/nbruin/riemann_surface to 6e689c2aa87ec16fa6ce335e171aa4bce1885120
 Resolution set to fixed
 Status changed from positive_review to closed
We'll likely be rewriting history on this branch.
New commits:
Creation of a vectorized gauss_legendre integrator
Change to node for degree=1, based on bugfix in mpmath:
some PEP8, reviewer comments; rationalized "degree" parameter for "nodes"
Initial checkin of RiemannSurface