##########################################################################
##Transportation.txt: Save this file as  Transportation.txt               #        
## To use it, stay in the                                                 #
##same directory, get into Maple (by typing: maple <Enter> )              #         
##and then type:  read Transportation.txt<Enter>                          #          
##Then follow the instructions given there                                #
##                                                                        #
##Written by Doron Zeilberger, Rutgers University ,                       #
#DoronZeil at gmail dot com                                               #
##########################################################################
 
#First Created: March 2019
#Extended Nov. 2023
 
print(`Created:  March 2019`):
print(`Extended:  Nov. 2023`):
print(` This is Transportation.txt `):
print(`It is a Maple package that accompany Doron Zeilberger's course Math 354: Introduction to Linear Optimization`):
print(``):
print(`It implementes the algorithm decribe in section 5.1 of the textbook of this class`):
print(``):
print(` "Elementary Linear Programming  with Applications" by Bernard Kolman and Robert E. Beck,2nd ed., Academic Press.`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AllPaths, ContiH, ContiV, Diff12, FindE, IsFe, Katan, OneStep, MaxC, MCR1, Pathsk, PickLine , TAe`):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: MCR, RandProb, TA,TAd, Vogel `):
 print(` `):


elif nops([args])=1 and op(1,[args])=AllPaths then
print(`AllPaths(v,S): inputs a lattice point v, and a set of lattice points S,  outputs all alternating vertical-horizontal`):
print(`paths that are self-avoding and that start and end at v. Try: `):
print(`AllPaths([3,4],{[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]});`):

elif nops([args])=1 and op(1,[args])=ContiH then
print(`ContiH(S,P): all the horizontal continuations of the path P starting at v via S.Try:`):
print(`ContiH({[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]},[[3,4],[1,4]]);`):

elif nops([args])=1 and op(1,[args])=ContiV then
print(`ContiV(S,P): all the vertical continuations of the path P starting at v via S.Try:`):
print(`ContiV({[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]},[[3,4]]);`):

elif nops([args])=1 and op(1,[args])=Diff12 then
print(`Diff12(L): the difference between the smallest and the next smallest in the list L. Try:`):
print(`Diff12([4,5,7]);`):

elif nops([args])=1 and op(1,[args])=FindE then
print(`FindE(C,s,d,M): inputs a cost matrix C, supply and demand vectors s and d and a current transportation tableau M`):
print(`finds the entering variable. it also returns the set of basic cells. Try:`):
print(`M:=MCR([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120])[1];`):
print(`FindE([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120],M);`):

elif nops([args])=1 and op(1,[args])=IsFe then
print(`IsFe(M,s,d): is the transportation tableu M  feasible where s and d are the supply and demand vectors?`):
print(`Try:`):
print(`IsFe([[100,0,0,20],[0,40,0,100],[0,20,80,0]],[120,140,100],[100,60,80,120] ); `);

elif nops([args])=1 and op(1,[args])=Katan then
print(`Katan(L,S): the smallest entry in the list L in the places S, also returns S minus that entry. Try:`):
print(`Katan([3,2,2,1,5,5],{1,2,3,5,6});`):

elif nops([args])=1 and op(1,[args])=MCR then
print(`MCR(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the`):
print(`initial basic feasible solution to the transportation problem using the Minimal Cost Rule. For example, to do Problem 3(a)`):

print(`It also returns the value`):

print(`of section 5.1 in the Kolman-Beck book "Introduction to Linear Programming" type:`):
print(`MCR([[6,6,3,9,2],[8,7,5,7,5],[3,4,8,2,7],[6,0,0,2,9]],[90,70,110,150],[100,40,100,60,120]);`):

print(`To do tableau 5.1 on p. 299 type:`):
print(`MCR([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120]);`):


elif nops([args])=1 and op(1,[args])=MCR1 then
print(`MCR1(C,s,d,M,i): inputs a trnasportation problem C,s,d, a partially filled tableau (up the first i-1 rows) and`):
print(`a row number i, performs the Minimal cost rule to row i.Try:`):
print(`MCR1([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120],[[0,0,0,0],[0,0,0,0],[0,0,0,0]],1);`):

elif nops([args])=1 and op(1,[args])=MaxC then
print(`MaxC(M,s,d,i,j): inputs a currently filled transportaion tableu M, with supply vector s, outputs`):
print(`and an entry [i,j] where it is currently 0, find the maximum allowed capacity`):
print(`the available load. Try:`):
print(`MaxC([[0,0,0,0],[0,0,0,0],[0,0,0,0]],[120,140,100],[100,60,80,120],1,1);`):
print(`MaxC([[0$5]$4],[90,70,110,150],[100,40,100,60,120],1,1);`):

elif nops([args])=1 and op(1,[args])=OneStep then
print(`OneStep(C,s,d,M): performs one step in the Transportation Problem algorithm. Try:`):
print(`M:=MCR([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120]);`):
print(`OneStep([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120],M);`):

elif nops([args])=1 and op(1,[args])=Pathsk then
print(`Pathsk(v,S,k): inputs a lattice point v, and a set of lattice points S,  outputs all alternating vertical-horizontal`):
print(`paths that are self-avoiding of`):
print(`length 2k-1 via S starting at v. Try:`):
print(`Pathsk([3,4],{[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]},2);`):

elif nops([args])=1 and op(1,[args])=PickEntry then
print(`PickEntry(C,RS,CS): inputs a cost matrix C, RS:a sublist of the set of rows, CS: a sublist of the set of columns`):
print(`outputs the entry,in Vogel's method to make non-zero. Try:`):
print(`PickEntry([[8,6,3,9],[2,6,1,4],[7,8,6,3]],{1,2,3},{1,2,3,4});`):

elif nops([args])=1 and op(1,[args])=PickLine then
print(`PickLine(C,RS,CS): inputs a cost matrix C, RS:a sublist of the set of rows, CS: a sublist of the set of columns`):
print(`outputs the row or column with the largest difference between smallest and second smallest, the output is a pair`):
print(`it is either [MemberOfRS,0] or  [MemberOfCS,1] . Try:`):
print(`PickLine([[8,6,3,9],[2,6,1,4],[7,8,6,3]],{1,2,3},{1,2,3,4});`):

elif nops([args])=1 and op(1,[args])=RandProb then
print(`RandProb(m,n,K,L): inputs a random m by n cost function, with entries from 3 to K, and random  supply vector of length m`):
print(`and a random demand function of length n whith entires between L and 2*L. Try:`):
print(`RandProb(3,3,10,100);`):

elif nops([args])=1 and op(1,[args])=TA then
print(`TA(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the`):
print(`optimal transportation tableau, followed by the cost.`):
print(`It uses the Minimal Cost Rule for the initial tableau.`):
print(`To do example 1 in setion 5.1 in the Kolman-Beck book, type:`):
print(``):
print(`TA([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120]);`):
print(``):
print(`To do example 2 in setion 5.1 in the Kolman-Beck book, type:`):
print(``):
print(`TA([[9,3,6,7,3],[7,5,2,10,6],[5,4,9,8,10]],[100,160,140],[90,60,80,100,70]);`):
print(``):
print(`To do example 3 in setion 5.1 in the Kolman-Beck book, (alias exercise 5) type:`):
print(``):
print(`TA([[9,3,6,7,3],[7,5,2,10,6],[5,4,9,8,10]],[100,160,100],[50,60,80,100,70]);`):
print(``):
print(`To do example 4 in setion 5.1 in the Kolman-Beck book, type:`):
print(``):
print(`TA([[8,6,3,9],[2,6,1,4],[7,8,6,3]],[120,140,100],[100,60,80,120]);`):
print(``):
print(`To do exercise 2 at the end of the section do`):
print(``):
print(`TA([[5,2,3,6],[2,7,7,4],[1,3,6,9]],[100,80,140],[40,60,80,120]);`):
print(``):
print(`To do exercise 7 at the end of the section do`):
print(``):
print(`TA([[6,6,3,9,2],[8,7,5,7,5],[3,4,8,2,7],[6,0,0,2,9]],[90,70,110,150],[100,40,100,60,120]);`):
print(``):
print(`To do exercise 4 at the end of the section do`):
print(``):
print(`TA([[4,6,7,5],[4,4,7,8],[5,3,6,5],[6,5,3,4]],[100,60,50,70],[40,70,120,50]);`):
print(``):
print(`To do exercise 8 at the end of the section do`):
print(``):
print(`TA([[2,5,6,3],[9,6,2,1],[7,7,2,4]],[100,90,130],[70,50,30,120]);`):
print(``):
print(`To do exercise 9 at the end of the section do`):
print(``):
print(`TA([[3,2,5,4],[6,5,7,8],[2,1,4,3],[4,3,5,2]],[80,60,50,100],[70,70,70,80]);`):
print(``):
print(`To do exercise 10 at the end of the section do`):
print(``):
print(`TA([[4,3,2,5,6],[8,3,4,5,7],[6,8,6,7,5],[4,3,5,2,4]],[70,80,60,120],[60,50,50,70,100]);`):
print(``):
print(``):
print(`To do exercise 11 at the end of the section do`):
print(``):
print(`TA([[4,2,9,7],[7,8,5,6],[3,3,4,1],[7,5,2,6]],[1,1,1,1],[1,1,1,1]);`):
print(``):
print(``):
print(`To do exercise 12 at the end of the section do`):
print(``):
print(`TA([[5,6,7,4],[2,9,7,5],[8,5,8,7]],[75,50,60],[45,50,25,50]);`):
print(``):
print(``):
#print(`To do exercise 13 at the end of the section do`):
#print(``):
#print(`TA([[6,4,3,5],[7,4,8,6],[8,3,2,5]],[100,60,50],[60,80,70,40]);`):
#print(``):

