Opened 2 years ago

Last modified 4 months ago

#29169 needs_review enhancement

Maximal and critical angles for cones

Reported by: mjo Owned by:
Priority: major Milestone: sage-9.7
Component: geometry Keywords:
Cc: jipilab, novoselt, tscrim, chapoton Merged in:
Authors: Michael Orlitzky Reviewers: Matthias Koeppe
Report Upstream: N/A Work issues:
Branch: u/mjo/ticket/29169 (Commits, GitHub, GitLab) Commit: fb55e7fc7e04ebc904f1b4c16f7716533be3e36d
Dependencies: #30605, #30654, #32877 Stopgaps:

Status badges

Description

My paper with two algorithms for finding the maximal/critical angles between polyhedral convex cones was accepted:

http://www.optimization-online.org/DB_HTML/2019/01/7048.html

I've already written these algorithms in sage, so it would be nice to add them to the library. I just need to clean up some of the documentation.

Change History (87)

comment:1 Changed 2 years ago by mjo

Work in progress on u/mjo/ticket/29169. Doesn't build or anything yet.

comment:2 Changed 2 years ago by jipilab

  • Cc jipilab added

comment:3 Changed 2 years ago by mjo

  • Dependencies set to 26623

comment:4 Changed 2 years ago by mkoeppe

  • Milestone changed from sage-9.1 to sage-9.2

Batch modifying tickets that will likely not be ready for 9.1, based on a review of the ticket title, branch/review status, and last modification date.

comment:5 Changed 2 years ago by jipilab

  • Branch set to u/mjo/ticket/29169
  • Commit set to 99d77f01268d0d9e381956ba8f0a228a6ab4c19f

New commits:

c9437f1Trac #26623: add global entries for pending cone catalog citations.
792ac2bTrac #26623: add cones.<whatever> shortcuts for common cones.
ddf922cTrac #26623: add sage.geometry.cone_catalog to the documentation index.
f4a4d34Trac #26623: update sage.geometry.cone doctests to use the cone catalog.
9c3bf30Trac #29169: add the Or2020 reference.
6b73284Trac #29169: add the geometry/cone_critical_angles module.
08dd8aaTrac #29169: add the critical_angles and max_angle methods.
99d77f0UNDO ME: add the critical angles module to the doc index.

comment:6 Changed 2 years ago by git

  • Commit changed from 99d77f01268d0d9e381956ba8f0a228a6ab4c19f to fc255ba01aa5ce3d2f176aa3f217600c76b7a199

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

1edea71Trac #29169: add the Or2020 reference.
a2a025fTrac #29169: add the geometry/cone_critical_angles module.
4d15305Trac #29169: add the critical_angles and max_angle methods.
fc255baUNDO ME: add the critical angles module to the doc index.

comment:7 Changed 2 years ago by git

  • Commit changed from fc255ba01aa5ce3d2f176aa3f217600c76b7a199 to 016564a88267480c9d5f0a684845053a6cfc2576

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

468103bTrac #29169: add the Or2020 reference.
0481039Trac #29169: add the geometry/cone_critical_angles module.
b5ac943Trac #29169: add the critical_angles and max_angle methods.
a3c4535Trac #29169: fix a broken random_cone() test.
016564aUNDO ME: add the critical angles module to the doc index.

comment:8 Changed 2 years ago by mjo

  • Status changed from new to needs_review

I can keep tweaking the documentation forever, but I think this is in good condition to get some third-party feedback.

comment:9 Changed 2 years ago by mjo

(You can check the docs with sage -docbuild --include-tests-blocks --underscore)

comment:10 Changed 2 years ago by chapoton

  • Dependencies changed from 26623 to #26623
  • dependency field was missing #
  • Author field is empty

comment:11 Changed 2 years ago by mjo

  • Authors set to Michael Orlitzky
  • Dependencies #26623 deleted

comment:12 Changed 2 years ago by chapoton

patchbot complains about coverage and unused import

comment:13 Changed 2 years ago by git

  • Commit changed from 016564a88267480c9d5f0a684845053a6cfc2576 to 9e917d3730c1783f9f4507039420383284e394e6

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

468103bTrac #29169: add the Or2020 reference.
7d57406Trac #29169: add the geometry/cone_critical_angles module.
9c4af89Trac #29169: add the critical_angles and max_angle methods.
e4683f4Trac #29169: fix a broken random_cone() test.
9e917d3UNDO ME: add the critical angles module to the doc index.

