DEP: special: deprecation warning for sph_harm + comments (scipy#22167)

* DEP: Add deprecation warning for sph_harm

* TST: Suppress deprecation warnings in sph_harm tests

* DOC: Add deprecation removal version to sph_harm doc

* TST: Add test for sph_harm deprecation warning

* TST: Add comments explaining why we're testing deprecated functions

Author: steppi

scipy/special/_special_ufuncs_docs.cpp

@@ -5887,8 +5887,8 @@ const char *sph_harm_doc = R"(
     where :math:`P_n^m` are the associated Legendre functions; see `lpmv`.
 
     .. deprecated:: 1.15.0
-        This function is deprecated and will be removed in a future version.
-        Use `scipy.special.sph_harm_y` instead.
+        This function is deprecated and will be removed in SciPy 1.17.0.
+        Please use `scipy.special.sph_harm_y` instead.
 
     Parameters
     ----------

scipy/special/tests/test_basic.py

@@ -4129,6 +4129,8 @@ def test_round(self):
         rndrl = (10,10,10,11)
         assert_array_equal(rnd,rndrl)
 
+# sph_harm is deprecated and is implemented as a shim around sph_harm_y.
+# The following two tests are maintained to verify the correctness of the shim.
 
 def test_sph_harm():
     # Tests derived from tables in
@@ -4139,35 +4141,39 @@ def test_sph_harm():
     sqrt = np.sqrt
     sin = np.sin
     cos = np.cos
-    assert_array_almost_equal(sh(0,0,0,0),
-           0.5/sqrt(pi))
-    assert_array_almost_equal(sh(-2,2,0.,pi/4),
-           0.25*sqrt(15./(2.*pi)) *
-           (sin(pi/4))**2.)
-    assert_array_almost_equal(sh(-2,2,0.,pi/2),
-           0.25*sqrt(15./(2.*pi)))
-    assert_array_almost_equal(sh(2,2,pi,pi/2),
-           0.25*sqrt(15/(2.*pi)) *
-           exp(0+2.*pi*1j)*sin(pi/2.)**2.)
-    assert_array_almost_equal(sh(2,4,pi/4.,pi/3.),
-           (3./8.)*sqrt(5./(2.*pi)) *
-           exp(0+2.*pi/4.*1j) *
-           sin(pi/3.)**2. *
-           (7.*cos(pi/3.)**2.-1))
-    assert_array_almost_equal(sh(4,4,pi/8.,pi/6.),
-           (3./16.)*sqrt(35./(2.*pi)) *
-           exp(0+4.*pi/8.*1j)*sin(pi/6.)**4.)
+    with suppress_warnings() as sup:
+        sup.filter(category=DeprecationWarning)
+        assert_array_almost_equal(sh(0,0,0,0),
+               0.5/sqrt(pi))
+        assert_array_almost_equal(sh(-2,2,0.,pi/4),
+               0.25*sqrt(15./(2.*pi)) *
+               (sin(pi/4))**2.)
+        assert_array_almost_equal(sh(-2,2,0.,pi/2),
+               0.25*sqrt(15./(2.*pi)))
+        assert_array_almost_equal(sh(2,2,pi,pi/2),
+               0.25*sqrt(15/(2.*pi)) *
+               exp(0+2.*pi*1j)*sin(pi/2.)**2.)
+        assert_array_almost_equal(sh(2,4,pi/4.,pi/3.),
+               (3./8.)*sqrt(5./(2.*pi)) *
+               exp(0+2.*pi/4.*1j) *
+               sin(pi/3.)**2. *
+               (7.*cos(pi/3.)**2.-1))
+        assert_array_almost_equal(sh(4,4,pi/8.,pi/6.),
+               (3./16.)*sqrt(35./(2.*pi)) *
+               exp(0+4.*pi/8.*1j)*sin(pi/6.)**4.)
 
 
 def test_sph_harm_ufunc_loop_selection():
     # see https://github.com/scipy/scipy/issues/4895
     dt = np.dtype(np.complex128)
-    assert_equal(special.sph_harm(0, 0, 0, 0).dtype, dt)
-    assert_equal(special.sph_harm([0], 0, 0, 0).dtype, dt)
-    assert_equal(special.sph_harm(0, [0], 0, 0).dtype, dt)
-    assert_equal(special.sph_harm(0, 0, [0], 0).dtype, dt)
-    assert_equal(special.sph_harm(0, 0, 0, [0]).dtype, dt)
-    assert_equal(special.sph_harm([0], [0], [0], [0]).dtype, dt)
+    with suppress_warnings() as sup:
+        sup.filter(category=DeprecationWarning)
+        assert_equal(special.sph_harm(0, 0, 0, 0).dtype, dt)
+        assert_equal(special.sph_harm([0], 0, 0, 0).dtype, dt)
+        assert_equal(special.sph_harm(0, [0], 0, 0).dtype, dt)
+        assert_equal(special.sph_harm(0, 0, [0], 0).dtype, dt)
+        assert_equal(special.sph_harm(0, 0, 0, [0]).dtype, dt)
+        assert_equal(special.sph_harm([0], [0], [0], [0]).dtype, dt)
 
 
 class TestStruve:
@@ -4667,3 +4673,8 @@ def test_warn_xlpmn(self, xlpmn):
         message = f"`scipy.special.{xlpmn.__name__}` is deprecated..."
         with pytest.deprecated_call(match=message):
             _ = xlpmn(1, 1, 0)
+
+    def test_warn_sph_harm(self):
+        msg = "`scipy.special.sph_harm` is deprecated..."
+        with pytest.deprecated_call(match=msg):
+            _ = special.sph_harm(1, 1, 0, 0)

scipy/special/tests/test_data.py

@@ -56,6 +56,15 @@ def data_local(func, dataname, *a, **kw):
     return FuncData(func, DATASETS_LOCAL[dataname], *a, **kw)
 
 
+# The functions lpn, lpmn, clpmn, and sph_harm appearing below are
+# deprecated in favor of legendre_p_all, assoc_legendre_p_all,
+# assoc_legendre_p_all (assoc_legendre_p_all covers lpmn and clpmn),
+# and sph_harm_y respectively. The deprecated functions listed above are
+# implemented as shims around their respective replacements. The replacements
+# are tested separately, but tests for the deprecated functions remain to
+# verify the correctness of the shims.
+
+
 def ellipk_(k):
     return ellipk(k*k)
 
@@ -191,7 +200,9 @@ def spherical_yn_(n, x):
     return spherical_yn(n.astype('l'), x)
 
 def sph_harm_(m, n, theta, phi):
-    y = sph_harm(m, n, theta, phi)
+    with suppress_warnings() as sup:
+        sup.filter(category=DeprecationWarning)
+        y = sph_harm(m, n, theta, phi)
     return (y.real, y.imag)
 
 def cexpm1(x, y):

scipy/special/tests/test_legendre.py

@@ -10,6 +10,14 @@
 from scipy.special import (legendre_p, legendre_p_all, assoc_legendre_p,
     assoc_legendre_p_all, sph_legendre_p, sph_legendre_p_all)
 
+# The functions lpn, lpmn, clpmn, appearing below are
+# deprecated in favor of legendre_p_all, assoc_legendre_p_all, and
+# assoc_legendre_p_all (assoc_legendre_p_all covers lpmn and clpmn)
+# respectively. The deprecated functions listed above are implemented as
+# shims around their respective replacements. The replacements are tested
+# separately, but tests for the deprecated functions remain to verify the
+# correctness of the shims.
+
 # Base polynomials come from Abrahmowitz and Stegan
 class TestLegendre:
     def test_legendre(self):

scipy/special/tests/test_mpmath.py

@@ -671,6 +671,14 @@ def test_lambertw_smallz():
 # Systematic tests
 # ------------------------------------------------------------------------------
 
+# The functions lpn, lpmn, clpmn, and sph_harm appearing below are
+# deprecated in favor of legendre_p_all, assoc_legendre_p_all,
+# assoc_legendre_p_all (assoc_legendre_p_all covers lpmn and clpmn),
+# and sph_harm_y respectively. The deprecated functions listed above are
+# implemented as shims around their respective replacements. The replacements
+# are tested separately, but tests for the deprecated functions remain to
+# verify the correctness of the shims.
+
 HYPERKW = dict(maxprec=200, maxterms=200)
 
 
@@ -2015,7 +2023,9 @@ def test_spherharm(self):
         def spherharm(l, m, theta, phi):
             if m > l:
                 return np.nan
-            return sc.sph_harm(m, l, phi, theta)
+            with suppress_warnings() as sup:
+                sup.filter(category=DeprecationWarning)
+                return sc.sph_harm(m, l, phi, theta)
         assert_mpmath_equal(
             spherharm,
             mpmath.spherharm,

scipy/special/tests/test_sph_harm.py

@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, suppress_warnings
 import scipy.special as sc
 
 class TestSphHarm:
@@ -53,6 +53,9 @@ def test_first_harmonics():
     # `phi` as the polar angle, and include the Condon-Shortley
     # phase.
 
+    # sph_harm is deprecated and is implemented as a shim around sph_harm_y.
+    # This test is maintained to verify the correctness of the shim.
+
     # Notation is Ymn
     def Y00(theta, phi):
         return 0.5*np.sqrt(1/np.pi)
@@ -75,7 +78,9 @@ def Y11(theta, phi):
     theta, phi = np.meshgrid(theta, phi)
 
     for harm, m, n in zip(harms, m, n):
-        assert_allclose(sc.sph_harm(m, n, theta, phi),
-                        harm(theta, phi),
-                        rtol=1e-15, atol=1e-15,
-                        err_msg=f"Y^{m}_{n} incorrect")
+        with suppress_warnings() as sup:
+            sup.filter(category=DeprecationWarning)
+            assert_allclose(sc.sph_harm(m, n, theta, phi),
+                            harm(theta, phi),
+                            rtol=1e-15, atol=1e-15,
+                            err_msg=f"Y^{m}_{n} incorrect")

scipy/special/xsf_special.h

@@ -13,6 +13,14 @@ namespace {
 
 template <typename T>
 std::complex<T> sph_harm(long long int m, long long int n, T theta, T phi) {
+    PyGILState_STATE gstate = PyGILState_Ensure();
+    PyErr_WarnEx(
+        PyExc_DeprecationWarning,
+        "`scipy.special.sph_harm` is deprecated as of SciPy 1.15.0 and will be "
+        "removed in SciPy 1.17.0. Please use `scipy.special.sph_harm_y` instead.",
+        1
+    );
+    PyGILState_Release(gstate);
     if (n < 0) {
         xsf::set_error("sph_harm", SF_ERROR_ARG, "n should not be negative");
         return std::numeric_limits<T>::quiet_NaN();