Constructions for common polyhedral cones

[orlitzky]

When working with graphs, we have a nice tab-completed list of common graphs that we can construct; for example,

sage: graphs.BuckyBall()
Bucky Ball: Graph on 60 vertices

More and more, I find myself wishing we had the same thing for convex cones, particularly in test cases. I have constructions for all of the following and use them regularly:

• The nonnegative orthant.
• The full ambient space.
• The schur cone that induces the majorization ordering.
• The permutation-invariant "rearrangement" cone.
• The trivial cone, since trivial_cone(n) looks a lot nicer than Cone([], ToricLattice(n)).

et cetera; I'm sure I can think of more. Each of these makes sense in any dimension n.

I'm wondering if you think this would be a useful feature to make available as e.g.

sage: cones.nonnegative_orthant(5)
5-d cone in 5-d lattice N

They should all probably take a "lattice" argument too now that I think about it.

CC: @novoselt

Component: geometry

Author: Michael Orlitzky

Branch/Commit: b2aab91

Reviewer: Jonathan Kliem

Issue created by migration from https://trac.sagemath.org/ticket/26623

[novoselt]

comment:1

They certainly should take a lattice argument! From a toric perspective it is probably more natural to have these as methods of lattices, i.e. something like

ToricLattice(n).trivial_cone()

although of course this makes it impossible to use ZZ^n. And perhaps even more natural is to have not cones, but associated toric varieties which are already available:

sage: toric_varieties.A(5).fan().generating_cone(0).rays()
N(1, 0, 0, 0, 0),
N(0, 1, 0, 0, 0),
N(0, 0, 1, 0, 0),
N(0, 0, 0, 1, 0),
N(0, 0, 0, 0, 1)
in 5-d lattice N

[orlitzky]

comment:2

Replying to @novoselt:

They certainly should take a lattice argument! From a toric perspective it is probably more natural to have these as methods of lattices, i.e. something like

ToricLattice(n).trivial_cone()

although of course this makes it impossible to use ZZ^n.

I think from a design standpoint, that's where the method belongs, but there are two things that bug me about it.

First: in many cases, the user won't care about the lattice and will want "the default." This is basically what happens when you do ToricLattice(n).nonnegative_orthant(), but it feels a little weird to specify a particular lattice when you don't care about it. By analogy, the Cone() function should also be a method on a ToricLattice object. A lot of the time you just don't care though, so it's nice to be able to get a sensible default without having to construct a generic lattice of the correct size first. On the other hand, I like that this eliminates the n and lattice parameters from the cone-creating functions. And it would ensure that (for example) if you only create one lattice object, then you can intersect any two cones created by methods on it. I'll have to think harder about this but it may be the way to go.

Second: tab completion on graphs.<tab> is so nice. Having a list of twenty or so cones interspersed in the tab-completion list for lattices would be bad for both people who want to see lattice methods and people who want to see predefined cones.

(Might we have something like lattice.cones.<tab> where lattice is some ToricLattice?)

And perhaps even more natural is to have not cones, but associated toric varieties which are already available:

sage: toric_varieties.A(5).fan().generating_cone(0).rays()
...

My main objection to this is that I'm not smart enough to implement it or remember which cone is associated with which toric variety =)

Having the code above be the implementation for nonnegative_orthant(5) would be fine with me though.

[novoselt]

comment:3

"A lot of the time you just don't care though" - I always care, we are just using cones for different things ;-) But you are right - just because cones are associated to lattices does not mean that they have to be methods of lattices. So I think cones.etc is the right interface, but there has to be an option to specify the lattice and the default should be the same as for Cone(...).

[orlitzky]

comment:4

I have another question. There are a few vector spaces associated with each cone, that can be thought of as cones themselves:

• The span of the cone K.
• Its orthogonal complement, (span(K))-perp.
• Its lineality space (the largest subspace it contains).

The span of the cone and its lineality space are already available as,

sage: K = Cone([(1,1),(0,1)])
sage: K.span()
Sublattice <N(1, 0), N(0, 1)>
sage: span(K)
Free module of degree 2 and rank 2 over Integer Ring
Echelon basis matrix:
[1 0]
[0 1]
sage: K.linear_subspace()
Vector space of degree 2 and dimension 0 over Rational Field
Basis matrix:
[]

