BUG: optimize: In _linprog_simplex, handle the case where the status returned by _solve_simplex is not 0 or 2.

[perimosocordiae]

The local variable message wasn't being set in all cases, causing runtime errors.

[argriffing]

I think this change looks good.

[coveralls]

Coverage remained the same when pulling 41394a6 on perimosocordiae:patch-6 into 11d9bd2 on scipy:master.

[argriffing]

A new test would be helpful.

[WarrenWeckesser]

Also, put a bit more context in the commit message, perhaps something like "BUG: optimize: In _linprog_simplex, handle the case where the status returned by _solve_simplex is not 0 and not 2."

[perimosocordiae]

I'm having trouble coming up with a minimal test case. The bug occurred when status = 3, which indicates that _solve_simplex was able to find a pivot column, but not a pivot row. This happened to coincide with the objective function being zero (abs(T[-1, -1]) < tol), which meant that the status was not reset to 2.

Unfortunately, I don't know the algorithm well enough to work backwards to input data that generate this behavior.

[argriffing]

Merged, thanks! I'll add a github issue in case anyone wants to try adding the test.

[argriffing]

I might have merged this too early. If phase 1 attempts to find a feasible solution, then shouldn't an error in this phase report that message and status, even if the attempt to find a feasible solution happens to use linear programming which failed for some other reason (e.g. unboundedness)? I agree that it's not good to have the message and status out of sync but maybe it should go the other way, in other words both the status and message should perhaps report feasibility failure rather than both reporting unboundedness, if the linear programming error occurred in phase 1?

[perimosocordiae]

That's fair. In that case, a fix might look like:

if status == 0 and abs(T[-1, -1]) < tol:
    # Remove the pseudo-objective row from the tableau
    # ...etc...
else:
    # Failure to find a feasible starting point
    status = 2

[pv]

linprog is not in 0.14.0, so the API can be changed freely.

[argriffing]

@perimosocordiae on the other hand, is it actually detecting that the original problem is unbounded during phase 1 when it is nominally looking for a feasible solution? If so, then I think the PR works. But if it is solving a different linear programming problem in the process of looking for a feasible solution to the original problem, then if this search fails because the different linear program was unbounded, then I think the user would want to see a message reflecting failure to find a feasible solution to the original problem rather than the details of the mode of failure of that search. If I understood the algorithm better then I'd know which it is.

[perimosocordiae]

I agree entirely. Maybe we should consult @robfalck, because I don't know which makes more sense either.

[argriffing]

I checked the test suite, and the case where the status is something other than 0 or 2 at the end of phase 1 is not covered by any test.

[robfalck]

Good find. The code to better handle non-convergence of Phase 1 needs to be cleaned up. I can do that.

On a similar note, I've found a case where an artificial variable remains basic at the end of Phase 1 (which is known to happen sometimes). I need to add the ability to handle the situation, since the failure to do so leads to incorrect results. This may be related to the issue found by @perimosocordiae.

[robfalck]

Also, when I wrote this code, I was under the assumption that phase one could not fail due to unboundedness (only due to infeasibility). Phase 1 is run if the origin is not a basic feasible solution. If the origin is not a basic feasible solution, then I don't see how the problem could be unbounded. If you have a counter-example, I'd be more than happy to look at this more.

[argriffing]

I'm having trouble coming up with a minimal test case.

Maybe someone at euroscipy sprints could try brutely sampling random linear programs to look for an example? Of course this will fail if examples do not occur with positive probability..

[perimosocordiae]

@robfalck The bug did present when I used (-inf, inf) bounds, which I believe adds an artificial variable. That might be a good starting point.

[argriffing]

I haven't been able to construct a counterexample, but it's possible that np.inf is interacting badly with linprog's usage of masked arrays. numpy.ma will unilaterally mask-out infs.

[argriffing]

Here's an example where np.inf causes things to go weirdly, but in a different way.

    c = [-1,8,4,-6]
    A_ub = [[-7,-7,6,9],
        [1,-1,-3,0],
        [10,-10,-7,7],
        [6,-1,3,4]]
    b_ub = [-3,np.inf,-6,6]
    A_eq = [[-10,1,1,-8]]
    b_eq = [-4]
    res = linprog(
            c,
            A_ub=A_ub,
            b_ub=b_ub,
            A_eq=A_eq,
            b_eq=b_eq,
            bounds=(-np.inf, np.inf))

_linprog.py", line 215, in _pivot_row
    return True, np.ma.where(q == q.min())[0][0]
IndexError: index 0 is out of bounds for axis 0 with size 0

[robfalck]

Thanks for the example, I'll work on a patch.

On Sun, Dec 21, 2014 at 11:10 PM, argriffing notifications@github.com
wrote:

Here's an example where np.inf causes things to go weirdly, but in a
different way.

c = [-1,8,4,-6]
A_ub = [[-7,-7,6,9],
    [1,-1,-3,0],
    [10,-10,-7,7],
    [6,-1,3,4]]
b_ub = [-3,np.inf,-6,6]
A_eq = [[-10,1,1,-8]]
b_eq = [-4]
res = linprog(
        c,
        A_ub=A_ub,
        b_ub=b_ub,
        A_eq=A_eq,
        b_eq=b_eq,
        bounds=(-np.inf, np.inf))

_linprog.py", line 215, in _pivot_row
return True, np.ma.where(q == q.min())[0][0]
IndexError: index 0 is out of bounds for axis 0 with size 0

—
Reply to this email directly or view it on GitHub
#3888 (comment).

• Rob Falck

[argriffing]

A recent message on the numpy mailing list does not help my confidence in numpy.ma. I try to avoid it in the same way that I avoid numpy.matrix when possible.

[rgommers]

@argriffing it's not that bad and, unlike np.matrix, has interesting use cases. Just don't expect to be able to do with it all the fancy stuff you can do with R's NA.

[argriffing]

On a similar note, I've found a case where an artificial variable remains basic at the end of Phase 1 (which is known to happen sometimes). I need to add the ability to handle the situation, since the failure to do so leads to incorrect results.

Yes, I think this is causing problems like #4746 and #4594. Here's the relevant bit:

    # if pseudo objective is zero, remove the last row from the tableau and
    # proceed to phase 2
    if abs(T[-1, -1]) < tol:
        # Remove the pseudo-objective row from the tableau
        T = T[:-1, :]
        # Remove the artificial variable columns from the tableau
        T = np.delete(T, np.s_[n+n_slack:n+n_slack+n_artificial], 1)
    else:
        # Failure to find a feasible starting point
        status = 2