CI: Add a check that Python.h is included first in any file using Python features.

.github/workflows/array_api.yml

@@ -94,4 +94,5 @@ jobs:
         python dev.py --no-build test -b all -s fft -- --durations 3 --timeout=60
         python dev.py --no-build test -b all -t scipy.special.tests.test_support_alternative_backends -- --durations 3 --timeout=60
         python dev.py --no-build test -b all -t scipy._lib.tests.test_array_api
+        python dev.py --no-build test -b all -t scipy._lib.tests.test__util -- --durations 3 --timeout=60
         python dev.py --no-build test -b all -t scipy.stats.tests.test_stats -- --durations 3 --timeout=60

.github/workflows/lint.yml

@@ -47,3 +47,8 @@ jobs:
         python tools/lint.py --diff-against origin/$GITHUB_BASE_REF
         python tools/unicode-check.py
         python tools/check_test_name.py
+
+    - name: Check that Python.h is first in any file including it.
+      shell: bash
+      run: |
+        python tools/check_python_h_first.py 

dev.py

@@ -921,6 +921,23 @@ def task_lint(fix):
         'doc': 'Lint only files modified since last commit (stricter rules)',
     }
 
+@task_params([])
+def task_check_python_h_first():
+    # Lint just the diff since branching off of main using a
+    # stricter configuration.
+    # emit_cmdstr(os.path.join('tools', 'lint.py') + ' --diff-against main')
+    cmd = "{!s} --diff-against=main".format(
+        Dirs().root / 'tools' / 'check_python_h_first.py'
+    )
+    return {
+        'basename': 'check_python_h_first',
+        'actions': [cmd],
+        'doc': (
+            'Check Python.h order only files modified since last commit '
+            '(stricter rules)'
+        ),
+    }
+
 
 def task_unicode_check():
     # emit_cmdstr(os.path.join('tools', 'unicode-check.py'))
@@ -955,6 +972,7 @@ def run(cls, fix):
             'lint': {'fix': fix},
             'unicode-check': {},
             'check-testname': {},
+            'check_python_h_first': {},
         })
 
 

doc/source/dev/api-dev/array_api.rst

@@ -100,10 +100,13 @@ variable is set:
 Support is provided in `scipy.special` for the following functions:
 `scipy.special.log_ndtr`, `scipy.special.ndtr`, `scipy.special.ndtri`,
 `scipy.special.erf`, `scipy.special.erfc`, `scipy.special.i0`,
-`scipy.special.i0e`, `scipy.special.i1`, `scipy.special.i1e`, 
+`scipy.special.i0e`, `scipy.special.i1`, `scipy.special.i1e`,
 `scipy.special.gammaln`, `scipy.special.gammainc`, `scipy.special.gammaincc`,
 `scipy.special.logit`, and `scipy.special.expit`.
 
+Support is provided in `scipy.stats` for the following functions:
+`scipy.stats.pearsonr` and `scipy.stats.moment`.
+
 
 Implementation notes
 --------------------
@@ -164,7 +167,7 @@ have to change is:
 
 Input array validation uses the following pattern::
 
-   xp = array_namespace(arr) # where arr is the input array 
+   xp = array_namespace(arr) # where arr is the input array
    # alternatively, if there are multiple array inputs, include them all:
    xp = array_namespace(arr1, arr2)
 
@@ -235,10 +238,22 @@ The following pytest markers are available:
   other than the default NumPy backend.
   ``@pytest.mark.usefixtures("skip_xp_backends")`` must be used alongside this
   marker for the skipping to apply.
+* ``skip_xp_invalid_arg`` is used to skip tests that use arguments which
+  are invalid when ``SCIPY_ARRAY_API`` is used. For instance, some tests of
+  `scipy.stats` functions pass masked arrays to the function being tested, but
+  masked arrays are incompatible with the array API. Use of the
+  ``skip_xp_invalid_arg`` decorator allows these tests to protect against
+  regressions when ``SCIPY_ARRAY_API`` is not used without resulting in failures
+  when ``SCIPY_ARRAY_API`` is used. In time, we will want these functions to emit
+  deprecation warnings when they receive array API invalid input, and this
+  decorator will check that the deprecation warning is emitted without it
+  causing the test to fail. When ``SCIPY_ARRAY_API=1`` behavior becomes the
+  default and only behavior, these tests (and the decorator itself) will be
+  removed.
 
 The following is an example using the markers::
 
-  from scipy.conftest import array_api_compatible
+  from scipy.conftest import array_api_compatible, skip_xp_invalid_arg
   ...
   @pytest.mark.skip_xp_backends(np_only=True,
                                  reasons=['skip reason'])
@@ -258,6 +273,13 @@ The following is an example using the markers::
       a = xp.asarray([1, 2, 3])
       b = xp.asarray([0, 2, 5])
       toto(a, b)
+  ...
+  # Do not run when SCIPY_ARRAY_API is used
+  @skip_xp_invalid_arg
+  def test_toto_masked_array(self):
+      a = np.ma.asarray([1, 2, 3])
+      b = np.ma.asarray([0, 2, 5])
+      toto(a, b)
 
 Passing a custom reason to ``reasons`` when ``cpu_only=True`` is unsupported
 since ``cpu_only=True`` can be used alongside passing ``backends``. Also,
@@ -269,7 +291,7 @@ compatibility, one can reduce verbosity by telling ``pytest`` to apply the
 markers to every test function using ``pytestmark``::
 
     from scipy.conftest import array_api_compatible
-    
+
     pytestmark = [array_api_compatible, pytest.mark.usefixtures("skip_xp_backends")]
     skip_xp_backends = pytest.mark.skip_xp_backends
     ...

scipy/_build_utils/_wrappers_common.py

@@ -38,8 +38,8 @@
 USE_OLD_ACCELERATE = ['lsame', 'dcabs1']
 
 C_PREAMBLE = """
-#include "fortran_defs.h"
 #include "npy_cblas.h"
+#include "fortran_defs.h"
 """
 
 LAPACK_DECLS = """
@@ -128,7 +128,8 @@ def get_blas_macro_and_name(name, accelerate):
         elif name == 'xerbla_array':
             return '', name + '__'
     if name in WRAPPED_FUNCS:
-        return '', name + 'wrp_'
+        name = name + 'wrp'
+        return 'F_FUNC', f'{name},{name.upper()}'
     return 'BLAS_FUNC', name
 
 

scipy/_build_utils/src/npy_cblas.h

@@ -26,17 +26,19 @@ enum CBLAS_SIDE {CblasLeft=141, CblasRight=142};
 
 #define CBLAS_INDEX size_t  /* this may vary between platforms */
 
-#ifdef ACCELERATE_NEW_LAPACK
-#define NO_APPEND_FORTRAN
-#define BLAS_SYMBOL_SUFFIX $NEWLAPACK
-#endif
-
 #ifdef NO_APPEND_FORTRAN
 #define BLAS_FORTRAN_SUFFIX
 #else
 #define BLAS_FORTRAN_SUFFIX _
 #endif
 
+// New Accelerate suffix is always $NEWLAPACK (no underscore)
+#ifdef ACCELERATE_NEW_LAPACK
+#undef BLAS_FORTRAN_SUFFIX
+#define BLAS_FORTRAN_SUFFIX
+#define BLAS_SYMBOL_SUFFIX $NEWLAPACK
+#endif
+
 #ifndef BLAS_SYMBOL_PREFIX
 #define BLAS_SYMBOL_PREFIX
 #endif

scipy/_build_utils/src/wrap_dummy_g77_abi.c

@@ -22,8 +22,8 @@ passing a pointer to a variable in which to store the computed result. Unlike
 return values, struct complex arguments work without segfaulting.
 */
 
-#include "fortran_defs.h"
 #include "npy_cblas.h"
+#include "fortran_defs.h"
 
 #ifdef __cplusplus
 extern "C" {

scipy/_build_utils/src/wrap_g77_abi.c

@@ -22,8 +22,8 @@ passing a pointer to a variable in which to store the computed result. Unlike
 return values, struct complex arguments work without segfaulting.
 */
 
-#include "fortran_defs.h"
 #include "npy_cblas.h"
+#include "fortran_defs.h"
 
 #ifdef __cplusplus
 extern "C" {

scipy/_lib/_util.py

@@ -15,7 +15,7 @@
 )
 
 import numpy as np
-from scipy._lib._array_api import array_namespace
+from scipy._lib._array_api import array_namespace, is_numpy
 
 
 AxisError: type[Exception]
@@ -708,23 +708,29 @@ def _nan_allsame(a, axis, keepdims=False):
 
 
 def _contains_nan(a, nan_policy='propagate', use_summation=True,
-                  policies=None):
-    if not isinstance(a, np.ndarray):
+                  policies=None, *, xp=None):
+    if xp is None:
+        xp = array_namespace(a)
+    not_numpy = not is_numpy(xp)
+
+    if not_numpy:
         use_summation = False  # some array_likes ignore nans (e.g. pandas)
     if policies is None:
         policies = ['propagate', 'raise', 'omit']
     if nan_policy not in policies:
         raise ValueError("nan_policy must be one of {%s}" %
                          ', '.join("'%s'" % s for s in policies))
 
-    if np.issubdtype(a.dtype, np.inexact):
-        # The summation method avoids creating a (potentially huge) array.
+    inexact = (xp.isdtype(a.dtype, "real floating")
+               or xp.isdtype(a.dtype, "complex floating"))
+    if inexact:
+        # The summation method avoids creating another (potentially huge) array
         if use_summation:
             with np.errstate(invalid='ignore', over='ignore'):
-                contains_nan = np.isnan(np.sum(a))
+                contains_nan = xp.isnan(xp.sum(a))
         else:
-            contains_nan = np.isnan(a).any()
-    elif np.issubdtype(a.dtype, object):
+            contains_nan = xp.any(xp.isnan(a))
+    elif is_numpy(xp) and np.issubdtype(a.dtype, object):
         contains_nan = False
         for el in a.ravel():
             # isnan doesn't work on non-numeric elements
@@ -738,6 +744,10 @@ def _contains_nan(a, nan_policy='propagate', use_summation=True,
     if contains_nan and nan_policy == 'raise':
         raise ValueError("The input contains nan values")
 
+    if not_numpy and contains_nan and nan_policy=='omit':
+        message = "`nan_policy='omit' is incompatible with non-NumPy arrays."
+        raise ValueError(message)
+
     return contains_nan, nan_policy
 
 

scipy/_lib/tests/test__util.py

@@ -10,9 +10,10 @@
 from pytest import raises as assert_raises
 import hypothesis.extra.numpy as npst
 from hypothesis import given, strategies, reproduce_failure  # noqa: F401
-from scipy.conftest import array_api_compatible
+from scipy.conftest import array_api_compatible, skip_xp_invalid_arg
 
-from scipy._lib._array_api import xp_assert_equal, xp_assert_close
+from scipy._lib._array_api import (xp_assert_equal, xp_assert_close, is_numpy,
+                                   copy as xp_copy)
 from scipy._lib._util import (_aligned_zeros, check_random_state, MapWrapper,
                               getfullargspec_no_self, FullArgSpec,
                               rng_integers, _validate_int, _rename_parameter,
@@ -302,7 +303,7 @@ def test_policy(self):
         with pytest.raises(ValueError, match=msg):
             _contains_nan(data, nan_policy="nan")
 
-    def test_contains_nan_1d(self):
+    def test_contains_nan(self):
         data1 = np.array([1, 2, 3])
         assert not _contains_nan(data1)[0]
 
@@ -312,17 +313,18 @@ def test_contains_nan_1d(self):
         data3 = np.array([np.nan, 2, 3, np.nan])
         assert _contains_nan(data3)[0]
 
-        data4 = np.array([1, 2, "3", np.nan])  # converted to string "nan"
+        data4 = np.array([[1, 2], [3, 4]])
         assert not _contains_nan(data4)[0]
 
-        data5 = np.array([1, 2, "3", np.nan], dtype='object')
+        data5 = np.array([[1, 2], [3, np.nan]])
         assert _contains_nan(data5)[0]
 
-    def test_contains_nan_2d(self):
-        data1 = np.array([[1, 2], [3, 4]])
+    @skip_xp_invalid_arg
+    def test_contains_nan_with_strings(self):
+        data1 = np.array([1, 2, "3", np.nan])  # converted to string "nan"
         assert not _contains_nan(data1)[0]
 
-        data2 = np.array([[1, 2], [3, np.nan]])
+        data2 = np.array([1, 2, "3", np.nan], dtype='object')
         assert _contains_nan(data2)[0]
 
         data3 = np.array([["1", 2], [3, np.nan]])  # converted to string "nan"
@@ -331,6 +333,33 @@ def test_contains_nan_2d(self):
         data4 = np.array([["1", 2], [3, np.nan]], dtype='object')
         assert _contains_nan(data4)[0]
 
+    @array_api_compatible
+    @pytest.mark.parametrize("nan_policy", ['propagate', 'omit', 'raise'])
+    def test_array_api(self, xp, nan_policy):
+        rng = np.random.default_rng(932347235892482)
+        x0 = rng.random(size=(2, 3, 4))
+        x = xp.asarray(x0)
+        x_nan = xp_copy(x, xp=xp)
+        x_nan[1, 2, 1] = np.nan
+
+        contains_nan, nan_policy_out = _contains_nan(x, nan_policy=nan_policy)
+        assert not contains_nan
+        assert nan_policy_out == nan_policy
+
+        if nan_policy == 'raise':
+            message = 'The input contains...'
+            with pytest.raises(ValueError, match=message):
+                _contains_nan(x_nan, nan_policy=nan_policy)
+        elif nan_policy == 'omit' and not is_numpy(xp):
+            message = "`nan_policy='omit' is incompatible..."
+            with pytest.raises(ValueError, match=message):
+                _contains_nan(x_nan, nan_policy=nan_policy)
+        elif nan_policy == 'propagate':
+            contains_nan, nan_policy_out = _contains_nan(
+                x_nan, nan_policy=nan_policy)
+            assert contains_nan
+            assert nan_policy_out == nan_policy
+
 
 def test__rng_html_rewrite():
     def mock_str():

scipy/_lib/tests/test_warnings.py

@@ -115,23 +115,3 @@ def test_warning_calls_filters(warning_calls):
             "found in:\n    {}".format(
                 "\n    ".join(bad_filters)))
 
-
-@pytest.mark.slow
-@pytest.mark.xfail(reason="stacklevels currently missing")
-def test_warning_calls_stacklevels(warning_calls):
-    bad_filters, bad_stacklevels = warning_calls
-
-    msg = ""
-
-    if bad_filters:
-        msg += ("warning ignore filter should not be used, instead, use\n"
-                "numpy.testing.suppress_warnings (in tests only);\n"
-                "found in:\n    {}".format("\n    ".join(bad_filters)))
-        msg += "\n\n"
-
-    if bad_stacklevels:
-        msg += "warnings should have an appropriate stacklevel:\n    {}".format(
-                "\n    ".join(bad_stacklevels))
-
-    if msg:
-        raise AssertionError(msg)

scipy/conftest.py

@@ -151,6 +151,10 @@ def check_fpu_mode(request):
 
 array_api_compatible = pytest.mark.parametrize("xp", xp_available_backends.values())
 
+skip_xp_invalid_arg = pytest.mark.skipif(SCIPY_ARRAY_API,
+    reason = ('Test involves masked arrays, object arrays, or other types '
+              'that are not valid input when `SCIPY_ARRAY_API` is used.'))
+
 
 @pytest.fixture
 def skip_xp_backends(xp, request):

scipy/io/_fast_matrix_market/src/_fmm_core.cpp

@@ -2,7 +2,8 @@
 // Use of this source code is governed by the BSD 2-clause license found in the LICENSE.txt file.
 // SPDX-License-Identifier: BSD-2-Clause
 
-#include <Python.h>
+#include "_fmm_core.hpp"
+
 #include <fast_matrix_market/types.hpp>
 #include <cstdint>
 namespace fast_matrix_market {
@@ -17,8 +18,6 @@ namespace fast_matrix_market {
 }
 #include <fast_matrix_market/fast_matrix_market.hpp>
 
-#include "_fmm_core.hpp"
-
 ////////////////////////////////////////////////
 //// Header methods
 ////////////////////////////////////////////////

scipy/linalg/_cythonized_array_utils.pyx

@@ -198,7 +198,8 @@ cdef inline (int, int) band_check_internal_c(const np_numeric_t[:, ::1]A) noexce
             if A[r, c] != zero:
                 upper_band = c - r
                 break
-        if upper_band == c:
+        # If existing band falls outside matrix; we are done
+        if r + 1 + upper_band > m - 1:
             break
 
     return lower_band, upper_band
@@ -229,7 +230,8 @@ cdef inline (int, int) band_check_internal_noncontig(const np_numeric_t[:, :]A)
             if A[r, c] != zero:
                 upper_band = c - r
                 break
-        if upper_band == c:
+        # If existing band falls outside matrix; we are done
+        if r + 1 + upper_band > m - 1:
             break
 
     return lower_band, upper_band