But none of those return cones directly. Having a cone is handy when you want to compute e.g. the intersection of one of these things with some other cone.

If we were to add similar methods that return cones, how would you prefer to do it? We could either...

1. Add methods to the cones themselves, for things like K.span_cone(), K.lineality_space(), K.orthogonal_complement()
2. Add them as part of the new cones.<whatever> interface, with something like cones.span_of(K).

The first one seems more correct to me, but risks polluting the namespace with a bunch of methods that all look the same but return different things.

The second one does make it obvious that you're going to get a cone back, but only if people know to look there for what should in principle be a method on the cone itself.

[orlitzky]

comment:5

Option 3:

For the span and the orthogonal complement, we could sensibly add a cone() method somewhere in the module hierarchy, like

sage: K.orthogonal_sublattice().cone()

which would be a wrapper around

sage: Cone( c*g for g in K.orthogonal_sublattice().gens() for c in [-1, 1] )

The problem with that is, when starting from a cone and passing through a vector space, the lattice information is lost. If I start with a cone in a lattice M, take its linear_subspace(), and then turn it into cone...

sage: K = Cone([(1,0,0),(-1,0,0),(0,1,0)], ToricLattice(3,'M'))
sage: Cone(K.linear_subspace().gens())
1-d cone in 3-d lattice N

it should live in the same lattice as the original cone. I don't see how to fix that immediately.

[novoselt]

comment:6

Just a reminder - orthogonality in the toric setting is handled via dual lattices:

sage: c = Cone([(1,2,3), (4,5,6)])
sage: c
2-d cone in 3-d lattice N
sage: c.orthogonal_sublattice()
Sublattice <M(1, -2, 1)>
sage: c.span()
Sublattice <N(1, 2, 3), N(0, 3, 6)>

I can't imagine use cases where I would like to treat all these spaces as cones. Can you please explain it in more details, i.e. what is wrong with treating them as, well, spaces? If you really do want to have cones associated to spaces, then you can easily attach a method to toric lattices and sublattices, or, alternatively or in addition, you can modify the Cone constructor to handle lattices and spaces. I just tried Cone(c.span()) for the above example and probably hit some infinite recursion that is in addition very slow - after a few minutes it is still doing something.

But I am still unclear why would one want to turn a space into a cone, what is there to gain from such a representation?

[orlitzky]

comment:7

Oof, it looks like I'm missing trac emails again.

It would just be nice to have the cone methods available on subspaces. For example, Corollary 1 in these notes,

http://michael.orlitzky.com/notes/on_z-operators_and_viability_theorems.pdf

is very easy to test as the intersection of a closed convex (dual) cone and a subspace. If span(x) is a cone, and if its orthogonal complement is a cone, then I can use the cone intersection method that we already have to test it.

Another example is the "maximal angle between two convex cones." This is a newish concept (from 2016), and I've implemented it as a method on cones. However, the "principal angle between subspaces" is a classical idea going back to the 1950s, and they turn out to be equivalent for cones that are subspaces. If subspaces can be treated as cones, then I can use the same method that I have on cones to compute principal angles between subspaces. And so on.

The reason Cone(c.span()) "hangs" is because it carefully enumerates every element of the vector space.

I think there's plenty of work to do just to get the trivial cone, nonnegative orthant, and a few others implemented and documented. I'll focus on that and worry about the other questions later.

[orlitzky]

Commit: eb072ce

[orlitzky]

comment:8

Here's my first attempt at this. I'm open to other suggestions if you don't like the module or class names.

---

New commits:

6645911	Trac #26623: add global entries for pending sage.geometry.cones citations.
84f5635	Trac #26623: add cones. shortcuts for common cones.
dd759dc	Trac #26623: add the sage.geometry.cones module to the documentation index.
eb072ce	Trac #26623: update doctests in cone.py to use new "cones" methods.

[orlitzky]

Author: Michael Orlitzky

[orlitzky]

Branch: u/mjo/ticket/26623

[videlec]

comment:9

I am fine with the class names but the objects would better be like graphs (and many other constructors): that is in caml case like Trivial, NonnegativeOrthant, etc

[orlitzky]

comment:10

