VoronoiDiagram returns empty regions

edit

This looks like a bug. I have reported it on trac: ticket 29176

Thanks for reporting this!

I think I have fixed the bug in ticket 29176.

Hopefully, it will be available from Sage version 9.1+.

In the meantime you can use the following workaround.

from sage.all import *
from sage.geometry.polyhedron.parent import Polyhedra_RDF_cdd
from sage.geometry.polyhedron.backend_cdd import Polyhedron_RDF_cdd
class Polyhedron_RDF_cdd_workaround(Polyhedron_RDF_cdd):
    def _init_from_cdd_output(self, cddout):
        Polyhedron_RDF_cdd._init_from_cdd_output(self, cddout)
        cddout = cddout.splitlines()

        def parse_indices(count, cdd_indices, cdd_indices_to_sage_indices=None):
            cdd_indices = [int(x) for x in cdd_indices]
            if cdd_indices_to_sage_indices is None:
                cdd_indices_to_sage_indices = {i: i-1 for i in cdd_indices}
            if count < 0:
                assert cdd_indices_to_sage_indices is not None, "Did not expect negative counts here"
                count = -count
                cdd_indices = list(set(cdd_indices_to_sage_indices.keys()) - set(cdd_indices))
                assert count in [len(cdd_indices), len(cdd_indices) - 1]
            assert count == len(cdd_indices)
            return [cdd_indices_to_sage_indices[i] for i in cdd_indices if cdd_indices_to_sage_indices[i] is not None]

        def parse_incidence_matrix(intro, data):
            if 'incidence_matrix' in self.__dict__:
                raise NotImplementedError("can not replace internal representation as this breaks caching")
            N = len(self._Hrepresentation)
            M = len(self._Vrepresentation)
            ret = matrix(ZZ, M, N, 0)
            data.pop(0)
            data.reverse()
            for adjacencies in data:
                assert adjacencies[2] == ':', "Not a line of adjacency data"
                cdd_vertex = int(adjacencies[0])
                count = int(adjacencies[1])
                if cdd_vertex not in self._cdd_V_to_sage_V:
                    self._cdd_V_to_sage_V[cdd_vertex] = None
                v = self._cdd_V_to_sage_V[cdd_vertex]
                if v is None:
                    continue
                for w in parse_indices(count, adjacencies[3:], self._cdd_H_to_sage_H):
                    if w is None:
                        continue
                    ret[v, w] = 1 
            ret.set_immutable()
            self.incidence_matrix.set_cache(ret)

        Polyhedron_RDF_cdd._parse_block(cddout, 'Vertex incidence', parse_incidence_matrix)

class Polyhedra_RDF_cdd_workaround(Polyhedra_RDF_cdd):
    Element = Polyhedron_RDF_cdd_workaround

class VoronoiDiagram_workaround(VoronoiDiagram):
    def __init__(self, points):
        VoronoiDiagram.__init__(self, points)
        self._P = {}
        if self._n > 0:
            e = [([sum(vector(i)[k] ** 2
                   for k in range(self._d))] +
                  [(-2) * vector(i)[l] for l in range(self._d)] + [1])
                 for i in self._points]
            e = [[self._base_ring(i) for i in k] for k in e]
            parent = Polyhedra_RDF_cdd_workaround(RDF, self._d + 1, 'cdd')
            p = parent(None, [e, []])
        for i in range(self._n):
            parent = Polyhedra_RDF_cdd_workaround(RDF, self._d, 'cdd')
            equ = p.Hrepresentation(i)
            pvert = [[u[k] for k in range(self._d)] for u in equ.incident()
                     if u.is_vertex()]
            prays = [[u[k] for k in range(self._d)] for u in equ.incident()
                     if u.is_ray()]
            pline = [[u[k] for k in range(self._d)] for u in equ.incident()
                     if u.is_line()]
            (self._P)[self._points[i]] = parent([pvert, prays, pline], None)

Then

P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]
V = VoronoiDiagram_workaround(P)
R = V.regions()
R

will produce correctly

{P(-2687.19000000000, -2088.53000000000): A 2-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex, 1 ray, 1 line,
 P(-2686.81000000000, -2084.19000000000): A 2-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex, 1 ray, 1 line}