Rivin's test for circumcribability and inscribability

[mo271]

Let's add the test by Rivin to check whether a polyhedral graph is the graph of a circumscribed polyhedron. (And also for inscribed polyhedra). It's a fairly straightforward LP.

I think it would be good to add this to graphs and also to polyhedra.

Depends on #24634

Component: graph theory

Keywords: polytopes, polyhedra, planar graphs

Author: Moritz Firsching

Branch/Commit: dbadf62

Reviewer: David Coudert

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

[mo271]

Branch: u/moritz/is_circumscribable

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

d3ff9df	typo
8c1b2ef	another typo
40a3d67	'is_circumscribable' and 'is_inscribable'
23b2986	added to list of methods in the beginning

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Commit: 23b2986

[dcoudert]

comment:4

Same remark than for #24634: is this definition/method valid for both graphs and digraphs ?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

ee3e900	move function from 'generic_graph' to 'graph'
c1ad1b6	'is_circumscribable' and 'is_inscribable'
09945f1	added to list of methods in the beginning
9406f03	move circumscribable and inscribable methods to graph

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 23b2986 to 9406f03

[mo271]

comment:7

I moved the methos to where they belong.

[dcoudert]

Reviewer: David Coudert

[dcoudert]

comment:8

A few comments:

• is circumscribed is all of -> is circumscribed if all of

• For error messages, I was told several times to be short, not starting with capital letter and not ending with .. So you should change raise NotImplementedError('%s is not polyhedral. This method only works for polyhedral graphs.' % str(self)) -> raise NotImplementedError('this method only works for polyhedral graphs')

• Do not hesitate to add some empty line to ease de readability

• You create an undirected copy of the graph. Does it means that the definition is valid for both graphs and digraphs ? It's just to know in which file should be the method. More importantly, is it useful ?

• I don't understand why you use vertices_dict to relabel the vertices. Why is it important to use integer vertex labels ?

• The faces are completely determinend -> The faces are completely determined

• set([]) -> set()

• if len(set(cycle)) > 2: Here the cycle is assumed to be simple and the first and last vertices are the same. So if len(cycle) > 3: is be better.

• Why are you forcing the solver to ppl ?

I may have other comments, but so far I don't understand the LP formulation. It is hard to read...

[mo271]

comment:9

Thanks for the comments, David!

Replying to @dcoudert:

• is circumscribed is all of -> is circumscribed if all of

done

• For error messages, I was told several times to be short, not starting with capital letter and not ending with .. So you should change raise NotImplementedError('%s is not polyhedral. This method only works for polyhedral graphs.' % str(self)) -> raise NotImplementedError('this method only works for polyhedral graphs')

ok, done.

• Do not hesitate to add some empty line to ease de readability

I've added a few.

• You create an undirected copy of the graph. Does it means that the definition is valid for both graphs and digraphs ? It's just to know in which file should be the method.

I create the undirected graph for one reason only: I want to iterate over all_simple_cycles of the (undirected) graph. There doesn't seem to be a method for that, but for digraphs, there is such a method; therefore I iterate over all simple cycles of the directed version of the graph and throw away the cycles of type [x,y,x], i.e. of length 3 and also I get every simple cycle of the original graph twice. This is why I collect the inequalities and equations as sets first, to avoid having each constraint twice.

This does not mean that the definition makes sense for directed graph; it is in the right place

More importantly, is it useful ?

I think this is an interesting property not so much for graphs, but for polyhedra. But it only depends on the graph of the polyhedra, which is why I put it in graph.py and not in the polyhedra section. The most useful call will be something like P.graph().is_inscribable() for a polyhedron P; (which could then be shortened to P.is_inscribable()). (Notice that there is already a method is_inscribed, so at least some people are interested in these kinds of questions)

• I don't understand why you use vertices_dict to relabel the vertices. Why is it important to use integer vertex labels ?

The reason is that the variables-dictionary given by M.new_variable() for a MILP M don't seem to be flexible enough to handle complicated keys.

• The faces are completely determinend -> The faces are completely determined

done

• set([]) -> set()

done

• if len(set(cycle)) > 2: Here the cycle is assumed to be simple and the first and last vertices are the same. So if len(cycle) > 3: is be better.

ok, sure.

• Why are you forcing the solver to ppl ?

I want the results to be accurate: Using a numerical solver might lead to incorrect results, especially because I am checking in end if the result is positive. Using exact arithmetic over the rationals and comparing against zero in the end will be slower but hopefully gives the correct result.

I may have other comments, but so far I don't understand the LP formulation. It is hard to read...

