#  September 5, 2008. MAPLE code to pivot on a dictionary.
#  Updated version. March 11, 2012 gives tableaux in ASWCM format.
#
#  Usage:
#     Set Up: 
#      We must first specify an m x n matrix, A, of constraint coeffs.
#                            an m x 1 matrix, b, of RHS values (assumed >= 0)
#                            a  1 x n matrix, c, of objective fn coeffs.
#            and an initial value for z (the objective value), usually zero.
#
# Here, we provide three MAPLE functions:
#
#     ratios: ratios(j) determines all upper bounds for entering variable xj.
#    
#     lppivot:  lppivot(i,j) modifies A,b,c,z by pivoting on eqn i, variable j 
#
#      dict:  displays the current A, b, c and z in dictionary format.
#
myeps := .1^8:


mm :=  linalg[rowdim](A);
nn :=  linalg[coldim](A);

origc := matrix(1,nn,[seq(c[1,hh],hh=1..nn)]):
 ### This is a bad idea, but let's assume it's initially in standard form.
 if nops(lhslist)=0 
 then 
   lhslist := [seq(h,h=nn+1-mm..nn )]; 
   basrow := [seq(h,h=1..mm) ];
 fi:


lppivot := proc(ext, ent)  # Manipulate global variables A,b,c,z by
                        # pivoting on row rw where exiting variable x_ext is basic 
                        # and column ent of current tableau.
local ii,jj,pv,pivcolval, newbasis, rw, cl, h;
global A, b, c, z, mm,nn, basrow, lhslist;

cl := ent;  # Entering variable gives the pivot column of matrix A

# Need to search for pivot row:
rw:=0;
for ii to mm
  do 
     if abs( A[ii,ext] - 1.) < myeps 
     then 
       rw := ii;  ii := mm;
     fi;
  od;
for ii to mm
  do 
     if ii <> rw and abs( A[ii,ext] ) > myeps 
     then 
       rw := 0;  ii := mm;
     fi;
  od;
if rw=0 then printf(`ERROR (lppivot): Chosen exiting variable x%d is NOT basic.\n`,ext);
else

pv := A[rw,cl]; # Store value of pivot before modifying it.
if pv=0 then printf(`ERROR (lppivot): Division by zero.\n`);
else

  ####
  ####  Perform the pivot on the data.
  ####

  # Scale the pivot row.
  for jj from 1 to nn
   do
       A[rw,jj]  := A[rw,jj]/pv:
   od:
       b[rw,1] := b[rw,1]/pv:

  # Add appropriate multiples of pivot row to all other rows
  for ii from 1 to mm
   do
    if ii <> rw
    then
     pivcolval := A[ii,cl]:
     for jj from 1 to nn
      do
          A[ii,jj]  := A[ii,jj] - pivcolval*A[rw,jj]:
      od:
          b[ii,1]   :=  b[ii,1] - pivcolval*b[rw,1]:
    fi:
   od:
   
  # Add appropriate multiple of pivot row to reduced costs and obj. value.
  pivcolval := c[1,cl]:
  for jj from 1 to nn
   do
          c[1,jj]  := c[1,jj] - pivcolval*A[rw,jj]:
   od:
          z := z + pivcolval*b[rw,1]:

  for h to mm   # Swap variable out of basis
   do
    if lhslist[h] = ext then lhslist[h] := ent;   fi
   od;

   newbasis := [];
   for h to mm
    do 
      newbasis := [ op(newbasis), x||(lhslist[h]) ];
    od;

   tableau( op( newbasis ) );

fi:  # If pivot non-zero.
fi:  # If pivot row exists.
     
end:
###############################################################################


ratios := proc(jj)  # Give ratios for entering variable x.jj
local ii;
global A, b, mm, lhslist, basrow;

printf(`\n`);
for ii to mm
 do
   if A[basrow[ii],jj] <= 0 
      then 
           printf(`    row %2d (x%d basic):  No upper bound \n`,ii, lhslist[ii]);
      else 
           printf(`    row %2d (x%d basic):  Upper bound = %8.4f \n`,
                                     ii, lhslist[ii], 
                                     evalf(b[ basrow[ii],1]/A[ basrow[ii],jj]) );
   fi:
 od:
printf(`\n`);
end:
###############################################################################

