Merge pull request #18779 from tylerjereddy/treddy_1_11_1_prep

MAINT: SciPy 1.11.1 backports

doc/source/release/1.11.1-notes.rst

@@ -5,19 +5,39 @@ SciPy 1.11.1 Release Notes
 .. contents::
 
 SciPy 1.11.1 is a bug-fix release with no new features
-compared to 1.11.0.
-
+compared to 1.11.0. In particular, a licensing issue
+discovered after the release of 1.11.0 has been addressed.
 
 
 Authors
 =======
+
 * Name (commits)
+* h-vetinari (1)
+* Robert Kern (1)
+* Ilhan Polat (4)
+* Tyler Reddy (8)
 
+A total of 4 people contributed to this release.
+People with a "+" by their names contributed a patch for the first time.
+This list of names is automatically generated, and may not be fully complete.
 
 Issues closed for 1.11.1
 ------------------------
 
+* `#18739 <https://github.com/scipy/scipy/issues/18739>`__: BUG: run method of scipy.odr.ODR class fails when delta0 parameter...
+* `#18751 <https://github.com/scipy/scipy/issues/18751>`__: BUG: segfault in \`scipy.linalg.lu\` on x86_64 windows and macos...
+* `#18753 <https://github.com/scipy/scipy/issues/18753>`__: BUG: factorial return type inconsistent for 0-dim arrays
+* `#18759 <https://github.com/scipy/scipy/issues/18759>`__: determinant of a 1x1 matrix returns an array, not a scalar
+* `#18765 <https://github.com/scipy/scipy/issues/18765>`__: Licensing concern
 
 
 Pull requests for 1.11.1
 ------------------------
+
+* `#18741 <https://github.com/scipy/scipy/pull/18741>`__: BUG: Fix work array construction for various weight shapes.
+* `#18747 <https://github.com/scipy/scipy/pull/18747>`__: REL, MAINT: prep for 1.11.1
+* `#18754 <https://github.com/scipy/scipy/pull/18754>`__: BUG: fix handling for \`factorial(..., exact=False)\` for 0-dim...
+* `#18762 <https://github.com/scipy/scipy/pull/18762>`__: FIX:linalg.lu:Guard against permute_l out of bound behavior
+* `#18763 <https://github.com/scipy/scipy/pull/18763>`__: MAINT:linalg.det:Return scalars for singleton inputs
+* `#18778 <https://github.com/scipy/scipy/pull/18778>`__: MAINT: fix unuran licensing

scipy/linalg/_basic.py

@@ -1001,7 +1001,8 @@ def det(a, overwrite_a=False, check_finite=True):
     det : (...) float or complex
         Determinant of `a`. For stacked arrays, a scalar is returned for each
         (m, m) slice in the last two dimensions of the input. For example, an
-        input of shape (p, q, m, m) will produce a result of shape (p, q).
+        input of shape (p, q, m, m) will produce a result of shape (p, q). If
+        all dimensions are 1 a scalar is returned regardless of ndim.
 
     Notes
     -----
@@ -1066,11 +1067,17 @@ def det(a, overwrite_a=False, check_finite=True):
 
     # Scalar case
     if a1.shape[-2:] == (1, 1):
-        if a1.dtype.char in 'dD':
-            return np.squeeze(a1)
+        # Either ndarray with spurious singletons or a single element
+        if max(*a1.shape) > 1:
+            temp = np.squeeze(a1)
+            if a1.dtype.char in 'dD':
+                return temp
+            else:
+                return (temp.astype('d') if a1.dtype.char == 'f' else
+                        temp.astype('D'))
         else:
-            return (np.squeeze(a1).astype('d') if a1.dtype.char == 'f' else
-                    np.squeeze(a1).astype('D'))
+            return (np.float64(a1.item()) if a1.dtype.char in 'fd' else
+                    np.complex128(a1.item()))
 
     # Then check overwrite permission
     if not _datacopied(a1, a):  # "a"  still alive through "a1"

scipy/linalg/_decomp_lu_cython.pyx

@@ -2,17 +2,17 @@
 
 import cython
 from cpython.mem cimport PyMem_Malloc, PyMem_Free
-from libc.string cimport memcpy
-
 from scipy.linalg.cython_lapack cimport sgetrf, dgetrf, cgetrf, zgetrf
-from scipy.linalg._cythonized_array_utils cimport lapack_t, swap_c_and_f_layout
+from scipy.linalg._cythonized_array_utils cimport swap_c_and_f_layout
 
 cimport numpy as cnp
 cnp.import_array()
 
