Incorrect basis inverse after solving small LP

[tkralphs]

I've encountered a small LP for which the the basis inverse (lp.basisInverse) automatically calculated by CyLP seems to be wrong (along with associated data attrictures lp.tableau and lp.rhs). However, if I calculate the basis inverse myself using getBinvACol, it is correct. Could this have something to do with scaling? Is the basis inverse computed by CyLP with respect to a scaled version of the matrix rather than the original?

import numpy as np
from cylp.cy import CyClpSimplex
from cylp.py.modeling import CyLPArray

A = [[ -8, 30], [ -2, -4], [-14, 8], [ 2, -36], [30, -8], [10, 10]]
b = [115,  -5, 1, -5, 191, 127]
c = [1, -1]

lp = CyClpSimplex()
A = np.matrix(A)
b = CyLPArray(b)
c = CyLPArray(c)

x = lp.addVariable('x', 2)

lp += x >= 0
lp += A * x <= b
lp.objective = -c * x if sense[0] == 'Max' else c * x
lp.primal(startFinishOptions = 'x')
Binv = np.zeros(shape = (lp.nConstraints, lp.nConstraints))
for i in range(lp.nVariables, lp.nVariables+lp.nConstraints):
   lp.getBInvACol(i, Binv[i-lp.nVariables,:])
print 'Correct basis inverse:')
print Binv
print 'Incorrect basis inverse:')
print lp.basisInverse
print 'Correct RHS:')
np.dot(lp.basisInverse, myb)
print 'Incorrect RHS:')
print lp.rhs

[mpy]

Thank's for the issue, that looks bad. Scaling is probably the reason. But let me figure that out...

[choidami]

Was this problem ever resolved?

[tkralphs]

No, but the workaround of just computing the basis inverse manually has been working fine for me. There are bigger fish to fry :).

[choidami]

Sorry if I'm being stupid, but If I understand correctly, your computation of Binv uses getBInvACol, while the code in lp.basisInverse does getBInvCol. If I do:

import numpy as np
from cylp.cy import CyClpSimplex
from cylp.py.modeling import CyLPArray

A = [[ -8, 30], [ -2, -4], [-14, 8], [ 2, -36], [30, -8], [10, 10]]
b = [115,  -5, 1, -5, 191, 127]
c = [1, -1]

lp = CyClpSimplex()
A = np.matrix(A)
b = CyLPArray(b)
c = CyLPArray(c)

x = lp.addVariable('x', 2)

lp += x >= 0
lp += A * x <= b
lp.objective = -c * x
lp.primal(startFinishOptions = 'x')

myBinv = np.zeros((lp.nConstraints, lp.nConstraints))
for i in range(lp.nConstraints):
  tmp = np.zeros(lp.nConstraints)
  lp.getBInvCol(i, tmp)
  myBinv[:, i] = tmp
print("myBinv:\n", myBinv)
print("lp.basisInverse:\n", lp.basisInverse)

Then the two results match. Note that the tmp variable makes a difference, otherwise, the result gets stored in the rows instead of the column.