Add PointCollection

[novoselt]

This patch adds support for a PointCollection class, which is not yet used for anything, but I propose switching cones and fans to it and plan to extensively use it in the rewrite of lattice polytopes.

Module level documentation explains the benefits, basically, it is an "enhanced but still fast tuple" which will allow uniform treatment of things like ray generators, lattice points, and normals. Instead of writing

sage: cone.rays()
sage: cone.ray_matrix()
sage: cone.ray_set()
sage: cone.lines()
sage: cone.line_matrix()
sage: cone.line_set()

which requires 6 different methods (and more for facet normals and vertices/points of polytopes), one will be able to use

sage: cone.rays()
sage: cone.rays().matrix()
sage: cone.rays().set()
sage: cone.lines()
sage: cone.lines().matrix()
sage: cone.lines().set()

which uses only special methods for cones (plus 2 for point collections).

Another good point will be elimination of IntegralRayCollection class which is currently the common base of cones and fans, which seemed to be convenient, but is definitely a bit weird.

To make sure that the transition is painless, this ticket will only add the new class, while transition will be done in a consecutive one, with all the proper deprecations etc. The doctests show how cone.rays() will eventually behave.

Regarding speed: there is Sequence class which does a similar thing, "Sagifying" lists and tuples. But it is very slow in construction (1000 times with checks on, 20 times with checks off compared to tuples, my class is only 2 times slower) and two times slower than tuples in access (mine is as fast as tuples):

sage: from sage.geometry.point_collection import PointCollection
sage: f = toric_varieties.dP6xdP6().fan()
sage: timeit("PointCollection(f.rays())", number=1000, repeat=100)
1000 loops, best of 100: 3.59 µs per loop
sage: timeit("tuple(f.rays())", number=1000, repeat=100)
1000 loops, best of 100: 1.62 µs per loop
sage: timeit("seq(f.rays())", number=100, repeat=10)
100 loops, best of 10: 1.49 ms per loop
sage: timeit("seq(f.rays(), f.lattice(), check=False)", number=1000, repeat=100)
1000 loops, best of 100: 39.3 µs per loop
sage: pc = PointCollection(f.rays())
sage: timeit("pc[5]", number=1000, repeat=100)
1000 loops, best of 100: 561 ns per loop
sage: pct = tuple(f.rays())
sage: timeit("pct[5]", number=1000, repeat=100)
1000 loops, best of 100: 547 ns per loop
sage: pcs = seq(f.rays(), f.lattice(), check=False)
sage: timeit("pcs[5]", number=1000, repeat=100)
1000 loops, best of 100: 1.16 µs per loop

While plain tuples are still the fastest, I think that benefits of the proposed class worth the associated penalty (which is not the case for sequences). In fact, I think that overall speed will improve due to better cache sharing.

CC: @vbraun @hivert

Component: geometry

Author: Andrey Novoseltsev

Reviewer: Volker Braun

Merged: sage-5.0.beta6

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

[hivert]

comment:2

Hi Andrey,

I got two questions:

• are those objects instance of a subclass of Element, do they have a parent ?
• are they immutable as tuples ?

If so, Its fairly possible that we would use them to represent element of cartesian products. They are currently represented as tuple so they don't have a parent.

[novoselt]

comment:3

Hi Florent,

Currently they derive from SageObject and don't have a parent. I am not sure what would be the parent - they are arbitrary collections of points of the same space, repetitions are allowed. So parent is the set of finite sequences of the underlying module. Do we have it in Sage?

Implementation-wise, there is an internal tuple and I have added those access methods which I definitely need. It is of course possible to add slicing etc, but definitely they will remain immutable.

Note also that these are collections of points of the same space and the whole class is written with this restriction in mind, which is probably not what you want for cartesian products. So maybe the best option for them would be to implement something similar but with necessary tweaks. As you can see from my patch, it is quite simple and short, I have spent most of the time figuring out how to do things in Cython (which I am still a beginner in) and running tests to make sure that the speed is close to tuples.

I thought of deriving from Element to use cache from Simon's new patches, but his work is not merged yet. I suppose we can derive from Element and try to split this new class into two, one for "arbitrary points" and one for points of the same space. I am a little worried about the overhead, but perhaps in Cython it will not get much worse.

Thanks for the interest!
Andrey

[novoselt]

comment:4

I'll switch to needs info in case we decide to do the splitting, but the patch should work fine and without any conflicts as it is a new module.

Volker, what do you think about the proposal to switch to these containers and deriving from Element?

[vbraun]

comment:5

Can we also abstract points, so that PointCollection : Point == Parent : Element?

If you need a really large number of points, then the coordinates of all points should be stored in a numpy array. So PointCollection would own the array and Point would just reference a row in the array. I don't think this is necessary at this point. But if we abstract Points then we can change its implementation it in the future without changing any code that builds on it.

[novoselt]

comment:6

I am a bit paranoid about adding yet another intermediate level... How would it perform on really simple examples? One of the most important cases for us are simplicial cones, so they should build and work as fast as possible.

I also think that points should not link to collections since the same point may belong to multiple collections (e.g. points of a face are also points of the whole thing).

And a natural parent of points seems to be the space in which they live, not the sequence in which they have happened to appear.

In a way this class already abstracts the way in which points are stored and say, instead of acting by some transformation matrix on single points or the matrix representation of a collection we may add direct action on the collection itself with whatever internal implementation (that's actually something that I was going to add but forgot). So if we use a numpy array, then we can live without any kind of stored objects to represent points and if someone really wants to access a single one - __getitem__ will construct a corresponding vector.

[novoselt]

comment:7

Volker: I am switching it back to needs review as adding new methods and actions is not going to happen in the near future and I was not yet going to change anything existing in this patch. So let me know if there are any picks that have to be fixed and let's leave the rest to the next ticket (switching cone.rays, cone.facet_normals, etc. to this class).

Florent: apart from wrapping tuple special methods, it still seems to me that goals of this class and elements of Cartesian products are quite orthogonal, and assumptions on components are completely different, in particular this class can be switched from tuples to numpy arrays, as Volker suggested, quite easily, because all components live in the same space. So I don't think it makes sense to unite these classes, but on the other hand this ticket shows that tuples can be wrapped efficiently. Maybe one can do something even better, but I didn't want to mess with memory allocation and access myself.