-
-ctypedef double complex doubleComplex
-ctypedef float complex floatComplex
+ctypedef fused lapack_t:
+    cnp.float32_t
+    cnp.float64_t
+    cnp.complex64_t
+    cnp.complex128_t
 
 
 @cython.nonecheck(False)
@@ -24,12 +24,12 @@ cdef void lu_decompose(cnp.ndarray[lapack_t, ndim=2] a,
                        int[::1] perm,
                        bint permute_l) noexcept:
     """LU decomposition and copy operations using ?getrf routines
-    
+
     This function overwrites inputs. For interfacing LAPACK,
     it creates a memory buffer and copies into with F-order
     then swaps back to C order hence no need for dealing with
     Fortran arrays which are inconvenient.
-    
+
     After the LU factorization, to minimize the amount of data
     copied, for rectangle arrays, and depending on the size,
     the smaller portion is copied out to U and the rest becomes
@@ -59,24 +59,24 @@ cdef void lu_decompose(cnp.ndarray[lapack_t, ndim=2] a,
     dims[0] = m
     dims[1] = n
 
-    if lapack_t is float:
+    if lapack_t is cnp.float32_t:
         b = cnp.PyArray_SimpleNew(2, dims, cnp.NPY_FLOAT32)
-        bb = <float *>cnp.PyArray_DATA(b)
+        bb = <cnp.float32_t *>cnp.PyArray_DATA(b)
         swap_c_and_f_layout(aa, bb, m, n)
         sgetrf(&m, &n, bb, &m, ipiv, &info)
-    elif lapack_t is double:
+    elif lapack_t is cnp.float64_t:
         b = cnp.PyArray_SimpleNew(2, dims, cnp.NPY_FLOAT64)
-        bb = <double *>cnp.PyArray_DATA(b)
+        bb = <cnp.float64_t *>cnp.PyArray_DATA(b)
         swap_c_and_f_layout(aa, bb, m, n)
         dgetrf(&m, &n, bb, &m, ipiv, &info)
-    elif lapack_t is floatcomplex:
+    elif lapack_t is cnp.complex64_t:
         b = cnp.PyArray_SimpleNew(2, dims, cnp.NPY_COMPLEX64)
-        bb = <floatComplex*>cnp.PyArray_DATA(b)
+        bb = <cnp.complex64_t *>cnp.PyArray_DATA(b)
         swap_c_and_f_layout(aa, bb, m, n)
         cgetrf(&m, &n, bb, &m, ipiv, &info)
     else:
         b = cnp.PyArray_SimpleNew(2, dims, cnp.NPY_COMPLEX128)
-        bb = <doubleComplex *>cnp.PyArray_DATA(b)
+        bb = <cnp.complex128_t *>cnp.PyArray_DATA(b)
         swap_c_and_f_layout(aa, bb, m, n)
         zgetrf(&m, &n, bb, &m, ipiv, &info)
 
@@ -103,7 +103,9 @@ cdef void lu_decompose(cnp.ndarray[lapack_t, ndim=2] a,
         # as final. Solution without argsort : ipiv[perm] = np.arange(m)
         for ind1 in range(m):
             ipiv[perm[ind1]] = ind1
-        memcpy(&perm[0], ipiv, m*sizeof(int))
+        for ind1 in range(m):
+            perm[ind1] = ipiv[ind1]
+
     finally:
         PyMem_Free(ipiv)
 
@@ -117,7 +119,7 @@ cdef void lu_decompose(cnp.ndarray[lapack_t, ndim=2] a,
         for ind1 in range(mn):
             a[ind1, ind1] = 1
             a[ind1, ind1+1:mn] = 0
-            
+
     else:  # square or fat, "a" holds bigger U
 
         lu[0, 0] = 1
@@ -128,32 +130,41 @@ cdef void lu_decompose(cnp.ndarray[lapack_t, ndim=2] a,
         for ind2 in range(mn - 1):  # cols
             for ind1 in range(ind2+1, m):  # rows
                 a[ind1, ind2] = 0
-    
+
     if permute_l:
         # b still exists -> use it as temp array
         # we copy everything to b and pick back
         # rows from b as dictated by perm
-        memcpy(bb, 
-               (&a[0, 0] if m > n else &lu[0, 0]),
-               m*n*sizeof(lapack_t)
-              )
-        for ind1 in range(m):
-            if perm[ind1] == ind1:  # Row is not permuted
-                continue
-            else:
-                memcpy((&a[ind1, 0] if m > n else &lu[ind1, 0]),
-                       bb + (perm[ind1]*mn),  # stride rows
-                       mn*sizeof(lapack_t))  # copy 1 row of mem
+
+        if m > n:
+            b[:, :] = a[:, :]
+            # memcpy(bb, &a[0, 0], m*mn*sizeof(lapack_t))
+            for ind1 in range(m):
+                if perm[ind1] == ind1:
+                    continue
+                else:
+                    a[ind1, :] = b[perm[ind1], :]
+
+        else:  # same but for lu array
+            b[:mn, :mn] = lu[:, :]
+            # memcpy(bb, &lu[0, 0], mn*n*sizeof(lapack_t))
+            for ind1 in range(mn):
+                if perm[ind1] == ind1:
+                    continue
+                else:
+                    lu[ind1, :] = b[perm[ind1], :mn]
 
 
 @cython.nonecheck(False)
 @cython.initializedcheck(False)
 def lu_dispatcher(a, u, piv, permute_l):
     if a.dtype.char == 'f':
-        lu_decompose[float](a, u, piv, permute_l)
+        lu_decompose[cnp.float32_t](a, u, piv, permute_l)
     elif a.dtype.char == 'd':
-        lu_decompose[double](a, u, piv, permute_l)
+        lu_decompose[cnp.float64_t](a, u, piv, permute_l)
     elif a.dtype.char == 'F':
-        lu_decompose[floatcomplex](a, u, piv, permute_l)
+        lu_decompose[cnp.complex64_t](a, u, piv, permute_l)
+    elif a.dtype.char == 'D':
+        lu_decompose[cnp.complex128_t](a, u, piv, permute_l)
     else:
-        lu_decompose[doublecomplex](a, u, piv, permute_l)
+        raise TypeError("Unsupported type given to lu_dispatcher")

scipy/linalg/tests/test_basic.py

@@ -932,6 +932,23 @@ class TestDet:
     def setup_method(self):
         self.rng = np.random.default_rng(1680305949878959)
 
+    def test_1x1_all_singleton_dims(self):
+        a = np.array([[1]])
+        deta = det(a)
+        assert deta.dtype.char == 'd'
+        assert np.isscalar(deta)
+        assert deta == 1.
+        a = np.array([[[[1]]]], dtype='f')
+        deta = det(a)
+        assert deta.dtype.char == 'd'
+        assert np.isscalar(deta)
+        assert deta == 1.
+        a = np.array([[[1 + 3.j]]], dtype=np.complex64)
+        deta = det(a)
+        assert deta.dtype.char == 'D'
+        assert np.isscalar(deta)
+        assert deta == 1.+3.j
+
     def test_1by1_stacked_input_output(self):
         a = self.rng.random([4, 5, 1, 1], dtype=np.float32)
         deta = det(a)

scipy/odr/_odrpack.py

@@ -899,10 +899,16 @@ def _gen_work(self):
         elif len(self.data.we.shape) == 3:
             ld2we, ldwe = self.data.we.shape[1:]
         else:
-            # Okay, this isn't precisely right, but for this calculation,
-            # it's fine
+            we = self.data.we
             ldwe = 1
-            ld2we = self.data.we.shape[1]
+            ld2we = 1
+            if we.ndim == 1 and q == 1:
+                ldwe = n
+            elif we.ndim == 2:
+                if we.shape == (q, q):
+                    ld2we = q
+                elif we.shape == (q, n):
+                    ldwe = n
 
         if self.job % 10 < 2:
             # ODR not OLS

scipy/odr/tests/test_odr.py

@@ -531,3 +531,32 @@ def func(b, x):
         p = Model(func)
         p.set_meta(name='Sample Model Meta', ref='ODRPACK')
         assert_equal(p.meta, {'name': 'Sample Model Meta', 'ref': 'ODRPACK'})
+
+    def test_work_array_del_init(self):
+        """
+        Verify fix for gh-18739 where del_init=1 fails.
+        """
+        def func(b, x):
+            return b[0] + b[1] * x
+
+        # generate some data
+        n_data = 4
+        x = np.arange(n_data)
+        y = np.where(x % 2, x + 0.1, x - 0.1)
+        x_err = np.full(n_data, 0.1)
+        y_err = np.full(n_data, 0.1)
+
+        linear_model = Model(func)
+        # Try various shapes of the `we` array from various `sy` and `covy`
+        rd0 = RealData(x, y, sx=x_err, sy=y_err)
+        rd1 = RealData(x, y, sx=x_err, sy=0.1)
+        rd2 = RealData(x, y, sx=x_err, sy=[0.1])
+        rd3 = RealData(x, y, sx=x_err, sy=np.full((1, n_data), 0.1))
+        rd4 = RealData(x, y, sx=x_err, covy=[[0.01]])
+        rd5 = RealData(x, y, sx=x_err, covy=np.full((1, 1, n_data), 0.01))
+        for rd in [rd0, rd1, rd2, rd3, rd4, rd5]:
+            odr_obj = ODR(rd, linear_model, beta0=[0.4, 0.4],
+                          delta0=np.full(n_data, -0.1))
+            odr_obj.set_job(fit_type=0, del_init=1)
+            # Just make sure that it runs without raising an exception.
+            odr_obj.run()

scipy/special/_basic.py

@@ -2930,7 +2930,11 @@ def factorial(n, exact=False):
         return _exact_factorialx_array(n)
     # we do not raise for non-integers with exact=True due to
     # historical reasons, though deprecation would be possible
-    return _ufuncs._factorial(n)
+    res = _ufuncs._factorial(n)
+    if isinstance(n, np.ndarray):
+        # _ufuncs._factorial does not maintain 0-dim arrays
+        return np.array(res)
+    return res
 
 
 def factorial2(n, exact=False):

scipy/special/tests/test_basic.py

@@ -2129,6 +2129,8 @@ def _check(n, expected):
     def test_factorial_array_corner_cases(self, content, dim, exact, dtype):
         if dtype == np.int64 and any(np.isnan(x) for x in content):
             pytest.skip("impossible combination")
+        # np.array(x, ndim=0) will not be 0-dim. unless x is too
+        content = content if (dim > 0 or len(content) != 1) else content[0]
         n = np.array(content, ndmin=dim, dtype=dtype)
         result = None
         if not content:
@@ -2151,14 +2153,22 @@ def test_factorial_array_corner_cases(self, content, dim, exact, dtype):
             # no error
             result = special.factorial(n, exact=exact)
 
+        # assert_equal does not distinguish scalars and 0-dim arrays of the same value, see
+        # https://github.com/numpy/numpy/issues/24050
+        def assert_really_equal(x, y):
+            assert type(x) == type(y), f"types not equal: {type(x)}, {type(y)}"
+            assert_equal(x, y)
+
         if result is not None:
             # expected result is empty if and only if n is empty,
             # and has the same dtype & dimension as n
             with suppress_warnings() as sup:
                 sup.filter(DeprecationWarning)
-                r = special.factorial(n.ravel(), exact=exact) if n.size else []
+                # keep 0-dim.; otherwise n.ravel().ndim==1, even if n.ndim==0
+                n_flat = n.ravel() if n.ndim else n
+                r = special.factorial(n_flat, exact=exact) if n.size else []
             expected = np.array(r, ndmin=dim, dtype=dtype)
-            assert_equal(result, expected)
+            assert_really_equal(result, expected)
 
     @pytest.mark.parametrize("exact", [True, False])
     @pytest.mark.parametrize("n", [1, 1.1, 2 + 2j, np.nan, None],
@@ -2211,6 +2221,8 @@ def test_factorial2_int_reference(self, n):
     def test_factorial2_array_corner_cases(self, content, dim, exact, dtype):
         if dtype == np.int64 and any(np.isnan(x) for x in content):
             pytest.skip("impossible combination")
+        # np.array(x, ndim=0) will not be 0-dim. unless x is too
+        content = content if (dim > 0 or len(content) != 1) else content[0]
         n = np.array(content, ndmin=dim, dtype=dtype)
         if np.issubdtype(n.dtype, np.integer) or (not content):
             # no error
@@ -2262,6 +2274,8 @@ def test_factorialk_int_reference(self, n, k):
     def test_factorialk_array_corner_cases(self, content, dim, dtype):
         if dtype == np.int64 and any(np.isnan(x) for x in content):
             pytest.skip("impossible combination")
+        # np.array(x, ndim=0) will not be 0-dim. unless x is too
+        content = content if (dim > 0 or len(content) != 1) else content[0]
         n = np.array(content, ndmin=dim, dtype=dtype)
         if np.issubdtype(n.dtype, np.integer) or (not content):
             # no error; expected result is identical to n

scipy/stats/_unuran/meson.build

@@ -119,7 +119,6 @@ unuran_sources = [
   '../../_lib/unuran/unuran/src/specfunct/cephes_polevl.c',
   '../../_lib/unuran/unuran/src/specfunct/cgamma.c',
   '../../_lib/unuran/unuran/src/specfunct/hypot.c',
-  '../../_lib/unuran/unuran/src/specfunct/log1p.c',
   '../../_lib/unuran/unuran/src/tests/chi2test.c',
   '../../_lib/unuran/unuran/src/tests/correlation.c',
   '../../_lib/unuran/unuran/src/tests/countpdf.c',