SphereIcosahedral & SphereHealpix: allow subdivisions that aren't pow…

…ers of 2, consistent with SphereCubed

doc/background/spherical_samplings.rst

@@ -148,12 +148,12 @@ Subdivision of the icosahedron (:class:`~pygsp.graphs.SphereIcosahedral`)
 -------------------------------------------------------------------------
 
 The sampling is made of an `icosahedron`_ (a `regular and convex polyhedron`_
-made of 12 vertices and 20 triangular faces) whose faces are recursively
-subdivided (into four triangles) to the desired resolution.
-The result is a :math:`\{3,5+\}_{2^m,0}` `geodesic polyhedron`_ (made of
-:math:`20⋅4^m` triangles, :math:`n=10⋅4^m+2` vertices) or its dual, a
-:math:`\{5+,3\}_{2^m,0}` `Goldberg polyhedron`_ (made of :math:`10⋅4^m-10`
-hexagons, 12 pentagons, :math:`n=20⋅4^m` vertices).
+made of 12 vertices and 20 triangular faces) whose faces are subdivided into
+:math:`m^2` triangles projected to the sphere.
+The result is a :math:`\{3,5+\}_{m,0}` `geodesic polyhedron`_ (made of
+:math:`20⋅m^2` triangles, :math:`n=10⋅m^2+2` vertices) or its dual, a
+:math:`\{5+,3\}_{m,0}` `Goldberg polyhedron`_ (made of :math:`10⋅m^2-10`
+hexagons, 12 pentagons, :math:`n=20⋅m^2` vertices).
 The resulting `polyhedral graph`_ (i.e., a 3-vertex-connected planar graph) is
 the 1-`skeleton`_ of the polyhedron.
 All have `icosahedral symmetry`_.
@@ -176,8 +176,8 @@ SHT.
 Subdivision of the cube (:class:`~pygsp.graphs.SphereCubed`)
 ------------------------------------------------------------
 
-The sampling is made of a cube whose faces are subdivided into finer
-quadrilaterals to the desired resolution.
+The sampling is made of a cube whose faces are subdivided into :math:`m^2`
+finer quadrilaterals projected to the sphere.
 The result is a convex polyhedron made of :math:`n=6⋅m^2` faces.
 The graph vertices represent the quadrilateral faces.
 
@@ -190,9 +190,9 @@ Subdivision of the rhombic dodecahedron (:class:`~pygsp.graphs.SphereHealpix`)
 ------------------------------------------------------------------------------
 
 The sampling is made of a `rhombic dodecahedron`_ (a convex polyhedron made of
-12 rhombic faces) whose faces are recursively subdivided (into four
-quadrilaterals) to the desired resolution.
-The result is a convex polyhedron made of :math:`n=12⋅4^m` faces.
+12 rhombic faces) whose faces are subdivided into :math:`m^2` finer
+quadrilaterals projected to the sphere.
+The result is a convex polyhedron made of :math:`n=12⋅m^2` faces.
 The graph vertices represent the quadrilateral faces.
 
 The Hierarchical Equal Area isoLatitude Pixelisation (`HEALPix`_) [Go]_ was

pygsp/graphs/nngraphs/spherehealpix.py

@@ -27,8 +27,9 @@ class SphereHealpix(NNGraph):
     Parameters
     ----------
     subdivisions : int
-        Number of recursive subdivisions. Known as ``order`` in HEALPix
-        terminology, with ``nside=2**order`` and ``npix=12*4**order``.
+        Number of edges the dodecahedron's edges are divided into (``nside``),
+        resulting in ``12*subdivisions**2`` vertices (``npix``).
+        It must be a power of 2 if ``nest=True`` (to be hierarchical).
     indexes : array_like of int
         Indexes of the pixels from which to build a graph. Useful to build a
         graph from a subset of the pixels, e.g., for partial sky observations.
@@ -66,7 +67,7 @@ class SphereHealpix(NNGraph):
     Examples
     --------
     >>> import matplotlib.pyplot as plt
-    >>> G = graphs.SphereHealpix(1, k=8)
+    >>> G = graphs.SphereHealpix()
     >>> fig = plt.figure()
     >>> ax1 = fig.add_subplot(131)
     >>> ax2 = fig.add_subplot(132, projection='3d')
