1import numpy as np
2import itertools
3from numpy.testing import (assert_equal,
4                           assert_almost_equal,
5                           assert_array_equal,
6                           assert_array_almost_equal,
7                           suppress_warnings)
8import pytest
9from pytest import raises as assert_raises
10from pytest import warns as assert_warns
11from scipy.spatial import SphericalVoronoi, distance
12from scipy.optimize import linear_sum_assignment
13from scipy.constants import golden as phi
14from scipy.special import gamma
15
16
17TOL = 1E-10
18
19
20def _generate_tetrahedron():
21    return np.array([[1, 1, 1], [1, -1, -1], [-1, 1, -1], [-1, -1, 1]])
22
23
24def _generate_cube():
25    return np.array(list(itertools.product([-1, 1.], repeat=3)))
26
27
28def _generate_octahedron():
29    return np.array([[-1, 0, 0], [+1, 0, 0], [0, -1, 0],
30                     [0, +1, 0], [0, 0, -1], [0, 0, +1]])
31
32
33def _generate_dodecahedron():
34
35    x1 = _generate_cube()
36    x2 = np.array([[0, -phi, -1 / phi],
37                   [0, -phi, +1 / phi],
38                   [0, +phi, -1 / phi],
39                   [0, +phi, +1 / phi]])
40    x3 = np.array([[-1 / phi, 0, -phi],
41                   [+1 / phi, 0, -phi],
42                   [-1 / phi, 0, +phi],
43                   [+1 / phi, 0, +phi]])
44    x4 = np.array([[-phi, -1 / phi, 0],
45                   [-phi, +1 / phi, 0],
46                   [+phi, -1 / phi, 0],
47                   [+phi, +1 / phi, 0]])
48    return np.concatenate((x1, x2, x3, x4))
49
50
51def _generate_icosahedron():
52    x = np.array([[0, -1, -phi],
53                  [0, -1, +phi],
54                  [0, +1, -phi],
55                  [0, +1, +phi]])
56    return np.concatenate([np.roll(x, i, axis=1) for i in range(3)])
57
58
59def _generate_polytope(name):
60    polygons = ["triangle", "square", "pentagon", "hexagon", "heptagon",
61                "octagon", "nonagon", "decagon", "undecagon", "dodecagon"]
62    polyhedra = ["tetrahedron", "cube", "octahedron", "dodecahedron",
63                 "icosahedron"]
64    if name not in polygons and name not in polyhedra:
65        raise ValueError("unrecognized polytope")
66
67    if name in polygons:
68        n = polygons.index(name) + 3
69        thetas = np.linspace(0, 2 * np.pi, n, endpoint=False)
70        p = np.vstack([np.cos(thetas), np.sin(thetas)]).T
71    elif name == "tetrahedron":
72        p = _generate_tetrahedron()
73    elif name == "cube":
74        p = _generate_cube()
75    elif name == "octahedron":
76        p = _generate_octahedron()
77    elif name == "dodecahedron":
78        p = _generate_dodecahedron()
79    elif name == "icosahedron":
80        p = _generate_icosahedron()
81
82    return p / np.linalg.norm(p, axis=1, keepdims=True)
83
84
85def _hypersphere_area(dim, radius):
86    # https://en.wikipedia.org/wiki/N-sphere#Closed_forms
87    return 2 * np.pi**(dim / 2) / gamma(dim / 2) * radius**(dim - 1)
88
89
90def _sample_sphere(n, dim, seed=None):
91    # Sample points uniformly at random from the hypersphere
92    rng = np.random.RandomState(seed=seed)
93    points = rng.randn(n, dim)
94    points /= np.linalg.norm(points, axis=1, keepdims=True)
95    return points
96
97
98class TestSphericalVoronoi:
99
100    def setup_method(self):
101        self.points = np.array([
102            [-0.78928481, -0.16341094, 0.59188373],
103            [-0.66839141, 0.73309634, 0.12578818],
104            [0.32535778, -0.92476944, -0.19734181],
105            [-0.90177102, -0.03785291, -0.43055335],
106            [0.71781344, 0.68428936, 0.12842096],
107            [-0.96064876, 0.23492353, -0.14820556],
108            [0.73181537, -0.22025898, -0.6449281],
109            [0.79979205, 0.54555747, 0.25039913]]
110        )
111
112    def test_constructor(self):
113        center = np.array([1, 2, 3])
114        radius = 2
115        s1 = SphericalVoronoi(self.points)
116        # user input checks in SphericalVoronoi now require
117        # the radius / center to match the generators so adjust
118        # accordingly here
119        s2 = SphericalVoronoi(self.points * radius, radius)
120        s3 = SphericalVoronoi(self.points + center, center=center)
121        s4 = SphericalVoronoi(self.points * radius + center, radius, center)
122        assert_array_equal(s1.center, np.array([0, 0, 0]))
123        assert_equal(s1.radius, 1)
124        assert_array_equal(s2.center, np.array([0, 0, 0]))
125        assert_equal(s2.radius, 2)
126        assert_array_equal(s3.center, center)
127        assert_equal(s3.radius, 1)
128        assert_array_equal(s4.center, center)
129        assert_equal(s4.radius, radius)
130
131    def test_vertices_regions_translation_invariance(self):
132        sv_origin = SphericalVoronoi(self.points)
133        center = np.array([1, 1, 1])
134        sv_translated = SphericalVoronoi(self.points + center, center=center)
135        assert_equal(sv_origin.regions, sv_translated.regions)
136        assert_array_almost_equal(sv_origin.vertices + center,
137                                  sv_translated.vertices)
138
139    def test_vertices_regions_scaling_invariance(self):
140        sv_unit = SphericalVoronoi(self.points)
141        sv_scaled = SphericalVoronoi(self.points * 2, 2)
142        assert_equal(sv_unit.regions, sv_scaled.regions)
143        assert_array_almost_equal(sv_unit.vertices * 2,
144                                  sv_scaled.vertices)
145
146    def test_old_radius_api(self):
147        sv_unit = SphericalVoronoi(self.points, radius=1)
148        with suppress_warnings() as sup:
149            sup.filter(DeprecationWarning, "`radius` is `None`")
150            sv = SphericalVoronoi(self.points, None)
151            assert_array_almost_equal(sv_unit.vertices, sv.vertices)
152
153    def test_old_radius_api_warning(self):
154        with assert_warns(DeprecationWarning):
155            SphericalVoronoi(self.points, None)
156
157    def test_sort_vertices_of_regions(self):
158        sv = SphericalVoronoi(self.points)
159        unsorted_regions = sv.regions
160        sv.sort_vertices_of_regions()
161        assert_equal(sorted(sv.regions), sorted(unsorted_regions))
162
163    def test_sort_vertices_of_regions_flattened(self):
164        expected = sorted([[0, 6, 5, 2, 3], [2, 3, 10, 11, 8, 7], [0, 6, 4, 1],
165                           [4, 8, 7, 5, 6], [9, 11, 10], [2, 7, 5],
166                           [1, 4, 8, 11, 9], [0, 3, 10, 9, 1]])
167        expected = list(itertools.chain(*sorted(expected)))  # type: ignore
168        sv = SphericalVoronoi(self.points)
169        sv.sort_vertices_of_regions()
170        actual = list(itertools.chain(*sorted(sv.regions)))
171        assert_array_equal(actual, expected)
172
173    def test_sort_vertices_of_regions_dimensionality(self):
174        points = np.array([[1, 0, 0, 0],
175                           [0, 1, 0, 0],
176                           [0, 0, 1, 0],
177                           [0, 0, 0, 1],
178                           [0.5, 0.5, 0.5, 0.5]])
179        with pytest.raises(TypeError, match="three-dimensional"):
180            sv = SphericalVoronoi(points)
181            sv.sort_vertices_of_regions()
182
183    def test_num_vertices(self):
184        # for any n >= 3, a spherical Voronoi diagram has 2n - 4
185        # vertices; this is a direct consequence of Euler's formula
186        # as explained by Dinis and Mamede (2010) Proceedings of the
187        # 2010 International Symposium on Voronoi Diagrams in Science
188        # and Engineering
189        sv = SphericalVoronoi(self.points)
190        expected = self.points.shape[0] * 2 - 4
191        actual = sv.vertices.shape[0]
192        assert_equal(actual, expected)
193
194    def test_voronoi_circles(self):
195        sv = SphericalVoronoi(self.points)
196        for vertex in sv.vertices:
197            distances = distance.cdist(sv.points, np.array([vertex]))
198            closest = np.array(sorted(distances)[0:3])
199            assert_almost_equal(closest[0], closest[1], 7, str(vertex))
200            assert_almost_equal(closest[0], closest[2], 7, str(vertex))
201
202    def test_duplicate_point_handling(self):
203        # an exception should be raised for degenerate generators
204        # related to Issue# 7046
205        self.degenerate = np.concatenate((self.points, self.points))
206        with assert_raises(ValueError):
207            SphericalVoronoi(self.degenerate)
208
209    def test_incorrect_radius_handling(self):
210        # an exception should be raised if the radius provided
211        # cannot possibly match the input generators
212        with assert_raises(ValueError):
213            SphericalVoronoi(self.points, radius=0.98)
214
215    def test_incorrect_center_handling(self):
216        # an exception should be raised if the center provided
217        # cannot possibly match the input generators
218        with assert_raises(ValueError):
219            SphericalVoronoi(self.points, center=[0.1, 0, 0])
220
221    @pytest.mark.parametrize("dim", range(2, 6))
222    @pytest.mark.parametrize("shift", [False, True])
223    def test_single_hemisphere_handling(self, dim, shift):
224        n = 10
225        points = _sample_sphere(n, dim, seed=0)
226        points[:, 0] = np.abs(points[:, 0])
227        center = (np.arange(dim) + 1) * shift
228        sv = SphericalVoronoi(points + center, center=center)
229        dots = np.einsum('ij,ij->i', sv.vertices - center,
230                                     sv.points[sv._simplices[:, 0]] - center)
231        circumradii = np.arccos(np.clip(dots, -1, 1))
232        assert np.max(circumradii) > np.pi / 2
233
234    @pytest.mark.parametrize("n", [1, 2, 10])
235    @pytest.mark.parametrize("dim", range(2, 6))
236    @pytest.mark.parametrize("shift", [False, True])
237    def test_rank_deficient(self, n, dim, shift):
238        center = (np.arange(dim) + 1) * shift
239        points = _sample_sphere(n, dim - 1, seed=0)
240        points = np.hstack([points, np.zeros((n, 1))])
241        with pytest.raises(ValueError, match="Rank of input points"):
242            SphericalVoronoi(points + center, center=center)
243
244    @pytest.mark.parametrize("dim", range(2, 6))
245    def test_higher_dimensions(self, dim):
246        n = 100
247        points = _sample_sphere(n, dim, seed=0)
248        sv = SphericalVoronoi(points)
249        assert sv.vertices.shape[1] == dim
250        assert len(sv.regions) == n
251
252        # verify Euler characteristic
253        cell_counts = []
254        simplices = np.sort(sv._simplices)
255        for i in range(1, dim + 1):
256            cells = []
257            for indices in itertools.combinations(range(dim), i):
258                cells.append(simplices[:, list(indices)])
259            cells = np.unique(np.concatenate(cells), axis=0)
260            cell_counts.append(len(cells))
261        expected_euler = 1 + (-1)**(dim-1)
262        actual_euler = sum([(-1)**i * e for i, e in enumerate(cell_counts)])
263        assert expected_euler == actual_euler
264
265    @pytest.mark.parametrize("dim", range(2, 6))
266    def test_cross_polytope_regions(self, dim):
267        # The hypercube is the dual of the cross-polytope, so the voronoi
268        # vertices of the cross-polytope lie on the points of the hypercube.
269
270        # generate points of the cross-polytope
271        points = np.concatenate((-np.eye(dim), np.eye(dim)))
272        sv = SphericalVoronoi(points)
273        assert all([len(e) == 2**(dim - 1) for e in sv.regions])
274
275        # generate points of the hypercube
276        expected = np.vstack(list(itertools.product([-1, 1], repeat=dim)))
277        expected = expected.astype(np.float64) / np.sqrt(dim)
278
279        # test that Voronoi vertices are correctly placed
280        dist = distance.cdist(sv.vertices, expected)
281        res = linear_sum_assignment(dist)
282        assert dist[res].sum() < TOL
283
284    @pytest.mark.parametrize("dim", range(2, 6))
285    def test_hypercube_regions(self, dim):
286        # The cross-polytope is the dual of the hypercube, so the voronoi
287        # vertices of the hypercube lie on the points of the cross-polytope.
288
289        # generate points of the hypercube
290        points = np.vstack(list(itertools.product([-1, 1], repeat=dim)))
291        points = points.astype(np.float64) / np.sqrt(dim)
292        sv = SphericalVoronoi(points)
293
294        # generate points of the cross-polytope
295        expected = np.concatenate((-np.eye(dim), np.eye(dim)))
296
297        # test that Voronoi vertices are correctly placed
298        dist = distance.cdist(sv.vertices, expected)
299        res = linear_sum_assignment(dist)
300        assert dist[res].sum() < TOL
301
302    @pytest.mark.parametrize("n", [10, 500])
303    @pytest.mark.parametrize("dim", [2, 3])
304    @pytest.mark.parametrize("radius", [0.5, 1, 2])
305    @pytest.mark.parametrize("shift", [False, True])
306    @pytest.mark.parametrize("single_hemisphere", [False, True])
307    def test_area_reconstitution(self, n, dim, radius, shift,
308                                 single_hemisphere):
309        points = _sample_sphere(n, dim, seed=0)
310
311        # move all points to one side of the sphere for single-hemisphere test
312        if single_hemisphere:
313            points[:, 0] = np.abs(points[:, 0])
314
315        center = (np.arange(dim) + 1) * shift
316        points = radius * points + center
317
318        sv = SphericalVoronoi(points, radius=radius, center=center)
319        areas = sv.calculate_areas()
320        assert_almost_equal(areas.sum(), _hypersphere_area(dim, radius))
321
322    @pytest.mark.parametrize("poly", ["triangle", "dodecagon",
323                                      "tetrahedron", "cube", "octahedron",
324                                      "dodecahedron", "icosahedron"])
325    def test_equal_area_reconstitution(self, poly):
326        points = _generate_polytope(poly)
327        n, dim = points.shape
328        sv = SphericalVoronoi(points)
329        areas = sv.calculate_areas()
330        assert_almost_equal(areas, _hypersphere_area(dim, 1) / n)
331
332    def test_area_unsupported_dimension(self):
333        dim = 4
334        points = np.concatenate((-np.eye(dim), np.eye(dim)))
335        sv = SphericalVoronoi(points)
336        with pytest.raises(TypeError, match="Only supported"):
337            sv.calculate_areas()
338
339    @pytest.mark.parametrize("radius", [1, 1.])
340    @pytest.mark.parametrize("center", [None, (1, 2, 3), (1., 2., 3.)])
341    def test_attribute_types(self, radius, center):
342        points = radius * self.points
343        if center is not None:
344            points += center
345
346        sv = SphericalVoronoi(points, radius=radius, center=center)
347        assert sv.points.dtype is np.dtype(np.float64)
348        assert sv.center.dtype is np.dtype(np.float64)
349        assert isinstance(sv.radius, float)
350
351    def test_region_types(self):
352        # Tests that region integer type does not change
353        # See Issue #13412
354        sv = SphericalVoronoi(self.points)
355        dtype = type(sv.regions[0][0])
356        sv.sort_vertices_of_regions()
357        assert type(sv.regions[0][0]) == dtype
358        sv.sort_vertices_of_regions()
359        assert type(sv.regions[0][0]) == dtype
360