Python: PyLinpro, a pure simplex tableau solver for linear programming

I have changed the title of this updated blog from "A buggy pure Python ...." since I believe that it is better to overwrite the old buggy code with the corrected current version. Hope the user of the code will understand. I rewrote the Python code three times and in the last iteration of the source code I resorted to the use of the "indent" program to properly indent the C/C++ clinpro.cpp program as a base for the final Python program.The translation was mostly mechanical (though I did it by hand, the conversion of programs in C/C++ to Python can be automated with a little handwringing).

The C/C++ Linpro code ( see previous blog) works fast enough in my laptop, that I thought I would first translate it into pure Python before I implement the Revised simplex method. I have googled for a suitable pseudocode for the Revised Simplex but I found nothing which is easy to understand and the descriptions made me confused.

For instance, some references claim that the Revised Simplex Method involves only using matrix- vector operations, some claim that it uses the inverse of the Basis all throughout and some describes the RSM as faster than the simple Tableau method which perforrms a lot of extraneous generated values. The best description maybe the use of matrix LU decomposition for inverting the basis if needed. But we need a tableau solver for comparison of the performance. With Python version, one can experiment with various elements, using rationals for example.

To browse/scan the Python code horizontally, place the mouse cursor inside the codebox and click it, then use the left and right arrow keys.

?View Code PYTHON
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
"""
file            pylinpro.py
author      ernesto.adorio@gmail.com
                UPEPP UP at Clark Field
                Pampanga, the PHilippines 
                zip code: 2009
source      translation  from BASIC to C/C++ to Python
                of G.E. linpro program.
revisions   2009.01.25  first version
                2009.01.26  debugged and updated.
"""
 
from matlib import matzero, matdim, matprint, matsetblock
 
def pylinpro (AA,M,N,L,  E,G, artifcode,opticode,Basis, disp_code, ZTOL = 1.0e-6):
  """
  Parameters
  AA           AA is a dynamically allocated matrix
  and should be dimensioned (M+2) * JLIM
  where
  JLIM = N + G + L + 1 if artifcode = 0
  and
  JLIM = N + G + M + 1 if artifcode = 1.
 
  B coefficients are stored in AA[ ][0]
  Z coefficients are stored in AA[M][]
  Phase 1 coefficients are stored in AA[M+1][]
 
  M             Number of constraint equations
  N             Number of variables
  L             Number of less thans
  E             Number of equal tos
  G             Number of greater thans
 
  artifcode
  0 - do not process artificial variables,
  1 - otherwise
 
  opticode
  0 - minimization
  1 - maximization
 
  Basis[]  Output current Basis column.
  Variables are numbered 1 to L+E+G.
 
  LPcode   Return LPcode code of LP routine:
  0 - Successful
  1 - Array too small
  2 - Inconsistent equations
  3 - B coefficients must be positive !
  4 - Infeasible, artificial variable in solution
  5 - No solution, unbounded solution
  6 - Excessive number of iterations
  7 - Failed to allocate space for VC
 
  disp_code
  0   no display
  1   for display of results only
  2   for display of Basis and objective function
  3   for display of tableau at each iteration
 
  Description
  This is a C++ translation of the G.E.(General Electric)
  LINPRO program.
  """
 
  def disp_Basis():
      print "Iteration No. %d " % LPIter
      print "Basis            Value "
      # print(" Variable ")
      for sub_i in range(M): # (sub_i = 0 sub_i < M sub_i++):
          print " %4d  %15.5lf" % (Basis[sub_i], AA[sub_i][0])
 
      print "The value of the Objective function is %11.5lf" % (-AA[M][0] if (opticode == 0) else AA[M][0])#? -AA[M][0] : AA[M][0])
 
  def disp_Tableau():
      print "Tableau at iteration : %-d " % LPIter
      print "Basis      Value        ",
      for sub_j in range(1, JLIM): # (sub_j = 1 sub_j < JLIM sub_j++):
          print "%-8d" % VC[sub_j],
 
      print ""
 
      for sub_i in range(M+2): #(sub_i = 0 sub_i < M + 2 sub_i++):
          if (sub_i == M):
              print "%4c" % ' ',
              for sub_j in range((JLIM+1)*8):  # (sub_j = 0 sub_j < (JLIM + 1) * 8 sub_j++):
                  print "%c" % '=',
              print ""
 
          if (sub_i < M):
              print "%8d" % Basis[sub_i],
 
          else:
              print "%8c" % ' ',
 
 
          print "%12.3lf" % AA[sub_i][0],
          for sub_j in range(1, JLIM):   #(sub_j = 1 sub_j < JLIM sub_j++):
              print "%8.3lf" % AA[sub_i][sub_j],
          print 
 
 
  MAXDOUBLE = 1.0e100
  print "M = %d" % M
  print "N = %d" % N
  print "L = %d" % L
  print "E = %d" % E
  print "G = %d" % G
 
  K = N + G + M + 1
  if (artifcode):
      JLIM = K
  else :
      JLIM = N + G + L + 1
 
 
 
  # Check for  size of AA, and equation for L, E and G 
  nrows, ncols  = matdim(AA)
 
  print "nrows = %d" % nrows
  print "ncols = %d" % ncols
  print "JLIM = %d"  % JLIM
  print "M = %d" % M
 
  if ((nrows < M + 2) or ((ncols < JLIM))):
      LPcode = 1		#  Array too small  
      return LPcode
 
  elif (L + E + G != M):
      LPcode = 2		#  Inconsistent equations. 
      return LPcode
 
 
  # Check for B 
  for I in range(M+1): #(I = 0 I <= M I++):
      if (AA[I][0] < 0.0):
	  LPcode = 3
	  return LPcode
 
 
 
  # Initialize the Basis 
  for I in range(M): # (I = 0 I < M I++):
      Basis[I] = I + 1 + N + G
 
  for I in range(M, M+2): #(I = M I <= M + 1 I++):
      Basis[I] = 0
 
 
 
  # Initialize the Variable Header  
  VC = range(JLIM) #(int *) (malloc (sizeof (int) * JLIM))
  if (VC is None):
      LPcode = 7
      return LPcode
 
 
  for J in range(1, JLIM):# (J = 1 J < JLIM J++):
      VC[J] = J
 
 
  # Initialize to zero the surplus and slack regions
  for I in range(M+2): #(I = 0 I <= M + 1 I++):
      for J in range(N+1, JLIM): #(J = N + 1 J < JLIM J++):
	  AA[I][J] = 0.0
 
 
 
  # The  surplus part
  J = N + 1
  for I in range(L + E, M): #(I = L + E I < M I++):
      AA[I][J] = -1.0
      J+= 1
 
 
  # The slack part of the array (include artificial
  # part if JLIM == K)
  J = N + G + 1
  for I in range(M): #(I = 0 I < M I++):
      AA[I][J] = 1.0
      J+= 1
      if (J >= JLIM):
	break
 
 
  # Negate the objective coefficients for
  # maximization problems
  if (opticode == 1):
      for J in range(1, JLIM):  #(J = 1 J < JLIM J++):
	  AA[M][J] = -AA[M][J]
 
 
 
  # Compute the last row for phase 1 
  do_phase1 = False if (G + E == 0) else True # 
  if (do_phase1):
      for J in range(JLIM): #(J = 0 J < JLIM J++):
	  RSum = 0.0
	  for I in range(L, M): #(I = L I < M I++):
	      RSum = RSum + AA[I][J]
 
	  AA[M + 1][J] = -RSum
 
 
  else :
      for J in range(JLIM): #(J = 0 J < JLIM J++):
	  AA[M + 1][J] = 0.0
 
 
  Test_Row = M+1 if (do_phase1) else M # ? M + 1 : M
  maxLPIter = 2 * K
  for LPIter in range(1, maxLPIter+1): #(LPIter = 1 LPIter <= maxLPIter LPIter++):
 
      # Step 1:  find pivot column 
      while (1):
	  if (disp_code == 2):
	    disp_Basis()
	  if (disp_code == 3):
	    disp_Tableau()
 
	  # find the most negative number in the test row 
	  tempR = 0.0
	  Piv_Col = 0
	  for J in range(JLIM -1, 0, -1): #(J = JLIM - 1 J >= 1 J--):
	      R = AA[Test_Row][J]
	      if (R < tempR):
		  if (do_phase1 or (AA[Test_Row + 1][J] == 0)):
		      tempR = R
		      Piv_Col = J
 
 
 
 
	  if ((Piv_Col == 0) and do_phase1):
	      Test_Row-= 1
	      do_phase1 = False
	      if (disp_code == 3):
		  print " End of Phase 1 "
 
 
	  else:
	      break
 
 
 
      if (Piv_Col == 0):
 
	  for I in range(M): #(I = 0 I < M I++):
	      if ((Basis[I] >= JLIM) and (abs (AA[I][0]) > ZTOL)):
		  if (disp_code >= 1):
		      print "Artificial variable in solution."
 
		  LPcode = 4
 
		  return LPcode
 
 
 
	  if (disp_code >= 1):
	      disp_Basis()
	      print "Succesful termination "
 
 
	  LPcode = 0		# successful termination 
 
	  return LPcode
 
 
      # Step 2: Find pivot row 
      R = MAXDOUBLE
      Piv_Row = -1
      Tie_Row = K + 2
 
      for I in range(M-1, -1, -1): # (I = M - 1 I >= 0 I--):
	  tempR = AA[I][Piv_Col]
	  if (tempR > ZTOL) :
	      Ratio = abs (AA[I][0] / tempR)
	      if (Ratio < R):
		  R = Ratio
		  Piv_Row = I
 
	      elif (Ratio == R): # break row ties 
		  if (Basis[Tie_Row] > Basis[I]):
		      Piv_Row = I
 
 
	      Tie_Row = Basis[I]
 
 
 
      if (Piv_Row == -1):
	  if (disp_code >= 1):
	      print "Unbounded Solution"
 
	  LPcode = 5
 
	  return LPcode
 
 
      if (disp_code >= 2):
	  print "Pivot Row: %d, Pivot Column: %d " % (Piv_Row, Piv_Col)
 
 
      # Step 3: Reduce the array 
      invPivot = 1.0 / AA[Piv_Row][Piv_Col]
      for I in range(Test_Row + 1): #(I = 0 I <= Test_Row I++):
	  R = AA[I][Piv_Col] * invPivot
	  if (I != Piv_Row):
	      for J in range( JLIM): # (J = 0 J < JLIM J++):
		  if (J != Piv_Col):
		      AA[I][J] = AA[I][J] - AA[Piv_Row][J] * R
		      if (abs (AA[I][J]) <= ZTOL):
			  AA[I][J] = 0.0
 
 
 
	      AA[I][Piv_Col] = 0.0
 
 
      for J in range(JLIM):   # (J = 0 J < JLIM J++):
	  AA[Piv_Row][J] = AA[Piv_Row][J] * invPivot
 
      AA[Piv_Row][Piv_Col] = 1.0
 
      # Step 4: replace old Basis variable 
      VC[Piv_Col] = Basis[Piv_Row]
      Basis[Piv_Row] = Piv_Col
 
 
  if (disp_code >= 1):
      print "Maximum iterations exceeded."
 
  LPcode = 6
 
  return LPcode
 
if __name__ == "__main__":
  M = 10
  N = 4
  L = 8
  E = 0
  G = 2
 
  artifcode = True
  opcode = 1
  if (artifcode):
      JLIM = N + G + M + 1
 
  else:
      JLIM = N + G + L + 1
  displaycode = 3
 
  Basis = range(JLIM)
  AA = matzero(M + 2, JLIM)
 
  matsetblock (AA, 0, 0, M, N,
	       # problem 1-7, pp.7-8 of Bronson[1982] 
	      [100000.0, 1.0, 1.0, 0.0, 0.0,
	       20000.0, 0.0, 0.0, 1.0, 1.0,
	       40000.0, 1.0, 0.0, 1.0, 0.0,
	       60000.0, 0.0, 1.0, 0.0, 1.0,
	       0.0, 1.0, -10.0, 0.0, 0.0,
	       0.0, 0.0, 0.0, 6.0, -5.0,
	       0.0, 2.0, -8.0, 0.0, 0.0,
	       0.0, 0.0, 0.0, 2.0, -8.0,
	       50000.0, 1.0, 1.0, 0.0, 0.0,
	       5000.0, 0.0, 0.0, 1.0, 1.0, 0.0, 
                 4.0, -3.0, 6.0, -1.0])
  print "Input matrix AA"
  matprint(AA)
 
  LPcode = pylinpro (AA, M, N,		# 
	   L, E, G,		# 
	   artifcode,		# process artificial vars?
	   opcode,		# maximization problem
	   Basis, displaycode)
 
  print " LP routine exit code:%d " % LPcode

