Add support for rational polyhedral fans

[novoselt]

This patch is a part of the following series adding support for cones/fans and toric varieties to Sage:

Prerequisites:

#8675 - Remove AmbientSpace._constructor and fix consequences

#8682 - Improve AlgebraicScheme_subscheme.__init__ and AmbientSpace._validate

#8694 - Improve schemes printing and LaTeXing

#8934 - Trivial bug in computing faces of non-full-dimensional lattice polytopes

#8936 - Expose facet inequalities for lattice polytopes

#8941 - _latex_ and origin for lattice polytopes

Main patches adding new modules:

#9062 - Add support for toric lattices

#8986 - Add support for convex rational polyhedral cones

#8987 - Add support for rational polyhedral fans (this ticket now also has a bug fix #9188 as a prerequisite)

#8988 - Add support for toric varieties

#8989 - Add support for Fano toric varieties

Everything was tested on sage.math using sage-4.4.2.rc0.

CC: @vbraun @loefflerd

Component: geometry

Author: Andrey Novoseltsev

Reviewer: Volker Braun

Merged: sage-4.5.2.alpha0

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

[novoselt]

Author: Andrey Novoseltsev

[novoselt]

comment:2

I will make some adjustments to this module following discussion here:
http://groups.google.com/group/sage-devel/browse_thread/thread/17743e17d67838ae

[novoselt]

Switch to toric lattices.

[novoselt]

comment:3

Attachment: trac_8987_add_support_for_rational_polyhedral_fans.patch.gz

New version of the patch using toric lattices from #9062. Switch of caching technique to allow efficient extension of class hierarchy is still pending.

[vbraun]

comment:4

Looks good so far. Adding the extra caching that was discussed on the mailinglist could be done later. Let me know when you want me to review the patches.

[novoselt]

comment:5

Hi Volker,

If you don't mind some delay of caching, then cones are ready for review.

In fans module I will add today or tomorrow derived classes EnhancedCone(Cone_of_fan) and EnhancedFan(RationalPolyhedralFan), with the syntax discussed on sage-devel before. For now they will make efficient copies of existing fans and cones, later they will be smoothly transitions to shared caching. Since these classes will be on top of the existing framework, the rest is ready to review. In fact, I can make a separate patch adding them in a new ticket, rather than modifying again patches here. Marking everything as ready!!!

Thank you!
Andrey

[novoselt]

comment:6

The second patch adds classes that will support cones and fans for domain and codomain of morphisms. While working on it, I moved the actual computation of the cone lattice to _compute_cone_lattice. I have also added a comparison function to cones that ignores types if both arguments are cones. I think that cone1 == cone2 should NOT take into account if these cones belong to fans, morphisms, etc. If their rays are the same - they are equal.

[novoselt]

Attachment: trac_8987_add_enhanced_cones_and_fans.patch.gz

Apply after the big patch.

[vbraun]

Reviewer: Volker Braun

[vbraun]

comment:7

Maybe I'm missing something, but is there a reason why comparison of cones depends on the ray ordering?

sage: Cone([(1,0), (0,1)]) == Cone([(1,0), (0,1)]) 
True
sage: Cone([(1,0), (0,1)]) == Cone([(0,1), (1,0)]) 
False

It seems more natural to compare ray_set() instead of rays().

[novoselt]

comment:8

Replying to @vbraun:

Maybe I'm missing something, but is there a reason why comparison of cones depends on the ray ordering?

sage: Cone([(1,0), (0,1)]) == Cone([(1,0), (0,1)]) 
True
sage: Cone([(1,0), (0,1)]) == Cone([(0,1), (1,0)]) 
False

It seems more natural to compare ray_set() instead of rays().

This is done like this:

sage: Cone([(1,0), (0,1)]).is_equivalent(Cone([(0,1), (1,0)]))
True

One reason for such behaviour of == is simplicity of comparison, which probably should be fast for sorting purposes (and I was implementing only __cmp__ for all classes, not __eq__). In the case of not strictly convex cones equality of ray sets is not necessary for equality of cones as sets of points (see the code of is_equivalent for cones).

Another reason is consistency with the implementation of fans. For fans checking mathematical equality is a bit more tedious (although still can be done quickly using sorted lists for rays and cones, see is_equivalent for fans), but the main reason is that I would prefer to be able to use plain integers to index divisors or charts. I was participating in the discussion how the equality of fans should be interpreted in Macaulay2, and in the documentation to one of the cohomology related functions there was a line about "canonical identification with ZZ^n". Well, such identifications can be canonical only if the order of things is fixed.