elif nops([args])=1 and op(1,[args])=TAd then
print(`TAd(C,s,d): like TA(C,s,d) but for the degenerate case. `):
print(`To do exercise 11 at the end of the section do`):
print(``):
print(`TAd[[4,2,9,7],[7,8,5,6],[3,3,4,1],[7,5,2,6]],[1,1,1,1],[1,1,1,1]);`):


elif nops([args])=1 and op(1,[args])=TAe then
print(`TAe(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the, such that supply equals demand`):
print(`optimal transportation tableau, followed by the cost.`):
print(`It uses the Minimal Cost Rule for the initial tableau.`):
print(`To do example 1 in setion 5.1 in the Kolman-Beck book, type:`):
print(``):

elif nops([args])=1 and op(1,[args])=Vogel then
print(`Vogel(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the`):
print(`initial basic feasible solution to the transportation problem using Vogel's method. For example, to do Problem 3(b)`):
print(`It also returns the value`):
print(`of section 5.1 in the Kolman-Beck book "Introduction to Linear Programming" type:`):
print(`Vogel([[6,6,3,9,2],[8,7,5,7,5],[3,4,8,2,7],[6,0,0,2,9]],[90,70,110,150],[100,40,100,60,120]);`):
print(`To do Example 4 in p. 320 type:`):
print(`Vogel([[8,6,3,9],[2,6,1,4],[7,8,6,3]],[120,140,100],[100,60,80,120]);`):


else
print(`There is no ezra for`,args):
fi:
 
end:



#MaxC(M,s,d,i,j): inputs a currently filled transportaion tableu M, with supply vector s, outputs
#and an entry [i,j] where it is currently 0, find the maximum allowed capacity
#the available load. Try:
#MaxC([[0$5]$4],[90,70,110,150],[100,40,100,60,120],1,1);
MaxC:=proc(M,s,d,i,j) local i1,j1:
if not (1<=i and i<=nops(M) and 1<=j and j<=nops(M[1])) then
 RETURN(FAIL):
fi:

if M[i][j]<>0 then
 RETURN(FAIL):
fi:

min(s[i]-add(M[i][j1],j1=1..nops(M[i])), d[j]-add(M[i1][j],i1=1..nops(M))):



end:


#IsFe(M,s,d): is the transportation tableu M  feasible where s and d are the supply and demand vectors?
IsFe:=proc(M,s,d) local i,j:

if not (nops(M)=nops(s) and nops(M[1])=nops(d)) then
print(`Bad input`):
 RETURN(FAIL):
fi:

if min(seq(seq(M[i][j],j=1..nops(M[i])),i=1..nops(M)))<0 then
 RETURN(false):
fi:

for i from 1 to nops(s) do
if add(M[i][j],j=1..nops(M[i])) >s[i] then
 RETURN(false):
fi:
od:

