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