Maybe these arguments are not very strong, but since there is a way to check equality in both senses, I think that it should be fine. The plan is also to have equivalence check relevant for toric varieties, i.e. up to GL(ZZ^n) action, but is_isomorphic functions throw NotImplementedError so far.

Would you like me to add some comments on this behaviour to the "main" documentation of cone and fan modules? It is already described in is_equivalent and is_isomorphic docstrings, but perhaps deserves more visibility.

[vbraun]

comment:9

1. I agree that we should test == without identifying GL(n,Z) images as it is expensive. But I'm afraid of the following, which is probably a common use case:

sage: diamond = lattice_polytope.octahedron(2) 
sage: P1xP1 = FaceFan(diamond)
sage: Cone([(0,1), (1,0)]) in P1xP1(dim=2)   
False
sage: Cone([(1,0), (0,1)]) in P1xP1(dim=2)   
True

Also, note that Fan() is implicitly sorting:

sage: fan1 = Fan([Cone([(1,0), (0,1)])])
sage: fan2 = Fan([Cone([(0,1), (1,0)])])
sage: fan1==fan2
True

I would argue that is is ok to treat fans and cones a bit different, one will compare cones often but fans only seldomly. I'm happy with comparison of fans to depend on the order of the rays and generating cones, otherwise all derived quantities (like cohomology generators of toric varieties) will not be equal (only isomorphic). But for cone1==cone2 to depend on the order of the rays is neither helpful nor intuitive.

2. Fan() should always error out if the cones are not generating cones. In particular, this one from the documentation

sage: cone1 = Cone([(1,0), (0,1)])
sage: cone2 = Cone([(0,1), (-1,-1)])
sage: cone3 = Cone([(-1,-1), (1,0)])
sage: P2 = Fan([cone1, cone2, cone2])
sage: P2.ngenerating_cones() 
2

should have raised an exception.

In normal use you'll never want to type in all the cones of the fan since there are so many. But its easy to get confused and add a cone that turns out to be not a generating cone, and its nice to catch this when generating the Fan.

There is certainly a need for a function that extracts the generating cones from a collection of cones, but I don't think it should be implicit in the Fan constructor.

3. Rename Cone_of_fan.fan_generating_cones() to star_generators().

There are similar methods to this one that we probably want to add later, so let me just throw out some names:

• Cone_of_fan.faces(): all subcones of the cone
• Cone_of_fan.facets(): subcones of one dimension lower
• Cone_of_fan.bounds(): supercones of one dimension higher
• Cone_of_fan.star(): the star of the cone
• Cone_of_fan.adjacent(): Adjacent cones (of the same dimension)

4. Do we need to know the set of ray indices, that is set(cone.fan_rays()), of a cone often? Right now there is the function cone_to_rays(cone_index) that is only used as a helper in cone_lattice(). If it is just a helper function, then it should be _cone_to_rays(). If you want to expose that functionality, it should be stored in the Cone_of_fan and retrieved via cone.fan_rays_set() or so.

I also think that cone.rays_idx() would be a better name than cone.fan_rays(), but if you disagree then I can live with the current name as well :-)

5. similarly to 4), ray_to_cones() is not very self-explanatory. We obviously need a way to find out which cones contain a given set of rays. I would prefer one (or all) of the following

• ray_to_cone(ray_index): the 1-cone spanned by the ray
• smallest_cone_containing(*Nlist): the smallest cone containing all points (which can be specified by a ray index or as a N-lattice element).

The functionality of ray_to_cones would then be provided by ray_to_cone(i).star_generators().

Let me know what you think...

[vbraun]

comment:10

I noticed that containing_cones(ray_indices) does essentially what I wanted in 5) in my previous comment. However I think it would be much more useful if it would return the actual cone object and not a set of ray indices. Generally speaking, I think the user should never need the ray indices. And if he still needs them, then they can be gotten easily enough from the Cone_of_fan object.

[novoselt]

comment:11

Replying to @vbraun:

1. I agree that we should test == without identifying GL(n,Z) images as it is expensive. But I'm afraid of the following, which is probably a common use case:

sage: diamond = lattice_polytope.octahedron(2) 
sage: P1xP1 = FaceFan(diamond)
sage: Cone([(0,1), (1,0)]) in P1xP1(dim=2)   
False
sage: Cone([(1,0), (0,1)]) in P1xP1(dim=2)   
True

You got a point here, I don't like this behaviour...

Also, note that Fan() is implicitly sorting:

sage: fan1 = Fan([Cone([(1,0), (0,1)])])
sage: fan2 = Fan([Cone([(0,1), (1,0)])])
sage: fan1==fan2
True

There was no sorting, rays were determined first as a set (union of ray sets of all cones) and then converted to a tuple. I am not sure if the same set will always lead to the same tuple or it depends on something. I definitely don't want to force sorting - if a user gives rays in a specific order, then (s)he probably wants them in that order.

