simplex.pl

00001 #! /usr/bin/perl
00002 # Implementation of Simplex Algorithm
00003 #  Mark Hutchinson 2004
00004 
00005 use CGI qw/:standard *table *center/;
00006 
00007 @tableau=();
00008 @maxmin=();
00009 @names=();
00010 @remove=();
00011 
00012 # Note: vertical counts from 0, horizontal from 1 (to allow space for
00013 # displayed calculations)
00014 
00015 sub data_error;
00016 
00017 # Alas, would that I knew how to import things...
00018 
00019 ######################### Start of fraction library #########################
00020 
00021 # Fraction routines
00022 #  Mark Hutchinson 2004
00023 
00024 sub f_split {
00025     ($_[0] =~ m!^([0-9\-]+)/([0-9\-]+)$!) or return ($_[0],'1');
00026     return ($1,$2);
00027 }
00028 
00029 sub f_join {
00030     return $_[0]."/".$_[1];
00031 }
00032 
00033 sub f_is_neg {
00034     my $neg;
00035     ($_[0] != 0) or return 0;
00036     $neg=$_[0]/abs($_[0]);
00037     $neg *= ($_[1]/abs($_[1]));
00038     return $neg;
00039 }
00040 
00041 sub f_toggle_neg {
00042     my $neg;
00043     my $num=$_[0];
00044     ($num != 0) or return '0/1';
00045     $neg=$num/abs($num);
00046     return ('-'.$num,$_[1]) if $neg==1;
00047     return (substr($num,1),$_[1]);
00048 }
00049 
00050 sub f_simplify {
00051     my @parts=&f_split($_[0]);
00052     &_f_simplify(@parts);
00053 }
00054 
00055 sub _f_simplify {
00056     my $neg,$a,$b,$c,$num,$denom;
00057     ($num,$denom)=@_;
00058     return "0/1" if $num==0;
00059     &data_error("Division by zero") if $denom==0;
00060     $neg=&f_is_neg(@_);
00061     $num=abs $num;
00062     $denom=abs $denom;
00063     if($num<$denom) {
00064         $a=$num;
00065         $b=$denom;
00066     } else {
00067         $a=$denom;
00068         $b=$num;
00069     }
00070     do {
00071         $c=-(int($a/$b)*$b-$a);
00072         if($c!=0){
00073             $a=$b;
00074             $b=$c;
00075         }
00076     } until ($c==0);
00077     $num/=$b;
00078     $denom/=$b;
00079     return &f_join($neg*$num,$denom);
00080 }
00081 
00082 sub f_proper_parts {
00083     my $num,$denom,$count,$neg,$proper;
00084     ($num,$denom)=&f_split($_[0]);
00085     return ('0','0/1') if $num==0;
00086     $neg=&f_is_neg($num,$denom);
00087     $num=abs $num;
00088     $denom= abs $denom;
00089     for($count=0;$num>=$denom;$count++){
00090         $num -= $denom;
00091     }
00092     $count *= $neg;
00093     return ('0',$_[0]) if $count==0;
00094     return ($count,'0/1') if $num==0;
00095     $proper=&f_join($num,$denom);
00096     return ($count,$proper);
00097 }
00098 
00099 sub f_make_proper {
00100     return &f_make_improper($_[0]) if $display_type eq 'Display as improper fractions';
00101     return &f_fract_to_num($_[0]) if $display_type eq 'Display as decimals';
00102     return &_f_make_proper($_[0]);
00103 }
00104 
00105 sub _f_make_proper {
00106     my $int,$fract,$num,$denom,$sign;
00107     if($display_sign eq 'on'){
00108         if(&f_lt($_[0],'0/1')){
00109             $sign='-';
00110         }else{
00111             $sign='+';
00112         }
00113     }
00114     ($int,$fract)=&f_proper_parts($_[0]);
00115     ($num,$denom)=&f_split($fract);
00116     return $int if $num==0;
00117     return '<sup>'.$num.'</sup>/<sub>'.$denom.'</sub>' if $int==0;
00118     return $int.$sign.'<sup>'.$num.'</sup>/<sub>'.$denom.'</sub>';
00119 }
00120 sub f_make_improper {
00121     my $num,$denom;
00122     ($num,$denom)=&f_split($_[0]);
00123     return '0' if ($num==0);
00124     return $num if ($denom==1);
00125     return '<sup>'.$num.'</sup>/<sub>'.$denom.'</sub>';
00126 }
00127 
00128 sub f_add {
00129     my @a,@b,$denom;
00130     @a=&f_split($_[0]);
00131     @b=&f_split($_[1]);
00132     $denom=$a[1]*$b[1];
00133     $a[0] *= $b[1];
00134     $b[0] *= $a[1];
00135     return &_f_simplify($a[0]+$b[0],$denom);
00136 }
00137 
00138 sub f_sub {
00139     my $bn,@b=&f_split($_[1]);
00140     @b=&f_toggle_neg(@b);
00141     $bn=&f_join(@b);
00142     return &f_add($_[0],$bn);
00143 }
00144 
00145 sub f_mult {
00146     my @a,@b;
00147     @a=&f_split($_[0]);
00148     @b=&f_split($_[1]);
00149     return &_f_simplify($a[0]*$b[0],$a[1]*$b[1]);
00150 }
00151 
00152 sub f_div {
00153     my @b=&f_split($_[1]);
00154     return &f_mult($_[0],&f_join($b[1],$b[0]));
00155 }
00156 
00157 sub f_fract_to_num {
00158     my @a=&f_split($_[0]);
00159     return $a[0]/$a[1];
00160 }
00161 
00162 sub f_eq {
00163     my $res=&f_sub($_[0],$_[1]);
00164     return ($res==0);
00165 }
00166 
00167 sub f_ne {
00168     my $res=&f_sub($_[0],$_[1]);
00169     return ($res!=0);
00170 }
00171 
00172 sub f_lt {
00173     my $res=&f_sub($_[0],$_[1]);
00174     return ($res<0);
00175 }
00176     
00177 sub f_le {
00178     my $res=&f_sub($_[0],$_[1]);
00179     return ($res<=0);
00180 }
00181     
00182 sub f_gt {
00183     my $res=&f_sub($_[0],$_[1]);
00184     return ($res>0);
00185 }
00186 
00187 sub f_ge {
00188     my $res=&f_sub($_[0],$_[1]);
00189     return ($res>=0);
00190 }
00191 
00192 sub f_num_to_fract {
00193     my $int,$no,$r,$s,$t,$num,$denom,$s,$a=$_[0];
00194     $int=int $a;
00195     $a -= $int;
00196     do {
00197         $no++;
00198         $r=$a*(10**$no);
00199         $s=int $r;
00200         $t=$r-$s;
00201     } until($t==0 or $no>=8);
00202     $num=$s+($int*(10**$no));
00203     $denom=10**$no;
00204     return &_f_simplify($num,$denom);
00205 }
00206 
00207 sub f_int_to_fract {
00208     return &f_join($_[0],1);
00209 }
00210 
00211 sub f_all_to_fract {
00212 # Handles integers, reals, and proper, improper and mixed fractions
00213     ($_[0] != "") or return "0/1";
00214     ($_[0] !~ m!^[0-9\-]+/[0-9\-]+$!) or return $_[0];
00215     ($_[0] !~ m/^[0-9\-]+$/) or return &f_int_to_fract($_[0]);
00216     
00217     if($_[0] !~ m!^([0-9\-]+)/([0-9\-]+)/([0-9\-]+)$!){
00218         ($_[0] =~ m/^[0-9\-]*\.[0-9\-]+$/) or &data_error("Fraction entered incorrectly");
00219         return &f_num_to_fract($_[0]);
00220     }
00221     return &f_add(&f_int_to_fract($1),&f_join($2,$3));
00222 }
00223 
00224 ########################## End of fraction library ##########################
00225 
00226 sub data_error {
00227     my $arg;
00228     $arg=$_[0];
00229     print header,
00230           start_html("Error: $arg"),
00231           h1("Error: $arg"),
00232           p(
00233           "There is a problem with the data supplied to the script ",
00234           i("simplex.pl"),
00235           "."),
00236           p(
00237           "Please ",
00238           a({-href=>"javascript:history.go(-1)"},"return to the previous page"),
00239           " and try again."),
00240           end_html;
00241           exit;
00242 }
00243 
00244 sub simplex_error {
00245     my $arg=$_[0];
00246     print h3("Error: $arg"),
00247           p(
00248           "There has been an error in applying the simplex algorithm: $arg."),
00249           end_html;
00250           exit;
00251 }
00252 
00253 sub print_tableau {
00254     my $title,$i,$j;
00255     $title=$_[0];
00256     print h3($title),
00257           start_center,
00258           start_table({-width=>'100%',-cols=>$vars+2,-border=>1,-cellspacing=>0,-cellpadding=>0}),
00259           '<tr align="center" valign="top">';
00260     for ($i=0;$i<=$vars;$i++){
00261         print td({-align=>'center'},$names[$i]);
00262     }
00263     print td({-align=>'center'},'RHS'),
00264           '</tr>';
00265     for ($j=0;$j<($objs+$cons);$j++){
00266         print '<tr align="center" valign="top"';
00267         print ' bgcolor="#DDDDDD"' if($j<$objs);
00268         print '>';
00269             print td{-align=>'center'},$tableau[$j][0];
00270         for ($i=1;$i<=$vars+1;$i++){
00271             print td({-align=>'center'},&f_make_proper($tableau[$j][$i]));
00272         }
00273         print '</tr>';
00274     }
00275     print end_table,
00276           end_center;
00277 }
00278 
00279 # Max/min-sensitive improvement function
00280 sub imp_poss {
00281     if($_[0]==1){
00282         return &f_lt($_[1],$_[2]);
00283     }
00284     return &f_gt($_[1],$_[2]);
00285 }
00286 
00287 sub display_solution {
00288     my $a,$one_one,$i,$j,$value,$row;
00289     for($i=0;$i<($objs+$cons);$i++){
00290         $tableau[$i][0]="R<sub>$i</sub>";
00291     }
00292     if($objs>1){
00293         &print_tableau("At end of stage $curr_o");
00294         return;
00295     } else {
00296         &print_tableau("Final tableau");
00297     }
00298     print h3("Solution"),
00299           p;
00300     for($i=1;$i<=$vars;$i++){
00301         $one_one=0;
00302         $value=0;
00303         $row=0;
00304         for($j=0;$j<($objs+$cons);$j++){
00305             $a=$tableau[$j][$i];
00306             ($a != -1) or ($a = -2);
00307             $one_one += abs($a);
00308             $row=$j if $a==1;
00309         }
00310         $value=$tableau[$row][$RHS] if $one_one==1;
00311         $value=&f_make_proper($value);
00312         print "$names[$i] = $value",
00313               br;
00314     }             
00315 }
00316 
00317                   
00318 ############################### Read in data ################################             
00319 &data_error("No data") if(param() eq "");
00320 
00321 ($vars,$objs,$cons)=(param('novars'),param('noobjs'),param('nocons'));
00322 foreach $a ($vars,$objs,$cons){
00323         &data_error("Incomplete data") if $a==0;
00324 }
00325 push @maxmin,"";
00326 for ($i=0;$i<$objs;$i++){
00327     $a=param('m'.$i);
00328     if($a eq 'Maximise'){
00329         push @maxmin,1;
00330     }elsif($a eq 'Minimise'){
00331         push @maxmin,0;
00332     }else{
00333         &data_error("Incomplete data");
00334     }
00335 }
00336 $display_type=param('fdisp');
00337 $display_sign=param('showsign');
00338 push @remove,-1;
00339 push @names,'&nbsp;';
00340 for ($i=1;$i<=$vars;$i++){
00341     $a=param('n'.$i);
00342     while ($a =~ s?\\?<sub>?){};
00343     while ($a =~ s?([^<])/?$1</sub>?){};
00344     push @names,$a;
00345     $a=param('r'.$i);
00346     if($a eq 'Never'){
00347         push @remove,-1;
00348     }else{
00349         push @remove,$a;
00350     }
00351 }    
00352 for ($i=0;$i<($objs+$cons);$i++){
00353     @temp=("R<sub>$i</sub>");
00354     for ($j=1;$j<=$vars+1;$j++){
00355         $temp=param('d'.$i.'_'.$j);
00356         $temp=&f_all_to_fract($temp);
00357         push @temp,$temp;
00358     }
00359     $tableau[$i]=[ @temp ];
00360 }
00361 print header,
00362       start_html('Simplex results'),
00363       h1('Simplex results'),
00364       p;
00365 
00366 ############################## Simplex algorithm ############################
00367 # need for loop for removal of objective rows after
00368 $copy_objs=$objs;
00369 
00370 outer_for: for($curr_o=1;$curr_o<=$copy_objs;$curr_o++,$objs--){
00371 print h2("Stage $curr_o"),
00372       p;
00373 &print_tableau("Initial tableau");
00374 $RHS=$vars+1;
00375 $pivot_no=1;
00376 
00377 outer_while: while(1) {
00378 $temp="0/1";
00379 for($i=$objs+1;$i<=$vars;$i++){
00380     if(&imp_poss($maxmin[$curr_o],$tableau[0][$i],$temp)){
00381         $temp=$tableau[0][$i];
00382         $pivotc=$i;
00383     }
00384 }
00385 if(&f_eq($temp,"0/1")){
00386     &display_solution;
00387     last outer_while;
00388 }
00389 
00390 if (&f_gt($tableau[$objs][$pivotc],"0/1") or (&f_lt($tableau[$objs][$pivotc],"0/1") and &f_lt($tableau[$objs][$RHS],"0/1"))){
00391     $pivotv = &f_div($tableau[$objs][$RHS],$tableau[$objs][$pivotc]);
00392 } else {
00393     $pivotv = '-1/1';
00394 }
00395 $pivotr=$objs;
00396 for($i=$objs+1;$i<($objs+$cons);$i++){
00397     next if(&f_eq($tableau[$i][$pivotc],'0/1'));
00398     $temp=&f_div($tableau[$i][$RHS],$tableau[$i][$pivotc]);
00399     if(&f_gt($temp,'0/1') and (&f_lt($pivotv,'0/1') or &f_lt($temp,$pivotv)))    {
00400         $pivotv=$temp;
00401         $pivotr=$i;
00402     }
00403 }
00404 ($pivotv != '-1/1') or &simplex_error("Unable to find pivot element (no feasible solution)");
00405 $value=$tableau[$pivotr][$pivotc];
00406 
00407 for($i=1;$i<=$RHS;$i++){
00408     $tableau[$pivotr][$i]=&f_div($tableau[$pivotr][$i],$value);
00409 }
00410 $tableau[$pivotr][0]="R<sub>$pivotr</sub>";
00411 for($j=0;$j<($objs+$cons);$j++){
00412     next if $j==$pivotr;
00413     $tosub=$tableau[$j][$pivotc];
00414     $char='- ';
00415     $char='+ ' if &f_lt($tosub,'0/1');
00416     if(not(&f_eq($tosub,'1/1') or &f_eq($tosub,'-1/1'))) {
00417         my $a=$tosub;
00418         $a=&f_join(&f_toggle_neg(&f_split($a))) if &f_lt($a,'0/1');
00419         $char=$char.&f_make_proper($a);
00420     }
00421     $char=$char."R<sub>$pivotr</sub>";
00422     $char="" if &f_eq($tosub,'0/1');
00423     $tableau[$j][0]="R<sub>$j</sub> ".$char;
00424     for($i=1;$i<=$RHS;$i++){
00425         $tableau[$j][$i]=&f_sub($tableau[$j][$i],&f_mult($tosub,$tableau[$pivotr][$i]));
00426     }
00427 }
00428 &print_tableau("Pivot $pivot_no: row R<sub>$pivotr</sub>, column $names[$pivotc].");
00429 $pivot_no++;
00430 }
00431 
00432 @new_remove=(-1);
00433 @new_names=("");
00434 @new_tableau=();
00435 $new_vars=$vars;
00436 
00437 last outer_for if($objs==1);
00438 
00439 if(($tableau[0][$RHS] != 0) and (@maxmin[$curr_o] == 0)){
00440     print p(
00441           "(Note: the first objective function has not reduced to 0.  This suggests that there is no feasible solution to the problem, and that the rest of this algorithm may be invalid.  Caution is necessary in interpreting the remaining stages.)");
00442 }
00443     
00444 
00445 for($i=1;$i<=$vars;$i++){
00446     $new_vars-- if($remove[$i] == $curr_o);
00447 }
00448 
00449 for($j=1;$j<($objs+$cons);$j++){
00450     $k=$j-1;
00451     @temp=("R<sub>$k</sub>");
00452     for($i=1;$i<=$vars+1;$i++){
00453         if($remove[$i] != $curr_o){
00454             push @temp,$tableau[$j][$i];
00455             push @new_names,$names[$i];
00456             push @new_remove,$remove[$i];
00457         }
00458     }
00459     @new_tableau[$j-1] = [ @temp ];
00460 }
00461 
00462 @tableau=@new_tableau;
00463 @remove=@new_remove;
00464 @names=@new_names;
00465 $vars=$new_vars;
00466 
00467 }
00468 $url=url(-relative=>1,-query=>1);
00469 $url =~ s/simplex\.pl/simplex_table\.pl/;
00470 
00471 print   p(
00472         a({-href=>"$url"},"Alter initial tableau and reprocess"),br,
00473         a({-href=>"../simplex.html"},"Return to the simplex page")),
00474         hr,
00475         font({-size=>'-1'},"This implementation of the Simplex algorithm is &copy; Mark Hutchinson 2004."),
00476         end_html;
00477         

Generated on Tue May 25 20:02:42 2004 for Simplex by doxygen 1.3.3