Revision history [back]
 | 1 | initial version |
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}
| | 2 | No.2 Revision |
I think I have fixed the bug in ticket 29176.
Hopefully 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}
