[system] / trunk / pg / macros / LinearProgramming.pl Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

# View of /trunk/pg/macros/LinearProgramming.pl

Wed Aug 6 18:59:26 2003 UTC (16 years, 6 months ago) by jj
File size: 9981 byte(s)
```New file with linear programming simplex method related macros.
```

```    1
3
4        LinearProgramming.pl
5
7
8        Macros related to the simplex method for Linear Programming.
9
10              lp_pivot_element()  - find the pivot element from a tableau
11              lp_solve()          - pivot until done
12              lp_current_value()  - given a tableau, find the value of a
13                                    requested variable
14              lp_display()        - display a tableau
15              lp_display_mm()     - display a tableau while in math mode
16              lp_pivot()          - perform one pivot on a tableau
17
18        Matrix display makes use of macros from PGmatrixmacros.pl, so it
19        must be included too.
20
22
23     These are macros for dealing with simplex tableau for linear
24     programming problems.  The tableau is a reference to an array
25     of arrays, which looks like, [[1,2,3], [4,5,6]].  The entries
26     can be real numbers or Fractions.
27
28     Tableaus are expected to be legal for the simplex method, such as
29
30     [[0,  3, -4,  5, 1, 0, 0, 28],
31      [0,  2,  0,  1, 0, 1, 0, 11],
32      [0,  1,  2,  3, 0, 0, 1,  3],
33      [1, -1,  2, -3, 0, 0, 0,  0]]
34
35     or something similar which arises after pivoting.
36
37 =cut
38
40
41 Take a tableau, the row and column number, and perform one pivot operation.
42 The tableau can be any matrix in the form reference to array of arrays, such
43 as [[1,2,3],[4,5,6]].  Row and column numbers start at 0.  An optional 4th
44 argument can be specified to indicate that the matrix has entries of type
45 Fraction (and then all entries should be of type Fraction).
46
47   \$m = [[1,2,3],[4,5,6]];
48   lp_pivot(\$m, 0,2);
49
50 This function is destructive - it changes the values of its matrix rather
51 than making a copy, and also returns the matrix, so
52
53   \$m = lp_pivot([[1,2,3],[4,5,6]], 0, 2);
54
55 will have the same result as the example above.
56
57 =cut
58
59
60 # perform a pivot operation
61 # lp_pivot([[1,2,3],...,[4,5,6]], row, col, fractionmode)
62 # row and col indecies start at 0
63 sub lp_pivot {
64   my \$a_ref = shift;
65   my \$row = shift;
66   my \$col = shift;
67   my \$fracmode = shift;
68   \$fracmode = 0 unless defined(\$fracmode);
69   my @matrows = @{\$a_ref};
70
71   if((\$fracmode and \$matrows[\$row][\$col]->scalar() == 0)
72      or \$matrows[\$row][\$col] == 0) {
73     warn "Pivoting a matrix on a zero element";
74     return(\$a_ref);
75   }
76   my (\$j, \$k, \$hld);
77   \$hld = \$matrows[\$row][\$col];
78   for \$j (1..scalar(@{\$matrows[0]})) {
79     \$matrows[\$row][\$j-1] = \$fracmode ? \$matrows[\$row][\$j-1]->divBy(\$hld) :
80       \$matrows[\$row][\$j-1]/\$hld;
81   }
82   for \$k (1..scalar(@matrows)) {
83     if(\$k-1 != \$row) {
84       \$hld = \$matrows[\$k-1][\$col];
85       for \$j (1..scalar(@{\$matrows[0]})) {
86         \$matrows[\$k-1][\$j-1] = \$fracmode ?
87           \$matrows[\$k-1][\$j-1]->minus(\$matrows[\$row][\$j-1]->times(\$hld)) :
88           \$matrows[\$k-1][\$j-1] - \$matrows[\$row][\$j-1]*\$hld;
89       }
90     }
91   }
92
93   return(\$a_ref);
94 }
95
96
98
99 Take a simplex tableau, and determine which element is the next pivot
100 element based on the algorithm in Mizrahi and Sullivan's Finite
101 Mathematics, section 4.2.  The tableau must represent a point in the
102 region of feasibility for a LP problem.  Otherwise, it can be any
103 matrix in the form reference to array of arrays, such as
104 [[1,2,3],[4,5,6]].  An optional 2nd argument can be specified to
105 indicate that the matrix has entries of type Fraction (and then all
106 entries should be of type Fraction).
107
108 It returns a pair [row, col], with the count starting at 0.  If there
109 is no legal pivot column (final tableau), it returns [-1,-1].  If
110 there is a column, but no pivot element (unbounded problem), it returns
111 [-1, col].
112
113 =cut
114
115 # Find pivot column for standard part
116
117 sub lp_pivot_element {
118   my \$a_ref = shift;
119   my \$fracmode = shift;
120   \$fracmode = 0 unless defined(\$fracmode);
121   my @m = @{\$a_ref};
122   my \$nrows = scalar(@m)-1; # really 1 less
123   my \$ncols = scalar(@{\$m[0]}) -1;  # really 1 less
124   # looking for minimum value
125   my \$minv = \$fracmode ? \$m[\$nrows][0]->scalar() : \$m[\$nrows][0];
126   my \$pcol=0;
127   my \$j;
128   for \$j (1..(\$ncols-1)) {
129     if((\$fracmode and \$m[\$nrows][\$j]->scalar()<\$minv) or \$m[\$nrows][\$j]<\$minv) {
130       \$minv = \$fracmode ? \$m[\$nrows][\$j]->scalar() :\$m[\$nrows][\$j];
131       \$pcol = \$j;
132     }
133   }
134   return (-1, -1) if (\$minv>=0); # This means we are done
135   # Now find the pivot row
136   my \$prow=-1;
137   for \$j (0..(\$nrows-1)) {
138     if(\$fracmode ? (\$m[\$j][\$pcol]->scalar() >0) : (\$m[\$j][\$pcol]>0)) { # found a candidate
139       if(\$prow == -1) {
140         \$prow = \$j;
141       } else { # Test to see if this is an improvement
142         if(\$fracmode ?
143             (\$m[\$prow][\$ncols]->scalar())/(\$m[\$prow][\$pcol]->scalar()) >
144                        (\$m[\$j][\$ncols]->scalar()/\$m[\$j][\$pcol]->scalar())
145            : (\$m[\$prow][\$ncols]/\$m[\$prow][\$pcol] >
146            \$m[\$j][\$ncols]/\$m[\$j][\$pcol])) {
147           \$prow = \$j;
148         }
149       }
150     }
151   }
152   return (\$prow, \$pcol);
153 }
154
156
157 Take a tableau, and perform simplex method pivoting until done.
158 The tableau can be any matrix in the form reference to array of arrays, such
159 as [[1,2,3],[4,5,6]], which represents a linear programming tableau at a
160 feasible point.  Options are specified in key/value pairs.
161
162    pivot_limit=> 10   (limit the number of pivots to at most 10 - default is 100)
163    fraction_mode=> 1  (entries are of type Fraction - defaults to 0, i.e., false)
164
165 This function is destructive - it changes the values of its matrix
166 rather than making a copy.  It returns a triple of the final tableau,
167 an endcode indicating the type of result, and the number of pivots
168 used.  The endcodes are 1 for success, 0 for unbounded.
169
170 Example:
171
172 \$m = [[0,  3, -4,  5, 1, 0, 0, 28],
173       [0,  2,  0,  1, 0, 1, 0, 11],
174       [0,  1,  2,  3, 0, 0, 1,  3],
175       [1, -1,  2, -3, 0, 0, 0,  0]];
176
177 (\$m, \$endcode, \$pivcount) = lp_solve(\$m, pivot_limit=>200);
178
179 =cut
180
181
182 # Solve a linear programming problem
183 # lp_solve([[1,2,3],[4,5,6]])
184 # It returns a triple of the final tableau, a code to say if we
185 #   succeeded, and the number of pivots
186 sub lp_solve {
187   my \$a_ref_orig = shift;
188   my %opts = @_;
189
190   set_default_options(\%opts,
191                       '_filter_name' => 'lp_solve',
192                       'pivot_limit' => 100,
193                       'fraction_mode' => 0,
194                       'allow_unknown_options'=> 0);
195
196   my (\$pcol, \$prow);
197   my \$a_ref;
198   # First we clone the matrix so that it isn't modified in place
199   for \$prow (1..scalar(@{\$a_ref_orig})) {
200     for \$pcol (1..scalar(@{\$a_ref_orig->[0]})) {
201       \$a_ref->[\$prow-1][\$pcol-1] = \$a_ref_orig->[\$prow-1][\$pcol-1];
202     }
203   }
204
205   my \$piv_count = 0;
206   my \$piv_limit = \$opts{'pivot_limit'}; # Just in case of cycling or bugs
207   my \$fracmode = \$opts{'fraction_mode'};
208   # First do alternate pivoting
209   # Now do regular pivots
210   do {
211     (\$prow, \$pcol) = lp_pivot_element(\$a_ref, \$fracmode);
212     if(\$prow>=0) {
213       \$a_ref = lp_pivot(\$a_ref, \$prow, \$pcol, \$fracmode);
214       \$piv_count++;
215     }
216   } until(\$prow<0 or \$piv_count>=\$piv_limit);
217   # code is 1 for success, 0 for unbounded
218   my \$endcode = 1;
219   \$endcode = 0 if (\$pcol>=0);
220   return(\$a_ref, \$endcode, \$piv_count);
221 }
222
224
225 Takes a simplex tableau and returns the value of a particular variable.
226 Variables are associated to column numbers which are indexed starting with
227 0.  So, usually this means that the objective function is 0, x_1 is 1, and
228 so on.  This can be used for slack variables too (assuming you know what
229 columns they are in).
230
231 =cut
232
233 # Get the current value of a variable from a tableau
234 # The variable is specified by column number, with 0 for P, 1 for x_1,
235 #  and so on
236
237 sub lp_current_value {
238   my \$col = shift;
239   my \$aref = shift;
240   my \$fractionmode = 0;
241   \$fractionmode =1 if(ref(\$aref->[0][0]) eq 'Fraction');
242
243   # Count how many ones there are.  If we hit non-zero/non-one, force count
244   # to fail
245   my (\$cnt,\$row,\$save) = (0,'',0);
246   for \$row (@{\$aref}) {
247     if(\$fractionmode) {
248       if(\$row->[\$col]->scalar() != 0) {
249   \$cnt += (\$row->[\$col]->scalar() == 1) ? 1 : 2;
250   \$save = \$row;
251       }
252     } else {
253       if(\$row->[\$col] != 0) {
254   \$cnt += (\$row->[\$col] == 1) ? 1 : 2;
255   \$save = \$row;
256       }
257     }
258   }
259   if(\$cnt != 1) {
260     return (\$fractionmode ? new Fraction(0) : 0);
261   }
262   \$cnt = scalar(@{\$save});
263   return \$save->[\$cnt-1];
264 }
265
267
268 Display a simplex tableau while in math mode.
269
270 \$m = [[0,  3, -4,  5, 1, 0, 0, 28],
271       [0,  2,  0,  1, 0, 1, 0, 11],
272       [0,  1,  2,  3, 0, 0, 1,  3],
273       [1, -1,  2, -3, 0, 0, 0,  0]];
274
275 \[ \{ lp_display_mm(\$m) \} \]
276
277 Accepts the same optional arguments as lp_display (see below), and
278 produces nicer looking results.  However, it cannot have answer rules
279 in the tableau (lp_display can have them for fill in the blank
280 tableaus).
281
282 =cut
283
284 # Display a tableau in math mode
285
286 sub lp_display_mm {
287   lp_display(@_, force_tex=>1);
288 }
289
290 # Make a copy of a tableau
291 sub lp_clone {
292   my \$a1_ref = shift;
293         my \$a_ref = []; # make a copy to modify
294         my \$nrows = scalar(@{\$a1_ref})-1; # really 1 less
295   my (\$j, \$k);
296   for \$j (0..\$nrows) {
297     if(\$a1_ref->[\$j] eq 'hline') {
298       \$a_ref->[\$j] = 'hline';
299     } else {
300       for \$k (0..(scalar(@{\$a1_ref->[\$j]}) -1)) {
301         \$a_ref->[\$j][\$k] = \$a1_ref->[\$j][\$k];
302       }
303     }
304   }
305   return(\$a_ref);
306 }
307
309
310 Display a simplex tableau while not in math mode.
311
312 \$m = [[0,  3, -4,  5, 1, 0, 0, 28],
313       [0,  2,  0,  1, 0, 1, 0, 11],
314       [0,  1,  2,  3, 0, 0, 1,  3],
315       [1, -1,  2, -3, 0, 0, 0,  0]];
316
317 \{ lp_display(\$m)\}
318
319 Takes the same optional arguments as display_matrix.  The default
320 for column alignment as "augmentation line" before the last column.
321 It also adds a horizontal line before the last row if it is not already
322 specified.
323
324 =cut
325
326
327 # Display a tableau
328 sub lp_display {
329   my \$a1_ref = shift;
330   my %opts = @_;
331   my \$nrows = scalar(@{\$a1_ref})-1; # really 1 less
332   my \$ncols = scalar(@{\$a1_ref->[0]}) -1;  # really 1 less
333
334   if(not defined(\$opts{'align'})) {
335     my \$align = "r" x \$ncols;
336     \$align .=  "|r";
337     \$opts{'align'} = \$align;
338   }
339   my \$a_ref = lp_clone(\$a1_ref);
340
341   if(\$a_ref->[\$nrows-1] ne 'hline') {
342     \$a_ref->[\$nrows+1] = \$a_ref->[\$nrows];
343     \$a_ref->[\$nrows] = 'hline';
344   }
345   display_matrix(\$a_ref, %opts);
346 }
347
348 # return 1 so that this file can be included with require
349 1
```