for j from 1 to nops(d) do
if add(M[i][j],i=1..nops(M))<d[j] then
 RETURN(false):
fi:
od:
true:
end:


#Katan(L,S): the smallest entry in the list L in the places S, also returns S minus that entry. Try:
#Katan([3,2,2,1,5,5],{1,2,3,5,6});
Katan:=proc(L,S) local aluf,si,i:

if L=[] then
 RETURN(FAIL):
fi:

if nops(S)>nops(L) then
 RETURN(FAIL):
fi:

aluf:=S[1]:
si:=L[S[1]]:

for i from 2 to nops(S) do
 if L[S[i]]<si then
   aluf:=S[i]:
   si:=L[S[i]]:
 fi:
od:
aluf,S minus {aluf}:
end:



#MCR1(C,s,d,M,i): inputs a trnasportation problem C,s,d, a partially filled tableau (up the first i-1 rows) and
#a row number i, performs the Minimal cost rule to row i.Try:
#MCR1([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120],[[0,0,0,0],[0,0,0,0],[0,0,0,0]],1);`):
MCR1:=proc(C,s,d,M,i) local S,i1,M1,gu,j1,ne:

 S:={seq(i1,i1=1..nops(C[i]))}:

M1:=M:

 while S<>{} do
 gu:=Katan(C[i],S):
  j1:=gu[1]:
   S:=gu[2]:
   ne:=MaxC(M1,s,d,i,j1):
  M1:=[op(1..i-1,M1),[op(1..j1-1,M1[i]),ne,op(j1+1..nops(M1[i]),M1[i] )],op(i+1..nops(M1),M1)]:
 od:
M1:
end:


#MCR(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the
#initial basic feasible solution to the transportation problem. For example, to do Probelm 3
#of section 5.1 in the Kolman-Beck book "Introduction to Linear Programming" type:
#MCR([[6,6,3,9,2],[8,7,5,7,5],[3,4,8,2,7],[6,0,0,0,2,9]],[90,70,110,150],[100,40,100,60,120]]);
MCR:=proc(C,s,d) local i,j,m,n,M1,opval:

if not (C<>[] and type(C,list))  then
 print(C,  `must have been a non-empty list`):
 RETURN(FAIL):
fi:
m:=nops(C):

if not {seq(type(C[i],list),i=1..m)}={true} then
 print(C, `must have been a list of lists`):
 RETURN(FAIL):
fi:
n:=nops(C[1]):

 if {seq(nops(C[i]),i=2..m)}<> {n} then
   print(C, `must be a list of lists all of the same length`):
 fi:

if not (type(s,list) and type(d,list) and nops(s)=m and nops(d)=n) then
 print(`Bad input`):
  RETURN(FAIL):
fi:

M1:=[[0$n]$m]:

for i from 1 to nops(C) do
 M1:=MCR1(C,s,d,M1,i):
od:
opval:=add(add(C[i][j]*M1[i][j],j=1..nops(C[i])),i=1..nops(C)):
M1,opval:
end:

#Diff12(L): the difference between the smallest and the next smallest in the list L. Try:
#Diff12([4,5,7]);
Diff12:=proc(L) local k1,k2,S,L1,i,j1,gu:
if nops(L)<2 then
 RETURN(FAIL):
fi:


gu:=Katan(L,{seq(i,i=1..nops(L))}):
j1:=gu[1]:
S:=gu[2]:
k1:=L[j1]:
L1:=[seq(L[i],i in S)]:
k2:=min(op(L1)):
k2-k1:
end:

#PickLine(C,RS,CS): inputs a cost matrix C, RS:a sublist of the set of rows, CS: a sublist of the set of columns
#outputs the row or column with the largest difference between smallest and second smallest, the output is a pair
#it is either [MemberOfRS,0] or  [MemberOfCS,1] . Try:
#PickLine([[8,6,3,9],[2,6,1,4],[7,8,6,3]],{1,2,3},{1,2,3,4});
PickLine:=proc(C,RS,CS) local i,j,R1,C1,si,halev,i1,j1,aluf:

R1:=[seq(C[RS[1]][CS[j]],j=1..nops(CS))]:
si:=Diff12(R1):
aluf:=[0,RS[1]]:

