Remove some unnecessary uses of scipy

.basedpyright/baseline.json

@@ -12606,21 +12606,13 @@
         ],
         "./sumpy/point_calculus.py": [
             {
-                "code": "reportUnknownMemberType",
+                "code": "reportAny",
                 "range": {
-                    "startColumn": 39,
-                    "endColumn": 66,
+                    "startColumn": 17,
+                    "endColumn": 60,
                     "lineCount": 1
                 }
             },
-            {
-                "code": "reportUnknownArgumentType",
-                "range": {
-                    "startColumn": 33,
-                    "endColumn": 30,
-                    "lineCount": 3
-                }
-            },
             {
                 "code": "reportUnknownMemberType",
                 "range": {
@@ -15479,14 +15471,6 @@
             }
         ],
         "./sumpy/test/geometries.py": [
-            {
-                "code": "reportAny",
-                "range": {
-                    "startColumn": 4,
-                    "endColumn": 8,
-                    "lineCount": 1
-                }
-            },
             {
                 "code": "reportUnknownArgumentType",
                 "range": {
@@ -21323,6 +21307,14 @@
                     "lineCount": 1
                 }
             },
+            {
+                "code": "reportArgumentType",
+                "range": {
+                    "startColumn": 39,
+                    "endColumn": 40,
+                    "lineCount": 1
+                }
+            },
             {
                 "code": "reportAny",
                 "range": {
@@ -21461,6 +21453,14 @@
                     "lineCount": 1
                 }
             },
+            {
+                "code": "reportArgumentType",
+                "range": {
+                    "startColumn": 43,
+                    "endColumn": 44,
+                    "lineCount": 1
+                }
+            },
             {
                 "code": "reportUnknownArgumentType",
                 "range": {
@@ -21488,8 +21488,8 @@
             {
                 "code": "reportUnknownMemberType",
                 "range": {
-                    "startColumn": 38,
-                    "endColumn": 48,
+                    "startColumn": 37,
+                    "endColumn": 47,
                     "lineCount": 1
                 }
             },
@@ -21504,8 +21504,8 @@
             {
                 "code": "reportAny",
                 "range": {
-                    "startColumn": 44,
-                    "endColumn": 60,
+                    "startColumn": 43,
+                    "endColumn": 59,
                     "lineCount": 1
                 }
             }

doc/misc.rst

@@ -111,7 +111,7 @@ The FAQ is maintained collaboratively on the
 Acknowledgments
 ===============
 
-Work on meshmode was supported in part by
+Work on sumpy was supported in part by
 
 * the US National Science Foundation under grant numbers DMS-1418961,
   DMS-1654756, SHF-1911019, and OAC-1931577.

sumpy/point_calculus.py

@@ -117,8 +117,9 @@ def __init__(self,
             weights_1d = None
 
         elif nodes == "legendre":
-            from scipy.special import legendre
-            points_1d, weights_1d, _ = legendre(npoints).weights.T
+            from numpy.polynomial.legendre import leggauss
+
+            points_1d, weights_1d = leggauss(npoints)
             points_1d = points_1d * (h/2)
             weights_1d = weights_1d * (h/2)
 
@@ -166,10 +167,13 @@ def basis(self) -> Sequence[
             a high-order interpolation basis on the :py:attr:`points`.
         """
 
-        from scipy.special import eval_chebyt
-
         from pytools import indices_in_shape
 
+        def eval_chebyt(n: int,
+                        x: Array1D[np.floating[Any]]) -> Array1D[np.floating[Any]]:
+            # T_n(x) = cos(n * arccos(x)), valid for x in [-1, 1]
+            return np.cos(n * np.arccos(x))
+
         def eval_basis(
                 ind: tuple[int, ...],
                 x: Array2D[np.floating[Any]]

sumpy/test/curve.py

@@ -23,36 +23,45 @@
 THE SOFTWARE.
 """
 
-from typing import TYPE_CHECKING, final
+from typing import TYPE_CHECKING, Any, final
 
 import numpy as np
-from scipy import fftpack
 
 
 if TYPE_CHECKING:
-    from numpy.typing import NDArray
+    import optype.numpy as onp
+
+
+def fftdiff(x: onp.Array1D[np.floating[Any]], *,
+            period: float = 1.0) -> onp.Array1D[np.floating[Any]]:
+    n = len(x)
+    return np.fft.ifft(
+        2j * np.pi * np.fft.fftfreq(n, d=period/n) * np.fft.fft(x)
+    ).real
 
 
 @final
 class CurveGrid:
-    pos: NDArray[np.floating]
-    mean_curvature: NDArray[np.floating]
-    normal: NDArray[np.floating]
+    pos: onp.Array2D[np.floating[Any]]
+    mean_curvature: onp.Array1D[np.floating[Any]]
+    normal: onp.Array2D[np.floating[Any]]
 
-    def __init__(self, x: NDArray[np.floating], y: NDArray[np.floating]):
+    def __init__(self,
+                 x: onp.Array1D[np.floating[Any]],
+                 y: onp.Array1D[np.floating[Any]]) -> None:
         self.pos = np.vstack([x, y]).copy()
-        xp = self.xp = fftpack.diff(x, period=1)
-        yp = self.yp = fftpack.diff(y, period=1)
-        xpp = self.xpp = fftpack.diff(xp, period=1)
-        ypp = self.ypp = fftpack.diff(yp, period=1)
+        xp = self.xp = fftdiff(x, period=1)
+        yp = self.yp = fftdiff(y, period=1)
+        xpp = self.xpp = fftdiff(xp, period=1)
+        ypp = self.ypp = fftdiff(yp, period=1)
         self.mean_curvature = (xp*ypp-yp*xpp)/((xp**2+yp**2)**(3/2))
 
         speed = self.speed = np.sqrt(xp**2+yp**2)
         self.normal = (np.vstack([yp, -xp])/speed).copy()
 
-    def __len__(self):
+    def __len__(self) -> int:
         return len(self.pos)
 
-    def plot(self):
+    def plot(self) -> None:
         import matplotlib.pyplot as pt
         pt.plot(self.pos[:, 0], self.pos[:, 1])

sumpy/test/geometries.py

@@ -24,27 +24,27 @@
 """
 
 from dataclasses import dataclass
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any
 
 import numpy as np
 import numpy.linalg as la
 
 
 if TYPE_CHECKING:
-    from numpy.typing import NDArray
+    import optype.numpy as onp
 
 
 @dataclass(frozen=True)
 class Geometry:
-    nodes: NDArray[np.floating]
-    normals: NDArray[np.floating]
-    weights: NDArray[np.floating]
-    area_elements: NDArray[np.floating]
+    nodes: onp.ArrayND[np.floating[Any]]
+    normals: onp.ArrayND[np.floating[Any]]
+    weights: onp.Array1D[np.floating[Any]]
+    area_elements: onp.Array1D[np.floating[Any]]
 
 
 def _geometry_from_curve_nodes(
-            x: NDArray[np.floating],
-            y: NDArray[np.floating]
+            x: onp.Array1D[np.floating[Any]],
+            y: onp.Array1D[np.floating[Any]]
         ):
     npts: int
     npts, = x.shape
@@ -95,12 +95,16 @@ def make_torus(
         r_minor * np.sin(v)
         ])
 
-    def diff2d(ary: NDArray[np.floating]):
-        from scipy import fftpack
-        return np.array([
-                            fftpack.diff(ary[idx])
-                            for idx in np.ndindex(ary.shape[:-1])
-                        ])
+    def fftdiff(x: onp.Array1D[np.floating[Any]], *,
+                period: float = 2.0 * np.pi) -> onp.Array1D[np.floating[Any]]:
+        n = len(x)
+        return np.fft.ifft(
+            2j * np.pi * np.fft.fftfreq(n, d=period/n) * np.fft.fft(x)
+        ).real
+
+    def diff2d(ary: onp.Array2D[np.floating[Any]]) -> onp.Array2D[np.floating[Any]]:
+        return np.array([fftdiff(ary[idx]) for idx in np.ndindex(ary.shape[:-1])])
+
     dnodes_du = np.array(
                     [diff2d(nodes[i].T).T for i in range(3)])
     dnodes_dv = np.array(

sumpy/test/test_l2l_coeffs.py

@@ -32,7 +32,6 @@
 
 import numpy as np
 import pytest
-import scipy.special as spsp
 
 from arraycontext import (
     pytest_generate_tests_for_array_contexts,
@@ -52,7 +51,7 @@
     LaplaceKernel,
     YukawaKernel,
 )
-from sumpy.tools import build_matrix
+from sumpy.tools import add_mi, build_matrix, mi_factorial, mi_power
 
 
 if TYPE_CHECKING:
@@ -153,12 +152,12 @@ def test_l2l_coefficient_differences(
     max_abs_error = 0.0
 
     for i, nu_i in enumerate(full_identifiers):
-        i_card = sum(np.array(nu_i))
+        i_card = sum(nu_i)
 
         # Compute error by formula
         error = 0.0 + 0.0j
         for k, nu_jk in enumerate(stored_identifiers):
-            jk_card = sum(np.array(nu_jk))
+            jk_card = sum(nu_jk)
             if jk_card >= i_card:
                 continue
 
@@ -167,24 +166,24 @@ def test_l2l_coefficient_differences(
 
             for q_idx in range(start_idx, end_idx):
                 nu_q = full_identifiers[q_idx]
-                nu_sum = tuple(map(sum, zip(nu_q, nu_jk, strict=True)))
+                nu_sum = add_mi(nu_q, nu_jk)
                 if nu_sum not in full_identifiers:
                     continue
 
                 deriv_idx = full_identifiers.index(nu_sum)
                 gamma_deriv = to_scalar(p2l_full.coeffs[deriv_idx])
-                h_pow = np.prod(h ** np.array(nu_q))
-                fact_nu_q = np.prod(spsp.factorial(nu_q))
+                h_pow = float(mi_power(h, nu_q))
+                fact_nu_q = mi_factorial(nu_q)
 
                 error += -M[i, k] * gamma_deriv * h_pow / fact_nu_q
 
-        error /= np.prod(spsp.factorial(nu_i))
+        error /= mi_factorial(nu_i)
         error *= global_const
 
         # Compute direct difference
         true_i_idx = lexpn_idx.index(nu_i)
         mu_full = to_scalar(p2l2l_full.coeffs[true_i_idx])
-        direct_diff = (mu_full - mu_c[i]) / np.prod(spsp.factorial(nu_i))
+        direct_diff = (mu_full - mu_c[i]) / mi_factorial(nu_i)
         direct_diff *= global_const
 
         # Compute errors

sumpy/test/test_m2m_coeffs.py

@@ -32,7 +32,6 @@
 
 import numpy as np
 import pytest
-import scipy.special as spsp
 
 from arraycontext import (
     pytest_generate_tests_for_array_contexts,
@@ -55,7 +54,7 @@
     LaplaceKernel,
     YukawaKernel,
 )
-from sumpy.tools import build_matrix
+from sumpy.tools import add_mi, build_matrix, mi_factorial, mi_power
 
 
 if TYPE_CHECKING:
@@ -191,7 +190,7 @@ def test_m2m_coefficient_differences(
 
         error = 0.0 + 0.0j
         for s, nu_js in enumerate(stored_identifiers):
-            nu_js_card = sum(np.array(nu_js))
+            nu_js_card = sum(nu_js)
             inner_sum = 0.0 + 0.0j
 
             if nu_js_card <= k_card:
@@ -200,14 +199,14 @@ def test_m2m_coefficient_differences(
 
                 for idx in range(start_idx, end_idx):
                     nu_l = full_identifiers[idx]
-                    nu_sum = tuple(a + b for a, b in zip(nu_l, nu_js, strict=True))
+                    nu_sum = add_mi(nu_l, nu_js)
 
                     if nu_sum not in full_identifiers:
                         continue
 
                     derivative_idx = full_identifiers.index(nu_sum)
-                    h_pow = np.prod(h ** np.array(nu_l))
-                    fact_nu_l = np.prod(spsp.factorial(nu_l))
+                    h_pow = float(mi_power(h, nu_l))
+                    fact_nu_l = mi_factorial(nu_l)
 
                     inner_sum += coeffs_full[derivative_idx] * h_pow / fact_nu_l
 
@@ -218,8 +217,8 @@ def test_m2m_coefficient_differences(
 
         if verbose:
             print(f"{k:3d} | {nu_k!s:>15s} | {k_card:6d} | "
-                  f"{error.real: .8e}{error.imag:+.8e}j | "
-                  f"{direct_diff.real: .8e}{direct_diff.imag:+.8e}j | "
+                  f"{error.real:.8e}{error.imag:+.8e}j | "
+                  f"{direct_diff.real:.8e}{direct_diff.imag:+.8e}j | "
                   f"{abs_err:9.2e}")
 
     if verbose: