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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5658 - (download) (as text) (annotate)
Sat May 3 17:43:29 2008 UTC (11 years, 8 months ago) by sh002i
File size: 11197 byte(s)
markup for ww-symbol-map

    1 ################################################################################
    2 # WeBWorK Online Homework Delivery System
    3 # Copyright  2000-2007 The WeBWorK Project, http://openwebwork.sf.net/
    4 # $CVSHeader: pg/macros/LinearProgramming.pl,v 1.2 2007/10/25 17:11:59 sh002i Exp $
    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