I would argue that is is ok to treat fans and cones a bit different, one will compare cones often but fans only seldomly. I'm happy with comparison of fans to depend on the order of the rays and generating cones, otherwise all derived quantities (like cohomology generators of toric varieties) will not be equal (only isomorphic). But for cone1==cone2 to depend on the order of the rays is neither helpful nor intuitive.

In general, I agree, but it may potentially lead to cone1 == cone2 while Fan([cone1]) != Fan([cone2]). Are you OK with this?

As I understand, it is not enough to just add __eq__ method to cones, since it will conflict wiht __cmp__ in the sense that it will be possible to have cone1 < cone2 and cone1 == cone2 at the same time, which is not cool. In fact, I do not see any good (and fast) way to compare cones in agreement with mathematical equivalence. It will require constructing Polyhedron objects for each cone to check strict convexity and deal with possible non-uniqueness of ray generators, this will make it highly undesirable to sort sequences of cones. It will also make hashing of cones either slow or inconsistent with equality check (which is accepted in Sage, but should not be the case unless necessary).

Yeah... I think my position is: I will do this change if you still want it after arguments above, but I am very against it, especially since a method performing mathematical equality check is provided. A compromise is to add __contains__ and contains to the Fan class so that one can do

sage: diamond = lattice_polytope.octahedron(2) 
sage: P1xP1 = FaceFan(diamond)
sage: Cone([(0,1), (1,0)]) in P1xP1   # no (dim=2)   
True
sage: Cone([(1,0), (0,1)]) in P1xP1   # no (dim=2)   
True

It would be also nice, perhaps, to be able to check points in such a way, i.e. if a point is in the support of a fan.

2. Fan() should always error out if the cones are not generating cones. In particular, this one from the documentation

sage: cone1 = Cone([(1,0), (0,1)])
sage: cone2 = Cone([(0,1), (-1,-1)])
sage: cone3 = Cone([(-1,-1), (1,0)])
sage: P2 = Fan([cone1, cone2, cone2])
sage: P2.ngenerating_cones() 
2

should have raised an exception.

In normal use you'll never want to type in all the cones of the fan since there are so many. But its easy to get confused and add a cone that turns out to be not a generating cone, and its nice to catch this when generating the Fan.

There is certainly a need for a function that extracts the generating cones from a collection of cones, but I don't think it should be implicit in the Fan constructor.

My intention was to write a fan constructor that will construct a fan if at all possible. (Similarly, the cone constructor by default will discard non-generating rays.) One may, for example, start with some fan and add more cones to it. In this case some of the old generating cones may become unnecessary. Or, if cones of a fan come one by one from a certain procedure (e.g. in computing GKZ decomposition) and it is not known in advance how to easily select generating cones (e.g. all generating cones are full-dimensional for complete fans), it would be nice if the fan constructor could automatically select necessary cones. If you really against silent discard, we can perhaps either add a warning about a non-generating cone present, or a default parameter generating_cones=True to Fan(...) and then throw an exception if there are such cones. If one (like me ;-)) wants to discard them, it will be possible to use generating_cones=False option.

3. Rename Cone_of_fan.fan_generating_cones() to star_generators().

OK.

There are similar methods to this one that we probably want to add later, so let me just throw out some names:

• Cone_of_fan.faces(): all subcones of the cone
• Cone_of_fan.facets(): subcones of one dimension lower

Note that these are inherited from plain cones. However, they are likely to be ConeFace objects and for fans it would make more sense to return other Cone_of_fan ones...

• Cone_of_fan.bounds(): supercones of one dimension higher

I don't like the name. I guess it means that self bounds those that are returned by this function. But it can also be interpreted as that it returns bounds of self in the sense its facets. How about facet_of for this one?

• Cone_of_fan.star(): the star of the cone

What exactly do you mean by star here? Do you want to get back the fan in the quotient lattice? Anything else attached to it? I think we can use coercion system to define the canonical map from the lattice of the original fan to this one.

• Cone_of_fan.adjacent(): Adjacent cones (of the same dimension)

Is this name standard? I see that it may be convenient to have such a function, but I am not sure I would guess from the name what it does.

4. Do we need to know the set of ray indices, that is set(cone.fan_rays()), of a cone often? Right now there is the function cone_to_rays(cone_index) that is only used as a helper in cone_lattice(). If it is just a helper function, then it should be _cone_to_rays(). If you want to expose that functionality, it should be stored in the Cone_of_fan and retrieved via cone.fan_rays_set() or so.

I think it is just a helper and can go to _cone_to_rays to clean up the namespace.

