[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 1453 - (download) (as text) (annotate)
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 
    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