#  August 28, 2017. MAPLE code to pivot on a dictionary for MA3231.
#  Use this -- or the equivalent MATLAB code -- to check your calculations.
#
#  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.
#
#     dratios: dratios(i) determines all upper bounds for leaving variable xi.
#
myeps := .1^8; # When working numerically, treat this as zero


 ### This is a bad idea, but let's assume it's initially in standard form.
 if nops(lhslist)=0 
 then 
   lhslist := [seq(h,h=linalg[coldim](A)+1-linalg[rowdim](A)..linalg[coldim](A))]; 
   basrow := [seq(h,h=1..linalg[rowdim](A))];
 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, basrow, lhslist;

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

# Need to search for pivot row:
rw:=0;
for ii to linalg[rowdim](A) 
  do 
     if abs( A[ii,ext] - 1.) < myeps 
     then 
       rw := ii;  ii := linalg[rowdim](A);
     fi;
  od;
for ii to linalg[rowdim](A) 
  do 
     if ii <> rw and abs( A[ii,ext] ) > myeps 
     then 
       rw := 0;  ii := linalg[rowdim](A);
     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 linalg[coldim](A) 
   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 linalg[rowdim](A)
   do
    if ii <> rw
    then
     pivcolval := A[ii,cl]:
     for jj from 1 to linalg[coldim](A) 
      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 linalg[coldim](A) 
   do
          c[1,jj]  := c[1,jj] - pivcolval*A[rw,jj]:
   od:
          z := z + pivcolval*b[rw,1]:

  for h to linalg[rowdim](A)  # Swap variable out of basis
   do
    if lhslist[h] = ext then lhslist[h] := ent;   fi
   od;

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

   dict( 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, lhslist, basrow;

printf("\n");
for ii to linalg[rowdim](A)
 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, lhslist, basrow;
 
   ####
   #### First figure out what basis the user wants
   ####
   #### (Note: this is an ugly bit of code.)
   if nargs <> linalg[rowdim](A)
   then
   printf("ERROR (dict): user must provide an ordered list of %d basic variables.\n", linalg[rowdim](A));
   else
      lhslist := [];
      for ii to nargs 
      do
       for jj to linalg[coldim](A)
       do
          if args[ii] = x||jj then lhslist := [op(lhslist),jj]; jj:=linalg[coldim](A); fi;
       od
      od;
      #### Now check if it really is a current basis:
      basrow := [];
      goodbasisflag := true;
      for ii to nargs
       do
          for jj to linalg[rowdim](A)
           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:=linalg[rowdim](A); 
                      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:=linalg[rowdim](A); 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 linalg[coldim](A)
       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 linalg[coldim](A)
       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 linalg[coldim](A)
         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;
###############################################################################


###############################################################################


dratios := proc(ii)  # Dual Simplex: Give ratios for exiting variable x.ii
local jj;
global A, c, lhslist, basrow;

printf("\n");
for jj to linalg[coldim](A)
 do
  if not member( jj , {op(lhslist)})
  then
   if A[basrow[ii],jj] >= 0 
      then 
           printf("    non-basic x%d:  No upper bound \n",jj);
      else 
           printf("    non-basic x%d:  Upper bound = %8.4f \n",
                                     jj, evalf(c[1,jj]/A[ basrow[ii],jj]) );
   fi:
  fi:
 od:
printf("\n");
end;
###############################################################################