Replying to @videlec:

I am fine with the class names but the objects would better be like graphs (and many other constructors): that is in caml case like Trivial, NonnegativeOrthant, etc

Doesn't that violate PEP8? They're not really constructors, even though they do mimic constructors (in the sense that they would be constructors if I didn't want the cones. prefix). However this is python where almost every function creates a new object and returns it, so the distinction seems a little arbitrary to me.

[videlec]

comment:11

Please have a look at: graphs., designs., polytopes., manifolds., ... All of them use caml case and most of them are not classes. I am not talking about PEP8 but consistency within sage.

[novoselt]

comment:12

It seems that when you do specify a lattice, you also have to specify its dimension and there is a check that you do not make a mistake. It would be more pleasant to provide only one input argument: either a number to get the cone in the default lattice of that dimension, or a lattice and get a cone in that lattice.

[orlitzky]

comment:13

Replying to @novoselt:

It seems that when you do specify a lattice, you also have to specify its dimension and there is a check that you do not make a mistake. It would be more pleasant to provide only one input argument: either a number to get the cone in the default lattice of that dimension, or a lattice and get a cone in that lattice.

But that's exactly how Cone() works, and Cone() is the function that underlies all of these convenience methods.

The additional check is not strictly required; it serves only to provide a more relevant error message. For example, if I were to delete the check and then pass a lattice of the wrong rank to Cone(), I would get

sage: M = ToricLattice(2)
sage: Cone([(1,-1,0),(0,1,-1)], lattice=M)
...
TypeError: cannot convert (1, -1, 0) to Vector space of dimension 2 over Rational Field!

If you know the generators of (say) the Schur cone, then you might be able to figure out what that error message means. But, if you don't know what those generators are off the top of your head, then an error message complaining about converting (1,-1,0) to a rational vector space could be confusing. That said, I'm not too attached to the additional sanity check and error message, and would be happy to remove them.

It would be more pleasant to provide only one input argument: either a number to get the cone in the default lattice of that dimension

I don't think this is workable, since in applications we often want to intersect the dual of some cone (which will live in a non-default lattice) with e.g. a nonnegative orthant.

only one input argument... a lattice, and get a cone in that lattice.

This is the best API, but as a user-interface is pretty annoying to use. The "pass a lattice, but only if you don't want the default one" approach of Cone() is IMO a good compromise. Having to do e.g.

sage: cones.trivial(ToricLattice(3))

is even worse than

sage: Cone([], ToricLattice(3))

In the common case, the user is left wondering why he has to type ToricLattice(n) all the time instead of just n. The default is usually fine until you start doing complicated stuff.

[novoselt]

comment:14

Cone works for arbitrary cones and it may be convenient to give input as integer vectors, especially as a user input, in which case it is handy to specify some lattice and cast all vectors to it. Just specifying a lattice does not make much sense for Cone except for a trivial one, perhaps. But in your constructors everything is determined by dimension, there are no explicit rays to give, so what I propose is to have two options: cones.trivial(3) and cones.trivial(my_lattice). I do not suggest that you eliminate the form when only the integer is given, rather do not require this integer in those cases, when the lattice was given explicitly!

[orlitzky]

comment:15

Replying to @novoselt:

I do not suggest that you eliminate the form when only the integer is given, rather do not require this integer in those cases, when the lattice was given explicitly!

Ah, sorry I misunderstood. This sounds like a good idea.

[sagetrac-git]

Changed commit from eb072ce to 5533bd2

[sagetrac-git]

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

53f7b42	Trac #26623: add global entries for pending sage.geometry.cones citations.
cab6be8	Trac #26623: add cones. shortcuts for common cones.
f48f779	Trac #26623: add the sage.geometry.cones module to the documentation index.
eda8548	Trac #26623: update doctests in cone.py to use new "cones" methods.
8122dc0	Trac #26623: factor out common "cones" argument processing.
5533bd2	Trac #26623: use Pascal case for "cones" method names.

[orlitzky]

comment:17

Each method will now make an attempt to determine the ambient dimension from the lattice and vice-versa. An error is raised if they are both left unspecified, or if they disagree. I don't think it's necessary to prohibit e.g. nonnegative_orthant(3,M) if the rank of M is 3.

I also added a commit on top of the others to switch to PascalCase names. My mother taught me that "everyone else is doing it" isn't a good excuse, so I'm not sold on the motivation, but it's there. I don't feel that strongly about it if I'm in the minority.

[sagetrac-git]

Changed commit from 5533bd2 to 59fd912

[sagetrac-git]

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

7ddd803	Trac #26623: add global entries for pending sage.geometry.cones citations.
0588289	Trac #26623: add cones. shortcuts for common cones.
a80699c	Trac #26623: add the sage.geometry.cones module to the documentation index.
1f8aeaf	Trac #26623: update doctests in cone.py to use new "cones" methods.
d266a33	Trac #26623: factor out common "cones" argument processing.
59fd912	Trac #26623: use Pascal case for "cones" method names.

[orlitzky]

comment:19

Force-pushed a rebase onto the "develop" branch to fix a conflict in src/doc/en/reference/references/index.rst.

[jplab]

comment:20

There seems to be a conflict with version 8.9.beta8.

[kliem]

comment:35

Thanks. Overall this will be a nice improvement.

A few comments for now.

• I would suggest using commas here:

+globally-available ``cones`` prefix, to create some common cones:
+
+- The nonnegative orthant,
+- the rearrangement cone of order ``p``,
+- the Schur cone,
+- the trivial cone.
-   globally-available ``cones`` prefix, to create some common cones:
-
-- The nonnegative orthant
-- The rearrangement cone of order ``p``
-- The Schur cone
-- The trivial cone

• The warning at the beginning should include mentioning, that any import will be imported at startup.
• If you lazy import Cone and random_cone than the bots will be fine with the number of startup modules. Also this should be done anyway, shouldn't it?

[sagetrac-git]

Changed commit from 42e9e2a to 282a694

[sagetrac-git]

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

fdfd1e4	Trac #26623: add cones. shortcuts for common cones.
ffa1e9b	Trac #26623: add sage.geometry.cone_catalog to the documentation index.
282a694	Trac #26623: update sage.geometry.cone doctests to use the cone catalog.

[orlitzky]

comment:37

Thanks, I added the commas, and expanded the warning at the beginning of the module.

As for the bots: I don't think they're complaining about Cone or random_cone:

1. Both of those are loaded by sage.geometry.all anyway.

2. Those import statements aren't executed until you run one of the functions in the cone_catalog module, so importing them is essentially "lazy" even if you're using sage as a library and not interactively.

Instead, I think the bots are pointing out that sage.geometry.cone_catalog is now loaded into the global namespace as cones, but that's almost unavoidable since that's the whole point. However, it is indeed slower than the lazy import of the class was, should you decide not to use it.

There may be some clever way to delay even that import, while still retaining the name "cones" along with its tab-completion list. For example, if I give a top-level variable name to the module, and then import that object lazily... e.g. this seems to work in cone_catalog.py

import sys

# A python object that can be imported lazily with lazy_import(...).
cones = sys.modules[__name__]

def __dir__():
    # Don't expose any global imports or variables to tab-completion.
    return ['nonnegative_orthant',
            'rearrangement',
            'schur',
            'trivial']

if I then do

lazy_import('sage.geometry.cone_catalog', 'cones')

in sage/geometry/all.py. But that's veering dangerously close to "too clever" to save a few microseconds in my opinion.

---

New commits:

fdfd1e4	Trac #26623: add cones. shortcuts for common cones.
ffa1e9b	Trac #26623: add sage.geometry.cone_catalog to the documentation index.
282a694	Trac #26623: update sage.geometry.cone doctests to use the cone catalog.

[kliem]

comment:38

One could just do

lazy_import("sage.geometry", "cone_catalog", "cones")

I don't know if it's faster, but I don't see a disadvantage. Tab completion and everything still works.

[sagetrac-git]

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

7567874	Trac #26623: add global entries for pending cone catalog citations.
39fd483	Trac #26623: add cones. shortcuts for common cones.
c6a1b3a	Trac #26623: add sage.geometry.cone_catalog to the documentation index.
e8fd707	Trac #26623: update sage.geometry.cone doctests to use the cone catalog.

[sagetrac-git]

Changed commit from 282a694 to e8fd707

[orlitzky]

comment:40

Replying to @kliem:

I don't know if it's faster, but I don't see a disadvantage. Tab completion and everything still works.

Thank you! I wasted a lot of time trying to figure out how to do that. I had even written a new lazy_import_module function to implement the hack in my comment:37, but this is much simpler and will make the bots happy.

[kliem]

comment:41

I failing doctest, due to the fact that there is a new submodule of geometry.

File "src/sage/misc/dev_tools.py", line 164, in sage.misc.dev_tools.load_submodules
Failed example:
    sage.misc.dev_tools.load_submodules(sage.geometry, "database$|lattice")
Expected:
    load sage.geometry.fan_isomorphism... succeeded
    load sage.geometry.hyperplane_arrangement.affine_subspace... succeeded
    ...
    load sage.geometry.riemannian_manifolds.surface3d_generators... succeeded
Got:
    load sage.geometry.cone_catalog... succeeded
    load sage.geometry.fan_isomorphism... succeeded
    load sage.geometry.hyperbolic_space.hyperbolic_coercion... succeeded
    load sage.geometry.hyperbolic_space.hyperbolic_constants... succeeded
    load sage.geometry.hyperbolic_space.hyperbolic_geodesic... succeeded
    load sage.geometry.hyperbolic_space.hyperbolic_interface... succeeded
    load sage.geometry.hyperplane_arrangement.affine_subspace... succeeded
    load sage.geometry.hyperplane_arrangement.arrangement... succeeded
    load sage.geometry.hyperplane_arrangement.check_freeness... succeeded
    load sage.geometry.hyperplane_arrangement.library... succeeded
    load sage.geometry.hyperplane_arrangement.plot... succeeded
    load sage.geometry.newton_polygon... succeeded
    load sage.geometry.polyhedron.cdd_file_format... succeeded
    load sage.geometry.polyhedron.combinatorial_polyhedron.base... succeeded
    load sage.geometry.polyhedron.double_description... succeeded
    load sage.geometry.polyhedron.double_description_inhomogeneous... succeeded
    load sage.geometry.polyhedron.face... succeeded
    load sage.geometry.polyhedron.library... succeeded
    load sage.geometry.polyhedron.plot... succeeded
    load sage.geometry.pseudolines... succeeded
    load sage.geometry.ribbon_graph... succeeded
    load sage.geometry.riemannian_manifolds.parametrized_surface3d... succeeded
    load sage.geometry.riemannian_manifolds.surface3d_generators... succeeded

[sagetrac-git]

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

bfd757f

Trac #26623: fix a doctest for sage.misc.dev_tools.load_submodules().

[sagetrac-git]

Changed commit from e8fd707 to bfd757f

[orlitzky]

comment:43

Replying to @kliem:

I failing doctest, due to the fact that there is a new submodule of geometry.

Oof, fixed. The switch back to a lazy import caused that one.

[kliem]

comment:44

Again many minor things.

• There is a tripple dash

+At the moment, only convex rational polyhedral cones are
+supported---specifically, those cones that can be built using the

• Instead of importing in each function, you could do something as

from sage.rings.all import ZZ as _ZZ

Don't know if that is better.

• I don't know, if this is going to link correctly in the documentation:

+    - ``lattice`` - a :class:`ToricLattice` in which the cone will live

• I think we are supposed to use double quotation marks for error messages:

+        raise ValueError('lattice rank=%d and ambient_dim=%d '
+                         'are incompatible' % (lattice.rank(), ambient_dim))

• Missing semicolons:

-    - ``ambient_dim`` -- (default: ``None``) the dimension of the
-      ambient space, a nonnegative integer
-
-    - ``lattice`` -- (default: ``None``) the toric lattice in which
-      the cone will live
+    - ``ambient_dim`` -- (default: ``None``); the dimension of the
+      ambient space, a nonnegative integer
+
+    - ``lattice`` -- (default: ``None``); the toric lattice in which
+      the cone will live

• Again I'm not sure that this will work in the documentation

+    A :class:`.ConvexRationalPolyhedralCone` living in ``lattice``

• I don't know if I would mention the errors after OUTPUT at all. Isn't this clear that this OUTPUT is given only if the INPUT is as required and otherwise an error is raised?

-    I = matrix.identity(ZZ,ambient_dim)
+    I = matrix.identity(ZZ, ambient_dim)

-    I = matrix.identity(ZZ,ambient_dim)
-    M = matrix.ones(ZZ,ambient_dim) - p*I
+    I = matrix.identity(ZZ, ambient_dim)
+    M = matrix.ones(ZZ, ambient_dim) - p*I

• The rearrangment cone will give nonsense, when p is not an integer. But I don't know if this is a problem.

[sagetrac-git]

Changed commit from bfd757f to ebb5a76

[sagetrac-git]

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

7660105	Trac #26623: use docstring conventions from the developer's guide.
f7e8af6	Trac #26623: add a few spaces to the cone catalog for readability.
e91d993	Trac #26623: make rearrangement cones of non-integer orders an error.
ebb5a76	Trac #26623: use double-quotes as first choice for string delimeters.

[orlitzky]

comment:46

Replying to @kliem:

Again many minor things.

• There is a tripple dash

This is actually the documented way to get sphinx to create an "em" dash in the HTML output. It's displaying correctly for me in the built docs.

• Instead of importing in each function, you could do something as

from sage.rings.all import ZZ as _ZZ

Don't know if that is better.

Clever =)