@@ -80,44 +81,43 @@ class SphereHealpix(NNGraph):
 
     >>> import matplotlib.pyplot as plt
     >>> fig, axes = plt.subplots(1, 2)
-    >>> graph = graphs.SphereHealpix(1, nest=False, k=8)
+    >>> graph = graphs.SphereHealpix(2, nest=False, k=8)
     >>> graph.set_coordinates('sphere', dim=2)
     >>> _ = graph.plot(indices=True, ax=axes[0], title='RING ordering')
-    >>> graph = graphs.SphereHealpix(1, nest=True, k=8)
+    >>> graph = graphs.SphereHealpix(2, nest=True, k=8)
     >>> graph.set_coordinates('sphere', dim=2)
     >>> _ = graph.plot(indices=True, ax=axes[1], title='NESTED ordering')
 
     """
 
-    def __init__(self, subdivisions=1, indexes=None, nest=False, **kwargs):
+    def __init__(self, subdivisions=2, indexes=None, nest=False, **kwargs):
         hp = _import_hp()
 
-        nside = hp.order2nside(subdivisions)
         self.subdivisions = subdivisions
         self.nest = nest
 
         if indexes is None:
-            npix = hp.nside2npix(nside)
+            npix = hp.nside2npix(subdivisions)
             indexes = np.arange(npix)
 
-        x, y, z = hp.pix2vec(nside, indexes, nest=nest)
+        x, y, z = hp.pix2vec(subdivisions, indexes, nest=nest)
         coords = np.stack([x, y, z], axis=1)
 
         k = kwargs.pop('k', 20)
         try:
             kernel_width = kwargs.pop('kernel_width')
         except KeyError:
             try:
-                kernel_width = _OPTIMAL_KERNEL_WIDTHS[k][nside]
+                kernel_width = _OPTIMAL_KERNEL_WIDTHS[k][subdivisions]
             except KeyError:
-                raise ValueError('No known optimal kernel width for {} '
-                                 'neighbors and nside={}.'.format(k, nside))
+                raise ValueError('No optimal kernel width for {} neighbors and'
+                                 ' {} subdivisions.'.format(k, subdivisions))
 
         super(SphereHealpix, self).__init__(coords, k=k,
                                             kernel_width=kernel_width,
                                             **kwargs)
 
-        lat, lon = hp.pix2ang(nside, indexes, nest=nest, lonlat=False)
+        lat, lon = hp.pix2ang(subdivisions, indexes, nest=nest, lonlat=False)
         self.signals['lat'] = np.pi/2 - lat  # colatitude to latitude
         self.signals['lon'] = lon
 
@@ -130,7 +130,7 @@ def _get_extra_repr(self):
         return attrs
 
 
-# TODO: find an interpolation between nside and k (#neighbors).
+# TODO: find an interpolation between subdivisions (nside) and k (#neighbors).
 _OPTIMAL_KERNEL_WIDTHS = {
     8: {
         1:    0.02500 * 32,  # extrapolated

pygsp/graphs/nngraphs/sphereicosahedral.py

@@ -26,7 +26,10 @@ class SphereIcosahedral(NNGraph):
     Parameters
     ----------
     subdivisions : int
-        Number of recursive subdivisions.
+        Number of edges the icosahedron's edges are divided into, resulting in
+        ``10*subdivisions**2+2`` vertices (or ``20*subdivisions**2`` if
+        ``dual=True``).
+        It must be a power of 2 in the current implementation.
     dual : bool
         Whether the graph vertices correspond to the vertices (``dual=False``)
         or the triangular faces (``dual=True``) of the subdivided icosahedron.
@@ -77,15 +80,15 @@ class SphereIcosahedral(NNGraph):
 
     >>> import matplotlib.pyplot as plt
     >>> fig, axes = plt.subplots(1, 2)
-    >>> graph = graphs.SphereIcosahedral(0, dual=False, k=5)
+    >>> graph = graphs.SphereIcosahedral(1, dual=False, k=5)
     >>> graph.set_coordinates('sphere', dim=2)
     >>> _ = graph.plot(indices=True, ax=axes[0], title='Icosahedron')
-    >>> graph = graphs.SphereIcosahedral(0, dual=True, k=3)
+    >>> graph = graphs.SphereIcosahedral(1, dual=True, k=3)
     >>> graph.set_coordinates('sphere', dim=2)
     >>> _ = graph.plot(indices=True, ax=axes[1], title='Dodecahedron')
 
     """
