[system] / trunk / pg / macros / LinearProgramming.pl Repository:
ViewVC logotype

Annotation of /trunk/pg/macros/LinearProgramming.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6058 - (view) (download) (as text)

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9