Opened 4 years ago
Closed 3 years ago
#18860 closed enhancement (fixed)
Faster Poyhedron.graph()
Reported by:  ncohen  Owned by:  

Priority:  major  Milestone:  sage6.8 
Component:  geometry  Keywords:  
Cc:  vdelecroix, dimpase, vbraun, novoselt  Merged in:  
Authors:  Nathann Cohen  Reviewers:  Dima Pasechnik 
Report Upstream:  N/A  Work issues:  
Branch:  5223cda (Commits)  Commit:  5223cda9f15b2a417b550039b978a07e126f9053 
Dependencies:  #18779  Stopgaps: 
Description (last modified by )
With this branch, polytopes.Gosset_3_21().graph()
takes 300ms, against 16s before.
The 'definition' used to compute adjacent vertices can be found in 10
And of course the following works:
sage: g=graphs.GossetGraph() sage: g.is_isomorphic(polytopes.Gosset_3_21().graph()) True
Change History (39)
comment:1 Changed 4 years ago by
 Branch set to public/18860
comment:2 Changed 4 years ago by
 Commit set to 3b5eabb9735a867678c2c6a13dc7700c24d82732
comment:3 Changed 4 years ago by
 Cc vbraun added
comment:4 followup: ↓ 5 Changed 4 years ago by
+ # Why 26 ?? I was expecting 'd1'. Are there some useless inequalities + # in the list?
well, you have 26 facets on each edge  none are more useless than the others. You can do the following (note that 0th and 1st vertices are adjacent).
sage: p=polytopes.Gosset_3_21() sage: M=p.incidence_matrix() sage: a=0; b=1; sage: vv=[i for i in [0..702] if M[a,i]*M[b,i]>0] sage: matrix([p.inequalities()[j].vector() for j in vv]).T.kernel() Free module of degree 9 and rank 1 over Integer Ring Echelon basis matrix: [0 1 0 0 0 0 0 0 0]
or just
sage: m=matrix([p.inequalities()[j].vector() for j in vv]) sage: m.rank() 8
which makes sense, as m
is of size 27x9. (27, as there is not fulldimensional).
General code for adjacency of two vertices a
and b
of a polytope in d
dimensional space
should check that m
has corank 1
.
(There could be special cases when one of the vertices is the origin, or some inequalities are homogeneous).
comment:5 in reply to: ↑ 4 ; followup: ↓ 6 Changed 4 years ago by
General code for adjacency of two vertices
a
andb
of a polytope ind
dimensional space should check thatm
has corank1
.
Would you know how to adapt this code for that? Hoping that it will remain faster than what we currently have?
Nathann
comment:6 in reply to: ↑ 5 ; followup: ↓ 7 Changed 4 years ago by
 Cc novoselt added