-    def __init__(self, subdivisions=1, dual=False, **kwargs):
+    def __init__(self, subdivisions=2, dual=False, **kwargs):
 
         self.subdivisions = subdivisions
         self.dual = dual
@@ -115,7 +118,11 @@ def normalize(vertices):
             """Project the vertices on the sphere."""
             vertices /= np.linalg.norm(vertices, axis=1)[:, None]
 
-        for _ in range(subdivisions):
+        if np.log2(subdivisions) % 1 != 0:
+            raise NotImplementedError('Only recursive subdivisions by two are '
+                                      'implemented. Choose a power of two.')
+
+        for _ in range(int(np.log2(subdivisions))):
             mesh = mesh.subdivide()
             # TODO: shall we project between subdivisions? Some do, some don't.
             # Projecting pushes points away from the 12 base vertices, which

pygsp/tests/test_graphs.py

@@ -556,21 +556,22 @@ def test_sphere_gausslegendre(self, nrings=11):
         self.assertRaises(NotImplementedError, graphs.SphereGaussLegendre,
                           reduced='glesp-equal-area')
 
-    def test_sphere_icosahedral(self, subdivisions=3):
-        n_faces = 20 * 4**subdivisions
-        n_edges = 30 * 4**subdivisions
+    def test_sphere_icosahedral(self, subdivisions=8):
+        n_faces = 20 * subdivisions**2
+        n_edges = 30 * subdivisions**2
         n_vertices = n_edges - n_faces + 2
         graph = graphs.SphereIcosahedral(subdivisions, kind='radius',
-                                         radius=1.6/2**subdivisions)
+                                         radius=1.6/subdivisions)
         self.assertEqual(graph.n_vertices, n_vertices)
         self.assertEqual(graph.n_edges, n_edges)
         self._test_sphere(graph)
-        self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*2**subdivisions)
+        self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*subdivisions)
         graph = graphs.SphereIcosahedral(subdivisions, dual=True, k=3)
         self.assertEqual(graph.n_vertices, n_faces)
         self.assertEqual(graph.n_edges, n_edges)
         self._test_sphere(graph)
-        self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*2**subdivisions)
+        self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*subdivisions)
+        self.assertRaises(NotImplementedError, graphs.SphereIcosahedral, 3)
 
     def test_sphere_cubed(self, subdivisions=9):
         for spacing in ['equiangular', 'equidistant']:
@@ -585,17 +586,20 @@ def test_sphere_cubed(self, subdivisions=9):
                   [-1, 0, 0], [0, -1, 0], [0, 0, -1]]
         np.testing.assert_allclose(graph.coords, coords)
 
-    def test_sphere_healpix(self, subdivisions=2):
-        nside = 2**subdivisions
+    def test_sphere_healpix(self, subdivisions=4):
         graph = graphs.SphereHealpix(subdivisions)
-        self.assertEqual(graph.n_vertices, 12*nside**2)
+        self.assertEqual(graph.n_vertices, 12*subdivisions**2)
         self._test_sphere(graph)
         graphs.SphereHealpix(subdivisions, k=20)
         graphs.SphereHealpix(subdivisions, k=21, kernel_width=0.1)
         self.assertRaises(ValueError, graphs.SphereHealpix, subdivisions, k=21)
         graphs.SphereHealpix(subdivisions, indexes=range(10), k=8)
         graphs.SphereHealpix(subdivisions, nest=True)
-        self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*nside)
+        self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*subdivisions)
+        # Subdivisions must be powers of 2 for NESTED ordering.
+        graphs.SphereHealpix(3, nest=False, kernel_width=1)
+        self.assertRaises(ValueError, graphs.SphereHealpix, 3, nest=True,
+                          kernel_width=1)
 
     def test_twomoons(self):
         graphs.TwoMoons(moontype='standard')