Repository: skypyproject/skypy
Pull request #643 (2025-01-27 15:45:49)

MAINT: Utility function to check if a value is a pole of the gamma function
## Description

The SkyPy implementation of the upper incomplete gamma function uses `scipy.special.gammasgn` to identify poles of the gamma function. However since SciPy 1.15 the behaviour of that function has changed :

```python
# SciPy 1.14.1
>>> scipy.special.gammasgn([-2, -1, 0])
array([0., 0., 0.])

# SciPy 1.15.0
>>> scipy.special.gammasgn([-2, -1, 0])
array([nan, nan,  1.])
```

This MR introduces a utility function `_is_gamma_pole` that correctly identifies poles of the gamma function using `gammasgn` for all supported versions of SciPy

## Checklist
- [x] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/main/CONTRIBUTING.rst)
- [x] Write unit tests
- [x] Write documentation strings
- [x] Assign someone from your working team to review this pull request
- [x] Assign someone from the infrastructure team to review this pull request


ntessore: Is this not actually a bug in `_gammaincc`? For `x = 0` I would now naively expect

```math
Q(a, 0) = \begin{cases}
  1 & a > 0 \\
  0 & a = 0 \\
  \infty & a < 0
\end{cases}
```

but I cannot remember where we needed the special case here.
rrjbca: Physically all distributions with x=0 and a<=0 would imply infinite source density right? So practically they are of no use to us. I don't remember if there was some mathematical reason to handle integer and non-integer a differently, I just went with the minimal change to maintain the old behaviour and unit tests with all versions of SciPy. 

Patch:
diff --git a/skypy/utils/special.py b/skypy/utils/special.py
index a2fae69..1766a05 100644
--- a/skypy/utils/special.py
+++ b/skypy/utils/special.py
@@ -10,10 +10,22 @@ This module computes useful special functions.
 
 """
 
+from importlib.metadata import version
+from packaging.version import parse
+
 import numpy as np
 import scipy.special as sc
 
 
+def _is_gamma_pole(a):
+    '''Utility function that returns True if a is a pole of the gamma function'''
+    s = sc.gammasgn(a)
+    if parse(version('scipy')) < parse('1.15'):
+        return s == 0
+    else:
+        return np.logical_or(a == 0, np.isnan(s))
+
+
 def _gammaincc(a, x):
     '''gammaincc for positive or negative indices and scalar inputs'''
     if x < 0:
@@ -21,10 +33,9 @@ def _gammaincc(a, x):
     if x == 0:
         if a > 0:
             return 1
-        s = sc.gammasgn(a)
-        if s == 0:
+        if _is_gamma_pole(a):
             return 0
-        return s*np.inf
+        return sc.gammasgn(a) * np.inf
     if np.isinf(x):
         return 0
     if a < 0:
@@ -95,7 +106,7 @@ def gammaincc(a, x):
         # do x = 0 for nonpositive a; depends on Gamma(a)
         i = i & (x == 0)
         s = sc.gammasgn(a[i])
-        g[i] = np.where(s == 0, 0, np.copysign(np.inf, s))
+        g[i] = np.where(_is_gamma_pole(a[i]), 0, np.copysign(np.inf, s))
         n[i] = 0
 
         # these are still to do


Test patch:
diff --git a/skypy/utils/tests/test_special.py b/skypy/utils/tests/test_special.py
index 92b7f10..c228d98 100644
--- a/skypy/utils/tests/test_special.py
+++ b/skypy/utils/tests/test_special.py
@@ -3,7 +3,7 @@ import numpy.testing as npt
 from astropy.utils.data import get_pkg_data_filename
 import pytest
 
-from skypy.utils.special import gammaincc
+from skypy.utils.special import gammaincc, _is_gamma_pole
 
 
 def test_gammaincc_scalar():
@@ -90,3 +90,25 @@ def test_gammaincc_neg_x_array():
     # negative x in array raises an exception
     with pytest.raises(ValueError):
         gammaincc(0.5, [3.0, 2.0, 1.0, 0.0, -1.0])
+
+
+gamma_poles = np.array([-2.0, -1.0, 0.0])
+not_gamma_poles = np.array([-2.5, 0.5, 1.0])
+
+
+@pytest.mark.parametrize('z', gamma_poles)
+def test_is_gamma_pole_scalar(z):
+    assert _is_gamma_pole(z)
+
+
+def test_is_gamma_pole_array():
+    assert np.all(_is_gamma_pole(gamma_poles))
+
+
+@pytest.mark.parametrize('z', not_gamma_poles)
+def test_is_not_gamma_pole_scalar(z):
+    assert not _is_gamma_pole(z)
+
+
+def test_is_not_gamma_pole_array():
+    assert not np.any(_is_gamma_pole(not_gamma_poles))