In [1]:
import sys,os
import numpy as np
import matplotlib as mpl
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import matplotlib.colors as colors
#this works apparently only for savefig stuff
mpl.rcParams['figure.figsize']=(6.0,4.0) #(6.0,4.0)
mpl.rcParams['font.size']=12 #10
mpl.rcParams['savefig.dpi']=100 #72
mpl.rcParams['figure.subplot.bottom']=.1 #.125
plt.rc('font', family='serif')
plt.rc('text', usetex=True)
In [2]:
#inline Shit
%matplotlib inline
%config InlineBackend.figure_format='svg'
%config InlineBackend.rc = {'figure.facecolor': 'white', 'figure.subplot.bottom': 0.125, 'figure.edgecolor': 'white', 'savefig.dpi': 300, 'figure.figsize': (8, 6), 'font.size': 10}
In [3]:
#GUi shit
%matplotlib tk
In [4]:
mpl.get_configdir()
Out[4]:
In [5]:
from matplotlib.table import Table
from matplotlib.table import Cell
class SimplexCalculator:
def __init__(self,A,b,c,indBasic, manual=True, minimize=True, maxSteps = 250):
self.m = A.shape[0];
self.n = A.shape[1];
self.A = A;
self.b = b;
self.c = c;
self.indBasic = indBasic;
self.manual = manual;
self.minimize = minimize;
self.maxSteps = maxSteps;
if(self.m != self.b.shape[0] and self.n != self.c.shape[0] and m != self.indBasic.shape[0]):
raise NameError("Wrong Problem size!")
indNonBasic = list(range(0,self.n)); # variables index 0,1,2,3,4,5...
for i in sorted(self.indBasic, reverse=True):
del indNonBasic[i]
self.indNonBasic = np.array(indNonBasic)
print("NonBasic Indices: " + str(self.indNonBasic))
print("Matrix A: \n" + str(self.A))
print("Vector b: \n" + str(self.b))
AbasisInv = self.A[:,indBasic].I;
self.AbasisInvA = AbasisInv * self.A;
self.xBasic = AbasisInv * self.b
print("Vector xBasis: \n" + str(self.xBasic))
self.cred = self.c.T - self.c[self.indBasic].T * self.AbasisInvA;
print("Vector cred: \n" + str(self.cred))
self.negopt = - self.c[self.indBasic].T * self.xBasic;
print("negOpt: \n" + str(self.negopt))
self.plotList = [];
plot = self.drawTableau(0)
self.plotList.append(plot)
def doSimplex(self, removeIndices=None,newc=None):
if(removeIndices!=None and newc!=None):
# Remove the indices and shrink matrix
for i in removeIndices:
if i in self.indBasic:
raise NameError('Index ' + str(i) + 'appears in the basis, it cannot be removed')
if(np.size(np.nonzero(removeIndices >= self.n)[0]) > 0):
raise NameError('RemoveIndices out of bound!')
# Change Matrix size:
self.AbasisInvA = np.delete(self.AbasisInvA,removeIndices,1)
print("AbasisInvA" + str(self.AbasisInvA))
self.n = self.n - np.size(removeIndices);
self.c = newc
self.cred = self.c.T - self.c[self.indBasic].T * self.AbasisInvA;
print("cred" + str(self.cred))
self.negopt = - self.c[self.indBasic].T * self.xBasic;
print("negOpt: \n" + str(self.negopt))
#Make new basis
self.indNonBasic = np.array([i for i in self.indNonBasic if i not in removeIndices])
print("indNonBasic" + str(self.indNonBasic))
plot = self.drawTableau(0)
#do simplex iterations:
self.step = 1;
abort = False;
print("=====================================================")
while(not abort and self.step <= self.maxSteps):
print(type(self.cred),self.cred.ravel())
if(self.minimize):
self.dualFeasible = np.size(np.nonzero(self.cred>=0)[1]) == self.n
else: # if maximize, the dual reduced cost needs to be negative to be dual feasible
self.dualFeasible = np.size(np.nonzero(self.cred.ravel()>=0)[1]) == self.n
self.primalFeasible = np.size(np.nonzero( self.xBasic>=0)[0]) == self.m
if( not self.primalFeasible and not self.dualFeasible):
raise NameError("The initial tableau is not primal feasible (x>=0, for min problem) and not dual feasible (reduced cost >= 0, for min problem ")
if(self.primalFeasible and self.dualFeasible):
print("Tableau is already at optimum!")
return;
print("Tableau step : " + str(self.step));
if(self.chooseIndex() == "abort"):
break;
if(self.getIndexLeaveEnter() == "abort"):
break;
self.doTableauStep()
print("indBasic: " + str(self.indBasic))
print("indNonBasic: " + str(self.indNonBasic))
#print("Vector cred: \n" + str(self.cred))
#print("Neg. Opt. value: \n" + str(self.negopt))
#print("Matrix AbasisInvA: \n" + str(self.AbasisInvA))
print("xBasic : \n" + str(self.xBasic))
#tableau[0,1::] = cred;
#print("Tableau: \n" + str(tableau))
#tableau[1::,0] = self.xBasic;
#tableau[1::,1::] = self.A;
#tableau = np.mat(np.zeros((self.m+1,self.n+1)));
#update last table with some infos
#entpacke parameter
self.updateTableau(*self.plotList[self.step-1], multFactors = np.concatenate( ( [self.multFactorCred] , self.multFactorA)),
leaveInd = self.leaveIndex, enterInd = self.enterIndex );
plot = self.drawTableau( self.step,self.leaveIndex,self.enterIndex)
#(fig,ax)
self.plotList.append( plot );
self.step+=1
print("=====================================================")
#endwhile
def switchBases(self,leaveIndex,enterIndex):
self.leaveIndex = np.nonzero(self.indBasic == leaveIndex)[0][0];
self.enterIndex = enterIndex;
self.doTableauStep()
plot = self.drawTableau( self.step,self.leaveIndex,self.enterIndex)
self.plotList.append( plot );
self.step+=1
def chooseIndex(self):
if(self.manual): # select the direction!
#Prepare index set to choose from
if(self.primalFeasible):
if(self.minimize):
indicesToChoose = np.nonzero(np.asarray(self.cred[0,self.indNonBasic]).squeeze() < 0)
else:
indicesToChoose = np.nonzero(np.asarray(self.cred[0,self.indNonBasic]).squeeze() > 0)
indicesToChoose = self.indNonBasic[indicesToChoose]
else:
indicesToChoose = np.nonzero(np.asarray(self.xBasic).squeeze() < 0)
indicesToChoose = self.indBasic[indicesToChoose]
if indicesToChoose.size == 0:
print("Optimum reached!")
return "abort"
if(self.primalFeasible):
print("Primal Step: Non Basic Directions to choose enterIndex: "+ str(indicesToChoose))
else:
print("Dual Step: Basic Directions to choose leaveIndex: " + str(indicesToChoose))
#Choose Index from Input
choosenIndex = 0;
choosenIndexOK = False;
while( not choosenIndexOK):
s = input();
if(s == "exit" or s =="end" or s=="stop" or s == "break"):
print("Abort iteration!")
return "abort"
choosenIndex = int(s)
if choosenIndex in indicesToChoose:
choosenIndexOK=True;
else:
print("You have choosen an index which is not valid to choose from, try again...")
#endwhile
if(self.primalFeasible):
self.enterIndex = choosenIndex;
else:
print(self.indBasic)
self.leaveIndex = np.nonzero(self.indBasic == choosenIndex)[0][0];
else: #non manual
if(self.primalFeasible):
if(self.minimize):
enterIndNonBasic = self.cred[0,self.indNonBasic].argmin(axis=1)[0,0]
self.enterIndex = self.indNonBasic[enterIndNonBasic]
#print(\"NonBasic Index entering: \" + str(enterIndex));
if( self.cred[0,self.enterIndex] >= 0.0 ):
print("Optimum reached!")
return "abort"
else:
enterIndNonBasic = self.cred[0,self.indNonBasic].argmax(axis=1)[0,0];
self.enterIndex = self.indNonBasic[enterIndNonBasic]
#print(\"NonBasic Index entering: \" + str(enterIndex));
if( self.cred[0,self.enterIndex] <= 0.0 ):
print("Optimum reached!")
return "abort"
else: #Dual Feasible
leaveIndBasic = self.xBasic.argmin(axis=1)[0,0];
self.leaveIndex = leaveIndBasic;
if(self.primalFeasible):
print("Enter index: " + str(self.enterIndex))
else:
print("Leave index: " + str(self.leaveIndex))
def getIndexLeaveEnter(self):
if(self.primalFeasible):
# determine smallest factor theta in this enter direction j;
self.leaveIndex = -1;
minTheta = 0;
for i in range(self.m):
# if u > 0 then:
if self.AbasisInvA[i,self.enterIndex] > 0.0:
value = self.xBasic[i,0] / self.AbasisInvA[i,self.enterIndex];
if(value < minTheta or self.leaveIndex == -1):
self.leaveIndex = i;
minTheta = value;
if(self.leaveIndex == -1):
print("Problem is unbounded!")
return "abort";
#=============================================================
else: # Dual
# determine smallest factor theta in this leave direction i;
self.enterIndex = -1;
minTheta = 0;
for j in self.indNonBasic:
if self.AbasisInvA[self.leaveIndex,j] < 0.0:
value = self.cred[0,j] / abs(self.AbasisInvA[self.leaveIndex,j]);
if(value < minTheta or self.enterIndex == -1):
self.enterIndex = j;
minTheta = value;
if(self.enterIndex == -1):
print("Problem is unbounded!")
return "abort";
#=============================================================
if(self.primalFeasible):
print("Leave index: " + str(self.leaveIndex))
else:
print("Enter index: " + str(self.enterIndex))
def doTableauStep(self):
# Calculate mult factor for matrix pivoting
self.multFactorA = np.array([0.0]*self.m);
for i in range(self.m):
if i!=self.leaveIndex:
self.multFactorA[i] = - self.AbasisInvA[i,self.enterIndex] / self.AbasisInvA[self.leaveIndex,self.enterIndex]
self.multFactorCred = -self.cred[0,self.enterIndex] / self.AbasisInvA[self.leaveIndex,self.enterIndex];
print("Multfactors A \n" + str(self.multFactorA));
print("Multfactors cred \n" + str(self.multFactorCred));
# Multiplying through the matrix AbasisInvA
# for all rows
for i in range(self.m):
if(i != self.leaveIndex):
self.AbasisInvA[i] += self.multFactorA[i] * self.AbasisInvA[self.leaveIndex];
self.xBasic[i] += self.multFactorA[i] * self.xBasic[self.leaveIndex]
# do the cred row
print(self.AbasisInvA[self.leaveIndex])
self.cred += self.multFactorCred * self.AbasisInvA[self.leaveIndex];
self.negopt += self.multFactorCred * self.xBasic[self.leaveIndex];
# do the leaveIndex Row
self.xBasic[self.leaveIndex] /= self.AbasisInvA[self.leaveIndex,self.enterIndex];
self.AbasisInvA[self.leaveIndex] /= self.AbasisInvA[self.leaveIndex,self.enterIndex];
#switch indexes;
self.indNonBasic[self.indNonBasic == self.enterIndex] = self.indBasic[self.leaveIndex];
self.indBasic[self.leaveIndex] = self.enterIndex;
def updateTableau(self,fig,ax,tb, width,height,offsetRow,offsetCol,fmt,multFactors,leaveInd,enterInd):
# Add cells multFactor
print("MultFactors: " + str(multFactors))
for i in range(np.shape(multFactors)[0]):
tb.add_cell(i+offsetRow, offsetCol+self.n+1, width/2, height, text=fmt.format(multFactors[i]), loc='center', facecolor='white')
#ax.plot([1.5*width +enterInd*width,1.5*width +enterInd*width + width],[0.1+ m*height - (leaveInd+1)*height,0.1+ m*height - (leaveInd+1)*height],linestyle='-',
# linewidth=2, color='k')
pivotCell = tb.get_celld()[(leaveInd+offsetRow+1,enterInd+offsetCol+1)]
if(pivotCell != None):
pivotCell.set_facecolor('#999900')
fig.canvas.draw()
def drawTableau(self,step,leaveInd=-1,enterInd=-1, fmt='{:.2f}', bkg_colors=['yellow', 'white']):
fig, ax = plt.subplots(figsize=(10,6))
#fig.set_size_inches(10,6)
ax.set_axis_on()
ax.autoscale(False)
tb = Table(ax, bbox=[0,0.1,1,1-0.1])
nRows = self.m+1+1;
nCols = self.n+1+1+1;
width, height = 1.0 / nCols, 1.0 / nRows
offsetRow = 1;
offsetCol = 1;
# Add label cells
for (i,j),x in np.ndenumerate(self.xBasic):
tb.add_cell(i+1+offsetRow, 0, width/2.0, height, text=r'\Huge $x_%i$' % self.indBasic[i], loc='center', facecolor='white')
for (i,j),x in np.ndenumerate(self.cred):
if j in self.indBasic:
color = 'blue'
else:
color = 'white'
tb.add_cell(0, j+1+offsetCol, width, height/2.0, text=r'\Huge $x_%i$' % j, loc='center', facecolor=color)
# Add value cells
for (i,j), val in np.ndenumerate(self.AbasisInvA):
y = j+1+offsetCol;
x = i+1+offsetRow;
tb.add_cell(x,y, width, height, text=fmt.format(val), loc='center', facecolor='green')
# Add cells xBasic
for (i,j), val in np.ndenumerate(self.xBasic):
tb.add_cell(i+1+offsetRow, j+offsetCol, width, height, text=fmt.format(val), loc='center', facecolor='red')
# Add cells cred
for (i,j), val in np.ndenumerate(self.cred):
color = 'cyan'
tb.add_cell(0+offsetRow, j+1+offsetCol, width, height, text=fmt.format(val), loc='center', facecolor=color)
# Add cells negopt
tb.add_cell(0+offsetRow, 0+offsetCol, width, height, text=fmt.format(self.negopt[0,0]), loc='center', facecolor='yellow')
ax.add_table(tb)
ax.set_title("Simplex Tableau at step %i" % step)
pad = 0.01
t1=plt.Text(x=width ,y=0.1-pad, text=r'\huge $\mathbf{x}_B=\mathbf{A}_B^{-1} \mathbf{b}$',va='top', ha='center')
t2=plt.Text(x=1.5*width + width*self.n/2.0, y=0.1-pad, text=r'\huge $\mathbf{A}_B^{-1} \mathbf{A}$', va='top',ha='center')
t3=plt.Text(x=width ,y=0.1+height*(self.m+1)+pad, text=r'\huge $-\mathbf{c}_B^{\top}\mathbf{x}_B$', va='bottom',ha='center')
ax.add_artist(t1)
ax.add_artist(t2)
ax.add_artist(t3)
fig.canvas.draw()
return (fig,ax,tb,width,height,offsetRow,offsetCol, fmt)
In [5]:
# Input Matrices
plt.close('all')
np.set_printoptions(precision=2)
# Feasible Problem
#A = np.matrix([[1,3,-1,0,0,1,0,0],[4,2,0,-1,0,0,1,0],[1,2,0,0,1,0,0,1]])
#b = np.matrix([9,16,20]).T
#c = np.matrix([0,0,0,0,0,1,1,1]).T
#indBasic = np.array([5,6,7]); # variables index 0,1,2,3,4,5...
#doSimplexAlgorithm(A,b,c,indBasic,False,True,2)
# Infeasible Problem
A = np.matrix([[-1,-3,1,0,0,-1],[-4,-2,0,1,0,-1],[1,2,0,0,1,-1]])
b = np.matrix([-9,-16,20]).T
c = np.matrix([0,0,0,0,0,1]).T
indBasic = np.array([2,3,4]); # variables index 0,1,2,3,4,5...
s = SimplexCalculator(A,b,c,indBasic,False,True,4);
s.doSimplex();
In [7]:
s.switchBases(3,1);
#s.switchBases(4,2);
c = np.matrix([-2,-4,0,0,0]).T
s.doSimplex(removeIndices=np.array([5]),newc=c);
In [6]:
A = np.matrix([[-1,-3,1,0,0],[-4,-2,0,1,0],[1,2,0,0,1]])
b = np.matrix([-9,-16,20]).T
c = np.matrix([0,0,0,0,0]).T
indBasic = np.array([2,3,4]); # variables index 0,1,2,3,4,5...
s = SimplexCalculator(A,b,c,indBasic,False,True,4);
s.doSimplex();
In [8]: