Let's consider the following integer program:
$\begin{array}{rccccl}\max&&&x_2&&\\&-5x_1&+&3x_2&\leq&0\\&5x_1&+&3x_2&\leq&5\\&x_1,&&x_2&\geq&0\\&x_1,&&x_2&&\text{integer}\end{array}$
{{{id=171| pts = []; for i in range(10): for j in range(10): pts.append( (i,j) ); P = list_plot( pts ); x = var('x'); P += plot( (5/3)*x, (x, 0, 20) ); P += plot( (-5/3)*x + 5, (x,0,20) ); P.set_axes_range( 0, 5, 0, 4 ); P ///
}}}
{{{id=2|
def Pivot(A, B, x, y):
#print "Pivoting on row %d and column %d." % (x,y);
m = len(A);
n = len(A[0]);
val = A[x][y];
for j in range(n):
A[x][j] = A[x][j] / val;
for i in range(m):
if i != x:
val = A[i][y];
for j in range(n):
A[i][j] = A[i][j] - val * A[x][j];
B[x] = y; # modify basis element
///
}}}
{{{id=141|
def Rebase(A, B):
# make sure the basis is basic
for i in range(1, len(B)):
Pivot(A, B, i, B[i]);
///
}}}
{{{id=170|
def PrimalSimplex(A, B):
for j in range(1,len(A[0])):
if A[0][j] < 0:
# find a pivot point
bestrow=-1;
bestratio = 10000;
for i in range(1, len(A)):
if A[i][j] > 0 and A[i][0] / A[i][j] < bestratio:
bestrow = i;
bestratio = A[i][0] / A[i][j];
if bestrow > 0:
Pivot(A, B, bestrow, j);
return PrimalSimplex(A, B);
else:
return None, None; # unbounded!
return A, B;
///
}}}
{{{id=130|
def DualSimplex(A, B):
for i in range(1,len(A)):
if A[i][0] < 0:
# find a pivot point
bestcol=-1;
bestratio = -10000;
for j in range(1, len(A[i])):
if A[i][j] < 0 and A[0][j] / A[i][j] > bestratio:
bestcol = j;
bestratio = A[0][j] / A[i][j];
if bestcol > 0:
Pivot(A, B, i, bestcol);
return DualSimplex(A, B);
else:
return None, None; # unbounded!
return A, B;
///
}}}
{{{id=3|
def PrintTableau(A,B):
str_length=4;
col_lengths = {};
for j in range(len(A[0])):
col_lengths[j] = 0;
for ii in range(len(A)):
if len(str(A[ii][j])) > col_lengths[j]:
col_lengths[j] = len(str(A[ii][j]));
for i in range(len(A)):
row_str = "[";
for j in range(len(A[i])):
s = str(A[i][j]);
for t in range(col_lengths[j] - len(s)):
row_str += " ";
row_str += s;
if j == 0:
row_str+="|";
else:
row_str+=" ";
print row_str + "]";
print;
///
}}}
{{{id=135|
def PrintActiveTableaus(Active):
for A, B in Active:
print "Current LP Lower Bound: ", -A[0][0], " approx: ", N(-A[0][0]);
PrintTableau(A, B);
print;
///
}}}
{{{id=136|
def BranchOnVariable(T, B, var, x):
# make two deep copies, with extra slack variable column
T1 = [ [ T[i][j] for j in range(len(T[i])) ] + [0] for i in range(len(T)) ];
T2 = [ [ T[i][j] for j in range(len(T[i])) ] + [0] for i in range(len(T)) ];
# make two new constraint rows
C1 = [ 0 ] * len(T[0]) + [1];
C1[0] = floor(x);
C1[var] = 1;
C2 = [ 0 ] * len(T[0]) + [-1];
C2[0] = (floor(x)+1);
C2[var] = 1;
T1.append(C1);
T2.append(C2);
# make two new bases
B1 = [ B[i] for i in range(len(B)) ] + [ len(T[0]) ];
B2 = [ B[i] for i in range(len(B)) ] + [ len(T[0]) ];
#Rebase(T1, B1);
#Rebase(T2, B2);
return (T1,B1), (T2,B2);
///
}}}
{{{id=163|
colors = [ "red", "green", "purple", "orange", "black", "blue" ];
fillcolors = [ (1,0.25,0.25), (0.25, 1, 0.25), (0.25, 0.25, 1), (1,1,0.5), (1,0.5,1), (0.5,1,1), (0.5,0.5,0.5) ];
def PlotActivePoints(OrigP, Active, cindex):
P = OrigP;
for A,B in Active:
inv = [ -1, -1];
for j in range(1,len(B)):
if B[j] == 1:
inv[0] = j;
if B[j] == 2:
inv[1] = j;
P += list_plot( [ (A[inv[0]][0], A[inv[1]][0])], size=50, color=colors[cindex]);
return P;
///
}}}
{{{id=137|
A = [\
[ 0, -1, -4, 0, 0, 0 ],\
[ 0, (5/3), -1, -1, 0, 0],\
[ 5, (5/3), 1, 0, 1, 0]];
B = [ 0, 3, 4 ];
PrintTableau(A,B);
///
[0| -1 -4 0 0 0 ]
[0|5/3 -1 -1 0 0 ]
[5|5/3 1 0 1 0 ]
}}}
{{{id=139|
Rebase(A,B);
PrintTableau(A,B);
PrimalSimplex(A, B);
PrintTableau(A,B);
///
[0| -1 -4 0 0 0 ]
[0|-5/3 1 1 0 0 ]
[5| 5/3 1 0 1 0 ]
[23/2|0 0 17/10 23/10 0 ]
[ 5/2|0 1 1/2 1/2 0 ]
[ 3/2|1 0 -3/10 3/10 0 ]
}}}
{{{id=145|
Active = [ (A,B) ];
PrintActiveTableaus(Active);
///
Current LP Lower Bound: -23/2 approx: -11.5000000000000
[23/2|0 0 17/10 23/10 0 ]
[ 5/2|0 1 1/2 1/2 0 ]
[ 3/2|1 0 -3/10 3/10 0 ]
}}}
{{{id=140|
pts = [];
for i in range(10):
for j in range(10):
pts.append( (i,j) );
P = list_plot( pts );
x = var('x');
P += plot( (5/3)*x, (x, 0, 20) );
P += plot( (-5/3)*x + 5, (x,0,20) );
P = PlotActivePoints(P, Active, 0);
P.set_axes_range( 0, 5, 0, 4 );
P
///
}}}
{{{id=143|
A0, B0 = Active[0];
Active.remove( (A0,B0) );
for X in BranchOnVariable(A0, B0, 1, 1.5):
Active.append(X);
A, B = X;
PrintTableau(A,B);
Rebase(A,B);
DualSimplex(A, B);
PrintActiveTableaus(Active);
///
[23/2|0 0 17/10 23/10 0 0 ]
[ 5/2|0 1 1/2 1/2 0 0 ]
[ 3/2|1 0 -3/10 3/10 0 0 ]
[ 1|1 0 0 0 0 1 ]
[23/2|0 0 17/10 23/10 0 0 ]
[ 5/2|0 1 1/2 1/2 0 0 ]
[ 3/2|1 0 -3/10 3/10 0 0 ]
[ 2|1 0 0 0 0 -1 ]
Current LP Lower Bound: -23/3 approx: -7.66666666666667
[23/3|0 0 4 0 0 23/3 ]
[ 5/3|0 1 1 0 0 5/3 ]
[ 1|1 0 0 0 0 1 ]
[ 5/3|0 0 -1 1 0 -10/3 ]
Current LP Lower Bound: -26/3 approx: -8.66666666666667
[26/3|0 0 0 4 0 17/3 ]
[ 5/3|0 1 0 1 0 5/3 ]
[ 2|1 0 0 0 0 -1 ]
[ 5/3|0 0 1 -1 0 -10/3 ]
}}}
{{{id=152|
P = polygon2d([[1.05,0], [1.95,0], [1.95,6], [1.05,6]], rgbcolor=fillcolors[0], alpha=0.5) + P;
P = PlotActivePoints(P, Active, 1);
P.set_axes_range(0,5,0,4);
P
///
}}}
{{{id=144|
A0, B0 = Active[0];
Active.remove( (A0,B0) );
for X in BranchOnVariable(A0, B0, 2, 1.5):
A, B = X;
PrintTableau(A, B);
Rebase(A,B);
A, B = DualSimplex(A, B);
if A is None:
print "Infeasible (Dual Unbounded)!\n";
else:
Active.append( (A, B) );
PrintActiveTableaus(Active);
///
[23/3|0 0 4 0 0 23/3 0 ]
[ 5/3|0 1 1 0 0 5/3 0 ]
[ 1|1 0 0 0 0 1 0 ]
[ 5/3|0 0 -1 1 0 -10/3 0 ]
[ 1|0 1 0 0 0 0 1 ]
[23/3|0 0 4 0 0 23/3 0 ]
[ 5/3|0 1 1 0 0 5/3 0 ]
[ 1|1 0 0 0 0 1 0 ]
[ 5/3|0 0 -1 1 0 -10/3 0 ]
[ 2|0 1 0 0 0 0 -1 ]
Infeasible (Dual Unbounded)!
Current LP Lower Bound: -26/3 approx: -8.66666666666667
[26/3|0 0 0 4 0 17/3 ]
[ 5/3|0 1 0 1 0 5/3 ]
[ 2|1 0 0 0 0 -1 ]
[ 5/3|0 0 1 -1 0 -10/3 ]
Current LP Lower Bound: -5 approx: -5.00000000000000
[ 5|0 0 0 0 0 1 4 ]
[ 1|0 1 0 0 0 0 1 ]
[ 1|1 0 0 0 0 1 0 ]
[7/3|0 0 0 1 0 -5/3 -1 ]
[2/3|0 0 1 0 0 5/3 -1 ]
}}}
{{{id=153|
P = polygon2d([[0,1.05],[0,1.95], [1.05,1.95], [1.05,1.05]], rgbcolor=(0.5,0.5,0.5), alpha=0.5) + P;
P = PlotActivePoints(P, Active, 2);
P.set_axes_range(0,5,0,4);
P
///
}}}
{{{id=146|
A0, B0 = Active[0];
Active.remove( (A0,B0) );
for X in BranchOnVariable(A0, B0, 2, 1.5):
A, B = X;
Rebase(A,B);
PrintTableau(A, B);
A, B = DualSimplex(A, B);
if A is None:
print "Infeasible (Dual Unbounded)!\n";
else:
Active.append( (A, B) );
PrintActiveTableaus(Active);
///
[26/3|0 0 0 4 0 17/3 0 ]
[ 5/3|0 1 0 1 0 5/3 0 ]
[ 2|1 0 0 0 0 -1 0 ]
[ 5/3|0 0 1 -1 0 -10/3 0 ]
[-2/3|0 0 0 -1 0 -5/3 1 ]
[26/3|0 0 0 4 0 17/3 0 ]
[ 5/3|0 1 0 1 0 5/3 0 ]
[ 2|1 0 0 0 0 -1 0 ]
[ 5/3|0 0 1 -1 0 -10/3 0 ]
[-1/3|0 0 0 1 0 5/3 1 ]
Infeasible (Dual Unbounded)!
Current LP Lower Bound: -5 approx: -5.00000000000000
[ 5|0 0 0 0 0 1 4 ]
[ 1|0 1 0 0 0 0 1 ]
[ 1|1 0 0 0 0 1 0 ]
[7/3|0 0 0 1 0 -5/3 -1 ]
[2/3|0 0 1 0 0 5/3 -1 ]
Current LP Lower Bound: -32/5 approx: -6.40000000000000
[32/5|0 0 0 3/5 0 0 17/5 ]
[ 1|0 1 0 0 0 0 1 ]
[12/5|1 0 0 3/5 0 0 -3/5 ]
[ 3|0 0 1 1 0 0 -2 ]
[ 2/5|0 0 0 3/5 0 1 -3/5 ]
}}}
{{{id=154|
P = polygon2d([[1.95,1.05],[1.95,1.95], [6,1.95], [6,1.05]], rgbcolor=(0.25,1,0.25), alpha=0.5) + P;
P = PlotActivePoints(P, Active, 2);
P.set_axes_range(0,5,0,4);
P
///
}}}
{{{id=147|
A0, B0 = Active[1];
Active.remove( (A0,B0) );
for X in BranchOnVariable(A0, B0, 1, 2.5):
A, B = X;
Rebase(A,B);
PrintTableau(A, B);
A, B = DualSimplex(A, B);
if A is None:
print "Infeasible (Dual Unbounded)!\n";
else:
Active.append( (A, B) );
PrintActiveTableaus(Active);
///
[32/5|0 0 0 3/5 0 0 17/5 0 ]
[ 1|0 1 0 0 0 0 1 0 ]
[12/5|1 0 0 3/5 0 0 -3/5 0 ]
[ 3|0 0 1 1 0 0 -2 0 ]
[ 2/5|0 0 0 3/5 0 1 -3/5 0 ]
[-2/5|0 0 0 -3/5 0 0 3/5 1 ]
[32/5|0 0 0 3/5 0 0 17/5 0 ]
[ 1|0 1 0 0 0 0 1 0 ]
[12/5|1 0 0 3/5 0 0 -3/5 0 ]
[ 3|0 0 1 1 0 0 -2 0 ]
[ 2/5|0 0 0 3/5 0 1 -3/5 0 ]
[-3/5|0 0 0 3/5 0 0 -3/5 1 ]
Current LP Lower Bound: -5 approx: -5.00000000000000
[ 5|0 0 0 0 0 1 4 ]
[ 1|0 1 0 0 0 0 1 ]
[ 1|1 0 0 0 0 1 0 ]
[7/3|0 0 0 1 0 -5/3 -1 ]
[2/3|0 0 1 0 0 5/3 -1 ]
Current LP Lower Bound: -6 approx: -6.00000000000000
[ 6|0 0 0 0 0 0 4 1 ]
[ 1|0 1 0 0 0 0 1 0 ]
[ 2|1 0 0 0 0 0 0 1 ]
[7/3|0 0 1 0 0 0 -1 5/3 ]
[ 0|0 0 0 0 0 1 0 1 ]
[2/3|0 0 0 1 0 0 -1 -5/3 ]
Current LP Lower Bound: -3 approx: -3.00000000000000
[3|0 0 0 4 0 0 0 17/3 ]
[0|0 1 0 1 0 0 0 5/3 ]
[3|1 0 0 0 0 0 0 -1 ]
[5|0 0 1 -1 0 0 0 -10/3 ]
[1|0 0 0 0 0 1 0 -1 ]
[1|0 0 0 -1 0 0 1 -5/3 ]
}}}
{{{id=148|
P = polygon2d([[2.05,0],[2.95,0], [2.95,1.05], [2.05,1.05] ], rgbcolor=(0.25,0.25,1), alpha=0.5) + P;
P = PlotActivePoints(P, Active, 2);
P.set_axes_range(0,5,0,4);
P
///
}}}
{{{id=158|
///
}}}
{{{id=172|
///
}}}