Replying to ncohen:
General code for adjacency of two vertices
a
andb
of a polytope ind
dimensional space should check thatm
has corank1
.Would you know how to adapt this code for that? Hoping that it will remain faster than what we currently have?
My understanding is that hacking hasse_diagram.py
to stop earlier is what one actually needs to do.
Except that it goes from bottom to top, which probably means that one gets adjacency of facets as the first step. But it should be possible to change it that goes the other way around.
We can write a paper about it ;).
comment:7 in reply to: ↑ 6 ; followup: ↓ 8 Changed 4 years ago by
My understanding is that hacking
hasse_diagram.py
to stop earlier is what one actually needs to do.
I do not know how to do that. If you know how to fix this implementation, at least it will be done. Otherwise chances are that nothing will happen.
Nathann
comment:8 in reply to: ↑ 7 Changed 4 years ago by
Replying to ncohen:
My understanding is that hacking
hasse_diagram.py
to stop earlier is what one actually needs to do.I do not know how to do that. If you know how to fix this implementation, at least it will be done. Otherwise chances are that nothing will happen
the whole Polyhedron code does a potentially very slow thing: enumerating the inequalities at construction time, instead of doing it optionally or/and lazily.
Indeed, e.g. the problem at hand, deciding if two vertices are adjacent, does not need precomputed inequalities. Namely, checking that v_1 and v_2 are adjacent can be done by solving the following LP:
variables: x_0..x_d,x_{d+1} objective: max x_0 constraints: 1 <= x_i <=1 for i in [1..d]; x_{d+1} >= 0; sum_{j=1}^d v_{i,j}*x_j = x_{d+1}, for i=1,2; sum_{j=1}^d v_{i,j}*x_j >= x_{d+1} + x_0, for i>2;
if the optimum exists and is strictly positive, they form an edge, and otherwise they don't.
So we are trying to find an inequality that is valid for the convex hull of v_k's and turns into an equation on v_1 and v_2, but is a strict inequality on the rest of v_k's.
(The LP above will not work if there are more v_j's on the line joining v_1 and v_2  but this case can be taken care of trivially.)
comment:9 Changed 4 years ago by
Oh. That's a nice formulation of the problem, thanks. I will not forget it.
Would you know a combinatorial way to compute it? It seems that we can already rely on cddlib (see #18861) in order to build the graph over RDF.
Nathann
comment:10 followup: ↓ 19 Changed 4 years ago by
Here is a new 'combinatorial' version. It uses an 'idea', which I am not sure is a proper definition of the object, but well... You tell me :P
1) For every inequality I
, build the set set(I)
of points for which it is an equality
2) For any pair of points i,j
, list all inequalities I_1, ..., I_k
that contain them both
3) Intersection set(I_1), ..., set(I_k)
. If this intersection has cardinality 2, then ij is an edge.
Works for Gosset (at least) :P
Nathann
comment:11 Changed 4 years ago by
 Commit changed from 3b5eabb9735a867678c2c6a13dc7700c24d82732 to d084b1f2d0a7fa39fd522079aba39494f1704526
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
d084b1f  trac #18860: Faster Poyhedron.graph()

comment:12 followup: ↓ 21 Changed 4 years ago by
 Description modified (diff)
A bit cleaner: I moved the code to the .graph()
, method for apparently the _vertex_adjacency_matrix
is not what it claims to be. Its dimension may not even be the number of vertices (if you have rays n your polyhedron) [1].
With this, all tests pass in the geometry
folder.
Nathann
[1] You can see what goes wrong by adding return self.graph().adjacency_matrix()
in the first line of _vertex_adjacency_matrix
.
comment:13 Changed 4 years ago by
 Status changed from new to needs_review
comment:14 Changed 4 years ago by
 Commit changed from d084b1f2d0a7fa39fd522079aba39494f1704526 to 74d8412010748142918edf6178636bffdb14e60a
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
74d8412  trac #18860: Faster Poyhedron.graph()