comment:14 Changed 2 years ago by mjo

The unused import was legit (and is now fixed), but the coverage warning is for the critical_angles and max_angle functions that are wrapped by the methods of the same name in the cone class. All of the documentation (and unit tests) is on the cone methods, because that's what people will use.

(I can add some trivial tests if a happy coverage report is worth it, but so far I have avoided it as it's a complete waste of CPU time.)

comment:15 Changed 2 years ago by chapoton

I think that coverage is mandatory, sorry. Very simple tests would do, and you can tag them with #indirect doctest if they do not use the method directly.

More annoying, the latest patchbot reports a Timed Out, which is not acceptable.

comment:16 Changed 2 years ago by git

  • Commit changed from 9e917d3730c1783f9f4507039420383284e394e6 to 017cf1ad3a571dfe69c566b3b73057f2931644f6

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

4037b4eTrac #29169: add the geometry/cone_critical_angles module.
dfd587eTrac #29169: add the critical_angles and max_angle methods.
90a1914Trac #29169: fix a broken random_cone() test.
017cf1aUNDO ME: add the critical angles module to the doc index.

comment:17 Changed 2 years ago by mjo

I think that will make patchbot happy. I added some trivial tests where they were absent, and fixed the timeout issue (fingers crossed) by adding # long time in a bunch of places and reducing the max test size in others.

comment:18 Changed 2 years ago by gh-kliem

I'm willing to review this, but I'd much rather go through it one by one. Is it possible to implement the functions one by one?

I didn't check, maybe it all depends on each other.

comment:19 Changed 2 years ago by mjo

They're largely intertwined, but I could still commit some of the helper functions individually if that would make reviewing easier.

The reason for all of the separate functions (and the size of the module) is that critical_angles and max_angle do a lot of the same stuff. The max angle is a critical angle -- so critical_angles will find it -- but if all we're looking for is the biggest one then we can bail out as soon as we find it. That makes max_angle faster, but naturally the two share a lot of code, all of which I've tried to factor out into sensibly-named helper functions.

Once I defined the helper functions, they needed docstrings, tests, and examples.... meaning that a lot of relatively trivial things shared between critical_angles and max_angle take up a page or two of code. The inverse_spd function, for example, is two lines of implementation and a page of docs. If only one of critical_angles or max_angle were present, all of those functions would be inlined. Some of them make sense out of context, but others like check_gevp_feasibility are 100% tied to the implementation of critical_angles and max_angle, and are separate only to avoid copy/pasting their code.

comment:20 Changed 2 years ago by gh-kliem

Ok. I think I can also just look at the helper functions and work my way up.

comment:21 Changed 23 months ago by gh-kliem

  • Totally minor: If your in this situation, your cone is probably the trivial cone or the ambient space for real. It sounds like you claim that it isn't:
-            raise ValueError("this cone cannot be trivial or the ambient space")
+            raise ValueError("the cone may not be trivial or the ambient space")

or just delete the "this"

Other than this, the changes to src/sage/geometry/cone.py seem fine to me.

comment:22 Changed 23 months ago by mjo

I see what you're saying, but in informal English I think "may not" suffers from the same problem. Likewise for "must not". How about "the cone should not be..."? I have been unable to trick my brain into reading that incorrectly =)

comment:23 Changed 23 months ago by gh-kliem

the cone should not be is fine. I wouldn't want to worry about it too much.

comment:24 follow-up: Changed 23 months ago by gh-kliem

  • inverse_spd: Wouldn't that fit better into matrix/matrix2.pyx. This is where cholesky is defined.
  • _normalize_gevp_solution: Fine.
  • This error message is a full sentence.
    +        ValueError: The ambient dimension must be a positive integer (there
    +        are no nontrivial cones in dimension 0).
    
    I think just the part in the parenthesis would be a fine error message.

otherwise _random_admissible_cone is fine

  • +    # The standard cone containment tests don't work with normalized
    +    # vectors in QQbar.
    +    u_is_in_P = all( g.inner_product(u) >= 0 for g in P.dual() )
    +    v_is_in_Q = all( h.inner_product(v) >= 0 for h in Q.dual() )
    
    I would propose fixing this.
    sage: 3/2*vector([1,-1,0]) in P                                                                                                                                                     
    True
    sage: QQbar(2).sqrt()*vector([1,-1,0]) in P  
    ... # very long traceback
    ValueError: Cannot coerce non-integral Algebraic Real 1.414213562373095? to Integer
    
    This is awful and inconsistent. There is no way you can convert the point [3/2, -3/2, 0] to be a lattice point. According to the documentation it should have returned False. I would propose fixing this as a dependency of this ticket:
    • outsource dual_rays or similar from dual
    • use dual_rays to improve _contains to behave coercion-like with your code
    • then you can just compute u in P.

If I ask whether a rational or real algebraic point is contained in a cone, I want coercion.

The second thing: all( z.inner_product(h) >= 0 for h in Q ) could then be reduced to z in Q.dual() (which is a bit wasteful, but not more than the current state). This could be exposed somehow as well as this is the general question whether Q satisfies some inequality, but that is probably not a priority.

  • gevp_licis not fast, but if the time of running this is not relevant, it should be fine.

comment:25 Changed 23 months ago by gh-kliem

In examples you define this G_index_sets = list(gevp_licis(G)) without ever using it or printing it.

comment:26 in reply to: ↑ 24 ; follow-up: Changed 23 months ago by mjo

  • Cc novoselt added

Replying to gh-kliem:

This is awful and inconsistent... If I ask whether a rational or real algebraic point is contained in a cone, I want coercion.

Me too, but I anticipate that Andrey will object. The cone module is fairly strict; if you want to do anything "dirty," you have to leave the lattice and move to your own vector space. You've probably seen this in other tickets of mine where I have to throw away the lattice structure of (say) a dual cone to be able to treat its rays as objects in the primal space. It's often annoying, but it's consistent within the cone module. And it's technically correct, even from the convex analysis/geometry perspective, so long as you use the traditional definition of "dual cone" where the elements lie in a dual vector space. So I've learned to live with it.

In either case, you are right that it should return False as documented.

comment:27 Changed 23 months ago by git

  • Commit changed from 017cf1ad3a571dfe69c566b3b73057f2931644f6 to 21f1cd2b9ce9838adb7f3f8cda6258dd62e407c2

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

c9cb3d5Trac #29169: add the Or2020 reference.
583967dTrac #29169: add the geometry/cone_critical_angles module.
61a46dcTrac #29169: add the critical_angles and max_angle methods.
c323228Trac #29169: fix a broken random_cone() test.
ec72863UNDO ME: add the critical angles module to the doc index.
f32eb63Trac #29169: tweak two max/critical angle error messages.
c8dd547Trac #29169: improve error message for _random_admissible_cone().
21f1cd2Trac #29169: remove two unused example lines.

comment:28 in reply to: ↑ 26 ; follow-up: Changed 23 months ago by gh-kliem

How about def contains(self, coerce=False) or similar. So Andrey can have his way, but we are still not creating copies of stuff at places where they do not belong.

It's still strange that it would accept a vector but not a point in the dual lattice as input.

In a way it makes sense:

sage: 2 in range(2,5)                                                                                                                                                                                                                                                                                                                                                      
True
sage: 2.5 in range(2,5)                                                                                                                                                                                                                                                                                                                                                    
False

This is completely natural and I can see that you want this behavior for a cone. However, the other question should still be accessible. As you cannot change the base ring of a cone (can you?) it would be nice for the cone to provide some answer for you.

Replying to mjo:

Replying to gh-kliem:

This is awful and inconsistent... If I ask whether a rational or real algebraic point is contained in a cone, I want coercion.

Me too, but I anticipate that Andrey will object. The cone module is fairly strict; if you want to do anything "dirty," you have to leave the lattice and move to your own vector space. You've probably seen this in other tickets of mine where I have to throw away the lattice structure of (say) a dual cone to be able to treat its rays as objects in the primal space. It's often annoying, but it's consistent within the cone module. And it's technically correct, even from the convex analysis/geometry perspective, so long as you use the traditional definition of "dual cone" where the elements lie in a dual vector space. So I've learned to live with it.

In either case, you are right that it should return False as documented.

comment:29 in reply to: ↑ 28 ; follow-up: Changed 23 months ago by novoselt

Replying to gh-kliem:

It's still strange that it would accept a vector but not a point in the dual lattice as input.

For any vector space with a basis it is quite natural to have conversion from the set of numbers, i.e. a vector from ZZ^n to that vector space. And the backward map is fairly natural as well. However a composition of such maps does not give a natural map between two vector spaces with chosen bases. On that ground dual lattice elements are not accepted as points in the primal lattice. And the reason why this particular case is explicitly prohibited is that it is a fairly common mistake in toric geometry to work with elements/cones of the wrong lattice, so it helps to catch mistakes. Volker's alternative code was actually having completely unrelated classes for a lattice and its dual to make sure that there is no conversion/coercion happening between them.

So I am strongly against changing this behaviour in any situation. Having explicit non-default way to mix lattices is perfectly fine, of course.

comment:30 follow-up: Changed 23 months ago by novoselt

Regarding checking rational/algebraic/float vectors for cone containment - it would be nice for it to work. The problem is, I guess, how to write down the check. For my purposes working with rational numbers was always sufficient and it has the advantage of exact computations, so no need to worry about precision issues. If you can add code to handle more cases - please do it and it has nothing to do with distinction between dual lattices!

comment:31 in reply to: ↑ 29 ; follow-up: Changed 23 months ago by mkoeppe

Replying to novoselt:

[...] dual lattice elements are not accepted as points in the primal lattice. And the reason why this particular case is explicitly prohibited is that it is a fairly common mistake in toric geometry to work with elements/cones of the wrong lattice, so it helps to catch mistakes. Volker's alternative code was actually having completely unrelated classes for a lattice and its dual to make sure that there is no conversion/coercion happening between them.

As a side remark, a weakness of the vector space / free module code in sage.modules is that dual spaces / modules are not implemented. The newer code in sage.tensor.modules fixes this, and it would be a great improvement if sage.geometry.cone and sage.geometry.toric_lattice could integrate with that.

comment:32 in reply to: ↑ 30 Changed 23 months ago by mjo

Replying to novoselt:

If you can add code to handle more cases - please do it and it has nothing to do with distinction between dual lattices!

I was mainly worried about being able to take, say, M(1,1) and (through coercion) come to the conclusion that it lives in the nonnegative quadrant of the lattice N. But thanks to some implementation details, that appears not to be a problem. I'll post something soon.

comment:33 Changed 23 months ago by mjo

  • Dependencies set to #30605

comment:34 in reply to: ↑ 31 Changed 23 months ago by novoselt

Replying to mkoeppe:

As a side remark, a weakness of the vector space / free module code in sage.modules is that dual spaces / modules are not implemented. The newer code in sage.tensor.modules fixes this, and it would be a great improvement if sage.geometry.cone and sage.geometry.toric_lattice could integrate with that.

Thank you so much for the pointer! I surely wished for this when I was implementing toric stuff.

comment:35 Changed 23 months ago by git

  • Commit changed from 21f1cd2b9ce9838adb7f3f8cda6258dd62e407c2 to adeb746e1e5628bc9bb8d9d548dd70889a881069

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

f8b8b37Trac #30605: update polyhedral cone containment documentation.
db132e0Trac #29169: add the Or2020 reference.
3ff8035Trac #29169: add the geometry/cone_critical_angles module.
68f339bTrac #29169: add the critical_angles and max_angle methods.
358aec4Trac #29169: fix a broken random_cone() test.
fd5ac83UNDO ME: add the critical angles module to the doc index.
f4fc161Trac #29169: tweak two max/critical angle error messages.
64c0784Trac #29169: improve error message for _random_admissible_cone().
31b47afTrac #29169: remove two unused example lines.
adeb746Trac #29169: take advantage of improvements in trac #30605.

comment:36 Changed 23 months ago by mjo

  • Dependencies changed from #30605 to #30605, #30654

comment:37 Changed 23 months ago by mjo

  • Status changed from needs_review to needs_work

I'll need to rebase this again once #30654 is reviewed.

comment:38 Changed 23 months ago by git

  • Commit changed from adeb746e1e5628bc9bb8d9d548dd70889a881069 to be2d48ec8a2a112e80178e00cb8614531365029b

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

aa506fcTrac #30654: mention inverse_positive_definite() in inverse().
5a834d2Trac #29169: add the Or2020 reference.
4341a39Trac #29169: add the geometry/cone_critical_angles module.
cfce9c1Trac #29169: add the critical_angles and max_angle methods.
970db5dTrac #29169: fix a broken random_cone() test.
714e3b7UNDO ME: add the critical angles module to the doc index.
f72376eTrac #29169: tweak two max/critical angle error messages.
511c360Trac #29169: improve error message for _random_admissible_cone().
4b452d5Trac #29169: remove two unused example lines.
be2d48eTrac #29169: take advantage of improvements in trac #30605.

comment:39 Changed 23 months ago by mjo

  • Status changed from needs_work to needs_review

Cherry-picked the two commits from #30654, and replaced the inverse_spd function with the inverse_positive_definite method.

comment:40 Changed 22 months ago by mkoeppe

  • Milestone changed from sage-9.2 to sage-9.3

comment:41 Changed 22 months ago by chapoton

  • Status changed from needs_review to needs_work

red branch => needs work

comment:42 Changed 21 months ago by git

  • Commit changed from be2d48ec8a2a112e80178e00cb8614531365029b to 470ee05f27f23e2a8f21fb5b6376b3e15320e238

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

dd51796Trac #29169: add the Or2020 reference.
7b24315Trac #29169: add the geometry/cone_critical_angles module.
3e35062Trac #29169: add the critical_angles and max_angle methods.
fb46bcbTrac #29169: fix a broken random_cone() test.
b9c0e4bUNDO ME: add the critical angles module to the doc index.
efca675Trac #29169: tweak two max/critical angle error messages.
115e90aTrac #29169: improve error message for _random_admissible_cone().
5fc86f2Trac #29169: remove two unused example lines.
470ee05Trac #29169: take advantage of improvements in trac #30605.

comment:43 Changed 21 months ago by mjo

  • Status changed from needs_work to needs_review

Needed a rebase, I think.

comment:44 Changed 17 months ago by mkoeppe

  • Milestone changed from sage-9.3 to sage-9.4

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

comment:45 Changed 13 months ago by mkoeppe

  • Milestone changed from sage-9.4 to sage-9.5

Setting a new milestone for this ticket based on a cursory review.

comment:46 Changed 8 months ago by git

  • Commit changed from 470ee05f27f23e2a8f21fb5b6376b3e15320e238 to 8d982c9457430ecbbcd060ea4de7facbc0feaa35

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

716cb50Trac #29169: add the Or2020 reference.
4b9d3d5Trac #29169: add the geometry/cone_critical_angles module.
3972994Trac #29169: add the critical_angles and max_angle methods.
6c76caeTrac #29169: fix a broken random_cone() test.
63fa030UNDO ME: add the critical angles module to the doc index.
5ea8a65Trac #29169: tweak two max/critical angle error messages.
8416c4dTrac #29169: improve error message for _random_admissible_cone().
3996831Trac #29169: remove two unused example lines.
2f9235cTrac #29169: take advantage of improvements in trac #30605.
8d982c9Trac #29169: remove set_random_seed() calls.

comment:47 Changed 8 months ago by mjo

  • Dependencies changed from #30605, #30654 to #30605, #30654, #32877

comment:48 Changed 8 months ago by mkoeppe

  • Milestone changed from sage-9.5 to sage-9.6

comment:49 Changed 5 months ago by mkoeppe

Typo: greather -> greater

comment:50 follow-up: Changed 5 months ago by mkoeppe

Wondering if there shouldn't be a variant of these methods that returns the tan of the angle rather the angle. Could one then expect that the output would be rational for rational input?

comment:51 Changed 5 months ago by mkoeppe

Another typo: constuct -> construct.

Run tox -e codespell?

comment:52 in reply to: ↑ 50 ; follow-up: Changed 5 months ago by mkoeppe

Replying to mkoeppe:

Wondering if there shouldn't be a variant of these methods that returns the tan of the angle rather the angle. Could one then expect that the output would be rational for rational input?

Or, in the methods that return a witness pair of angles, just do not output the angle, as it is trivial to recover?

This would eliminate the dependency on sage.symbolic (for #32432)

comment:53 Changed 5 months ago by mkoeppe

Also there seems to be no doctest for exact=False

comment:54 follow-up: Changed 5 months ago by mkoeppe

I'd suggest to get rid of _lists_equivalent and just write set(l1) == set(l2)

comment:55 Changed 5 months ago by mkoeppe

Perhaps _solve_gevp_naive and solve_gevp_nonzero should be generators?

comment:56 Changed 5 months ago by git

  • Commit changed from 8d982c9457430ecbbcd060ea4de7facbc0feaa35 to 70cd301036ce8c958aae6bfe89b65d7e91c40024

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

545900fTrac #29169: add the critical_angles and max_angle methods.
50a0aadTrac #29169: fix a broken random_cone() test.
2d69edcUNDO ME: add the critical angles module to the doc index.
b61a34eTrac #29169: tweak two max/critical angle error messages.
1e6dd4dTrac #29169: improve error message for _random_admissible_cone().
4a5ec4dTrac #29169: remove two unused example lines.
a462418Trac #29169: take advantage of improvements in trac #30605.
1e17f39Trac #29169: remove set_random_seed() calls.
208a528Trac #29169: fix a few typos.
70cd301Trac #29169: convert solve_gevp* functions to generators.

comment:57 in reply to: ↑ 54 ; follow-up: Changed 5 months ago by mjo

Replying to mkoeppe:

I'd suggest to get rid of _lists_equivalent and just write set(l1) == set(l2)

Doesn't work when the lists contain (mutable, not-hashable) vectors unfortunately. I also tried sage's Set() in case I had overlooked that originally but it's not checking equality element-wise.

I fixed the typos and converted those methods to generators. I'll add tests for exact=False when I get a chance.

comment:58 in reply to: ↑ 57 ; follow-up: Changed 5 months ago by mkoeppe

Replying to mjo:

Doesn't work when the lists contain (mutable, not-hashable) vectors unfortunately.

Perhaps just change the code so it creates immutable vectors?

comment:59 in reply to: ↑ 58 ; follow-up: Changed 5 months ago by mjo

Replying to mkoeppe:

Perhaps just change the code so it creates immutable vectors?

That still doesn't work because

sage: set([AA(0)]) == set([ZZ(0)])
False

The _lists_equivalent function is used to compare eigenvalue/vector results coming from three different places, two of which use extend=True and can "silently" jump from e.g. the rationals to the algebraic field. I could probably get them all to return eigenvalues in the same field, but it really doesn't matter for the tests. It's a lot easier to compare them with == than it would be to make the tests or the functions themselves more complicated.

comment:60 Changed 5 months ago by git

  • Commit changed from 70cd301036ce8c958aae6bfe89b65d7e91c40024 to 6db717c066b84a482be12bfc7873a38e556f684a

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

e6223c8Trac #29169: convert solve_gevp* functions to generators.
823d5cbTrac #29169: add a warning about inexact critical angle arithmetic.
6db717cTrac #29169: add an exact=False example for cone max_angle().

comment:61 in reply to: ↑ 52 Changed 5 months ago by mjo

Replying to mkoeppe:

Replying to mkoeppe:

Wondering if there shouldn't be a variant of these methods that returns the tan of the angle rather the angle. Could one then expect that the output would be rational for rational input?

Or, in the methods that return a witness pair of angles, just do not output the angle, as it is trivial to recover?

This would eliminate the dependency on sage.symbolic (for #32432)

This makes some engineering sense, but I'd prefer to keep the nice interface (e.g max_angle versus min_cos_of_angle). The internal functions could easily be reworked, but in the user-facing portion, I think we should give people what they're actually looking for and not make them compute it themselves (which would require sage.symbolic anyway).

comment:62 follow-up: Changed 5 months ago by mkoeppe

+            raise ValueError('eigenspace of dimension %d > 1 '
+                             'corresponding to eigenvalue %s'
+                              % (mult, cos_theta))

Shouldn't this be a NotImplementedError?

comment:63 Changed 5 months ago by mkoeppe

gevp_licis can be replaced by sage.matroids.linear_matroid.LinearMatroid(G).independent_sets()

comment:64 Changed 5 months ago by mkoeppe

+    if exact:
+        ring = QQbar

Do you really need QQbar or is AA enough?

comment:65 Changed 5 months ago by git

  • Commit changed from 6db717c066b84a482be12bfc7873a38e556f684a to 9f28b6469ea7482f572b1e70dec9966d424d594b

Branch pushed to git repo; I updated commit sha1. New commits:

790423eTrac #29169: clean up two import lists in cone_critical_angles.py.
b5026c3Trac #29169: use matroid implementation for gevp_licis().
9f28b64Trac #29169: switch from QQbar to AA for cone critical angles.

comment:66 in reply to: ↑ 62 Changed 5 months ago by mjo

Replying to mkoeppe:

Shouldn't this be a NotImplementedError?

IIRC I checked the python docs for NotImplementedError and thought some other error would be more accurate. It's not that we've chosen not to implement this; rather, we encountered a value that we don't know how to handle while performing some other computation. If you prefer NotImplementedError though I don't care all that much.

gevp_licis can be replaced by sage.matroids.linear_matroid.LinearMatroid(G).independent_sets()

Not entirely, but the latter gives a much faster implementation, thanks! I also added some docs stating that our subsets need to be ordered similarly. It's in the paper but wasn't in the code.

Do you really need QQbar or is AA enough?

I don't remember why I wrote it that way but I can't find a reason for QQbar now, so I've changed it.

comment:67 follow-up: Changed 5 months ago by mkoeppe

+    for I in G_index_sets:
+        G_I = G[range(n),I]
+        I_complement = [ i for i in range(P.nrays()) if not i in I ]
+        G_I_c_T = G[range(n),I_complement].transpose()
+
+        for J in H_index_sets:
+            J_complement = [ j for j in range(Q.nrays()) if not j in J ]

I think this code would be slightly cleaner if I and J were just frozenset. Then you can write G_I = G[range(n),sorted(I)] and I_complement = range(P.nrays()) - I

comment:68 in reply to: ↑ 67 ; follow-up: Changed 5 months ago by mjo

Replying to mkoeppe:

I think this code would be slightly cleaner if I and J were just frozenset. Then you can write ... I_complement = range(P.nrays()) - I

Should that work? I get unsupported operand type(s) for -: 'range' and 'frozenset'.

If you meant frozenset(range(P.nrays())), then that would work; still, we'd have to sort the whole thing again to make sure that I_complement has the same order as I itself (and that we can use it as a matrix index).

comment:69 in reply to: ↑ 68 Changed 5 months ago by mkoeppe

Replying to mjo:

If you meant frozenset(range(P.nrays()))

I guess that's what I meant

comment:70 in reply to: ↑ 59 Changed 5 months ago by mkoeppe

Replying to mjo:

sage: set([AA(0)]) == set([ZZ(0)])
False

The _lists_equivalent function is used to compare eigenvalue/vector results coming from three different places, two of which use extend=True and can "silently" jump from e.g. the rationals to the algebraic field.

The doctests could be changed to just compare with the results in the correct field, so that you don't run into the hashing contract violation hash(AA(0)) == hash(ZZ(0)) => False.

Last edited 5 months ago by mkoeppe (previous) (diff)

comment:71 Changed 5 months ago by git

  • Commit changed from 9f28b6469ea7482f572b1e70dec9966d424d594b to 4d9af99fece20a7feea9cfc7b2a1f8956900680c

Branch pushed to git repo; I updated commit sha1. New commits:

7979f25Trac #29169: use base ring's zero for the zero-eigenvalue subproblem.
c5303a9Trac #29169: return immutable vectors from eigenvalue problems.
4d9af99Trac #29169: replace _lists_equivalent() with set equality.

comment:72 Changed 5 months ago by mjo

That should do it. Even if the input vectors are rational we convert them to AA, and then the resulting matrix is over AA, so there's actually not too many options for the eigenvalue field (the other is RDF).

comment:73 Changed 5 months ago by mkoeppe

Thanks for making this change!

comment:74 follow-up: Changed 5 months ago by mkoeppe

+Finding the maximal (or equivalently, the minimal) angle between two

How is this equivalent?

comment:75 in reply to: ↑ 74 Changed 5 months ago by mjo

Replying to mkoeppe:

How is this equivalent?

The minimal angle between P and Q is pi minus the maximal angle between P and -Q, akin to how minimizing some real f is the same as maximizing -f and flipping the result.

comment:76 Changed 5 months ago by mkoeppe

Thanks.

comment:77 follow-up: Changed 5 months ago by mkoeppe

Does the change to src/sage/misc/randstate.pyx belong here on this ticket?

comment:78 Changed 5 months ago by mkoeppe

  • Cc tscrim chapoton added

Looks good to me now, but probably the PEP 8 police may want to have a word

comment:79 Changed 5 months ago by mkoeppe

  • Reviewers set to Matthias Koeppe

comment:80 Changed 5 months ago by git

  • Commit changed from 4d9af99fece20a7feea9cfc7b2a1f8956900680c to 9a8fb8c979dc27cd47a20d58f45100f4c2485a2f

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

60d5dcbTrac #29169: fix a few typos.
f87ea0dTrac #29169: convert solve_gevp* functions to generators.
d01643dTrac #29169: add a warning about inexact critical angle arithmetic.
33e073bTrac #29169: add an exact=False example for cone max_angle().
2b6e4ebTrac #29169: clean up two import lists in cone_critical_angles.py.
91017baTrac #29169: use matroid implementation for gevp_licis().
3a6304dTrac #29169: switch from QQbar to AA for cone critical angles.
1ed1b43Trac #29169: use base ring's zero for the zero-eigenvalue subproblem.
4491553Trac #29169: return immutable vectors from eigenvalue problems.
9a8fb8cTrac #29169: replace _lists_equivalent() with set equality.

comment:81 in reply to: ↑ 77 ; follow-up: Changed 5 months ago by mjo

Replying to mkoeppe:

Does the change to src/sage/misc/randstate.pyx belong here on this ticket?

Nope, dropped. There's also a commit titled UNDO ME that should be dropped before merging, but for now it allows the docs to build.

comment:82 Changed 5 months ago by tscrim

You don't need the \ here:

            if (other.lattice_dim() != self.lattice_dim()):
                raise ValueError("lattice dimensions of self and other "\
                                 "must agree")
            if other.is_trivial() or other.is_full_space():
                raise ValueError("other cone cannot be trivial or the "\
                                 "ambient space")

(The parentheses take care of it.)

No period/full-stop at the end of INPUT: lines.

In _normalize_gevp_solution, do you need to pass the full tuple? I guess this makes the code easier. If so, I would replace irrelevant with ignored.

No space around the list; e.g.:

-G_index_sets = [ s for s in gevp_licis(G) if not len(s) == n ]
+G_index_sets = [s for s in gevp_licis(G) if not len(s) == n]

Some of your # long time are misaligned compared to others in their block (it isn't a big issue, it was just visually jarring for me).

comment:83 Changed 5 months ago by git

  • Commit changed from 9a8fb8c979dc27cd47a20d58f45100f4c2485a2f to 911140a7eac21b73ff75ba2c576164d62eed3edb

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

f0679b3Trac #29169: add a warning about inexact critical angle arithmetic.
a68f088Trac #29169: add an exact=False example for cone max_angle().
a89ab18Trac #29169: clean up two import lists in cone_critical_angles.py.
9e47e6bTrac #29169: use matroid implementation for gevp_licis().
571d626Trac #29169: switch from QQbar to AA for cone critical angles.
ca3a18bTrac #29169: use base ring's zero for the zero-eigenvalue subproblem.
bef682dTrac #29169: return immutable vectors from eigenvalue problems.
ba237dbTrac #29169: replace _lists_equivalent() with set equality.
0926a6fTrac #29169: kill periods at the end of INPUT/OUTPUT list items.
911140aTrac #29169: kill spaces around [ list comprehensions ].

comment:84 Changed 5 months ago by git

  • Commit changed from 911140a7eac21b73ff75ba2c576164d62eed3edb to fb55e7fc7e04ebc904f1b4c16f7716533be3e36d

Branch pushed to git repo; I updated commit sha1. New commits:

fb55e7fTrac #29169: change "irrelevant" to "ignored" in a docstring.

comment:85 in reply to: ↑ 81 ; follow-up: Changed 5 months ago by mkoeppe

Replying to mjo:

Nope, dropped. There's also a commit titled UNDO ME that should be dropped before merging, but for now it allows the docs to build.

Yes, what's up with that?

comment:86 in reply to: ↑ 85 Changed 5 months ago by mjo

Replying to mkoeppe:

Replying to mjo:

Nope, dropped. There's also a commit titled UNDO ME that should be dropped before merging, but for now it allows the docs to build.

Yes, what's up with that?

Without it, the docs aren't built for cone_critical_angles.py, so there's no way to tell if the syntax/formatting is totally broken. I can remove it if everyone's happy with the rest.

comment:87 Changed 4 months ago by mkoeppe

  • Milestone changed from sage-9.6 to sage-9.7
Note: See TracTickets for help on using tickets.