for i from 2 to nops(RS) do
 R1:=[seq(C[RS[i]][CS[j1]],j1=1..nops(CS))]:
 halev:=Diff12(R1):
  if halev>si then
   aluf:=[0,RS[i]]:
  fi:
od:

for j from 1 to nops(CS) do
 C1:=[seq(C[RS[i1]][CS[j]],i1=1..nops(RS))]:
 halev:=Diff12(C1):
  if halev>si then
   aluf:=[1,CS[j]]:
  fi:
od:

aluf:

end:


#PickEntry(C,RS,CS): inputs a cost matrix C, RS:a sublist of the set of rows, CS: a sublist of the set of columns
#outputs the entry,in Vogel's method to make non-zero. Try:
#PickEntry([[8,6,3,9],[2,6,1,4],[7,8,6,3]],{1,2,3},{1,2,3,4});
PickEntry:=proc(C,RS,CS) local gu,i,j,R1,C1,i1,mu:

gu:=PickLine(C,RS,CS):

if gu[1]=0 then
 i:=gu[2]:
  R1:=C[i]:
 mu:=Katan(R1,CS)[1]:
RETURN([i,mu]):
elif gu[1]=1 then
  j:=gu[2]:
  C1:=[seq(C[i1][j],i1=1..nops(C))]:
   mu:=Katan(C1,RS)[1]:
  RETURN([mu,j]):
else
 RETURN(FAIL):
fi:
end:

  

#Vogel(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the
#initial basic feasible solution to the transportation problem. For example, to do Probelm 3
#of section 5.1 in the Kolman-Beck book "Introduction to Linear Programming" type:
#Vogel([[6,6,3,9,2],[8,7,5,7,5],[3,4,8,2,7],[6,0,0,0,2,9]],[90,70,110,150],[100,40,100,60,120]]);
Vogel:=proc(C,s,d) local RS,CS,NZset,T,i,j,Ts,Td,lu,ku,i1,j1,T1,M1,opval:
RS:={seq(i,i=1..nops(C))}:
CS:={seq(j,j=1..nops(C[1]))}:

NZset:={}:

for i from 1 to nops(s) do
 Ts[i]:=s[i]:
od:

for j from 1 to nops(d) do
 Td[j]:=d[j]:
od:

while nops(RS)>=2 and nops(CS)>=2 do
 lu:=PickEntry(C,RS,CS):
  NZset:=NZset union {lu}:
  i:=lu[1]:
  j:=lu[2]:
  ku:=min(Ts[i],Td[j]):
  T[lu]:=ku:
  Ts[i]:=Ts[i]-ku:
  Td[j]:=Td[j]-ku:
  if Ts[i]=0 then
   RS:=RS minus {i}:
  elif Td[j]=0 then
   CS:=CS minus {j}:
  fi:
od:



if nops(CS)=1 then
 NZset:=NZset union {seq([RS[i1],CS[1]],i1=1..nops(RS))}:
 for i1 from 1 to nops(RS) do
 T[[RS[i1],CS[1]]]:=Ts[RS[i1]]:
 od:
else
 NZset:=NZset union {seq([RS[1],CS[j1]],j1=1..nops(CS))}:
  for j1 from 1 to nops(CS) do
 T[[RS[1],CS[j1]]]:=Td[CS[j1]]:
 od:
fi:

for i from 1 to nops(C) do
 for j from 1 to nops(C[i]) do
 if member([i,j],NZset) then
  T1[[i,j]]:=T[[i,j]]:
 else
  T1[[i,j]]:=0:
 fi:
 od:
od:

M1:=[seq([seq(T1[[i1,j1]],j1=1..nops(C[i1]))],i1=1..nops(C))]:
opval:=add(add(C[i][j]*M1[i][j],j=1..nops(C[i])),i=1..nops(C)):
M1,opval:
end:

    
  

#FindE(C,s,d,M): inputs a cost matrix C, supply and demand vectors s and d and a current transportation tableau M
#finds the cell of the entering variable, it also returns the set of basic cells.  Try:
#print(`M:=MCR([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120])[1];`):
#FindE([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120],M);`):
FindE:=proc(C,s,d,M) local m,n,i,j,B,NB,v,w,var,eq,b,gad:
m:=nops(C):
n:=nops(C[1]):

if not (nops(s)=m and nops(d)=n) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

B:={}:
NB:={}:

for i from 1 to m do
 for j from 1 to n do
   if M[i][j]>0 then
    B:=B union {[i,j]}:
   else
     NB:=NB union {[i,j]}:
   fi:
 od:
od:

if nops(B)<m+n-1 then
#  print(`degenerate case not yet implemented`):
  RETURN(FAIL):
elif nops(B)>m+n-1 then
  print(`bad input`):
  RETURN(FAIL):
fi:

var:={seq(v[i],i=1..m), seq(w[j],j=1..n)}:
eq:={}:

for b in B do
 eq:=eq union {v[b[1]]+w[b[2]]-C[b[1]][b[2]]}:
od:

eq:=eq union {v[1]=0}:

var:=solve(eq,var):

v:=subs(var,[seq(v[i],i=1..m)]):
w:=subs(var,[seq(w[j],j=1..n)]):

gad:=max(seq(v[b[1]]+w[b[2]]-C[b[1]][b[2]], b in NB)):
gad:

if gad<=0 then
 RETURN(M):
fi:

for b in NB do
 if v[b[1]]+w[b[2]]-C[b[1]][b[2]]=gad then
  RETURN(b,B):
 fi:
od:
FAIL:
end:


#ContiV(S,P): all the vertical continuations of the path P starting at v via S.Try:
#ContiV([3,4],{[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]},[[3,4]]);
ContiV:=proc(S,P) local S1,s1,sof,gu:
S1:=S minus {op(P)}:
gu:={}:
sof:=P[nops(P)]:

for s1 in S1 do
 if s1[2]=sof[2] then
  gu:=gu union {[op(P),s1]}:
 fi:
od:
gu:
end:

#ContiH(S,P): all the horizontal continuations of the path P starting at v via S.Try:
#ContiH({[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]},[[3,4],[1,4]]);
ContiH:=proc(S,P) local S1,s1,sof,gu:
S1:=S minus {op(P)}:
gu:={}:
sof:=P[nops(P)]:

for s1 in S1 do
 if s1[1]=sof[1] then
  gu:=gu union {[op(P),s1]}:
 fi:
od:
gu:
end:


#Pathsk(v,S,k): inputs a lattice point v, and a set of lattice points S,  outputs all alternating vertical-horizontal
#paths that are self-avoding of
#length 2k-1 via S starting at v. Try
#Pathsk([3,4],{[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]},2);
Pathsk:=proc(v,S,k) local gu,gu1:
option remember:
if k=1 then
 RETURN(ContiV(S,[v])):
fi:
gu:=Pathsk(v,S,k-1):

gu:={seq(op(ContiH(S,gu1)),gu1 in gu)}:

gu:={seq(op(ContiV(S,gu1)),gu1 in gu)}:
gu:

end:

#AllPaths(v,S): inputs a lattice point v, and a set of lattice points S,  outputs all alternating vertical-horizontal
#paths that are self-avoding and that start and end at v. Try
#AllPaths([3,4],{[1,1],[1,4],[2,2],[2,4],[3,2],[3,3]});

AllPaths:=proc(v,S) local k,gu,lu,lu1:
gu:={}:

for k from 2 while lu<>{} do
 lu:=Pathsk(v,S,k):
  for lu1 in lu do
    if lu1[nops(lu1)][1]=v[1] then
      gu:=gu union {[op(lu1),v]}:
    fi:
 od:
od:
gu:
end:

#OneStep(C,s,d,M): performs one step in the Transportation Problem algorithm. Try:
#M:=MCR([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120]);
#OneStep([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120],M);

OneStep:=proc(C,s,d,M) local gu,lu,k,T,i,j:
gu:=FindE(C,s,d,M):
if gu=FAIL then
 RETURN(FAIL):
fi:


if gu=M then
 RETURN(M):
fi:

gu:=AllPaths(gu):

if gu={} then
 RETURN(FAIL):
fi:

gu:=gu[1]:

k:=(nops(gu)-1)/2:

lu:=min(seq(M[gu[2*i][1]][gu[2*i][2]],i=1..k)):

for i from 1 to nops(M) do
for j from 1 to nops(M[i]) do
 T[i,j]:=M[i][j]:
od:
od:

