Merge pull request #36 from mcmtroffaes/expose_incidence

Expose adjacency and incidence.

CHANGELOG.rst

@@ -1,6 +1,9 @@
 Version 2.1.1 (in development)
 ------------------------------
 
+* exposed adjacency and incidence (see issues #33, #34, and #36,
+  contributed by bobmyhill)
+
 Version 2.1.0 (15 October 2018)
 -------------------------------
 

cdd.pyx

@@ -26,7 +26,7 @@ cimport libc.stdlib
 from fractions import Fraction
 import numbers
 
-__version__ = "2.1.1a0"
+__version__ = "2.1.1a1"
 
 # also need time_t
 cdef extern from "time.h":
@@ -150,6 +150,36 @@ cdef _set_set(set_type set_, pset):
         else:
             set_delelem(set_, elem)
 
+cdef _get_dd_setfam(dd_SetFamilyPtr setfam):
+    """Create tuple of Python frozensets from dd_SetFamilyPtr, and
+    free the pointer. The indexing of the sets start at 0, unlike the
+    string output from cddlib, which starts at 1.
+    """
+    cdef long elem
+    if setfam == NULL:
+        raise ValueError("failed to get set family")
+    result = tuple(frozenset([elem - 1
+                              for elem from 1 <= elem <= setfam.setsize
+                              if set_member(elem, setfam.set[i])])
+                   for i in range(setfam.famsize))
+    dd_FreeSetFamily(setfam)
+    return result
+
+cdef _get_ddf_setfam(ddf_SetFamilyPtr setfam):
+    """Create tuple of Python frozensets from ddf_SetFamilyPtr, and
+    free the pointer. The indexing of the sets start at 0, unlike the
+    string output from cddlib, which starts at 1.
+    """
+    cdef long elem
+    if setfam == NULL:
+        raise ValueError("failed to get set family")
+    result = tuple(frozenset([elem - 1
+                              for elem from 1 <= elem <= setfam.setsize
+                              if set_member(elem, setfam.set[i])])
+                   for i in range(setfam.famsize))
+    ddf_FreeSetFamily(setfam)
+    return result
+
 cdef _raise_error(dd_ErrorType error, msg):
     """Convert error into string and raise it."""
     cdef libc.stdio.FILE *pfile
@@ -844,6 +874,30 @@ cdef class Polyhedron(NumberTypeable):
         else:
             return _make_ddf_matrix(ddf_CopyGenerators(self.ddf_poly))
 
+    def get_adjacency(self):
+        if self.dd_poly:
+            return _get_dd_setfam(dd_CopyAdjacency(self.dd_poly))
+        else:
+            return _get_ddf_setfam(ddf_CopyAdjacency(self.ddf_poly))
+
+    def get_input_adjacency(self):
+        if self.dd_poly:
+            return _get_dd_setfam(dd_CopyInputAdjacency(self.dd_poly))
+        else:
+            return _get_ddf_setfam(ddf_CopyInputAdjacency(self.ddf_poly))
+
+    def get_incidence(self):
+        if self.dd_poly:
+            return _get_dd_setfam(dd_CopyIncidence(self.dd_poly))
+        else:
+            return _get_ddf_setfam(ddf_CopyIncidence(self.ddf_poly))
+
+    def get_input_incidence(self):
+        if self.dd_poly:
+            return _get_dd_setfam(dd_CopyInputIncidence(self.dd_poly))
+        else:
+            return _get_ddf_setfam(ddf_CopyInputIncidence(self.ddf_poly))
+
 # module initialization code comes here
 # initialize module constants
 dd_set_global_constants()

docs/polyhedron.rst

@@ -41,6 +41,46 @@ Methods and Attributes
         column vector with `n` ones followed by `s` zeroes, and `V` is the
         stacked matrix of `n` vertex row vectors on top of `s` ray row vectors.
 
+.. method:: Polyhedron.get_adjacency()
+
+        Get the adjacencies.
+
+        :returns: Adjacency list.
+        :rtype: :class:`tuple`
+
+        H-representation: For each vertex, list adjacent vertices.
+        V-representation: For each face, list adjacent faces.
+
+.. method:: Polyhedron.get_input_adjacency()
+
+        Get the input adjacencies.
+
+        :returns: Input adjacency list.
+        :rtype: :class:`tuple`
+
+        H-representation: For each face, list adjacent faces.
+        V-representation: For each vertex, list adjacent vertices.
+
+.. method:: Polyhedron.get_incidence()
+
+        Get the incidences.
+
+        :returns: Incidence list.
+        :rtype: :class:`tuple`
+
+        H-representation: For each vertex, list adjacent faces.
+        V-representation: For each face, list adjacent vertices.
+
+.. method:: Polyhedron.get_input_incidence()
+
+        Get the input incidences.
+
+        :returns: Input incidence list.
+        :rtype: :class:`tuple`
+
+        H-representation: For each face, list adjacent vertices.
+        V-representation: For each vertex, list adjacent faces.
+
 .. attribute:: Polyhedron.rep_type
 
         Representation (see :class:`~cdd.RepType`).
@@ -50,8 +90,8 @@ Methods and Attributes
     The H-representation and/or V-representation are not guaranteed to
     be minimal, that is, they can still contain redundancy.
 
-Example
--------
+Examples
+--------
 
 This is the sampleh1.ine example that comes with cddlib.
 
@@ -79,3 +119,108 @@ begin
 end
 >>> print(list(ext.lin_set)) # note: first row is 0, so fourth row is 3
 [3]
+
+
+The following example illustrates how to get adjacencies and incidences.
+
+>>> import cdd
+>>> # We start with the H-representation for a square
+>>> # 0 <= 1 + x1 (face 0)
+>>> # 0 <= 1 + x2 (face 1)
+>>> # 0 <= 1 - x1 (face 2)
+>>> # 0 <= 1 - x2 (face 3)
+>>> mat = cdd.Matrix([[1, 1, 0], [1, 0, 1], [1, -1, 0], [1, 0, -1]])
+>>> mat.rep_type = cdd.RepType.INEQUALITY
+>>> poly = cdd.Polyhedron(mat)
+>>> # The V-representation can be printed in the usual way:
+>>> gen = poly.get_generators()
+>>> print(gen)
+V-representation
+begin
+ 4 3 rational
+ 1 1 -1
+ 1 1 1
+ 1 -1 1
+ 1 -1 -1
+end
+>>> # graphical depiction of vertices and faces:
+>>> #
+>>> #   2---(3)---1
+>>> #   |         |
+>>> #   |         |
+>>> #  (0)       (2)
+>>> #   |         |
+>>> #   |         |
+>>> #   3---(1)---0
+>>> #
+>>> # vertex 0 is adjacent to vertices 1 and 3
+>>> # vertex 1 is adjacent to vertices 0 and 2
+>>> # vertex 2 is adjacent to vertices 1 and 3
+>>> # vertex 3 is adjacent to vertices 0 and 2
+>>> print([list(x) for x in poly.get_adjacency()])
+[[1, 3], [0, 2], [1, 3], [0, 2]]
+>>> # vertex 0 is the intersection of faces (1) and (2)
+>>> # vertex 1 is the intersection of faces (2) and (3)
+>>> # vertex 2 is the intersection of faces (0) and (3)
+>>> # vertex 3 is the intersection of faces (0) and (1)
+>>> print([list(x) for x in poly.get_incidence()])
+[[1, 2], [2, 3], [0, 3], [0, 1]]
+>>> # face (0) is adjacent to faces (1) and (3)
+>>> # face (1) is adjacent to faces (0) and (2)
+>>> # face (2) is adjacent to faces (1) and (3)
+>>> # face (3) is adjacent to faces (0) and (2)
+>>> print([list(x) for x in poly.get_input_adjacency()])
+[[1, 3], [0, 2], [1, 3], [0, 2], []]
+>>> # face (0) intersects with vertices 2 and 3
+>>> # face (1) intersects with vertices 0 and 3
+>>> # face (2) intersects with vertices 0 and 1
+>>> # face (3) intersects with vertices 1 and 2
+>>> print([list(x) for x in poly.get_input_incidence()])
+[[2, 3], [0, 3], [0, 1], [1, 2], []]
+>>> # add a vertex, and construct new polyhedron
+>>> gen.extend([[1, 0, 2]])
+>>> vpoly = cdd.Polyhedron(gen)
+>>> print(vpoly.get_inequalities())
+H-representation
+begin
+ 5 3 rational
+ 1 0 1
+ 2 1 -1
+ 1 1 0
+ 2 -1 -1
+ 1 -1 0
+end
+>>> # so now we have:
+>>> # 0 <= 1 + x2
+>>> # 0 <= 2 + x1 - x2
+>>> # 0 <= 1 + x1
+>>> # 0 <= 2 - x1 - x2
+>>> # 0 <= 1 - x1
+>>> #
+>>> # graphical depiction of vertices and faces:
+>>> #
+>>> #        4
+>>> #       / \
+>>> #      /   \
+>>> #    (1)   (3)
+>>> #    /       \
+>>> #   2         1
+>>> #   |         |
+>>> #   |         |
+>>> #  (2)       (4)
+>>> #   |         |
+>>> #   |         |
+>>> #   3---(0)---0
+>>> #
+>>> # for each face, list adjacent faces
+>>> print([list(x) for x in vpoly.get_adjacency()])
+[[2, 4], [2, 3], [0, 1], [1, 4], [0, 3]]
+>>> # for each face, list adjacent vertices
+>>> print([list(x) for x in vpoly.get_incidence()])
+[[0, 3], [2, 4], [2, 3], [1, 4], [0, 1]]
+>>> # for each vertex, list adjacent vertices
+>>> print([list(x) for x in vpoly.get_input_adjacency()])
+[[1, 3], [0, 4], [3, 4], [0, 2], [1, 2]]
+>>> # for each vertex, list adjacent faces
+>>> print([list(x) for x in vpoly.get_input_incidence()])
+[[0, 4], [3, 4], [1, 2], [0, 2], [1, 3]]