comment:15 followup: ↓ 16 Changed 4 years ago by
The adjacency matrix for unbound polyhedra is essentially the projectvization; See also the adjacency matrix docs for a more downtoearth definition.
comment:16 in reply to: ↑ 15 Changed 4 years ago by
The adjacency matrix for unbound polyhedra is essentially the projectvization; See also the adjacency matrix docs for a more downtoearth definition.
I was not talking about the .adjacency_matrix
method. I was talking about ._vertex_adjacency_matrix
, which should (from its name) have as many rows/cols as there are vertices.
Nathann
comment:17 followup: ↓ 18 Changed 4 years ago by
Again its really the Vrepresentation adjacency matrix; The noncompact case is just the projectivization. So in a sense its really the vertex adjacency matrix; I know its somewhat confusing but most people are interested in the compact case so being overly nitpicky about wording is going to be confusing too.
comment:18 in reply to: ↑ 17 Changed 4 years ago by
Again its really the Vrepresentation adjacency matrix; The noncompact case is just the projectivization. So in a sense its really the vertex adjacency matrix; I know its somewhat confusing but most people are interested in the compact case so being overly nitpicky about wording is going to be confusing too.
So was this line incorrect return Graph(self.vertex_adjacency_matrix(), loops=False)
?
I expect the .vertex_graph()
to give me a graph whose vertices are the vertices of the polyhedron. And the incidence matrix of that graph is *not* the output of the _vertex_adjacency_matrix
whose rows are not the vertices only. Do we agree?
Nathann
comment:19 in reply to: ↑ 10 ; followup: ↓ 20 Changed 4 years ago by
Replying to ncohen:
Here is a new 'combinatorial' version. It uses an 'idea', which I am not sure is a proper definition of the object, but well... You tell me
:P
1) For every inequality
I
, build the setset(I)
of points for which it is an equality 2) For any pair of pointsi,j
, list all inequalitiesI_1, ..., I_k
that contain them both 3) Intersectionset(I_1), ..., set(I_k)
. If this intersection has cardinality 2, then ij is an edge.
yes, that's correct for the polytopes  provided there is no redundancy. E.g. if a polytope is not fulldimensional, then adding an equation satisfied by all the vertices to an inequality is a "copy" of this inequality, vertexvise.
Works for Gosset (at least)
:P
Nathann
comment:20 in reply to: ↑ 19 Changed 4 years ago by
yes, that's correct for the polytopes  provided there is no redundancy.
Excellent news! Then we have a very nice speedup ahead :PPPP
As for redundancy, it is true that this algorithm will not behave very well if there are several vertices at the same place which are actually the same. Redundancy of inequalities, however, is not a problem.
E.g. if a polytope is not fulldimensional, then adding an equation satisfied by all the vertices to an inequality is a "copy" of this inequality, vertexvise.
Yeah but in terms of sets it produces the very same one. So it's cool.
Nathann
comment:21 in reply to: ↑ 12 ; followups: ↓ 22 ↓ 23 Changed 4 years ago by
Replying to ncohen:
A bit cleaner: I moved the code to the
.graph()
, method for apparently the_vertex_adjacency_matrix
is not what it claims to be. Its dimension may not even be the number of vertices (if you have rays n your polyhedron) [1].With this, all tests pass in the
geometry
folder.
Aren't we changing the behaviour of some functions for the noncompact case?
comment:22 in reply to: ↑ 21 Changed 4 years ago by
Aren't we changing the behaviour of some functions for the noncompact case?
This patch only touches the function .vertex_graph
. Its behaviour changes indeed for the noncompact case, but unless you can defend that the previous behavior was not a bug (there were more vertices in the graph than vertices in the polyhedron) I would say that this is not a problem.
Nathann
comment:23 in reply to: ↑ 21 ; followup: ↓ 24 Changed 4 years ago by
Aren't we changing the behaviour of some functions for the noncompact case?
also note that not a single doctest changed. But to be honest, I still believe that the ._vertex_adjacency_matrix
should be either fixed or renamed. It is not what it claims to be.
Nathann
comment:24 in reply to: ↑ 23 ; followup: ↓ 25 Changed 4 years ago by
Replying to ncohen:
Aren't we changing the behaviour of some functions for the noncompact case?
also note that not a single doctest changed. But to be honest, I still believe that the
._vertex_adjacency_matrix
should be either fixed or renamed. It is not what it claims to be.
presently we have
adjacency_matrix = vertex_adjacency_matrix
in geometry/polhedron/base.py
. I'd think we should rename ._vertex_adjacency_matrix
to something like ._Vrep_adjacency_matrix
, and, dually ._facet_adjacency_matrix
to ._Hrep_adjacency_matrix
.
And then set adjacency_matrix = _Vrep_adjacency_matrix
, and provide a meaningful implementation of vertex_adjacency_matrix
, as you propose here.
comment:25 in reply to: ↑ 24 ; followup: ↓ 26 Changed 4 years ago by
Seems reasonable. Though I would very much like to see the definition of a V_rep adjacency matrix
.
Nathann
comment:26 in reply to: ↑ 25 Changed 4 years ago by
Replying to ncohen:
Seems reasonable. Though I would very much like to see the definition of a
V_rep adjacency matrix
.
well, check out sage.geometry.polyhedron.constructor?
. In short, any polyhedron is the Minkowski sum of a polytope P (spanned by vertices), a pointed polyhedral cone (spanned by rays), and a subspace S (spanned by lines).
 If S is empty or of dimension bigger than 1, it can be ignored.
 Otherwise S is a line, and it is adjacent to all the vertices, and to no ray; and, counterintuitively, no vertices are adjacent to each other.
 Rays are not adjacent to each other, and might be adjacent to some vertices.
How these general adjacencies are used, I don't know.
comment:27 Changed 4 years ago by
Oh. Then it seems that it will be possible to rewrite that method using the algorithm implemented here, then. And everything will be faster.
Nathann
comment:28 followup: ↓ 29 Changed 3 years ago by
Are we still in needs review
here? I struggle to remember what was the sticking point.
comment:29 in reply to: ↑ 28 Changed 3 years ago by
Are we still in
needs review
here? I struggle to remember what was the sticking point.
I do not think that there was any. We had discussions on the name of other functions, but this new implementation only touches the "vertex graph" and clearly improves the running time.
Nathann
comment:30 Changed 3 years ago by
 Reviewers set to Dima Pasechnik
 Status changed from needs_review to positive_review
OK. Merges in 6.8, no problem...
comment:31 Changed 3 years ago by
Thanks!
comment:32 Changed 3 years ago by
 Status changed from positive_review to needs_work
sage t long src/sage/graphs/generators/random.py ********************************************************************** File "src/sage/graphs/generators/random.py", line 806, in sage.graphs.generators.random.RandomTriangulation Failed example: for i in range(10): g = graphs.RandomTriangulation(30,embed=True) assert g.is_planar() and g.size() == 3*g.order()6 Exception raised: Traceback (most recent call last): File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 496, in _run self.compile_and_execute(example, compiler, test.globs) File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 858, in compile_and_execute exec(compiled, globs) File "<doctest sage.graphs.generators.random.RandomTriangulation[3]>", line 3, in <module> assert g.is_planar() and g.size() == Integer(3)*g.order()Integer(6) File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/sitepackages/sage/graphs/generic_graph.py", line 4035, in is_planar planar = is_planar(G,kuratowski=kuratowski,set_pos=set_pos,set_embedding=set_embedding) File "sage/graphs/planarity.pyx", line 99, in sage.graphs.planarity.is_planar (/mnt/highperf/buildbot/slave/sage_git/build/src/build/cythonized/sage/graphs/planarity.c:1800) g.relabel(ffrom) File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/sitepackages/sage/graphs/generic_graph.py", line 18006, in relabel new_attr[perm[v]] = value KeyError: 0 ********************************************************************** 1 item had failures: 1 of 5 in sage.graphs.generators.random.RandomTriangulation [71 tests, 1 failure, 5.56 s]
Would have been easy to catch if you had checked the patchbot output...
comment:33 Changed 3 years ago by
 Commit changed from 74d8412010748142918edf6178636bffdb14e60a to 5223cda9f15b2a417b550039b978a07e126f9053
comment:34 followup: ↓ 35 Changed 3 years ago by
 Status changed from needs_work to needs_review
comment:35 in reply to: ↑ 34 Changed 3 years ago by
Replying to ncohen:
Could you explain what happens in the last commit? (E.g.. why did it work earlier?)
comment:36 Changed 3 years ago by
Oh, no problem: there is a difference between the previous output or .vertex_graph
and the new version: in the former, the vertices of the graph were labelled on integers. Now, the vertices of that graph are the 'vertex' objects contained in the polyhedron.
I adapted the code, for it expected the keys of the position dictionary (the vertices) to be integers, instead of the new points.
Actually, it convinced me that the previous version of this function was incorrect in theory. It was assumed that the vertices of the polyhedron appeared in the very same order as how they were first given to the polyhedron constructor. I don't know if it is true in practice, but for sure there is no reason to assume that a priori. So you can see this as a bugfix ;)
Nathann
comment:37 Changed 3 years ago by
Vertices are definitely reordered in Polyhedron
.
comment:38 Changed 3 years ago by
 Status changed from needs_review to positive_review
OK, looks good.
comment:39 Changed 3 years ago by
 Branch changed from public/18860 to 5223cda9f15b2a417b550039b978a07e126f9053
 Resolution set to fixed
 Status changed from positive_review to closed
Branch pushed to git repo; I updated commit sha1. New commits:
trac #18779: polytopes.gosset_3_21 and graphs.GossetGraph
trac #18779: Reviewer's (dead right) comments
capitalise Gosset
trac #18860: Faster Poyhedron.graph()