I also think that cone.rays_idx() would be a better name than cone.fan_rays(), but if you disagree then I can live with the current name as well :-)

I disagree, because idx is cryptic (indices?) and fan_rays clearly tells one that it has something to do with the fan, even if it is not obvious that it will return indices rather than rays. (But what else can one expect? Rays themselves have no explicit relation to the fan.)

5. similarly to 4), ray_to_cones() is not very self-explanatory.

Agreed.

We obviously need a way to find out which cones contain a given set of rays. I would prefer one (or all) of the following

• ray_to_cone(ray_index): the 1-cone spanned by the ray

I was tempted to say that

sage: fan(1)[i]

will do exactly this, but it will not since fan(1) purposefully returns the list of rays, rather than one-dimensional cones... Did you notice this convention? What do you think of it?

Well, I think at least

sage: fan.cones(1)[i]

should do this work, even if it may fail to do so now because there is no particular sorting of 1-dimensional cones in the lattice. This actually seriously bugs me, since it was very annoying in LatticePolytope until I fixed it. Fixing it for cones and fans was on my todo list for the future, but I guess the best time to do it is now. Assuming that this will be done, I propose to get rid of ray_to_cones (or perhaps rename it to _ray_to_cones) and do not add ray_to_cone.

• smallest_cone_containing(*Nlist): the smallest cone containing all points (which can be specified by a ray index or as a N-lattice element).

Can I remove smallest_? I think since cone_containing promises to return a single cone, it is quite reasonable to expect that it will be smallest (and of course it will be written in the documentation).

The functionality of ray_to_cones would then be provided by ray_to_cone(i).star_generators().

Or, as I propose

fan.cones(1)[i].star_generators()

Let me know what you think now!

I am switching this ticket back to needs work. Regarding the last two tickets in the sequence - while I am working on this one, you can probably go ahead and review them too in the sense of comments. They will likely need some minor changes after updating this patch, but I currently don't plan to do anything else unless you request it.

[novoselt]

comment:12

Replying to @vbraun:

Generally speaking, I think the user should never need the ray indices. And if he still needs them, then they can be gotten easily enough from the Cone_of_fan object.

Yes. I tried to make cones and fans better compared to LatticePolytopes where everything is based on indices, but there is still room for improvement. Will work on it.

[vbraun]

comment:13

Replying to @novoselt:

There was no sorting, rays were determined first as a set (union of ray sets of all cones) and then converted to a tuple.

I think the integer vectors get sorted lexicographically by their entries, and I think the rays inherit that. So the result should be deterministic :-)

In general, I agree, but it may potentially lead to cone1 == cone2 while Fan([cone1]) != Fan([cone2]). Are you OK with this?

I think that right now Fan( cone_list_1 ) == Fan( cone_list_2 ) as long as the cones are in the same order even if the order of the rays in each cone is different. I think that would be good enough :-). If somebody by hand specifies different orders of the rays of the fan then its only fair to treat the fans as different.

The "mathematical" comparison of the cones is only slow if they are not strictly convex, so I don't think that would be a big problem. Also, in the non-strict case one could work with facet normals which might be cached already... I'm not necessarily insisting that we do it that way, but just throwing out some options.

2. I think the warning would be the fine, too.

3. I only want to traverse the face lattice with those methods. By "star", I mean without quotient (in the triangulation sense, not how star is used for fans usually). I'm open to other names :) And adjacent works in the same way as "adjacency matrix", but I don't have a reference at hand right now that uses that.

I was tempted to say that

sage: fan(1)[i]

will do exactly this, but it will not since fan(1) purposefully returns the list of rays, rather than one-dimensional cones... Did you notice this convention? What do you think of it?

Oww I didn't notice that... Can we have it return the 1-cones in the same order as the rays of the fan?

• smallest_cone_containing(*Nlist): the smallest cone containing all points (which can be specified by a ray index or as a N-lattice element).

Can I remove smallest_?

Fine by me!

[vbraun]

comment:45

I think the patch(es) are in great shape now. Positive review!

[novoselt]

comment:46

Thank you!

For the release manager: apply all patches in the order listed in the ticket.

[novoselt]

Attachment: trac_8987_cmp_fix.patch.gz

[novoselt]

comment:47

I made a systematic mistake in doctests of __cmp__ methods of all patches, discovered in #9062. So now I am posting these one-line patches to fix them, please review!

[vbraun]

comment:49

Looks good!

[loefflerd]

comment:50

This patch won't apply without #8986, which is currently at "needs work", so it can't be merged at present.

[qed777]

Merged: sage-4.5.2.alpha0

[qed777]

comment:54

One or more of #8986, #8987, and #9062 may cause the doctest failures listed at #9590.