dict := proc()
local ii, jj, goodbasisflag;
global A, b, c, z, mm,nn, lhslist, basrow;
 
   ####
   #### First figure out what basis the user wants
   ####
   #### (Note: this is an ugly bit of code.)
   if nargs <> mm
   then
   printf(`ERROR (dict): user must provide an ordered list of %d basic variables.\n`, mm);
   else
      lhslist := [];
      for ii to nargs 
      do
       for jj to nn
       do
          if args[ii] = x||jj then lhslist := [op(lhslist),jj]; jj:=nn;  fi;
       od
      od;
      #### Now check if it really is a current basis:
      basrow := [];
      goodbasisflag := true;
      for ii to nargs
       do
          for jj to mm
           do
              if abs( A[jj,lhslist[ii]] - 1.) < myeps
               then
                  if nops( basrow ) < ii 
                  then  # good! we've found the entry we want.
                     basrow := [ op(basrow), jj];
                  else # Uh oh! We've already got one non-zero entry in that column!
                      goodbasisflag := false; 
                      jj:=mm;
                      ii:=nargs;
                  fi;
               else # If there is any other nonzero entry in this column, we're dead.
                  if abs( A[jj, lhslist[ii]] ) > myeps
                  then goodbasisflag := false; 
                         jj:=mm; ii:=nargs;
                  fi;
               fi;  # if entry is 1
           od;  # for jj
       od; # for ii
   
     if goodbasisflag
     then
      ####
      ####  Display the current dictionary
      ####
    
      printf(`\n`);
      # Display objective function.
      printf(`     z   = %8.2f `,evalf( z ) );
      
      for jj to nn
       do
          if not member( jj , {op(lhslist)})
          then
              if c[1,jj] > myeps 
              then
                printf(` + %6.2f x%d `,evalf(c[1,jj]), jj);
              else
                if c[1,jj] < -myeps 
                then
                  printf(` - %6.2f x%d `,evalf(-c[1,jj]), jj);
                else
                  printf(`              `);
                fi       
              fi       
          fi
       od;

      ## Dashed line between objective function and constraints
      printf(`\n    -----------------`);
      for jj to nn
       do
          if not member( jj , {op(lhslist)})
          then
                  printf(`--------------`);
          fi
       od;
       printf(`\n`);
   
   
      # Display current A and b values 
      # in dictionary format
   
      for ii to nargs
      do
        # Display i-th row of dictionary
        printf(`     x%d  = %8.2f `,lhslist[ii], evalf( b[ basrow[ii] , 1] ) );
        
        for jj to nn
         do
            if not member( jj , {op(lhslist)})
            then
                if A[ basrow[ii] ,jj] > myeps 
                then
                  printf(` - %6.2f x%d `,evalf(A[basrow[ii] ,jj]), jj);
                else
                  if A[basrow[ii] ,jj] < -myeps 
                  then
                    printf(` + %6.2f x%d `,evalf(-A[basrow[ii] ,jj]), jj);
                  else
                    printf(`             `);
                  fi       
                fi       
            fi
         od; # for jj
         printf(`\n`);
      od;  # for ii
      printf(`\n`);

    fi; # if good basis
  
   fi; # nargs correct
     
end:
###############################################################################

tableau := proc()
local ii, jj, goodbasisflag;
global A, b, c, origc, z, mm, nn, lhslist, basrow;
 
   ####
   #### First figure out what basis the user wants
   ####
   #### (Note: this is an ugly bit of code.)
   if nargs <> mm
   then
   printf(`ERROR (dict): user must provide an ordered list of %d basic variables.\n`, mm);
   else
      lhslist := [];
      for ii to nargs 
      do
       for jj to nn
       do
          if args[ii] = x||jj then lhslist := [op(lhslist),jj]; jj:=nn; fi;
       od
      od;
      #### Now check if it really is a current basis:
      basrow := [];
      goodbasisflag := true;
      for ii to nargs
       do
          for jj to mm
           do
              if abs( A[jj,lhslist[ii]] - 1.) < myeps
               then
                  if nops( basrow ) < ii 
                  then  # good! we've found the entry we want.
                     basrow := [ op(basrow), jj];
                  else # Uh oh! We've already got one non-zero entry in that column!
                      goodbasisflag := false; 
                      jj:=mm;
                      ii:=nargs;
                  fi;
               else # If there is any other nonzero entry in this column, we're dead.
                  if abs( A[jj, lhslist[ii]] ) > myeps
                  then goodbasisflag := false; 
                         jj:=mm;  ii:=nargs;
                  fi;
               fi;  # if entry is 1
           od;  # for jj
       od; # for ii
   
     if goodbasisflag
     then
      ####
      ####  Display the current tableau
      ####
    
      printf(`\n               |`);
      # List variables
      for jj to nn  do printf(`   x%-2d `,jj); od;
      printf(`  |\n   Basis   c_B |`);
      for jj to nn  do printf(` %4d  `,origc[1,jj]); od;
      printf(`  |   RHS \n  _____________|`);
      for jj to nn  do printf(`_______`); od;
      printf(`__|__________\n`);
      for ii to mm 
       do 
         printf(`    x%-2d  %5.2f |  `,lhslist[ii], origc[1, lhslist[ii] ]);
         for jj to nn do printf(` % 4.2f `,A[ii,jj]); od;
         printf(`| %6.2f \n`,b[ii,1]);
       od;
      printf(`  _____________|`);
      for jj to nn  do printf(`_______`); od;
      printf(`__|__________\n       zj      |  `);
      for jj to nn do printf(` % 4.2f `,add( origc[1,lhslist[hh]]*A[hh,jj],hh=1..mm) ); od;
      printf(`|%8.2f \n`, add( origc[1,lhslist[hh]]*b[hh,1],hh=1..mm) );
      printf(`    cj - zj    |  `);
      for jj to nn do printf(` % 4.2f `,origc[1,jj]-add( origc[1,lhslist[hh]]*A[hh,jj],hh=1..mm) ); od;
      printf(`|\n\n`);
       
       
      

    fi; # if good basis
  
   fi; # nargs correct
     
end:
###############################################################################