The matsetblock() function copies the contents of a list to a rectangular part ot the matrix with positions bounded by
the (upper,left) and (lower, right) indices. Here is the code.

?View Code PYTHON
1
2
3
4
5
6
7
8
9
10
11
12
13
14
def matsetblock(AA, top, left, right, bottom, A):
    """
    copies contents of A to the rectangular region of AA
    with indices (top,left)-(right,bottom)
    """
    maxelts = len(A)
    k = 0
    for i in range(left, right+1):
        for j in range(top, bottom+1):
            AA[i][j] = A[k]
            k += 1
            if k >= maxelts:
               return 1
    return 0

When the program is run by typing "python pylinpro.py" on the command line it will output

To browse/scan the genetated output horizontally, place the mouse cursor inside the codebox and click it, then use the left and right
arrow keys.


"""
Tableau at iteration : 9
Basis Value 11 13 12 16 5 8 7 8 9 10 9 12 13 14 15 16
7 61363.636 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.568 -1.250 0.000 0.000 0.114 0.125 0.000 0.000 0.000
6 15000.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -1.000
11 46363.636 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.682 1.500 0.000 1.000 -0.136 -1.250 0.000 0.000 0.000
10 41363.636 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.432 -0.250 1.000 0.000 0.114 0.125 0.000 0.000 0.000
1 30909.091 1.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.455 1.000 0.000 0.000 -0.091 0.000 0.000 0.000 0.000
3 9090.909 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.455 0.000 0.000 0.000 0.091 0.000 0.000 0.000 0.000
2 7727.273 0.000 1.000 0.000 0.000 0.000 0.000 0.000 -0.114 0.250 0.000 0.000 -0.023 -0.125 0.000 0.000 0.000
14 69090.909 0.000 0.000 0.000 0.000 0.000 0.000 0.000 3.455 0.000 0.000 0.000 -0.909 0.000 1.000 0.000 0.000
15 11363.636 0.000 0.000 0.000 0.000 -1.000 0.000 0.000 0.568 -1.250 0.000 0.000 0.114 0.125 0.000 1.000 0.000
4 10909.091 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.545 0.000 0.000 0.000 -0.091 0.000 0.000 0.000 0.000
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
144090.909 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.705 3.250 0.000 0.000 0.341 0.375 0.000 0.000 0.000
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Iteration No. 9
Basis Value
7 61363.63636
6 15000.00000
11 46363.63636
10 41363.63636
1 30909.09091
3 9090.90909
2 7727.27273
14 69090.90909
15 11363.63636
4 10909.09091
The value of the Objective function is 144090.90909
Succesful termination
LP routine exit code:0
"""

Seems to match the objective value obtained from the C/.C++ version.

Hope you find this pure Python linear programming code useful in the educational and research setting.

The new download link (complete with supporting files) is http://adorio-research.org/download/Python/pylinpro.tgz. Download the tarred file, and extract using the command tar xzvvf pylinpro.tgz..

Sep. 28, 2009: I have now 'codeboxed' the code so it is easy to download and browse.

--Ernie

2 Responses to “Python: PyLinpro, a pure simplex tableau solver for linear programming”

  1. government grant scams Says:

    government grant information...

    Hello. You have a good site. You need to visit my free government grants to start a business webpage....

  2. Dave Says:

    Hey...

    Where did you learn to do html?...