simplex.pl

#! /usr/bin/perl
# Implementation of Simplex Algorithm
# © Mark Hutchinson 2004

use CGI qw/:standard *table *center/;

@tableau=();
@maxmin=();
@names=();
@remove=();

# Note: vertical counts from 0, horizontal from 1 (to allow space for
# displayed calculations)

sub data_error;

# Alas, would that I knew how to import things...

######################### Start of fraction library #########################

# Fraction routines
# © Mark Hutchinson 2004

sub f_split {
    ($_[0] =~ m!^([0-9\-]+)/([0-9\-]+)$!) or return ($_[0],'1');
    return ($1,$2);
}

sub f_join {
    return $_[0]."/".$_[1];
}

sub f_is_neg {
    my $neg;
    ($_[0] != 0) or return 0;
    $neg=$_[0]/abs($_[0]);
    $neg *= ($_[1]/abs($_[1]));
    return $neg;
}

sub f_toggle_neg {
    my $neg;
    my $num=$_[0];
    ($num != 0) or return '0/1';
    $neg=$num/abs($num);
    return ('-'.$num,$_[1]) if $neg==1;
    return (substr($num,1),$_[1]);
}

sub f_simplify {
    my @parts=&f_split($_[0]);
    &_f_simplify(@parts);
}

sub _f_simplify {
    my $neg,$a,$b,$c,$num,$denom;
    ($num,$denom)=@_;
    return "0/1" if $num==0;
    &data_error("Division by zero") if $denom==0;
    $neg=&f_is_neg(@_);
    $num=abs $num;
    $denom=abs $denom;
    if($num<$denom) {
        $a=$num;
        $b=$denom;
    } else {
        $a=$denom;
        $b=$num;
    }
    do {
        $c=-(int($a/$b)*$b-$a);
        if($c!=0){
            $a=$b;
            $b=$c;
        }
    } until ($c==0);
    $num/=$b;
    $denom/=$b;
    return &f_join($neg*$num,$denom);
}

sub f_proper_parts {
    my $num,$denom,$count,$neg,$proper;
    ($num,$denom)=&f_split($_[0]);
    return ('0','0/1') if $num==0;
    $neg=&f_is_neg($num,$denom);
    $num=abs $num;
    $denom= abs $denom;
    for($count=0;$num>=$denom;$count++){
        $num -= $denom;
    }
    $count *= $neg;
    return ('0',$_[0]) if $count==0;
    return ($count,'0/1') if $num==0;
    $proper=&f_join($num,$denom);
    return ($count,$proper);
}

sub f_make_proper {
    return &f_make_improper($_[0]) if $display_type eq 'Display as improper fractions';
    return &f_fract_to_num($_[0]) if $display_type eq 'Display as decimals';
    return &_f_make_proper($_[0]);
}

sub _f_make_proper {
    my $int,$fract,$num,$denom,$sign;
    if($display_sign eq 'on'){
        if(&f_lt($_[0],'0/1')){
            $sign='-';
        }else{
            $sign='+';
        }
    }
    ($int,$fract)=&f_proper_parts($_[0]);
    ($num,$denom)=&f_split($fract);
    return $int if $num==0;
    return '<sup>'.$num.'</sup>/<sub>'.$denom.'</sub>' if $int==0;
    return $int.$sign.'<sup>'.$num.'</sup>/<sub>'.$denom.'</sub>';
}
sub f_make_improper {
    my $num,$denom;
    ($num,$denom)=&f_split($_[0]);
    return '0' if ($num==0);
    return $num if ($denom==1);
    return '<sup>'.$num.'</sup>/<sub>'.$denom.'</sub>';
}

sub f_add {
    my @a,@b,$denom;
    @a=&f_split($_[0]);
    @b=&f_split($_[1]);
    $denom=$a[1]*$b[1];
    $a[0] *= $b[1];
    $b[0] *= $a[1];
    return &_f_simplify($a[0]+$b[0],$denom);
}

sub f_sub {
    my $bn,@b=&f_split($_[1]);
    @b=&f_toggle_neg(@b);
    $bn=&f_join(@b);
    return &f_add($_[0],$bn);
}

sub f_mult {
    my @a,@b;
    @a=&f_split($_[0]);
    @b=&f_split($_[1]);
    return &_f_simplify($a[0]*$b[0],$a[1]*$b[1]);
}

sub f_div {
    my @b=&f_split($_[1]);
    return &f_mult($_[0],&f_join($b[1],$b[0]));
}

sub f_fract_to_num {
    my @a=&f_split($_[0]);
    return $a[0]/$a[1];
}

sub f_eq {
    my $res=&f_sub($_[0],$_[1]);
    return ($res==0);
}

sub f_ne {
    my $res=&f_sub($_[0],$_[1]);
    return ($res!=0);
}

sub f_lt {
    my $res=&f_sub($_[0],$_[1]);
    return ($res<0);
}
    
sub f_le {
    my $res=&f_sub($_[0],$_[1]);
    return ($res<=0);
}
    
sub f_gt {
    my $res=&f_sub($_[0],$_[1]);
    return ($res>0);
}

sub f_ge {
    my $res=&f_sub($_[0],$_[1]);
    return ($res>=0);
}

sub f_num_to_fract {
    my $int,$no,$r,$s,$t,$num,$denom,$s,$a=$_[0];
    $int=int $a;
    $a -= $int;
    do {
        $no++;
        $r=$a*(10**$no);
        $s=int $r;
        $t=$r-$s;
    } until($t==0 or $no>=8);
    $num=$s+($int*(10**$no));
    $denom=10**$no;
    return &_f_simplify($num,$denom);
}

sub f_int_to_fract {
    return &f_join($_[0],1);
}

sub f_all_to_fract {
# Handles integers, reals, and proper, improper and mixed fractions
    ($_[0] != "") or return "0/1";
    ($_[0] !~ m!^[0-9\-]+/[0-9\-]+$!) or return $_[0];
    ($_[0] !~ m/^[0-9\-]+$/) or return &f_int_to_fract($_[0]);
    
    if($_[0] !~ m!^([0-9\-]+)/([0-9\-]+)/([0-9\-]+)$!){
        ($_[0] =~ m/^[0-9\-]*\.[0-9\-]+$/) or &data_error("Fraction entered incorrectly");
        return &f_num_to_fract($_[0]);
    }
    return &f_add(&f_int_to_fract($1),&f_join($2,$3));
}

########################## End of fraction library ##########################

sub data_error {
    my $arg;
    $arg=$_[0];
    print header,
          start_html("Error: $arg"),
          h1("Error: $arg"),
          p(
          "There is a problem with the data supplied to the script ",
          i("simplex.pl"),
          "."),
          p(
          "Please ",
          a({-href=>"javascript:history.go(-1)"},"return to the previous page"),
          " and try again."),
          end_html;
          exit;
}

sub simplex_error {
    my $arg=$_[0];
    print h3("Error: $arg"),
          p(
          "There has been an error in applying the simplex algorithm: $arg."),
          end_html;
          exit;
}

sub print_tableau {
    my $title,$i,$j;
    $title=$_[0];
    print h3($title),
          start_center,
          start_table({-width=>'100%',-cols=>$vars+2,-border=>1,-cellspacing=>0,-cellpadding=>0}),
          '<tr align="center" valign="top">';
    for ($i=0;$i<=$vars;$i++){
        print td({-align=>'center'},$names[$i]);
    }
    print td({-align=>'center'},'RHS'),
          '</tr>';
    for ($j=0;$j<($objs+$cons);$j++){
        print '<tr align="center" valign="top"';
        print ' bgcolor="#DDDDDD"' if($j<$objs);
        print '>';
            print td{-align=>'center'},$tableau[$j][0];
        for ($i=1;$i<=$vars+1;$i++){
            print td({-align=>'center'},&f_make_proper($tableau[$j][$i]));
        }
        print '</tr>';
    }
    print end_table,
          end_center;
}

# Max/min-sensitive improvement function
sub imp_poss {
    if($_[0]==1){
        return &f_lt($_[1],$_[2]);
    }
    return &f_gt($_[1],$_[2]);
}

sub display_solution {
    my $a,$one_one,$i,$j,$value,$row;
    for($i=0;$i<($objs+$cons);$i++){
        $tableau[$i][0]="R<sub>$i</sub>";
    }
    if($objs>1){
        &print_tableau("At end of stage $curr_o");
        return;
    } else {
        &print_tableau("Final tableau");
    }
    print h3("Solution"),
          p;
    for($i=1;$i<=$vars;$i++){
        $one_one=0;
        $value=0;
        $row=0;
        for($j=0;$j<($objs+$cons);$j++){
            $a=$tableau[$j][$i];
            ($a != -1) or ($a = -2);
            $one_one += abs($a);
            $row=$j if $a==1;
        }
        $value=$tableau[$row][$RHS] if $one_one==1;
        $value=&f_make_proper($value);
        print "$names[$i] = $value",
              br;
    }             
}

                  
############################### Read in data ################################             
&data_error("No data") if(param() eq "");

($vars,$objs,$cons)=(param('novars'),param('noobjs'),param('nocons'));
foreach $a ($vars,$objs,$cons){
        &data_error("Incomplete data") if $a==0;
}
push @maxmin,"";
for ($i=0;$i<$objs;$i++){
    $a=param('m'.$i);
    if($a eq 'Maximise'){
        push @maxmin,1;
    }elsif($a eq 'Minimise'){
        push @maxmin,0;
    }else{
        &data_error("Incomplete data");
    }
}
$display_type=param('fdisp');
$display_sign=param('showsign');
push @remove,-1;
push @names,'&nbsp;';
for ($i=1;$i<=$vars;$i++){
    $a=param('n'.$i);
    while ($a =~ s?\\?<sub>?){};
    while ($a =~ s?([^<])/?$1</sub>?){};
    push @names,$a;
    $a=param('r'.$i);
    if($a eq 'Never'){
        push @remove,-1;
    }else{
        push @remove,$a;
    }
}    
for ($i=0;$i<($objs+$cons);$i++){
    @temp=("R<sub>$i</sub>");
    for ($j=1;$j<=$vars+1;$j++){
        $temp=param('d'.$i.'_'.$j);
        $temp=&f_all_to_fract($temp);
        push @temp,$temp;
    }
    $tableau[$i]=[ @temp ];
}
print header,
      start_html('Simplex results'),
      h1('Simplex results'),
      p;