But, the extra imports are "free", so I don't think it's worth changing them around again at this point.

• I don't know, if this is going to link correctly in the documentation:

+    - ``lattice`` - a :class:`ToricLattice` in which the cone will live

You are right... ToricLattice isn't a class.

• I think we are supposed to use double quotation marks for error messages:

Is that written down anywhere (for my future reference)? I don't care either way, I changed them to double quotes.

• Missing semicolons:

Hmm.. someone has updated the developer's guide without fixing the existing docstrings or mentioning it to anyone. I tried to conform to the new examples in the guide.

• Again I'm not sure that this will work in the documentation

+    A :class:`.ConvexRationalPolyhedralCone` living in ``lattice``

This one works!

• I don't know if I would mention the errors after OUTPUT at all. Isn't this clear that this OUTPUT is given only if the INPUT is as required and otherwise an error is raised?

Other languages/projects have a standard location to document the exceptions that can be raised, but as far as I know, we don't. We probably should. I think it's nice to mention that a function can fail in the OUTPUT docs, just to be on the safe side in case the user skips the INPUT docs. But to me it feels like the relationship between the inputs belongs in the INPUT block. My mind could be changed on this, if everyone could agree to do it consistently.

-    I = matrix.identity(ZZ,ambient_dim)
+    I = matrix.identity(ZZ, ambient_dim)

