Repository: scipy/scipy
Pull request: TST: optimize.linprog/milp: add tests for various bug reports
#### Reference issue
Closes gh-18907  # already fixed by HiGHS 1.6.0, it seems
gh-20116
gh-20336

#### What does this implement/fix?
There were a few open bug reports that will may be closed after upgrading to the most recent HiGHS. In case we want regression tests, here they are. If not (since these are all upstream fixes), consider this a meta-issue of bugs for which we're looking for upstream fixes.

Problem statement:
BUG: `presolve` option in `milp` causing feasible problem to be reported as infeasible. 
### Describe your issue.

When setting `presolve=True` via `options` dictionary in `scipy.optimize.milp` a feasible problem is being reported as infeasible. The exact same problem is "made" feasible by simply swapping to `presolve=False`. 

Alternatively changing `c` also seems to make the problem feasible. Changing `c` should not change whether a problem is feasible or not. 

### Reproducing Code Example

```python
from scipy.optimize import milp
import numpy as np
inf= np.inf

# set up problem
c = np.array([-0.85850509, -0.82892676, -0.80026454, -0.63015535, -0.5099006 , -0.50077193, -0.4894404 , -0.47285865, -0.39867774, -0.38069646, -0.36733012, -0.36733012, -0.35820411, -0.31576141, -0.20626091, -0.12466144, -0.10679516, -0.1061887 , -0.1061887 , -0.1061887 , -0. , -0. , -0. , -0. , 0. , 0. , 0. , 0. ])

constraints=(np.array([[ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 1., 0., 0., 0., -25., -0., -0., -0.], [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., 0., 0., 0., 0., -1., 0., 0., 0., 2., 0., 0., 0.], [ 0., 0., 0., 0., 1., 1., 1., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0., -25., -0., -0.], [ 0., 0., 0., 0., -1., -1., -1., -1., 0., -1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0., 0., -0., -0., -25., -0.], [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., -1., -1., 0., 0., 0., 0., 0., 0., 2., 0.], [ 1., 1., 1., 1., 0., 0., 0., 0., 1., 0., 1., 1., 1., 1., 0., 1., 1., 0., 0., 0., 0., 1., 1., 1., -0., -0., -0., -25.], [ -1., -1., -1., -1., 0., 0., 0., 0., -1., 0., -1., -1., -1., -1., 0., -1., -1., 0., 0., 0., 0., -1., -1., -1., 0., 0., 0., 2.]]),
 np.array([-inf, -inf, -inf, -inf, -inf, -inf, -inf, -inf, -inf]), np.array([100., 0., 0., 0., 0., 0., 0., 0., 0.]))

integrality = 1

bounds=(0, 1)

# solve problem
# attempt 1
r1 = milp(c=c, constraints=constraints, integrality=integrality, bounds=bounds, options={'presolve': True}) # res.message = "The problem is infeasible. (HiGHS Status 8: model_status is Infeasible; primal_status is b'Basic')"

# attempt 2
r2 = milp(c=c, constraints=constraints, integrality=integrality, bounds=bounds, options={'presolve': False}) # res.message = 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'

# update c
c_updated = np.array([-0.85850509, -0.82892676, -0.80026454, -0.63015535, -0.5099006 , -0.50077193, -0.4894404 , -0.47285865, -0.39867774, -0.38069646, -0.36733012, -0.36733012, -0.35820411, -0.31576141, -0.20626091, -0.12466144, -0.10679516, -0.1061887 , -0.1061887 , -0.1061887 , -0.5 , -0. , -0. , -0. , 0. , 0. , 0. , 0. ])

# attempt 3
r3 = milp(c=c_updated, constraints=constraints, integrality=integrality, bounds=bounds, options={'presolve': True}) # res.message = 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'

print(r1)
print(r2)
print(r3)
```


### Error message

```shell
N/A
```


### SciPy/NumPy/Python version and system information

```shell
1.9.0 1.22.4 sys.version_info(major=3, minor=9, micro=5, releaselevel='final', serial=0) lapack_mkl_info: NOT AVAILABLE openblas_lapack_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib'] lapack_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib'] blas_mkl_info: NOT AVAILABLE blis_info: NOT AVAILABLE openblas_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib'] blas_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib']
```

