Repository: scipy/scipy
Pull request: ENH: linprog HiGHS marginals/sensitivy analysis
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->

- Closes scipygh-8836
- Closes gh-11848

#### What does this implement/fix?
<!--Please explain your changes.-->

- For HiGHS linprog methods (`highs-ipm`, `highs-ds`)
- Returns `ineqlin`, `eqlin`, `lower`, and `upper` results fields containing marginals and solution information with the following fields:
    - lower
        - marginals
        - freedom (`x - lb`)
    - upper
        - marginals
        - freedom (`ub - x`)
    - eqlin
        - marginals
        - slack (`b_ub - A_ub @ x`)
    - ineqlin
        - marginals
        - violation (`b_eq - A_eq @ x`)

#### Additional information
<!--Any additional information you think is important.-->

This PR implements long-desired functionality from the SciPy LP solver to provide some sensitivity analysis in the form of marginals, also commonly referred to as "Lagrange multipliers", "dual values", or "shadow prices".  The name `marginals`, `freedom`, `slack`, and `violation` were settled on after discussion with HiGHS maintainers and consideration for clarity and future sensitivity analysis extensions.

See original draft PR https://github.com/mckib2/scipy/pull/10 for original work with @mdhaber and discussion as well as additional follow-on development for ranging capabilities.

Problem statement:
scipy.optimize.linprog(method='simplex') needs to return duals
I have a huge linprog problem of almost 1k variables and restrictions. I can calculate the solution with scipy.optimize.linprog(method='simplex') but I need shadow prices (or opportunity costs / duals) of ~100 inequalities.

I'm able to calculate them by adding 1 to the right side of the inequality and then solving that problem. Then I get the shadow price substracting the objective functions values for both solutions: shadow_price_i = f_max_original - f_max_i. Then repeat 100 times. This method works but it's painfully slow (1h). Note: I could also pose the dual problem and solve it 

It's a shame that scipy does not return duals, so I'm opening a feature request. I'm new contributing to open projects, but I'm a CS graduate and have experience with both python and the simplex method, so I may be able to do so with some guidance.




How to get Lagrange / lambda multipliers out of 'linprog' optimize subroutine in scipy module ?
It is required to obtain Lagrangian multipliers while optimizing a linear programming problem.  Kindly help.

Best regards,
Bharat

Hints from issue discussion:
@j6e I'd be happy to help with this. The simplex method could use some maintenance, and it sounds like you have the background for it.

BTW I wrote the `interior-point` method but didn't return duals because the solution is, in general, on an optimal face - not necessarily at a vertex. The duals are not useful in this case for shadow prices, right? And  I haven't implemented a crossover routine to find an optimal basis - or implemented the revised simplex method, for that matter - because SciPy's primary sparse linear programming solver, SuperLU, doesn't support rank-one updates and can't easily be hacked to do so (#7700), so it would not be efficient.

We could start by getting your development environment set up. What's your OS?
> We could start by getting your development environment set up. What's your OS?

I have w10 at home, but I usually use the ubuntu 14.04 image for python development.

>BTW I wrote the interior-point method but didn't return duals because the solution is, in general, on an optimal face - not necessarily at a vertex. The duals are not useful in this case for shadow prices, right?

I was aware that the IPM did find it's solutions on an optimal plane, but I don't get why the duals are not usefull for shadow prices. I'm by no means an expert in the matter so who knows.

Since you're the writer of the IPM, maybe you want to know that it categorices my problem as infeasible (simplex solves it). I'm tinkering a bit, because I think it's some sort of numerical inestability (simplex got me troubles with this also). I'll open a new issue if I can't work it out.
My understanding is that at an optimal vertex you have a lot of active inequality constraints, and so then the duals tell you how modifying the constraint would change the optimal objective. On an optimal face, not as many equality constraints are active, so changing constraints doesn't change the objective, and thus the duals aren't as useful. I, too, may be totally wrong.

Are you able to send your problem? I definitely want to know why `interior-point` is failing on your problem.