-    I = matrix.identity(ZZ,ambient_dim)
-    M = matrix.ones(ZZ,ambient_dim) - p*I
+    I = matrix.identity(ZZ, ambient_dim)
+    M = matrix.ones(ZZ, ambient_dim) - p*I

Added spacing throughout.

• The rearrangement cone will give nonsense, when p is not an integer. But I don't know if this is a problem.

Thanks, I don't like that the function would silently return nonsense if you pass in e.g. p=1.5 I've added another error check to that cone.

[kliem]

Reviewer: Jonathan Kliem

[kliem]

comment:47

Missing comma for the if ... then statement:

+    Jeong's Corollary 5.2.4 [Jeong2017]_ states that if `p = n-1` in
+    an `n`-dimensional ambient space then the Lyapunov rank of the

Otherwise, you can set it on positive review on my behalf.

[sagetrac-git]

Changed commit from ebb5a76 to b2aab91

[sagetrac-git]

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

24d97c0	Trac #26623: add cones. shortcuts for common cones.
1f454d1	Trac #26623: add sage.geometry.cone_catalog to the documentation index.
12ee26c	Trac #26623: update sage.geometry.cone doctests to use the cone catalog.
b2aab91	Trac #26623: fix a doctest for sage.misc.dev_tools.load_submodules().

[orlitzky]

comment:49

Replying to @kliem:

Missing comma for the if ... then statement:

...

Otherwise, you can set it on positive review on my behalf.

Thanks again. I squashed this and the other recent commits into the main one. It will be nice to have this off my to-do list.

[vbraun]

Changed branch from u/mjo/ticket/26623-ng to b2aab91