Repository: scipy/scipy
Pull request: BUG: Wrong status returned by _check_result. See #18106. optimize
#### Reference issue
Closes #18106 

#### What does this implement/fix?
The _check_result function that would parse the results from highs didn't take into account integrality of the variables. The status would say that it encountered a problem when in fact it was correct and it had converged.

#### Additional information
Checks what are the variables that are out of bounds, then if they are 0. If any of them aren't 0 then the solver ran into problems. If they are all 0 then check if they have integrality > 1 (Semi-continuous or  Semi-integer).
Added a couple of unit tests (checks the status) including the exact code mentioned.


Problem statement:
BUG: Linprog reports failure despite success convergence, given a problem involving semi-continuous variables
### Describe your issue.

Given the following problem involving a semi-continuous variable, `linprog(method='highs')` reports a failure despite successful convergence to a solution.

The problem is the following:

```
maximize x
s.t
    x <= 0.5
    x is semi-continuous, bounds=[1;5]
```

Linprog reports the correct solution: `x=[0.]` and `slack=[0.5]`, but the `success` flag is `False` and there is an error message.

### Reproducing Code Example

```python
import numpy as np
from scipy.optimize import linprog

c = np.array([-1.])
A_ub = np.array([[1.]])
b_ub = np.array([0.5])
A_eq = np.zeros((0, 1), dtype=float)
b_eq = np.array([], dtype=float)
bounds = np.array([[1., 5.]])
integrality = np.array([2])

ans = linprog(
    c,
    A_ub,
    b_ub,
    A_eq,
    b_eq,
    bounds,
    integrality=integrality,
    method='highs'
)
print(f'ans.success: {ans.success}')
print(f'ans.x: {ans.x}')
print(f'ans.slack: {ans.slack}')
print(f'ans.message: {ans.message}')
```


### Error message

```shell
ans.success: False
ans.x: [0.]
ans.slack: [0.5]
ans.message: The solution does not satisfy the constraints within the required tolerance of 3.16E-04, yet no errors were raised and there is no certificate of infeasibility or unboundedness. Check whether the slack and constraint residuals are acceptable; if not, consider enabling presolve, adjusting the tolerance option(s), and/or using a different method. Please consider submitting a bug report.
```


### SciPy/NumPy/Python version and system information

```shell
1.10.1 1.24.2 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.2.1",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "0.29.33",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.2.1",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.2.1",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.12.1",
      "include directory": "/tmp/pip-build-env-xob7rtip/overlay/lib/python3.10/site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "linux"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "linux"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "OpenBLAS",
      "found": true,
      "version": "0.3.18",
      "detection method": "cmake",
      "include directory": "unknown",
      "lib directory": "unknown",
      "openblas configuration": "unknown",
      "pc file directory": "unknown"
    },
    "lapack": {
      "name": "OpenBLAS",
      "found": true,
      "version": "0.3.18",
      "detection method": "cmake",
      "include directory": "unknown",
      "lib directory": "unknown",
      "openblas configuration": "unknown",
      "pc file directory": "unknown"
    }
  },
  "Python Information": {
    "path": "/opt/python/cp310-cp310/bin/python",
    "version": "3.10"
  }
}
```

Hints from issue discussion:
I looked into the implementation of `linprog`. The bounds checking in `_check_results` doesn't make sense for semi-continuous variables:

https://github.com/scipy/scipy/blob/4a6f2d1e9a66a9447dd6a69a44be7170f2567b31/scipy/optimize/_linprog_util.py#L1491

Semi-continuous variables can always take value 0 (same for semi-integer variables), regardless of the bounds. At least that's what the documentation of `linprog` is saying for the `integrality` flag:
> 2 : Semi-continuous variable; decision variable must be within bounds or take value 0.
> 3 : Semi-integer variable; decision variable must be an integer within bounds or take value 0.

Thanks. Those checks were probably written before the semi-continuous bounds were added. We'll adjust that.

Commit message: BUG: Wrong status returned by _check_result. See #18106. optimize

Patch to review:
diff --git a/scipy/optimize/_linprog.py b/scipy/optimize/_linprog.py
index a11041f099e9..18b68bda1897 100644
--- a/scipy/optimize/_linprog.py
+++ b/scipy/optimize/_linprog.py
@@ -640,7 +640,7 @@ def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None,
                              **solver_options)
         sol['status'], sol['message'] = (
             _check_result(sol['x'], sol['fun'], sol['status'], sol['slack'],
-                          sol['con'], lp.bounds, tol, sol['message']))
+                          sol['con'], lp.bounds, tol, sol['message'], integrality))
         sol['success'] = sol['status'] == 0
         return OptimizeResult(sol)
 