for i from 1 to 2*k do
 if i mod 2=1 then
   T[op(gu[i])]:=T[op(gu[i])]+lu:
 else
   T[op(gu[i])]:=T[op(gu[i])]-lu:
 fi:
od:

[seq([seq(T[i,j],j=1..nops(M[i]))],i=1..nops(M))]:
end:



#TAe(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the
#optimal transportation tableau, followed by the cost.
#It uses the Minimal Cost Rule for the initial tableau.
#To do example 1 in setion 5.1 in the Kolman-Beck book, type:
#TAe([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120]]);
#To do example 2 in setion 5.1 in the Kolman-Beck book, type:
#TAe([[9,3,6,7,3],[7,5,2,10,6],[5,4,9,8,10]],[100,160,140],[90,60,80,100,70]]);
#To do example 3 in setion 5.1 in the Kolman-Beck book, (alias exercise 5) type:
#TAe([[9,3,6,7,3],[7,5,2,10,6],[5,4,9,8,10]],[100,160,100],[50,60,80,100,70]]);
#To do example 4 in setion 5.1 in the Kolman-Beck book, type:
#TAe([[8,6,3,9],[2,6,1,4],[7,8,6,3]],[120,140,100],[100,60,80,120]]);
#To do exercise 2 at the end of the section do
#TAe([[5,2,3,6],[2,7,7,4],[1,3,6,9]],[100,80,140],[40,60,80,120]]);
#To do exercise 3 at the end of the section do
#TAe([[6,6,3,9,2],[8,7,5,7,5],[3,4,8,2,7],[6,0,0,0,2,9]],[90,70,110,150],[100,40,100,60,120]]);
#To do exercise 4 at the end of the section do
#TAe([[4,6,7,5],[4,4,7,8],[5,3,6,5],[6,5,3,4]],[100,60,50,70],[40,70,120,50]]);

TAe:=proc(C,s,d) local i,j,M1,M2,M3:

if convert(s,`+`)<convert(d,`+`) then
 print(`The total supply must be >= the total demand`):
 RETURN(FAIL):
fi:

if convert(s,`+`)>convert(d,`+`) then
 print(`total supply must equal total demand`):
 RETURN(FAIL):
fi:

M1:=MCR(C,s,d)[1]:
if M1=FAIL  then
 RETURN(TAd(C,s,d)):
fi:

M2:=OneStep(C,s,d,M1):
if M2=FAIL then
 RETURN(TAd(C,s,d)):
fi:
while M2<>M1 do
 M3:=OneStep(C,s,d,M2):
 M1:=M2:
 M2:=M3:
od:

M2,add(add(C[i][j]*M2[i][j],j=1..nops(M2[i])),i=1..nops(M2)):
end:



#TA(C,s,d): given a cost matrix C, a supply vector s and a demand vector d, outputs the
#optimal transportation tableau, followed by the cost.
#It uses the Minimal Cost Rule for the initial tableau.
#To do example 1 in setion 5.1 in the Kolman-Beck book, type:
#TA([[5,7,9,6],[6,7,10,5],[7,6,8,1]],[120,140,100],[100,60,80,120]]);
#To do example 2 in setion 5.1 in the Kolman-Beck book, type:
#TA([[9,3,6,7,3],[7,5,2,10,6],[5,4,9,8,10]],[100,160,140],[90,60,80,100,70]]);
#To do example 3 in setion 5.1 in the Kolman-Beck book, (alias exercise 5) type:
#TA([[9,3,6,7,3],[7,5,2,10,6],[5,4,9,8,10]],[100,160,100],[50,60,80,100,70]]);
#To do example 4 in setion 5.1 in the Kolman-Beck book, type:
#TA([[8,6,3,9],[2,6,1,4],[7,8,6,3]],[120,140,100],[100,60,80,120]]);
#To do exercise 2 at the end of the section do
#TA([[5,2,3,6],[2,7,7,4],[1,3,6,9]],[100,80,140],[40,60,80,120]]);
#To do exercise 3 at the end of the section do
#TA([[6,6,3,9,2],[8,7,5,7,5],[3,4,8,2,7],[6,0,0,0,2,9]],[90,70,110,150],[100,40,100,60,120]]);
#To do exercise 4 at the end of the section do
#TA([[4,6,7,5],[4,4,7,8],[5,3,6,5],[6,5,3,4]],[100,60,50,70],[40,70,120,50]]);

