feat: remove simple usages of scipy

Author: alexfikl

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(