Ask Your Question
3

VoronoiDiagram returns empty regions

asked 2020-02-02 18:09:33 +0100

muradtuk gravatar image

updated 2020-02-02 21:33:17 +0100

vdelecroix gravatar image

Dear all,

Given the set of 2 points

P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]

I would like to compute the voronoi diagram of $P$. When doing so i saw that the regions of the voronoi are:

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

In my code, I do need that those regions are not empty! Why is this happening and how can i fix it?

Please advise and thanks in advance.

P.s. Note that when I make the points a little I get that the voronoi is not empty.

Full Code:

P = [[-2687.19, -2088.53], [-2686.81, -2084.19]]
V = VoronoiDiagram(P)
R = V.regions()
edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
2

answered 2020-03-19 16:46:15 +0100

Jonathan Kliem gravatar image

updated 2020-03-19 16:46:33 +0100

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}
edit flag offensive delete link more
1

answered 2020-02-10 13:47:43 +0100

jipilab gravatar image

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

Thanks for reporting this!

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2020-02-02 17:43:31 +0100

Seen: 678 times

Last updated: Mar 19 '20