test/test_adjacency_list.py

@@ -0,0 +1,71 @@
+import cdd
+import nose
+
+
+def _make_vertex_adjacency_list(number_type):
+
+    # The following lines test that poly.get_adjacency_list()
+    # returns the correct adjacencies.
+
+    # We start with the H-representation for a cube
+    mat = cdd.Matrix([[1, 1, 0 ,0],
+                      [1, 0, 1, 0],
+                      [1, 0, 0, 1],
+                      [1, -1, 0, 0],
+                      [1, 0, -1, 0],
+                      [1, 0, 0, -1]],
+                     number_type=number_type)
+    mat.rep_type = cdd.RepType.INEQUALITY
+    poly = cdd.Polyhedron(mat)
+    adjacency_list = poly.get_adjacency()
+
+    # Family size should equal the number of vertices of the cube (8)
+    nose.tools.assert_equal(len(adjacency_list), 8)
+
+    # All the vertices of the cube should be connected by three other vertices
+    nose.tools.assert_equal([len(adj) for adj in adjacency_list], [3]*8)
+
+    # The vertices must be numbered consistently
+    # The first vertex is adjacent to the second, fourth and eighth
+    # (note the conversion to a pythonic numbering system)
+    adjacencies = [[1, 3, 7],
+                   [0, 2, 6],
+                   [1, 3, 4],
+                   [0, 2, 5],
+                   [2, 5, 6],
+                   [3, 4, 7],
+                   [1, 4, 7],
+                   [0, 5, 6]]
+    for i in range(8):
+        nose.tools.assert_equal(list(adjacency_list[i]), adjacencies[i])
+
+
+def _make_facet_adjacency_list(number_type):
+    # This matrix is the same as in vtest_vo.ine
+    mat = cdd.Matrix([[0, 0, 0, 1],
+                      [5, -4, -2, 1],
+                      [5, -2, -4, 1],
+                      [16, -8, 0, 1],
+                      [16, 0, -8, 1],
+                      [32, -8, -8, 1]], number_type=number_type)
+
+    mat.rep_type = cdd.RepType.INEQUALITY
+    poly = cdd.Polyhedron(mat)
+
+    adjacencies = [[1, 2, 3, 4, 6],
+                   [0, 2, 3, 5],
+                   [0, 1, 4, 5],
+                   [0, 1, 5, 6],
+                   [0, 2, 5, 6],
+                   [1, 2, 3, 4, 6],
+                   [0, 3, 4, 5]]
+
+    adjacency_list = poly.get_input_adjacency()
+    for i in range(7):
+        nose.tools.assert_equal(list(adjacency_list[i]), adjacencies[i])
+
+def test_all():
+    _make_vertex_adjacency_list("fraction")
+    _make_vertex_adjacency_list("float")
+    _make_facet_adjacency_list("fraction")
+    _make_facet_adjacency_list("float")