Base sage.geometry.cone on the Parma Polyhedra Library (PPL)

[vbraun]

As a first useful application of the PPL Cython library interface I have changed sage.geometry.cone.Cone to use the PPL wrapper instead of cddlib. Here is a quick benchmark with a fan that was somewhat challenging:

sage: from sage.schemes.generic.toric_variety_library import toric_varieties_rays_cones
sage: rays, cones = toric_varieties_rays_cones['BCdlOG']
sage: timeit('Fan(cones,rays)')
5 loops, best of 3: 1.95 s per loop

With the old Polyhedron/cddlib interface, I got instead

5 loops, best of 3: 42.1 s per loop

Apply

1. attachment: trac_10140_sublattice_intersection.patch
2. attachment: trac_10140_base_cone_on_ppl_original.patch
3. attachment: trac_10140_reviewer.patch
4. attachment: trac_10140_switch_point_containment_to_PPL.patch

Apply trac_10140_sublattice_intersection.patch, trac_10140_base_cone_on_ppl_original.patch, trac_10140_reviewer.patch, trac_10140_switch_point_containment_to_PPL.patch

Depends on #10039

CC: @novoselt @nthiery

Component: geometry

Keywords: ppl

Author: Volker Braun, Andrey Novoseltsev

Reviewer: Andrey Novoseltsev, Volker Braun

Merged: sage-4.7.1.alpha1

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

[vbraun]

comment:1

In the current patch I fixed all doctest errors in cone.py and fan.py. There are some more broken doctests in other modules due to a different output ray ordering, and I will fix these once the ppl spkg and wrapper are reviewed.

[novoselt]

comment:2

Wonderful speed gain! However, when do you expect PPL to become a standard library? It is my impression that currently it is difficult to add new standard packages and it is more likely to become an optional one for a while. In which case there should be an option to use PPL, but its presence should not affect work of toric geometry.

Little picks at the patch:

1. I don't see why it is necessary to remove the possibility of Cone construction from a Polyhedron.
2. I think that new possible input C_Polyhedron should be documented in the INPUT section rather than a note in the end.
3. Is it possible to make C_Polyhedrons immutable (on demand)?
4. I prefer all line generators to be put in the end of ray list for non-strictly convex cones, if they were determined automatically. Because if one works with faces of such a cone, then all these generators are the same and, therefore, others are more "important". That's very minor, of course.
5. I am trying to follow PEP8 style guide http://www.python.org/dev/peps/pep-0008/ which says: When raising an exception, use "raise ValueError('message')" instead of the older form "raise ValueError, 'message'". New line 1167 tries to change this recommended form to an old one ;-)
6. Looking at what you have done in facet normals, maybe we should change cone.dual_lattice() method to always return ZZⁿ if there is no dual method for the lattice? Thinking of the class group situation, it seems like a more sensible default.
7. I think cones should not be printed explicitly (i.e. using "print") when ValueError is raised due to attempt to intersect cones in lattices of different dimension. The reason is that maybe some users want to intercept this exception and deal with it in their own way without showing any output. Also, I think that the check should be done not for equality of lattice dimensions, but for equality of lattices themselves, because cones of different lattices should not be intersected.

[vbraun]

comment:3

1. Well you can always write Cone(Polyhedron.rays()). I took out the possibility to directly pass a Polyhedron to make sure that I caught all instances of old code.

2. PPL stuff should not be visible to the user, I think. If anything, I'll probably move the note to a comment in the code.

3. In the C++ library there is no particular support for that. Of course, in C++ you just declare things const if you don't want them mutable. So I guess we should emulate that by get/set_immutable() methods.

4,5,6. Sounds good.

7. Oops left-over print() from debugging things :-)

[novoselt]

comment:4