TA:=proc(C,s,d) local i,C1,M,lu,d1,gu:

if convert(s,`+`)<convert(d,`+`) then
 print(`The total supply must be >= the total demand`):
 RETURN(FAIL):
fi:

if convert(s,`+`)=convert(d,`+`) then
 RETURN(TAe(C,s,d)):
fi:

lu:=convert(s,`+`)-convert(d,`+`):

C1:=[seq([op(C[i]),0],i=1..nops(C))]:
d1:=[op(d),lu]:

gu:=TAe(C1,s,d1):
lu:=gu[2]:
M:=gu[1]:
[seq([op(1..nops(M[i])-1,M[i])],i=1..nops(M))]:
M,lu:
end:



#TAeV(C,s,d): Like TAe(C,d,s) but using Vogel's method for the intial tableau
TAeV:=proc(C,s,d) local i,j,M1,M2,M3:

if convert(s,`+`)<convert(d,`+`) then
 print(`The total supply must be >= the total demand`):
 RETURN(FAIL):
fi:

if convert(s,`+`)>convert(d,`+`) then
 print(`total supply must equal total demand`):
 RETURN(FAIL):
fi:

M1:=Vogel(C,s,d)[1]:

if M1=FAIL  then
 RETURN(FAIL):
fi:


M2:=OneStep(C,s,d,M1):
if M2=FAIL then
 RETURN(FAIL):
fi:
while M2<>M1 do
 M3:=OneStep(C,s,d,M2):
 M1:=M2:
 M2:=M3:
od:

M2,add(add(C[i][j]*M2[i][j],j=1..nops(M2[i])),i=1..nops(M2)):
end:

TAv:=proc(C,s,d) local i,C1,M,lu,d1,gu:

if convert(s,`+`)<convert(d,`+`) then
 print(`The total supply must be >= the total demand`):
 RETURN(FAIL):
fi:

if convert(s,`+`)=convert(d,`+`) then
 RETURN(TAeV(C,s,d)):
fi:

lu:=convert(s,`+`)-convert(d,`+`):



C1:=[seq([op(C[i]),0],i=1..nops(C))]:
d1:=[op(d),lu]:

gu:=TAeV(C1,s,d1):


lu:=gu[2]:
M:=gu[1]:
[seq([op(1..nops(M[i])-1,M[i])],i=1..nops(M))]:
M,lu:
end:


#RandKatan(n,K): a random vector of length n that adds up to 0
RandKatan:=proc(n,K) local ra,lu,i:
ra:=rand(K..2*K):
lu:=[seq(1/ra(),i=1..n-1)]:
lu:=[op(lu),-convert(lu,`+`)]:
end:


TAd:=proc(C,s,d) local gu,s1,d1,co,i,j:
gu:=TAv(C,s,d):
if gu<>FAIL then
 RETURN(gu):
fi:

s1:=s+RandKatan(nops(s),10^7);
d1:=d+RandKatan(nops(d),10^7);


gu:=TAv(C,s1,d1):


if gu=FAIL then
 RETURN(FAIL):
fi:

co:=gu[2]:
gu:=gu[1]:

[seq([seq(round(gu[i][j]),j=1..nops(gu[i]))],i=1..nops(gu))],round(co):

end:

#RandProb(m,n,K,L): inputs a random m by n cost function, with entries from 3 to K, and random  supply vector of length m
#and a random demand function of length n whith entires between L and 2*L. Try
#RandProb(3,3,10,100);
RandProb:=proc(m,n,K,L) local ra1,ra2,i,j,C,s,d,Ns,Nd:
ra1:=rand(3..K):
ra2:=rand(L..2*L):

C:=[seq([seq(ra1(),j=1..n)],i=1..m)]:
s:=[seq(ra2(),i=1..m-1)]:
d:=[seq(ra2(),j=1..n-1)]:
Ns:=convert(s,`+`):
Nd:=convert(d,`+`):
if Ns>Nd then
 s:=[op(s),0]:
 d:=[op(d),Ns-Nd]:
elif Ns<Nd then
 s:=[op(s),Nd-Ns]:
 d:=[op(d),0]:
fi:

s:=s+[11$m]:
d:=d+[11$n]:

C,s,d:
end: