Rmsh

.github/workflows/workflow.yml

@@ -244,7 +244,6 @@ jobs:
       run: |
         conda info
         conda list
-        conda install pyshtools  # debug depends issues
         pip install -e ".[dev]"  # install aspire
         pip freeze
     - name: Test

environment-accelerate.yml

@@ -7,7 +7,6 @@ channels:
 dependencies:
   - pip
   - python=3.8
-  - pyshtools
   - numpy=1.24.1
   - scipy=1.10.1
   - scikit-learn

pyproject.toml

@@ -43,7 +43,6 @@ dependencies = [
     "pillow",
     "psutil",
     "pymanopt",
-    "pyshtools<=4.10.4",  # 4.11.7 might have a packaging bug
     "PyWavelets",
     "ray >= 2.9.2",
     "scipy >= 1.10.0",

src/aspire/basis/basis_utils.py

@@ -4,13 +4,11 @@
 """
 
 import logging
-import warnings
 
 import numpy as np
 from numpy import diff, exp, log, pi
 from numpy.polynomial.legendre import leggauss
-from pyshtools.expand import spharm_lm
-from scipy.special import jn, jv, sph_harm
+from scipy.special import jn, jv
 
 from aspire.utils import grid_2d, grid_3d
 
@@ -158,6 +156,33 @@ def norm_assoc_legendre(j, m, x):
     return px
 
 
+def sph_harm(j, m, theta, phi):
+    """
+    Compute spherical harmonics.
+
+    Note call signature convention may be different from other packages.
+
+    :param m: Order |m| <= j
+    :param j: Harmonic degree, j>=0
+    :param theta: latitude coordinate [0, pi]
+    :param phi: longitude coordinate [0, 2*pi]
+    :return: Complex array of evaluated spherical harmonics.
+    """
+
+    # Compute sph_harm for positive `abs(m)`
+    y = (
+        norm_assoc_legendre(j, abs(m), np.cos(theta))
+        * np.exp(1j * abs(m) * phi)
+        * np.sqrt(0.5 / np.pi)
+    )
+
+    # Use identity for negative `m`
+    if m < 0:
+        y = (-1) ** (m % 2) * np.conj(y)
+
+    return y
+
+
 def real_sph_harmonic(j, m, theta, phi):
     """
     Evaluate a real spherical harmonic
@@ -172,28 +197,8 @@ def real_sph_harmonic(j, m, theta, phi):
     """
     abs_m = abs(m)
 
-    # The `scipy` sph_harm implementation is much faster,
-    #   but incorrectly returns NaN for high orders.
-    # For higher order use `pyshtools`.
-    if j < 86:
-        y = sph_harm(abs_m, j, phi, theta)
-    else:
-        warnings.warn(
-            "Computing higher order spherical harmonics is slow."
-            "  Consider using `FFBBasis3D` or decreasing volume size.",
-            stacklevel=1,
-        )
-
-        y = spharm_lm(
-            j,
-            abs_m,
-            theta,
-            phi,
-            kind="complex",
-            degrees=False,
-            csphase=-1,
-            normalization="ortho",
-        )
+    # Note the calling convention here may not match other `sph_harm` packages
+    y = sph_harm(j, abs_m, theta, phi)
 
     if m < 0:
         y = np.sqrt(2) * np.imag(y)

tests/test_basis_utils.py

@@ -1,6 +1,7 @@
 from unittest import TestCase
 
 import numpy as np
+from scipy.special import sph_harm as sp_sph_harm
 
 from aspire.basis.basis_utils import (
     all_besselj_zeros,
@@ -9,10 +10,61 @@
     norm_assoc_legendre,
     real_sph_harmonic,
     sph_bessel,
+    sph_harm,
     unique_coords_nd,
 )
 
 
+def test_sph_harm_low_order():
+    """
+    Test the `sph_harm` implementation matches `scipy` at lower orders.
+    """
+    m = 3
+    j = 5
+    x = np.linspace(0, np.pi, 42)
+    y = np.linspace(0, 2 * np.pi, 42)
+
+    ref = sp_sph_harm(m, j, y, x)  # Note calling convention is different
+    np.testing.assert_allclose(sph_harm(j, m, x, y), ref)
+
+    # negative m
+    m *= -1
+    ref = sp_sph_harm(m, j, y, x)  # Note calling convention is different
+    np.testing.assert_allclose(sph_harm(j, m, x, y), ref)
+
+
+def test_sph_harm_high_order():
+    """
+    Test we remain finite at higher orders where `scipy.special.sph_harm` overflows.
+    """
+    m = 87
+    j = 87
+    x = 0.12345
+    y = 0.56789
+
+    # If scipy fixed their implementation for higher orders in the future,
+    # this check should fail and we can reconsider that package.
+    ref = sp_sph_harm(m, j, y, x)  # Note calling convention is different
+    assert not np.isfinite(ref)
+
+    # Can manually check against pyshtools,
+    # but we are avoiding that package dependency.
+    # Leaving this here intentionally for future developers.
+    # y = spharm_lm(
+    #     j,
+    #     abs_m,
+    #     theta,
+    #     phi,
+    #     kind="complex",
+    #     degrees=False,
+    #     csphase=-1,
+    #     normalization="ortho",
+    # )
+
+    # Check we are finite.
+    assert np.isfinite(sph_harm(j, m, x, y))
+
+
 class BesselTestCase(TestCase):
     def setUp(self):
         pass