SciPy 1.15.0 support (#376)

lmo/_poly.py

@@ -113,7 +113,8 @@ def eval_sh_jacobi(  # noqa: C901
         ) / 6
 
     # don't use `eval_sh_jacobi`: https://github.com/scipy/scipy/issues/18988
-    return sps.eval_jacobi(n, a, b, 2 * x - 1)
+    out = sps.eval_jacobi(n, a, b, 2 * x - 1)
+    return out.item() if isinstance(out, np.generic) or np.isscalar(x) else out
 
 
 def peaks_jacobi(n: int, a: float, b: float) -> onp.ArrayND[_Float]:

lmo/distributions/_genlambda.py

@@ -1,6 +1,4 @@
-"""
-Probability distributions, compatible with [`scipy.stats`][scipy.stats].
-"""
+"""Probability distributions, compatible with [`scipy.stats`][scipy.stats]."""
 
 from __future__ import annotations
 
@@ -62,6 +60,9 @@ def _genlambda_qdf(q: _XT, b: float, d: float, f: float) -> _XT:
     return (1 + f) * q ** (b - 1) + (1 - f) * (1 - q) ** (d - 1)  # pyright: ignore[reportReturnType]
 
 
+# pyright: reportIncompatibleMethodOverride=false
+
+
 def _genlambda_cdf0(  # noqa: C901
     x: float,
     b: float,
@@ -133,7 +134,7 @@ def _genlambda_cdf0(  # noqa: C901
 
 class genlambda_gen(rv_continuous):
     @override
-    def _argcheck(self, /, b: _F8, d: _F8, f: _F8) -> np.bool_:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _argcheck(self, /, b: _F8, d: _F8, f: _F8) -> np.bool_:
         return np.isfinite(b) & np.isfinite(d) & (f >= -1) & (f <= 1)
 
     @override
@@ -153,23 +154,24 @@ def _fitstart(
         /,
         data: _ArrF8,
         args: tuple[float, float, float] | None = None,
-    ) -> tuple[float, float, float, float, float]:
+    ) -> tuple[_F8, _F8, _F8, _F8, _F8]:
         #  Arbitrary, but the default f=1 is a bad start
-        return super()._fitstart(data, args or (1.0, 1.0, 0.0))  # pyright: ignore[reportReturnType]
+        loc, scale, b, d, f = super()._fitstart(data, args or (1.0, 1.0, 0.0))
+        return loc, scale, b, d, f
 
     @override
-    def _pdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _pdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:
         return 1 / self._qdf(self._cdf(x, b, d, f), b, d, f)  # pyright: ignore[reportReturnType]
 
     @override
-    def _cdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _cdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:
         return _genlambda_cdf(x, b, d, f)
 
     def _qdf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:
         return _genlambda_qdf(q, b, d, f)
 
     @override
-    def _ppf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _ppf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:
         return _genlambda_ppf(q, b, d, f)
 
     @override

lmo/distributions/_kumaraswamy.py

@@ -5,7 +5,6 @@
 from typing import Any, Final, TypeAlias, final
 
 import numpy as np
-import numpy.typing as npt
 import optype.numpy as onp
 import scipy.special as sps
 
@@ -81,5 +80,11 @@ def _entropy(self, a: float, b: float) -> float:
         return (1 - 1 / b) + (1 - 1 / a) * harmonic(b) - math.log(a * b)
 
     @override
-    def _munp(self, /, n: int | npt.NDArray[np.intp], a: float, b: float) -> _FloatND:
+    def _munp(
+        self,
+        /,
+        n: int | onp.ArrayND[np.intp],
+        a: float,
+        b: float,
+    ) -> _Float | _FloatND:
         return b * sps.beta(1 + n / a, b)

lmo/distributions/_wakeby.py

@@ -173,9 +173,12 @@ def _wakeby_sf0(x: _F8, /, b: _F8, d: _F8, f: _F8) -> _F8:  # noqa: C901
 _wakeby_lm: Final = get_lm_func("wakeby")
 
 
+# pyright: reportIncompatibleMethodOverride = false
+
+
 class wakeby_gen(lmt.rv_continuous):
     @override
-    def _argcheck(self, /, b: _F8, d: _F8, f: _F8) -> bool | np.bool_:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _argcheck(self, /, b: _F8, d: _F8, f: _F8) -> bool | np.bool_:
         return (
             np.isfinite(b)
             & np.isfinite(d)
@@ -207,28 +210,29 @@ def _fitstart(
         /,
         data: _ArrF8,
         args: tuple[float, float, float] | None = None,
-    ) -> tuple[float, float, float, float, float]:
+    ) -> tuple[_F8, _F8, _F8, _F8, _F8]:
         #  Arbitrary, but the default f=1 is a bad start
-        return super()._fitstart(data, args or (1.0, 1.0, 0.5))  # pyright: ignore[reportReturnType]
+        loc, scale, b, d, f = super()._fitstart(data, args or (1.0, 1.0, 0.5))
+        return loc, scale, b, d, f
 
     @override
-    def _pdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _pdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:
         # application of the inverse function theorem
         return 1 / self._qdf(self._cdf(x, b, d, f), b, d, f)  # pyright: ignore[reportReturnType]
 
     @override
-    def _cdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _cdf(self, /, x: _XT, b: float, d: float, f: float) -> _XT:
         return 1 - _wakeby_sf(x, b, d, f)
 
     def _qdf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:
         return _wakeby_qdf(q, b, d, f)
 
     @override
-    def _ppf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _ppf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:
         return _wakeby_isf(1 - q, b, d, f)
 
     @override
-    def _isf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:  # pyright: ignore[reportIncompatibleMethodOverride]
+    def _isf(self, /, q: _XT, b: float, d: float, f: float) -> _XT:
         return _wakeby_isf(q, b, d, f)
 
     @override

lmo/ostats.py

@@ -145,5 +145,5 @@ def from_cdf(
         msg = "F must lie between 0 and 1"
         raise ValueError(msg)
 
-    out = betainc(i + 1, n - i, p)
+    out = betainc(float(i + 1), float(n - i), p)
     return out.item() if out.ndim == 0 and np.isscalar(F) else out

lmo/special.py

@@ -120,17 +120,19 @@ def gamma2(
         - [`scipy.special.gammaincc`][scipy.special.gammaincc] for the
           regularized gamma function \( Q(a,\ x) \).
     """
+    x_ = np.asarray(x)
     if a == 0:
-        return sps.expm1(x, out=out)
+        out_ = np.expm1(x_, out=out)
+        return out_.item() if x_.ndim == 0 and np.isscalar(x) else out_
 
     g = sps.gamma(a)
 
     if out is not None:
-        out = sps.gammaincc(a, x, out=out)
+        sps.gammaincc(a, x_, out=out)
         np.multiply(out, g, out=out)
         return out
 
-    res = sps.gammaincc(a, x)
+    res = sps.gammaincc(a, x_)
     res *= g
     return res
 

lmo/theoretical/_utils.py

@@ -13,11 +13,13 @@
 
     import lmo.typing as lmt
 
-__all__ = ("ALPHA", "QUAD_LIMIT", "l_coef_factor", "l_const", "tighten_cdf_support")
+__all__ = "ALPHA", "QUAD_LIMIT", "l_coef_factor", "l_const", "tighten_cdf_support"
 
+###
 
 _Tss = ParamSpec("_Tss")
 
+###
 
 ALPHA: Final = 0.1
 QUAD_LIMIT: Final = 100
@@ -36,17 +38,14 @@ def l_const(r: int, s: float, t: float, k: int = 0) -> float:
             return 1 / (r - 1)
 
     # math.lgamma is faster than scipy.special.loggamma.
-    if r + s + t <= 20:
-        v = gamma(r + s + t + 1) / (gamma(r + s) * gamma(r + t))
-    elif r + s + t <= 128:
-        v = exp(lgamma(r + s + t + 1) - lgamma(r + s) - lgamma(r + t))
+    rst = r + s + t
+    if rst <= 20:
+        v = gamma(rst + 1) / (gamma(r + s) * gamma(r + t))
+    elif rst <= 128:
+        v = exp(lgamma(rst + 1) - lgamma(r + s) - lgamma(r + t))
     else:
         return exp(
-            lgamma(r + s + t + 1)
-            - lgamma(r + s)
-            - lgamma(r + t)
-            + lgamma(r - k)
-            - log(r)
+            lgamma(rst + 1) - lgamma(r + s) - lgamma(r + t) + lgamma(r - k) - log(r)
         )
 
     return factorial(r - 1 - k) / r * v
@@ -57,15 +56,18 @@ def l_coef_factor(
     s: float = 0,
     t: float = 0,
 ) -> onp.ArrayND[npc.floating]:
+    r_ = r if isinstance(r, np.ndarray) else int(r)
     if s == t == 0:
-        return np.sqrt(2 * r - 1)
+        return np.sqrt(2 * r_ - 1)
 
     assert s + t > -1
 
     import scipy.special as sps
 
-    rst = r + s + t
-    return np.sqrt((rst + r - 1) * sps.beta(r + s, r + t) / sps.beta(r, rst)) * r / rst
+    rs = r_ + s
+    rt = r_ + t
+    rst = rs + t
+    return np.sqrt((rs + rt - 1) * sps.beta(rs, rt) / sps.beta(r_, rst)) * r_ / rst
 
 
 def tighten_cdf_support(
@@ -75,16 +77,12 @@ def tighten_cdf_support(
     """Attempt to tighten the support by checking some common bounds."""
     a, b = (-np.inf, np.inf) if support is None else map(float, support)
 
-    # attempt to tighten the default support by checking some common bounds
-    if cdf(0) == 0:
-        # left-bounded at 0 (e.g. weibull)
+    if cdf(0) == 0:  # left-bounded at 0 (e.g. weibull)
         a = 0
 
-        if (u1 := cdf(1)) == 0:
-            # left-bounded at 1 (e.g. pareto)
+        if (u1 := cdf(1)) == 0:  # left-bounded at 1 (e.g. pareto)
             a = 1
-        elif u1 == 1:
-            # right-bounded at 1 (e.g. beta)
+        elif u1 == 1:  # right-bounded at 1 (e.g. beta)
             b = 1
 
     return a, b