############################## Simplex algorithm ############################
# need for loop for removal of objective rows after
$copy_objs=$objs;

outer_for: for($curr_o=1;$curr_o<=$copy_objs;$curr_o++,$objs--){
print h2("Stage $curr_o"),
      p;
&print_tableau("Initial tableau");
$RHS=$vars+1;
$pivot_no=1;

outer_while: while(1) {
$temp="0/1";
for($i=$objs+1;$i<=$vars;$i++){
    if(&imp_poss($maxmin[$curr_o],$tableau[0][$i],$temp)){
        $temp=$tableau[0][$i];
        $pivotc=$i;
    }
}
if(&f_eq($temp,"0/1")){
    &display_solution;
    last outer_while;
}

if (&f_gt($tableau[$objs][$pivotc],"0/1") or (&f_lt($tableau[$objs][$pivotc],"0/1") and &f_lt($tableau[$objs][$RHS],"0/1"))){
    $pivotv = &f_div($tableau[$objs][$RHS],$tableau[$objs][$pivotc]);
} else {
    $pivotv = '-1/1';
}
$pivotr=$objs;
for($i=$objs+1;$i<($objs+$cons);$i++){
    next if(&f_eq($tableau[$i][$pivotc],'0/1'));
    $temp=&f_div($tableau[$i][$RHS],$tableau[$i][$pivotc]);
    if(&f_gt($temp,'0/1') and (&f_lt($pivotv,'0/1') or &f_lt($temp,$pivotv)))    {
        $pivotv=$temp;
        $pivotr=$i;
    }
}
($pivotv != '-1/1') or &simplex_error("Unable to find pivot element (no feasible solution)");
$value=$tableau[$pivotr][$pivotc];

for($i=1;$i<=$RHS;$i++){
    $tableau[$pivotr][$i]=&f_div($tableau[$pivotr][$i],$value);
}
$tableau[$pivotr][0]="R<sub>$pivotr</sub>";
for($j=0;$j<($objs+$cons);$j++){
    next if $j==$pivotr;
    $tosub=$tableau[$j][$pivotc];
    $char='- ';
    $char='+ ' if &f_lt($tosub,'0/1');
    if(not(&f_eq($tosub,'1/1') or &f_eq($tosub,'-1/1'))) {
        my $a=$tosub;
        $a=&f_join(&f_toggle_neg(&f_split($a))) if &f_lt($a,'0/1');
        $char=$char.&f_make_proper($a);
    }
    $char=$char."R<sub>$pivotr</sub>";
    $char="" if &f_eq($tosub,'0/1');
    $tableau[$j][0]="R<sub>$j</sub> ".$char;
    for($i=1;$i<=$RHS;$i++){
        $tableau[$j][$i]=&f_sub($tableau[$j][$i],&f_mult($tosub,$tableau[$pivotr][$i]));
    }
}
&print_tableau("Pivot $pivot_no: row R<sub>$pivotr</sub>, column $names[$pivotc].");
$pivot_no++;
}

@new_remove=(-1);
@new_names=("");
@new_tableau=();
$new_vars=$vars;

last outer_for if($objs==1);

if(($tableau[0][$RHS] != 0) and (@maxmin[$curr_o] == 0)){
    print p(
          "(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.)");
}
    

for($i=1;$i<=$vars;$i++){
    $new_vars-- if($remove[$i] == $curr_o);
}

for($j=1;$j<($objs+$cons);$j++){
    $k=$j-1;
    @temp=("R<sub>$k</sub>");
    for($i=1;$i<=$vars+1;$i++){
        if($remove[$i] != $curr_o){
            push @temp,$tableau[$j][$i];
            push @new_names,$names[$i];
            push @new_remove,$remove[$i];
        }
    }
    @new_tableau[$j-1] = [ @temp ];
}

@tableau=@new_tableau;
@remove=@new_remove;
@names=@new_names;
$vars=$new_vars;

}
$url=url(-relative=>1,-query=>1);
$url =~ s/simplex\.pl/simplex_table\.pl/;

print   p(
        a({-href=>"$url"},"Alter initial tableau and reprocess"),br,
        a({-href=>"../simplex.html"},"Return to the simplex page")),
        hr,
        font({-size=>'-1'},"This implementation of the Simplex algorithm is &copy; Mark Hutchinson 2004."),
        end_html;
        

---