Hints from issue discussion:
@jajhall I confirmed the report. Is this of interest?
Yes, thanks. We know that there are bugs in the MIP presolve - as exposed by various open issues on GitHub. Simple instances like this that expose bugs are very useful. I'll see whether I can reproduce the issue when setting up the model from C++ (allowing me to debug).
Unfortunately, I reproduce the optimal objective function values (suggesting strongly that I've got the right MIP) but not the bug. This is with the latest working version of HiGHS. 

What version of HiGHS is being used if `SciPy/NumPy/Python version and system information` is

`1.9.0 1.22.4 sys.version_info(major=3, minor=9, micro=5, releaselevel='final', serial=0) lapack_mkl_info: NOT AVAILABLE openblas_lapack_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib'] lapack_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib'] blas_mkl_info: NOT AVAILABLE blis_info: NOT AVAILABLE openblas_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib'] blas_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] runtime_library_dirs = ['/usr/local/lib']`
@jajhall SciPy's version of HiGHS is at https://github.com/scipy/HiGHS/commits/4a122958a82e67e725d08153e099efe4dad099a2.
It looks like the closes HiGHS ancestor is https://github.com/scipy/HiGHS/commit/c0b554410c79db42c18c0e96ef45fca0fa56d28d. The only difference seems to be some merge commits.
Thanks. It would be good to verify that this exposes a bug in presolve that is  now fixed
Yes, with `c0b5544` infeasibility is reported if presolve is on, and I recognise the error messages printed by the MIP solver that went away when we fixed the bug.
OK, so we'll want to update HiGHS to fix this. @mckib2 would you be willing to submit a PR?
> OK, so we'll want to update HiGHS to fix this. @mckib2 would you be willing to submit a PR?

@mdhaber This PR is open: gh-18642.  It's waiting for review, it's a fairly recent revision of HiGHS and adds some other long standing improvements.  I'll make sure to update again before it goes in.  I haven't had a lot of time to champion the review cycle, so if you have any time to look at it and give feedback that would be helpful.  All tests were passing the last time time I ran it (I believe the circleci failure was not due to the changes in the PR)
gh-18642 is all the low-level stuff, so if it passes all the higher-level tests, I'm not much help there. 
I'm seeing the same issue with a different MIP problem. The following problem is optimized correctly if presolve is disabled. If presolve is enabled, the problem is incorrectly reported as infeasible.

```
import scipy.optimize

bounds = scipy.optimize.Bounds(lb=0, ub=1)
integrality = [1, 1, 0, 0]
c = [10, 9.52380952, -1000, -952.38095238]
constraints = [
    scipy.optimize.LinearConstraint([1, 1, 0, 0], ub=1),
    scipy.optimize.LinearConstraint([0, 0, 1, 1], ub=1),
    scipy.optimize.LinearConstraint([200, 0, 0, 0], ub=200),
    scipy.optimize.LinearConstraint([0, 200, 0, 0], ub=200),
    scipy.optimize.LinearConstraint([0, 0, 2000, 0], ub=1000),
    scipy.optimize.LinearConstraint([0, 0, 0, 2000], ub=1000),
    scipy.optimize.LinearConstraint([-1, 0, 1, 0], ub=0),
    scipy.optimize.LinearConstraint([-1, -1, 0, 1], ub=0),
]
result = scipy.optimize.milp(
    c=c,
    constraints=constraints,
    bounds=bounds,
    integrality=integrality,
    # Disable presolve as that is causing issues in some of the tests
    options={"presolve": False},
)
```

SciPy/NumPy/Python version and system information:
```
1.13.0 1.23.5 sys.version_info(major=3, minor=11, micro=7, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-ilw0s1wu/overlay/lib/python3.11/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp311-cp311/bin/python
  version: '3.11'
```

Commit message: Merge branch 'main' into test_bugs_linprog_milp

Patch to review:
diff --git a/scipy/optimize/tests/test_linprog.py b/scipy/optimize/tests/test_linprog.py
index cf4612a80d5c..664450e622b5 100644
--- a/scipy/optimize/tests/test_linprog.py
+++ b/scipy/optimize/tests/test_linprog.py
@@ -1863,7 +1863,7 @@ def test_complementary_slackness(self):
         # non-zero RHS
         assert np.allclose(res.ineqlin.marginals @ (b_ub - A_ub @ res.x), 0)
 
-    @pytest.mark.xfail(reason='Upstream / Wrapper issue, see gh-20589')
+
     def test_bug_20336(self):
         """
         Test that `linprog` now solves a poorly-scaled problem
diff --git a/scipy/optimize/tests/test_milp.py b/scipy/optimize/tests/test_milp.py
index a5809095d621..0fa7406f67e9 100644
--- a/scipy/optimize/tests/test_milp.py
+++ b/scipy/optimize/tests/test_milp.py
@@ -387,7 +387,7 @@ def test_mip_rel_gap_passdown():
     # monotonically with the mip_rel_gap parameter.
     assert np.all(np.diff(sol_mip_gaps) < 0)
 
-@pytest.mark.xfail(reason='Upstream / Wrapper issue, see gh-20589')
+
 def test_large_numbers_gh20116():
     h = 10 ** 12
     A = np.array([[100.4534, h], [100.4534, -h]])
@@ -400,9 +400,6 @@ def test_large_numbers_gh20116():
     assert np.all(A @ res.x < b)
 
 
-@pytest.mark.skipif(_IS_32BIT,
-    reason="Fails on 32 bit.",
-)
 def test_presolve_gh18907():
     from scipy.optimize import milp
     import numpy as np


Code review comments:
diff_hunk: @@ -387,7 +387,7 @@ def test_mip_rel_gap_passdown():
     # monotonically with the mip_rel_gap parameter.
     assert np.all(np.diff(sol_mip_gaps) < 0)
 
-@pytest.mark.xfail(reason='Upstream / Wrapper issue, see gh-20589')
+
path: scipy/optimize/tests/test_milp.py
text: @author:
```suggestion
@user1.mark.xfail(reason='Upstream / Wrapper issue, see gh-20116')
```

diff_hunk: @@ -1863,7 +1863,7 @@ def test_complementary_slackness(self):
         # non-zero RHS
         assert np.allclose(res.ineqlin.marginals @ (b_ub - A_ub @ res.x), 0)
 
-    @pytest.mark.xfail(reason='Upstream / Wrapper issue, see gh-20589')
+
path: scipy/optimize/tests/test_linprog.py
text: @author:
```suggestion
    @user1.mark.xfail(reason='Upstream / Wrapper issue, see gh-20589')
```