Parent Directory
| 
Revision Log
New file with linear programming simplex method related macros.
1 2 =head1 NAME 3 4 LinearProgramming.pl 5 6 =head1 SYNPOSIS 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 21 =head1 DESCRIPTION 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 39 =head2 lp_pivot 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 97 =head2 lp_pivot_element 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 155 =head2 lp_solve 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 223 =head2 lp_current_value 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 266 =head2 lp_display_mm 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 308 =head2 lp_display 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
| aubreyja at gmail dot com | ViewVC Help |
| Powered by ViewVC 1.0.9 |