Parent Directory
| 
Revision Log
test development branch
1 ################################################################################ 2 # WeBWorK Online Homework Delivery System 3 # Copyright © 2000-2007 The WeBWorK Project, http://openwebwork.sf.net/ 4 # $CVSHeader$ 5 # 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 17 =head1 NAME 18 19 LinearProgramming.pl - Macros for the simplex tableau for linear programming 20 problems. 21 22 =head1 SYNPOSIS 23 24 Macros related to the simplex method for Linear Programming. 25 26 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 33 Matrix display makes use of macros from PGmatrixmacros.pl, so it must be 34 included too. 35 36 =head1 DESCRIPTION 37 38 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 42 Tableaus are expected to be legal for the simplex method, such as 43 44 [[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 49 or something similar which arises after pivoting. 50 51 =head1 MACROS 52 53 =head2 lp_pivot 54 55 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 61 $m = [[1,2,3],[4,5,6]]; 62 lp_pivot($m, 0,2); 63 64 This function is destructive - it changes the values of its matrix rather than 65 making a copy, and also returns the matrix, so 66 67 $m = lp_pivot([[1,2,3],[4,5,6]], 0, 2); 68 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 # ^function lp_pivot 77 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 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 121 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 125 =cut 126 127 # Find pivot column for standard part 128 129 # ^function lp_pivot_element 130 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 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 175 =over 176 177 =item C<S<< pivot_limit => 10 >>> 178 179 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 Example: 191 192 $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 197 ($m, $endcode, $pivcount) = lp_solve($m, pivot_limit=>200); 198 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 # ^function lp_solve 206 # ^uses set_default_options 207 # ^uses lp_pivot_element 208 # ^uses lp_pivot 209 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 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 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 # ^function lp_current_value 261 # ^uses Fraction::new 262 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 $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 304 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 308 =cut 309 310 # Display a tableau in math mode 311 # ^function lp_display_mm 312 # ^uses lp_display 313 sub lp_display_mm { 314 lp_display(@_, force_tex=>1); 315 } 316 317 # Make a copy of a tableau 318 # ^function lp_clone 319 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 $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 349 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 353 =cut 354 355 # Display a tableau 356 # ^function lp_display 357 # ^uses lp_clone 358 # ^uses display_matrix 359 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 1;
| aubreyja at gmail dot com | ViewVC Help |
| Powered by ViewVC 1.0.9 |