1. My concern was that it removes publicly exposed feature that has been in Sage for several months (and is likely to stay there a few months more). So I vote for resurrecting this possibility, which is actually very easy when things are PPL based - you just replace Polyhedron with Polyhedron.rays(), as you said. More importantly, if PPL will become something that users have to install separately, then ALL of the old code related to polyhedron -> cone should remain untouched, but there should be some new path to go around it. By the way, I am not quite familiar with packages - is it possible to change Sage library during package installation? If not and all PPL-related code should be in the library no matter whether PPL is installed or not, then "doctest-fixing" may be a bad approach here. Maybe we should change doctests to make them deterministic. Although I am not quite sure how.

   2. I see your point, but having non-documented input rubs me the wrong way. If you don't want to mention PPL in the documentation (which is quite reasonable), maybe we can add a function like _Cone_from_PPL for such input? For the actual class there is no need in this - it is already not supposed to be called directly by users.

   3. Yes, that's what I had in mind. My concern is that when you use PPL representation, you may need to actually change it somehow to get the result. So far in your code you were constructing a new object based on existing one, so that the original is intact, but it seems to me that this step is quite easy to forget and it would be nice to have something forcing you to do so.

   4. That explains it ;-)

[jdemeyer]

comment:5

This is quite possibly a stupid question, but does it mean we can get rid of cddlib if this gets merged?

[sagetrac-mhampton]

comment:6

No! Currently cddlib is required by Gfan. It may be possible to convince Anders Jensen, the author of Gfan, to switch to PPL as well.

[vbraun]

comment:7

The updated version now cleans up all doctests.

[vbraun]

comment:8

For the patch bot:

Depends on #10233, #10039

Apply trac_10140_base_cone_on_ppl.patch, trac_10140_fix_toric_variety_doctests.patch

[vbraun]

comment:9

For the patch bot:

Depends on #9972, #10233, #10039.

Apply trac_10140_base_cone_on_ppl.patch, trac_10140_fix_toric_variety_doctests.patch

[kiwifb]

comment:10

I have trouble applying trac_10140_base_cone_on_ppl.patch. The following section doesn't apply at least:

@@ -1578,7 +1667,7 @@
         
         Let's take a 3-d cone on 4 rays::
 
-            sage: c = Cone([(1,0,1), (0,1,1), (-1,0,1), (0,-1,1)])
+            sage: c = Cone([(1,0,1), (0,1,1), (-1,0,1), (0,-1,1)], check=False)
             
         Then any ray generates a 1-d face of this cone, but if you construct
         such a face directly, it will not "sit" inside the cone::

Is it due to #9972? I can see that particular section in #9972 but slightly shifted in line number.

[novoselt]

comment:11

It is quite likely, since #9972 was invasive and written mostly by me, i.e. it was not coordinated with this one. What exactly do you mean by "have trouble"? I just tried to apply it on top of my queue and it succeeded with two fuzzy hunks.

For the record: I am going to review this ticket but I am waiting for PPL package to get approved ;-)

[kiwifb]

comment:12

Patching failed miserably. But I hadn't apply any of #9972. Since #9972 seems to going in 4.6.2, I will only apply the present patches from a 4.6.2.

[novoselt]

Work Issues: rebasement

[novoselt]

comment:13

This is on a fresh installation of 4.6.2.alpha2 without any patches applied:

applying trac_10140_base_cone_on_ppl.patch
patching file sage/geometry/cone.py
Hunk #30 succeeded at 1981 with fuzz 1 (offset 10 lines).
Hunk #39 succeeded at 2332 with fuzz 2 (offset 29 lines).
now at: trac_10140_base_cone_on_ppl.patch

[vbraun]

Attachment: trac_10140_fix_toric_variety_doctests.patch.gz

Rebased patch

[vbraun]

comment:14

I've uploaded my rebased patch.

[novoselt]

comment:15

François, does the new version work for you? (Now I don't get any fuzz for both patches.)

[novoselt]

Changed work issues from rebasement to none

[novoselt]

Changed work issues from remove/review dependencies to none

[novoselt]

comment:47

Attachment: trac_10140_sublattice_intersection.patch.gz

OK, the new patches apply cleanly to sage-4.7.alpha4 and pass all long tests, documentation builds without warnings. The first patch now includes some adjustments to lattice treatment, moved here from #10882 (I'll update the patch there shortly.)

Volker, if you are fine with my patches, please switch it to positive review!

[novoselt]

comment:48

Attachment: trac_10140_switch_point_containment_to_PPL.patch.gz

One more patch: I was learning how to use PPL and trying to switch point containment check in cones so that it does not call polyhedron method. In the process I have discovered a bug with constructing cones without rays, i.e. like Cone([], lattice=N): the PPL representation in this case didn't know the ambient space of this cone leading to mistakes. It is fixed in the second hunk of the last patch.

Regarding original goal, here are the timings before

sage: timeit("c = Cone([(1,0), (0,1)]); (1,1) in c", number=1000)
1000 loops, best of 3: 27.8 ms per loop
sage: c = Cone([(1,0), (0,1)])
sage: timeit("(1,1) in c", number=1000)
1000 loops, best of 3: 729 µs per loop

and after

sage: timeit("c = Cone([(1,0), (0,1)]); (1,1) in c", number=1000)
1000 loops, best of 3: 2.3 ms per loop
sage: c = Cone([(1,0), (0,1)])
sage: timeit("(1,1) in c", number=1000)
1000 loops, best of 3: 290 µs per loop

As we see, even on very simple example we get 10x speedup for "single checks" when most of the time is spent constructing different representations of the cone. When everything is precached and we count only actual containment, we have 3x speed up.

A more complicated example before:

sage: c = Cone([(4, 0, 1, 12, 1), (-2, -2, -73, 1, 0), (1, -2, 0, 1, -3), (5, -1, -1,
sage: 1, 1), (-4, 5, 1, 2, 3), (0, -1, -23, 0, 1), (-2, 1, -1, 1, -1), (0, 2,
sage: -3, 1, 0), (-1, 3, 2, -1, 0), (0, 2, -1, 4, 15)])
sage: c.ray_matrix()
[ -4  -2  -2  -1   0   5   0   1]
[  5  -2   1   3  -1  -1   2  -2]
[  1 -73  -1   2 -23  -1  -1   0]
[  2   1   1  -1   0   1   4   1]
[  3   0  -1   0   1   1  15  -3]
sage: timeit("(1,2,3,4,5) in c", number=1000)
1000 loops, best of 3: 4.52 ms per loop

and after

sage: c = Cone([(4, 0, 1, 12, 1), (-2, -2, -73, 1, 0), (1, -2, 0, 1, -3), (5, -1, -1,
sage: 1, 1), (-4, 5, 1, 2, 3), (0, -1, -23, 0, 1), (-2, 1, -1, 1, -1), (0, 2,
sage: -3, 1, 0), (-1, 3, 2, -1, 0), (0, 2, -1, 4, 15)])
sage: timeit("(1,2,3,4,5) in c", number=1000)
1000 loops, best of 3: 1.3 ms per loop

(I didn't bother with fresh start timing here.)

Conclusion: Volker's PPL wrapper rocks!

[vbraun]

comment:49

I'm happy with the reviewer patch, so positive review alltogether ;-)

[jdemeyer]

comment:51

@ Andrey Novoseltsev: Can you upload your reviewer patch again using hg export tip? The user and date fields are missing.

[novoselt]

comment:52

Attachment: trac_10140_reviewer.patch.gz

Done!

[jdemeyer]

Merged: sage-4.7.1.alpha1

[vbraun]

comment:54

I'm adding an unformatted "Apply" section in the ticket description to help the patch buildbot figure out the correct patches.

[vbraun]

comment:55

Apply trac_10140_sublattice_intersection.patch, trac_10140_base_cone_on_ppl_original.patch, trac_10140_reviewer.patch, trac_10140_switch_point_containment_to_PPL.patch

[dimpase]

comment:56

How hard would it be to set up PPL to do Vrepresentation and Hrepresentation of Polyhedron?

[novoselt]

comment:57

Replying to @dimpase:

How hard would it be to set up PPL to do Vrepresentation and Hrepresentation of Polyhedron?

See #11634. Not easy, I guess, but Volker has done it.