SciPy 1.13 breaks something with sparse matrices

[moorepants]

============================= test session starts ==============================
platform linux -- Python 3.10.14, pytest-8.2.2, pluggy-1.5.0
rootdir: /home/runner/work/cyipopt/cyipopt
configfile: pytest.ini
testpaths: cyipopt/tests
collected 107 items

cyipopt/tests/integration/test_hs071.py ......                           [  5%]
cyipopt/tests/unit/test_deprecations.py .....                            [ 10%]
cyipopt/tests/unit/test_deriv_errors.py ..........                       [ 19%]
cyipopt/tests/unit/test_ipopt_funcs.py ....ss....                        [ 28%]
cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py ...................... [ 49%]
....................X...                                                 [ 71%]
cyipopt/tests/unit/test_scipy_optional.py s.......................FF.... [100%]

=================================== FAILURES ===================================
___________________ test_minimize_ipopt_sparse_jac_if_scipy ____________________

    @pytest.mark.skipif("scipy" not in sys.modules,
                        reason="Test only valid of Scipy available")
    def test_minimize_ipopt_sparse_jac_if_scipy():
        """ `minimize_ipopt` works with objective gradient, and sparse
            constraint jacobian. Solves
            Hock & Schittkowski's test problem 71:

            min x0*x3*(x0+x1+x2)+x2
            s.t. x0**2 + x1**2 + x2**2 + x3**2 - 40  = 0
cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py::TestTrustRegionConstr::test_list_of_problems[prob2]
cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py::TestTrustRegionConstr::test_list_of_problems[prob7]
cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py::TestTrustRegionConstr::test_args
cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py::TestEmptyConstraint::test_empty_constraint
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:673: OptimizeWarning: Constraint options `finite_diff_jac_sparsity`, `finite_diff_rel_step`, `keep_feasible`, and `hess`are ignored by this method.
    constraints = optimize._minimize.standardize_constraints(constraints, x0,

cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py::TestTrustRegionConstr::test_list_of_problems[prob6]
cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py::TestTrustRegionConstr::test_list_of_problems[prob7]
  /usr/share/miniconda3/envs/test-environment/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:231: UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.
    self.H.update(self.x - self.x_prev, self.g - self.g_prev)

cyipopt/tests/unit/test_scipy_ipopt_from_scipy.py: 25 warnings
cyipopt/tests/unit/test_scipy_optional.py: 2 warnings
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:615: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
    x, info = nlp.solve(x0)

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_input_validation
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method a newt cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[cg]
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method cg cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[bfgs]
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method bfgs cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[newton-cg]
cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[newton-cg]
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method newton-cg cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[trust-constr]
  /usr/share/miniconda3/envs/test-environment/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:504: UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.
    self.H.update(delta_x, delta_g)

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[dogleg]
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method dogleg cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[trust-ncg]
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method trust-ncg cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[trust-exact]
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method trust-exact cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_jac_with_scipy_methods[trust-krylov]
  /home/runner/work/cyipopt/cyipopt/cyipopt/scipy_interface.py:556: RuntimeWarning: Method trust-krylov cannot handle bounds.
    res = optimize.minimize(fun, x0, args, method, jac, hess, hessp,

-- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html
=========================== short test summary info ============================
FAILED cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_sparse_jac_if_scipy - IndexError: tuple index out of range
FAILED cyipopt/tests/unit/test_scipy_optional.py::test_minimize_ipopt_sparse_and_dense_jac_if_scipy - IndexError: tuple index out of range
======= 2 failed, 101 passed, 3 skipped, 1 xpassed, 46 warnings in 5.10s =======

[moorepants]

A call to scipy.sparse.vstack() fails:

> /home/moorepants/src/cyipopt/cyipopt/scipy_interface.py(288)_get_sparse_jacobian_structure()
    286             jacobians.append(coo_array(np.ones((con_val.size, x0.size))))
    287             con_jac_is_sparse.append(False)
--> 288     J = scipy.sparse.vstack(jacobians)
    289     return con_jac_is_sparse, J.row, J.col
    290 

ipdb> jacobians
[<4 sparse array of type '<class 'numpy.float64'>'
    with 4 stored elements in COOrdinate format>, <4 sparse array of type '<class 'numpy.float64'>'
    with 4 stored elements in COOrdinate format>]

[moorepants]

In SciPy 1.12 it works:

ipdb> jacobians
[<1x4 sparse array of type '<class 'numpy.float64'>'
    with 4 stored elements in COOrdinate format>, <1x4 sparse array of type '<class 'numpy.float64'>'
    with 4 stored elements in COOrdinate format>]
ipdb> scipy.sparse.vstack(jacobians)
<2x4 sparse array of type '<class 'numpy.float64'>'
    with 8 stored elements in COOrdinate format>

[moorepants]

The list of two sparse column vectors is:

ipdb> jacobians[0].toarray()
array([[ 2., 10., 10.,  2.]])
ipdb> jacobians[1].toarray()
array([[25.,  5.,  5., 25.]])

[moorepants]

In SciPy 1.13.1 the shape of the jacobians is different:

ipdb> jacobians[0].toarray()
array([ 2., 10., 10.,  2.])
ipdb> jacobians[1].toarray()
array([25.,  5.,  5., 25.])
ipdb> jacobians[0].shape
(4,)
ipdb> jacobians[1].shape
(4,)

[moorepants]

Maybe the calls to coo_array() are returning different dimensions than before:

def _get_sparse_jacobian_structure(constraints, x0):
    con_jac_is_sparse = []
    jacobians = []
    x0 = np.asarray(x0)
    if isinstance(constraints, dict):
        constraints = (constraints, )
    if len(constraints) == 0:
        return [], [], []
    for con in constraints:
        con_jac = con.get('jac', False)
        if con_jac:
            if isinstance(con_jac, bool):
                _, jac_val = con['fun'](x0, *con.get('args', []),
                                        **con.get('kwargs', {}))
            else:
                jac_val = con_jac(x0, *con.get('args', []),
                                  **con.get('kwargs', {}))
            # check if dense or sparse
            if isinstance(jac_val, coo_array):
                jacobians.append(jac_val)
                con_jac_is_sparse.append(True)
            else:
                # Creating the coo_array from jac_val would yield to
                # wrong dimensions if some values in jac_val are zero,
                # so we assume all values in jac_val are nonzero
                jacobians.append(coo_array(np.ones_like(np.atleast_2d(jac_val))))
                con_jac_is_sparse.append(False)
        else:
            # we approximate this jacobian later (=dense)
            con_val = np.atleast_1d(con['fun'](x0, *con.get('args', []),
                                               **con.get('kwargs', {})))
            jacobians.append(coo_array(np.ones((con_val.size, x0.size))))
            con_jac_is_sparse.append(False)
    J = scipy.sparse.vstack(jacobians)
    return con_jac_is_sparse, J.row, J.col

[moorepants]

This seems to be the issue: https://github.com/scipy/scipy/issues/20415

[moorepants]

Fixed in #257 , the tests were simply updated to work with scipy 1.13 when calling coo_array() such that 2d arrays are passed in. No actual issue in the internal handling in cyipopt but users will have to give proper input (which changed without warning in SciPy 1.13).