@@ -692,7 +692,7 @@ def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None,
 
     x, fun, slack, con = _postsolve(x, postsolve_args, complete)
 
-    status, message = _check_result(x, fun, status, slack, con, lp_o.bounds, tol, message)
+    status, message = _check_result(x, fun, status, slack, con, lp_o.bounds, tol, message, integrality)
 
     if disp:
         _display_summary(message, status, fun, iteration)
diff --git a/scipy/optimize/_linprog_util.py b/scipy/optimize/_linprog_util.py
index 717735a59ee9..65ce0106d87d 100644
--- a/scipy/optimize/_linprog_util.py
+++ b/scipy/optimize/_linprog_util.py
@@ -1416,7 +1416,7 @@ def _postsolve(x, postsolve_args, complete=False):
     return x, fun, slack, con
 
 
-def _check_result(x, fun, status, slack, con, bounds, tol, message):
+def _check_result(x, fun, status, slack, con, bounds, tol, message, integrality=None):
     """
     Check the validity of the provided solution.
 
@@ -1489,6 +1489,21 @@ def _check_result(x, fun, status, slack, con, bounds, tol, message):
         is_feasible = False
     else:
         invalid_bounds = (x < bounds[:, 0] - tol).any() or (x > bounds[:, 1] + tol).any()
+
+
+        # in case integrality is 2 or 3, there needs to be an aditional check to see if
+        # the values that are out of bounds are both 0 and with integrality more than 1
+        if np.any(integrality) and invalid_bounds and (integrality > 1).any():
+            # Check if the vals that are out of bounds are all 0 (within tolerance),
+            # if they arent then its done -> invalid_bounds==True. If they are all 0 then
+            # check if those values have integrality > 1
+            out_of_bounds = np.logical_or(x < bounds[:, 0] - tol, x > bounds[:, 1] + tol)
+            if (np.abs(x[out_of_bounds]) < tol).all():
+                invalid_idx = np.where(out_of_bounds)
+                invalid_bounds = not np.isin(invalid_idx,
+                                             np.where(integrality > 1)).all()
+
+
         invalid_slack = status != 3 and (slack < -tol).any()
         invalid_con = status != 3 and (np.abs(con) > tol).any()
         is_feasible = not (invalid_bounds or invalid_slack or invalid_con)
diff --git a/scipy/optimize/tests/test_linprog.py b/scipy/optimize/tests/test_linprog.py
index a0c205c36e33..0bd37a0b3bff 100644
--- a/scipy/optimize/tests/test_linprog.py
+++ b/scipy/optimize/tests/test_linprog.py
@@ -2366,6 +2366,42 @@ def test_mip_rel_gap_passdown(self):
         assert not np.all(gap_diffs == 0)
 
 
+    def test_semi_continuous(self):
+        # See issue #18106. This test if the integrality parameter
+        # is being checked correctly (status returns 0),
+        # when its > 1 it allows values to be 0 even if its out of bounds.
+
+        c = np.array([-1.])
+        A_ub = np.array([[1.]])
+        b_ub = np.array([0.5])
+        A_eq = np.zeros((0, 1), dtype=float)
+        b_eq = np.array([], dtype=float)
+        bounds = np.array([[1., 5.]])
+        integrality = np.array([2])
+
+        res = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds,
+                      integrality=integrality, method='highs')
+
+        assert np.all(res.x == 0)
+        assert res.status == 0
+
+        c = np.array([-1., -1., 1])
+        A_ub = np.array([[1., 1., 0], [1.,1., 0], [0,0,1]])
+        b_ub = np.array([0.5, 0.5, 1.])
+        A_eq = np.zeros((0, 3), dtype=float)
+        b_eq = np.array([], dtype=float)
+        bounds = np.array([[1., 5.], [1., 5.],  [0.5, 5]])
+        integrality = np.array([2, 2, 0])
+
+        res = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds,
+                      integrality=integrality, method='highs')
+
+        np.testing.assert_allclose(res.x, [0, 0, 0.5])
+        assert res.status == 0
+
+
+
+
 ###########################
 # Autoscale-Specific Tests#
 ###########################


Code review comments:
diff_hunk: @@ -1489,6 +1489,21 @@ def _check_result(x, fun, status, slack, con, bounds, tol, message):
         is_feasible = False
     else:
         invalid_bounds = (x < bounds[:, 0] - tol).any() or (x > bounds[:, 1] + tol).any()
+
path: scipy/optimize/_linprog_util.py
text: @user1:
```suggestion
```
PEP8 doesn't allow for multiple blank lines within a function.