I tried to explain it a bit better and also gave a link to a write-up by David Eppstein. (https://www.ics.uci.edu/~eppstein/junkyard/uninscribable/) I hope this helps; I am looking forward to your other comments.. Thanks again!

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

81d87e1	added 'Predefined Classes'
74a45df	'Graph' not 'GenericGraph
f329b20	'is_circumscribable' and 'is_inscribable'
ade7514	added to list of methods in the beginning
9db46ad	move circumscribable and inscribable methods to graph
307b007	improvements suggested by dcoudert

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 9406f03 to 307b007

[dcoudert]

comment:11

At the beginning of method is_circumscribable, you have G = self.to_undirected() but you don't use G afterward, and self is undirected. I suspect that you want to remove loops and multiple edges. If so, you should use G in the following.

Concerning vertices_dict, the variables-dictionary given by M.new_variable() is flexible enough to accept any hashable type. So you don't need this dictionary.

I suggest to rewrite the code as:

    M = MixedIntegerLinearProgram(maximization=True, solver="ppl")
    e = M.new_variable(nonnegative=True)
    c = M.new_variable()
    M.set_min(c[0], -1)
    M.set_max(c[0], 1)
    M.set_objective(c[0])

    for u,v in self.edge_iterator(labels=0):
        if u > v:
            u,v = v,u
        M.set_max(e[u,v], ZZ(1)/ZZ(2))
        M.add_constraint(e[u,v] - c[0], min=0)
        M.add_constraint(e[u,v] + c[0], max=ZZ(1)/ZZ(2))

    from sage.misc.flatten import flatten
    # The faces are completely determined by the graph structure:
    # for polyhedral graph, there is only one way to choose the faces.
    efaces = self.faces()
    vfaces = [flatten(face) for face in efaces]

    # In order to generate all simple cycles of G, we use the "all_simple_cycles"
    # method of directed graphs, generating each cycle twice (in both directions)
    # The two sets below make sure only one direction gives rise to an (in)equality
    D = self.to_directed()
    equality_constraints = set()
    inequality_constraints = set()
    for cycle in D.all_simple_cycles():
        if len(cycle) > 3:
            edges = (tuple(sorted([cycle[i], cycle[i+1]])) for i in range(len(cycle)-1))
            scycle = set(cycle)
            if any(scycle.issubset(_) for _ in vfaces):
                equality_constraints.add(edges)
            else:
                inequality_constraints.add(edges)

    for eq in equality_constraints:
        M.add_constraint(M.sum(e[_] for _ in eq) == 1)
    for ieq in inequality_constraints:
        M.add_constraint(M.sum(e[_] for _ in ieq) - c[0] >= 1)

This code is working (I tried), and I find it easier to read.

Do we agree that

• the equality constraints are exactly the faces
• any simple cycle that is not a face induces an inequality constraint

If so, why using issubset ?

Otherwise, we could do something like:

    vfaces = set(Set(flatten(face)) for face in efaces)

    # In order to generate all simple cycles of G, we use the "all_simple_cycles"
    # method of directed graphs, generating each cycle twice (in both directions)
    # The two sets below make sure only one direction gives rise to an (in)equality
    D = self.to_directed()
    equality_constraints = set()
    inequality_constraints = set()
    for cycle in D.all_simple_cycles():
        if len(cycle) > 3:
            edges = (tuple(sorted([cycle[i], cycle[i+1]])) for i in range(len(cycle)-1))
            if not Set(cycle) in vfaces:
                equality_constraints.add(edges)
            else:
                inequality_constraints.add(edges)

and we can even save set equality_constraints since it is duplicate of efaces.

[mo271]

comment:12

Replying to @dcoudert:

At the beginning of method is_circumscribable, you have G = self.to_undirected() but you don't use G afterward, and self is undirected. I suspect that you want to remove loops and multiple edges. If so, you should use G in the following.

No, I simply forgot to delete that file when moving the method from generic_graph to graph. Graph with loops and multiple edges are caught by the is_polyhedral test before.

Concerning vertices_dict, the variables-dictionary given by M.new_variable() is flexible enough to accept any hashable type. So you don't need this dictionary.

Great! I don't know anymore why I was thinking these dictionary weren't flexible enough.

I suggest to rewrite the code as:

I adopted all your changes, they do make the code more readable! Not all doctest worked with what you had, but after added set to the following line

vfaces = [set(flatten(face)) for face in efaces]

and modifying a few things, everything now works again. Also, your version generated all the inequalities for non-facial cycles twice; I added a frozenset to remedy this. (I had tuple(sorted(edges)) before.)

Do we agree that

• the equality constraints are exactly the faces
• any simple cycle that is not a face induces an inequality constraint

Yes!

If so, why using issubset ?

Not using it anymore..

Otherwise, we could do something like: and we can even save set equality_constraints since it is duplicate of efaces.

done, the equalities are now treated directly from efaces before

Thank you: this certainly improved the code!

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

3b7916b

cleaner rewrite suggested by dcoudert

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 307b007 to 3b7916b

[dcoudert]

comment:14

Please add an empty line after M.add_constraint(M.sum(e[tuple(sorted(_))] for _ in edges) == 1).

vfaces should be a set, not a list. Testing if an item is in a list takes linear time in the size of the list, and much less in a set. See the example below.

sage: L = list(range(1000))
sage: shuffle(L)
sage: S = set(L)
sage: %time all(i in L for i in L)
CPU times: user 7.08 ms, sys: 423 µs, total: 7.5 ms
Wall time: 7.19 ms
True
sage: %time all(i in S for i in L)
CPU times: user 304 µs, sys: 182 µs, total: 486 µs
Wall time: 329 µs
True

So I propose to use

vfaces = set(frozenset(flatten(face)) for face in efaces))

and

            if len(cycle) > 3:
                scycle = frozenset(cycle)
                if scycle not in vfaces:
                    edges = (tuple(sorted([cycle[i], cycle[i+1]])) for i in range(len(cycle)-1))
                    inequality_constraints.add(frozenset(edges))

We need to use frozenset (or Set) since this type is hashable, and so we can have a set of frozensets.

Observe that since cycle is a list, we can do scycle = frozenset(cycle). No need to flatten the cycle first.

In the description of the method, you could add after by solving a linear program

that assigns weights between 0 and 1/2 on each edge of the polyhedron, so that the weights on any face add to exactly one and the weights on any non-facial cycle add to more than one. If this can be done, the polyhedron can be circumscribed.

[mo271]

comment:15

Replying to @dcoudert:

Please add an empty line after M.add_constraint(M.sum(e[tuple(sorted(_))] for _ in edges) == 1).

done

vfaces should be a set, not a list.

sure, that makes sense!

Observe that since cycle is a list, we can do scycle = frozenset(cycle). No need to flatten the cycle first.

Nice! Now I avoid flattening lists also when generating the list
vfaces = set(frozenset([_[0] for _ in face]) for face in efaces)' Usingflatteninstead of[[0] for in face]` leads to errors when the vertices of the Graph are lists themselves (as often happens when looking at the planar dual). Now it works.

In the description of the method, you could add after by solving a linear program

... done!

Thank you for the further suggestions!

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

c2549b5

make vfaces set of frozensets

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 3b7916b to c2549b5

[dcoudert]

comment:17

This is much better, but I have further comments (sorry):

• please ensure that the comments are written in 80 columns mode. The line A polyhedron is circumscribed if all of its facets are tangent to a sphere. exceeds the 80th column.

• In the examples, you could add

sage: O = graphs.OctahedralGraph()
sage: 0.is_circumscribable()
True

• in the example of is_inscribable:

sage: C = graphs.CubeGraph(3)
sage: v = next(C.vertex_iterator())
sage: triangle = [_ + v for _ in C.neighbors(v)]

Honestly, I don't understand the transformation, but it's not important.

May be it would be more informative to have:

sage: C = graphs.CubeGraph(3)
sage: C.is_inscribable()
True
sage: C.planar_dual().is_inscribable()
True
sage: H = graphs.HerschelGraph()
sage: H.is_inscribable()               # long time (> 1 sec)
False
sage: H.planar_dual().is_inscribable() # long time (> 1 sec)
True

[mo271]

comment:18

Replying to @dcoudert:

This is much better, but I have further comments (sorry):

That's great, I enjoy your comments very much!

• please ensure that the comments are written in 80 columns mode.

done.

• In the examples, you could add ...

I added the octahedron.

Honestly, I don't understand the transformation, but it's not important.

I added an explanation what it does (it is cutting off a vertex from the cube)

May be it would be more informative to have...

I added some more examples that I find pretty instructive. They include the smallest examples of non-(in/circum)scribable graphs. (Another would be dual to the truncated tetrahedron: the stacked tetrahedron, a.k.a. Triakis tetradhedron (https://en.wikipedia.org/wiki/Triakis_tetrahedron), which is not inscribed. But I don't know how to get the graph easily, except polytopes.truncated_tetrahedron().graph().planar_dual() or polytopes.truncated_tetrahedron().polar().graph() of course..)

Thank you for looking at the code so closely!

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c2549b5 to dacc52a

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

dacc52a

more examples and some pep8 improvements

[dcoudert]

comment:20

Perfect.

Now let's check the references.

• .. [HRS1993] ] Craig D. Hodgson, Igor Rivin and Warren D. Smith, \ -> .. [HRS1993] \C. D. Hodgson, I. Rivin and W. D. Smith.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from dacc52a to dbadf62

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

dbadf62

fixed reference

[mo271]

comment:22

fixed, thanks!

---

New commits:

dbadf62

fixed reference

[dcoudert]

comment:23

No more comments.

[vbraun]

Changed branch from u/moritz/is_circumscribable to dbadf62