- Have you confirmed that all the constraints are indeed satisfied by the solution returned by `simplex`? Sometimes the results it gives don't satisfy the constraints (e.g. #8174).
- What is the message given by `interior-point`? Did it say it was infeasible in presolve or did the algorithm run to completion?
[Here's](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.5719&rep=rep1&type=pdf
) some information about sensitivity analysis from interior-point solutions.
But it does not appear to be as simple as when you have a basic optimal solution.
> Have you confirmed that all the constraints are indeed satisfied by the solution returned by simplex? Sometimes the results it gives don't satisfy the constraints (e.g. #8174).

It was the case, that instance of the problem was in fact infeasible.

> My understanding is that at an optimal vertex you have a lot of active inequality constraints, and so then the duals tell you how modifying the constraint would change the optimal objective. On an optimal face, not as many equality constraints are active, so changing constraints doesn't change the objective, and thus the duals aren't as useful. I, too, may be totally wrong.

That being true, the equality constraints must hold, and therefore their shadow prices are relevant (In the report you mentioned there is an example of this on page 33). I checked on my own problem, and while this does not prove anything, the extra slack given by IPM in the zero slack inequalities of simplex is in the order of 0.01 (literally a cent in a problem in the range of millions) and their dual prices vary in the range of 0.1. 

I read somewhere that using the solution of IPM as a starting point of simplex makes it converge quickly to an optimal vertex, maybe it's what we are looking for using returning dual values both in IPM/Simplex... (or not). But I can't recall an LP optimizer that did not return duals, didn't matter if it was IPM or simplex. 
@j6e 
> I read somewhere that using the solution of IPM as a starting point of simplex makes it converge quickly to an optimal vertex, maybe it's what we are looking for using returning dual values both in IPM/Simplex... (or not).

There are "crossover" routines for getting from the interior point solution to an optimal vertex, but I didn't implement one because SciPy's primary sparse linear programming solver, SuperLU, doesn't support rank-one updates and can't easily be hacked to do so (#7700), so I couldn't implement it efficiently.

>But I can't recall an LP optimizer that did not return duals, didn't matter if it was IPM or simplex.

And I wanted to return the duals, but I didn't consider it to be a primary objective, so it didn't get done in the original push.  The reason it was not trivial is that before the interior point method sees the problem there is preprocessing (to convert the problem to standard form, eliminate redundant equality constraints, tighten bounds, eliminate fixed variables, etc...) and so the duals calculated by the algorithm need postprocessing to get the duals for the problem as originally posed. I would need help untangling all this. If you think this is something you could help with, I'd definitely be interested in working with you on it!
@mdhaber I will be thrilled to help!
@j6e Great! I think the first thing to do is to familiarize yourself (at a high level) with the code and the algorithm available in [this paper](https://link.springer.com/chapter/10.1007/978-1-4757-3216-0_8).  I tried to align the code as clearly as possible with the equations in the paper, but one notable difference is the use of the variable `z` in the code for what is _s_ in the paper.

What I need help with first is understanding how the most fundamental preprocessing step - converting the problem to standard form - affects the duals we're after. Specifically, you see in the `_get_Abc` function of `_linprog_ip.py`:

```
    Given a linear programming problem of the form:

    minimize:     c^T * x

    subject to:   A_ub * x <= b_ub
                  A_eq * x == b_eq
                  bounds[i][0] < x_i < bounds[i][1]

    return the problem in standard form:
    minimize:     c'^T * x'

    subject to:   A * x' == b
                  0 < x' < oo

    by adding slack variables and making variable substitutions as necessary.
```

At the end of the function `_ip_hsd` we have our final `x_hat` which is the solution to the standard form problem; `_postprocess` uses this `x` to get the solution to the original primary problem. How do we need to postprocess the `y` and `z` from the standard form problem to get the duals for the original problem? (Initially, let's ignore the fact that `_presolve` might have eliminated redundant constraints, etc...)

On a side note, you might also take a look at reference [4] from the MOSEK paper, which is about computing the optimal basis. Note that the LU factorization needs to be updated as the basis changes, but I haven't been able to hack SuperLU to update LU factorizations due to #7700, and I think computing the LU from scratch every time would make the procedure too slow to be worth the time coding.


Hi @j6e just thought I'd check in and see if you had any questions about the code?
> There are "crossover" routines for getting from the interior point solution to an optimal vertex, but I didn't implement one because SciPy's primary sparse linear programming solver, SuperLU, doesn't support rank-one updates and can't easily be hacked to do so (#7700), so I couldn't implement it efficiently.

Ref #11759 : HiGHS does have a crossover routine
Hi @bharat29, thanks for making the request. This is a duplicate of gh-8836, though, so I'm closing. Hopefully we'll add this with gh-11759.
Sure.  Thank you 😊

Commit message: TST: add complementary slackness test for HiGHS marginals

Patch to review:
diff --git a/scipy/optimize/_highs/cython/src/HConst.pxd b/scipy/optimize/_highs/cython/src/HConst.pxd
index 6cdad205f705..316c208a1421 100644
--- a/scipy/optimize/_highs/cython/src/HConst.pxd
+++ b/scipy/optimize/_highs/cython/src/HConst.pxd
@@ -35,14 +35,14 @@ cdef extern from "HConst.h" nogil:
 
 
     cdef enum HighsBasisStatus:
-        LOWER "HighsBasisStatus::LOWER" = 0, # (slack) variable is at its lower bound [including fixed variables]
-        BASIC "HighsBasisStatus::BASIC" # (slack) variable is basic
-        UPPER "HighsBasisStatus::UPPER" # (slack) variable is at its upper bound
-        ZERO "HighsBasisStatus::ZERO" # free variable is non-basic and set to zero
-        NONBASIC "HighsBasisStatus::NONBASIC" # nonbasic with no specific bound information - useful for users and postsolve
-        SUPER "HighsBasisStatus::SUPER" # Super-basic variable: non-basic and either free and
-                                        # nonzero or not at a bound. No SCIP equivalent
-
+        HighsBasisStatusLOWER "HighsBasisStatus::LOWER" = 0, # (slack) variable is at its lower bound [including fixed variables]
+        HighsBasisStatusBASIC "HighsBasisStatus::BASIC" # (slack) variable is basic
+        HighsBasisStatusUPPER "HighsBasisStatus::UPPER" # (slack) variable is at its upper bound
+        HighsBasisStatusZERO "HighsBasisStatus::ZERO" # free variable is non-basic and set to zero
+        HighsBasisStatusNONBASIC "HighsBasisStatus::NONBASIC" # nonbasic with no specific bound information - useful for users and postsolve
+        HighsBasisStatusSUPER "HighsBasisStatus::SUPER" # Super-basic variable: non-basic and either free and
+                                                        # nonzero or not at a bound. No SCIP equivalent
+                                        
     cdef enum SolverOption:
         SOLVER_OPTION_SIMPLEX "SolverOption::SOLVER_OPTION_SIMPLEX" = -1
         SOLVER_OPTION_CHOOSE "SolverOption::SOLVER_OPTION_CHOOSE"
diff --git a/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx b/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx
index 0436474bc86a..579f7e057ac1 100644
--- a/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx
+++ b/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx
@@ -1,7 +1,8 @@
 # distutils: language=c++
 # cython: language_level=3
 
-
+import numpy as np
+cimport numpy as np
 from scipy.optimize import OptimizeWarning
 from warnings import warn
 import numbers
@@ -35,6 +36,9 @@ from .HConst cimport (
     HighsOptionTypeINT,
     HighsOptionTypeDOUBLE,
     HighsOptionTypeSTRING,
+
+    HighsBasisStatusLOWER,
+    HighsBasisStatusUPPER,
 )
 from .Highs cimport Highs
 from .HighsStatus cimport (
@@ -47,6 +51,7 @@ from .HighsStatus cimport (
 from .HighsLp cimport (
     HighsLp,
     HighsSolution,
+    HighsBasis,
 )
 from .HighsInfo cimport HighsInfo
 from .HighsOptions cimport (
@@ -661,6 +666,9 @@ def _highs_wrapper(
     # We might need an info object if we can look up the solution and a place to put solution
     cdef HighsInfo info = highs.getHighsInfo() # it should always be safe to get the info object
     cdef HighsSolution solution
+    cdef HighsBasis basis
+    cdef double[:, ::1] marg_bnds = np.zeros((2, numcol))  # marg_bnds[0, :]: lower
+                                                           # marg_bnds[1, :]: upper
 
     # If the status is bad, don't look up the solution
     if model_status != HighsModelStatusOPTIMAL:
@@ -677,6 +685,15 @@ def _highs_wrapper(
     else:
         # Should be safe to read the solution:
         solution = highs.getSolution()
+        basis = highs.getBasis()
+
+        # lagrangians for bounds based on column statuses
+        for ii in range(numcol):
+            if HighsBasisStatusLOWER == basis.col_status[ii]:
+                marg_bnds[0, ii] = solution.col_dual[ii]
+            elif HighsBasisStatusUPPER == basis.col_status[ii]:
+                marg_bnds[1, ii] = solution.col_dual[ii]
+
         return {
             'status': <int> model_status,
             'message': highs.highsModelStatusToString(model_status).decode(),
@@ -689,6 +706,11 @@ def _highs_wrapper(
             # Note: this is for all constraints (A_ub and A_eq)
             'slack': [rhs[ii] - solution.row_value[ii] for ii in range(numrow)],
 
+            # slacks in HiGHS appear as Ax - s, not Ax + s, so lambda is negated;
+            # lambda are the lagrange multipliers associated with Ax=b
+            'lambda': [-1*solution.row_dual[ii] for ii in range(numrow)],
+            'marg_bnds': marg_bnds,
+
             'fun': info.objective_function_value,
             'simplex_nit': info.simplex_iteration_count,
             'ipm_nit': info.ipm_iteration_count,
diff --git a/scipy/optimize/_linprog_highs.py b/scipy/optimize/_linprog_highs.py
index f416093d4e79..aded503f5168 100644
--- a/scipy/optimize/_linprog_highs.py
+++ b/scipy/optimize/_linprog_highs.py
@@ -15,7 +15,7 @@
 
 import inspect
 import numpy as np
-from .optimize import _check_unknown_options, OptimizeWarning
+from .optimize import _check_unknown_options, OptimizeWarning, OptimizeResult
 from warnings import warn
 from ._highs._highs_wrapper import _highs_wrapper
 from ._highs._highs_constants import (
@@ -195,6 +195,25 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True,
                 The number of primal/dual pushes performed during the
                 crossover routine for ``solver='ipm'``.  This is ``0``
                 for ``solver='simplex'``.
+            marginals : OptimizeResult
+                The sensitivity (partial derivative) of the objective function
+                with respect to the right-hand side of each constraint.
+                A dictionary consisting of the fields:
+
+                ineqlin : np.ndarray
+                    The sensitivity with respect to the right-hand side
+                    of the inequality constraints, `b_ub`.
+
+                eqlin : np.ndarray
+                    The sensitivity with respect to the right-hand side
+                    of the equality constraints, `b_eq`.
+
+                lower, upper : np.ndarray
+                    The sensitivity with respect to the lower and upper
+                    `bounds`.
+
+                These partial derivatives are also referred to as
+                "Lagrange multipliers", "dual values", and "shadow prices".
 
     References
     ----------
@@ -334,6 +353,17 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True,
     else:
         slack, con = None, None
 
+    # lagrange multipliers for equalities/inequalities and upper/lower bounds
+    if 'lambda' in res:
+        lamda = res['lambda']
+        marg_ineqlin = np.array(lamda[:len(b_ub)])
+        marg_eqlin = np.array(lamda[len(b_ub):])
+        marg_upper = res['marg_bnds'][1, :]
+        marg_lower = res['marg_bnds'][0, :]
+    else:
+        marg_ineqlin, marg_eqlin = None, None
+        marg_upper, marg_lower = None, None
+
     # this needs to be updated if we start choosing the solver intelligently
     solvers = {"ipm": "highs-ipm", "simplex": "highs-ds", None: "highs-ds"}
 
@@ -352,6 +382,12 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True,
 
     sol = {'x': np.array(res['x']) if 'x' in res else None,
            'slack': slack,
+           'marginals': OptimizeResult({
+               'ineqlin': marg_ineqlin,
+               'eqlin': marg_eqlin,
+               'upper': marg_upper,
+               'lower': marg_lower,
+           }),
            'fun': res.get('fun'),
            'con': con,
            'status': status,
diff --git a/scipy/optimize/tests/test_linprog.py b/scipy/optimize/tests/test_linprog.py
index acb0c71e8b6f..112c76cfd4c3 100644
--- a/scipy/optimize/tests/test_linprog.py
+++ b/scipy/optimize/tests/test_linprog.py
@@ -8,6 +8,7 @@
                            assert_array_less, assert_warns, suppress_warnings)
 from pytest import raises as assert_raises
 from scipy.optimize import linprog, OptimizeWarning
+from scipy.optimize._numdiff import approx_derivative
 from scipy.sparse.linalg import MatrixRankWarning
 from scipy.linalg import LinAlgWarning
 import scipy.sparse
@@ -1675,6 +1676,68 @@ def test_crossover(self):
         # there should be nonzero crossover iterations for IPM (only)
         assert_equal(res.crossover_nit == 0, self.method != "highs-ipm")
 
+    def test_marginals(self):
+        # Ensure lagrange multipliers are correct by comparing the derivative
+        # w.r.t. b_ub/b_eq/ub/lb to the reported duals.
+        c, A_ub, b_ub, A_eq, b_eq, bounds = very_random_gen(seed=0)
+        res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
+                      bounds=bounds, method=self.method, options=self.options)
+        lb, ub = bounds.T
+
+        # sensitivity w.r.t. b_ub
+        def f_bub(x):
+            return linprog(c, A_ub, x, A_eq, b_eq, bounds,
+                           method=self.method).fun
+
+        dfdbub = approx_derivative(f_bub, b_ub, method='3-point', f0=res.fun)
+        assert_allclose(res.marginals.ineqlin, dfdbub)
+
+        # sensitivity w.r.t. b_eq
+        def f_beq(x):
+            return linprog(c, A_ub, b_ub, A_eq, x, bounds,
+                           method=self.method).fun
+
+        dfdbeq = approx_derivative(f_beq, b_eq, method='3-point', f0=res.fun)
+        assert_allclose(res.marginals.eqlin, dfdbeq)
+
+        # sensitivity w.r.t. lb
+        def f_lb(x):
+            bounds = np.array([x, ub]).T
+            return linprog(c, A_ub, b_ub, A_eq, b_eq, bounds,
+                           method=self.method).fun
+
+        with np.errstate(invalid='ignore'):
+            # approx_derivative has trouble where lb is infinite
+            dfdlb = approx_derivative(f_lb, lb, method='3-point', f0=res.fun)
+            dfdlb[~np.isfinite(lb)] = 0
+
+        assert_allclose(res.marginals.lower, dfdlb)
+
+        # sensitivity w.r.t. ub
+        def f_ub(x):
+            bounds = np.array([lb, x]).T
+            return linprog(c, A_ub, b_ub, A_eq, b_eq, bounds,
+                           method=self.method).fun
+
+        with np.errstate(invalid='ignore'):
+            dfdub = approx_derivative(f_ub, ub, method='3-point', f0=res.fun)
+            dfdub[~np.isfinite(ub)] = 0
+
+        assert_allclose(res.marginals.upper, dfdub)
+
+    def test_complementary_slackness(self):
+        # Ensure that the complementary slackness condition is satisfied.
+        c, A_ub, b_ub, A_eq, b_eq, bounds = very_random_gen(seed=42)
+        res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
+                      bounds=bounds, method=self.method, options=self.options)
+
+        resid = (-c + A_ub.T @ res.marginals.ineqlin +
+                 A_eq.T @ res.marginals.eqlin +
+                 res.marginals.upper +
+                 res.marginals.lower)
+        # assert_allclose has trouble at 0, so add 1
+        assert_allclose(resid + 1, 1)
+
 
 ################################
 # Simplex Option-Specific Tests#


Code review comments:
diff_hunk: @@ -1675,6 +1676,68 @@ def test_crossover(self):
         # there should be nonzero crossover iterations for IPM (only)
         assert_equal(res.crossover_nit == 0, self.method != "highs-ipm")
 
+    def test_marginals(self):
+        # Ensure lagrange multipliers are correct by comparing the derivative
+        # w.r.t. b_ub/b_eq/ub/lb to the reported duals.
+        c, A_ub, b_ub, A_eq, b_eq, bounds = very_random_gen(seed=0)
+        res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
+                      bounds=bounds, method=self.method, options=self.options)
+        lb, ub = bounds.T
+
+        # sensitivity w.r.t. b_ub
+        def f_bub(x):
+            return linprog(c, A_ub, x, A_eq, b_eq, bounds,
+                           method=self.method).fun
+
+        dfdbub = approx_derivative(f_bub, b_ub, method='3-point', f0=res.fun)
+        assert_allclose(res.marginals.ineqlin, dfdbub)
+
+        # sensitivity w.r.t. b_eq
+        def f_beq(x):
+            return linprog(c, A_ub, b_ub, A_eq, x, bounds,
+                           method=self.method).fun
+
+        dfdbeq = approx_derivative(f_beq, b_eq, method='3-point', f0=res.fun)
+        assert_allclose(res.marginals.eqlin, dfdbeq)
+
+        # sensitivity w.r.t. lb
+        def f_lb(x):
+            bounds = np.array([x, ub]).T
+            return linprog(c, A_ub, b_ub, A_eq, b_eq, bounds,
+                           method=self.method).fun
+
+        with np.errstate(invalid='ignore'):
+            # approx_derivative has trouble where lb is infinite
+            dfdlb = approx_derivative(f_lb, lb, method='3-point', f0=res.fun)
+            dfdlb[~np.isfinite(lb)] = 0
+
+        assert_allclose(res.marginals.lower, dfdlb)
+
+        # sensitivity w.r.t. ub
+        def f_ub(x):
+            bounds = np.array([lb, x]).T
+            return linprog(c, A_ub, b_ub, A_eq, b_eq, bounds,
+                           method=self.method).fun
+
+        with np.errstate(invalid='ignore'):
+            dfdub = approx_derivative(f_ub, ub, method='3-point', f0=res.fun)
+            dfdub[~np.isfinite(ub)] = 0
+
+        assert_allclose(res.marginals.upper, dfdub)
+
+    def test_complementary_slackness(self):
+        # Ensure that the complementary slackness condition is satisfied.
+        c, A_ub, b_ub, A_eq, b_eq, bounds = very_random_gen(seed=42)
+        res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
+                      bounds=bounds, method=self.method, options=self.options)
+
+        resid = (-c + A_ub.T @ res.marginals.ineqlin +
+                 A_eq.T @ res.marginals.eqlin +
+                 res.marginals.upper +
+                 res.marginals.lower)
+        # assert_allclose has trouble at 0, so add 1
+        assert_allclose(resid + 1, 1)
path: scipy/optimize/tests/test_linprog.py
text: @user1:
Set `atol` instead of adding 1.

@author:
Yep, that worked!  Removed `+1` and added `atol=1e-12`

diff_hunk: @@ -195,6 +195,25 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True,
                 The number of primal/dual pushes performed during the
                 crossover routine for ``solver='ipm'``.  This is ``0``
                 for ``solver='simplex'``.
+            marginals : OptimizeResult
path: scipy/optimize/_linprog_highs.py
text: @user1:
This needs to be moved/copied over to `_linprog_doc.py` to appear in the documentation. 

@author:
Copied over to `_linprog_doc.py`

@user1:
This went into the "Options" of  `method='highs'` only. It should be in"Returns" of  `method='highs'`, `method='highs-ipm'`, and `method='simplex'`. 

@author:
Yep, I did this incorrectly, thanks for seeing that.  I'll have time to fix this and review your comment below in the next few days

@author:
Latest push moves the entry from Options -> Results in `method='highs'` and adds the same entry to `method='highs-ipm` and `method='highs-ds`