######################################################################
## DiGpaths.txt: Save this file as  DiGpaths.txt to use it,          #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
#then type: read `DiGpaths.txt`: <Enter>                             #
## Then follow the instructions given there                          #
##                                                                   #
## Written by George Spahn and Doron Zeilberger, Rutgers University ,#
##  DoronZeil@gmail.com      .                                       # 
######################################################################

with(combinat):

Digits:=100:

print(`First Written: Oct. 2022: tested for Maple 20  `):
print():
print(`This is DiGpaths.txt, A Maple package`):
print(`accompanying George Spahn and Doron Zeilberger's  article: `):
print(` Variations on the Missionaries and Cannibals Puzzle `):
print(`-------------------------------`):

print(`-------------------------------`):
print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/DiGpaths.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`-------------------------------`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):

print(`-------------------------------`):
print(`For a list of the supporting functions type: ezra1();`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print():
 print(`-------------------------------`):


ezra1:=proc()
if args=NULL then

print(` The supporting procedures are:   `):
print(` AM, AMs, AMs1, CMg1, coin1, IsMC, KidsP, NuPaths, Paths,  PathsSA, WtPaths, WtPaths1 `):

else
 ezra(args):
fi:
end:

ezra:=proc()

if args=NULL then

print(` DiGpaths.txt: A Maple package for finding paths in Directed Graphs from the start vertex to the end vertex, with applications to Missionaries and Cannibals riddles`):
print(`The main procedures are :  MCgraph, PathsG, RG, SeqrBd, SolveMC, SSp,SSpLA `):


elif nargs=1 and args[1]=AM then
print(`AM(G): The adjacency matrix of the directed graph G`):

elif nargs=1 and args[1]=AMs then
print(`AMs(G,a): The symbolic adjacency matrix of the directed graph G with a[i,j] meaning an edge from i to j. Try:`):
print(`AMs([{1,2},{2,3},{1,2}],a);`):

elif nargs=1 and args[1]=AMs1 then
print(`AMs1(G,a,r): The symbolic adjacency matrix of the directed graph G with a[i,j,r], meaning an edge from i to j at step r. Try:`):
print(`AMs([{1,2},{2,3},{1,2}],a,4);`):

elif nargs=1 and args[1]=CMg1 then
print(`CMg1(M,C,B,d): inputs posiitve integers M, C, B,d and finds all legal vertices and edges. Try:`):
print(`CMg1(3,3,2,0);`):

elif nargs=1 and args[1]=coin1 then
print(`coin1(p): Inputs a positive rational number p between 0 and 1 and outputs 1 with prob. p and 0 with prob. 1-p. Try:`):
print(`coin1(1/3);`):


elif nargs=1 and args[1]=IsMC then
print(`IsMC(p,M,C,d): Is the pair p=[i,j] a legal vertex in a Missinaries and Cannibals puzzle?. Try:`):
print(`IsMC([1,1],3,3,0);`):

elif nargs=1 and args[1]=KidsP then
print(`KidsP(G,p): Given a directed graph given as a list of sets where L[i] are the vertices reachable by an edge from i, and a path,`):
print(`finds the set of continuations. Try:`):
print(`KidsP([{2,3},{3},{1,4},{1}],[1,2]);`):


elif nargs=1 and args[1]=MCgraph then
print(`MCgraph(M,C,B,d): The Missionaries and Cannibals graph with M in our data structure followed by the code. Try:`):
print(`The first vertex is the initial position [M,C,1] and the last one is the final position, [0,0,0]`):
print(`MCgraph(3,3,2,0);`):

elif nargs=1 and args[1]=NuPaths then
print(`NuPaths(G,i,j,k): Inputs a directed graph G vertices i and j (between 1 and nops(G)) and a positive inteter k, outputs the`):
print(`number of paths from i to j with k steps. Try:`):
print(`G:=RG(10,1/2); NuPaths(G,1,10,4);`):

elif nargs=1 and args[1]=Paths then
print(`Paths(G,i,k): inputs a directed graph G an edge i and a positive integer k, all the paths of length k that start at i. Try:`):
print(`G:=RG(10,1/2): Paths(G,1,5);`):

elif nargs=1 and args[1]=PathsG then
print(`PathsG(G): Inputs a directed graph G and outputs all the shortest paths from vertex 1 to vertex nops(G). Try:`):
print(`G:=RG(30,1/10); PathsG(G);`):
print(`G:=MCgraph(3,3,2,0)[1]; PathsG(G);`):

elif nargs=1 and args[1]=PathsSA then
print(`PathsSA(G,i,k): inputs a directed graph G an edge i and a positive integer k, all the Self-avoiding paths of length k that start at i. Try:`):
print(`G:=RG(10,1/2): PathsSA(G,1,5);`):

elif nargs=1 and args[1]=RG then
print(`RG(n,p): A random directed graph with prob. p for an edge. Try:`):
print(`RG(30,1/3);`):



elif nargs=1 and args[1]=SeqrBd then
print(`SeqrBd(r,B,d,K): The sequence`):
print(`[SSp(MCgraph(i+r,i,B,d)[2],i=1..K)]`):
print(`In other words the number of solutions with i+r missionaries, i cannibals, Boat size B, and safety margin d for i from 1 to K. Try:`):
print(`SeqrBd(5,3,1,5);`):

elif nargs=1 and args[1]=SolveMC then
print(`SolveMC(M,C,B,d): Gives ONE solution of the Missionaries and Cannibals problem with M missionaries, C cannibals, boat capacity B, and safety margin d. Try:`):
print(`SolveMC(3,3,2,0);`):

elif nargs=1 and args[1]=SSp then
print(`SSp(G): inputs a directed graph G, uses explicit construction to find the length of the shortest path from 1 to nops(G) and how many there are. Try`):
print(`SSp(MCgraph(3,3,2,0)[1]);`):


elif nargs=1 and args[1]=SSpLA then
print(`SSpLA(G): inputs a directed graph G, uses Linear Algebrato find the length of the shortest path from 1 to nops(G) and how many there are. `):
print(`Same output as SSp(G). Usually much slower, just for fun. Try:`):
print(`SSpLA(MCgraph(3,3,2,0)[1]):`):


elif nargs=1 and args[1]=WtPaths then
print(`WtPaths(G,i,j,k,a): Inputs a directed graph G vertices i and j (between 1 and nops(G)) and a positive inteter k, and a symbol a, outputs the`):
print(`weight-enumerator of the number of paths from i to j with k steps, where the weight of an edge from i to j is a[i,j]. Try:`):
print(`G:=RG(10,1/2); WtPaths(G,1,10,4,a);`):

elif nargs=1 and args[1]=WtPaths1 then
print(`WtPaths1(G,i,j,k,a): Inputs a directed graph G vertices i and j (between 1 and nops(G)) and a positive inteter k, and a symbol a, outputs the`):
print(`weight-enumerator of the number of paths from i to j with k steps, where the weight of an edge from i to j at step r is a[i,j,r]. Try:`):
print(`G:=RG(10,1/2); WtPaths1(G,1,10,4,a);`):

else
 print(`There is no such thing as`, args):

fi:



end:

with(linalg):

with(combinat):

#KidsP(G,p): Given a directed graph given as a list of sets where L[i] are the vertices reachable by an edge from i, and a path,
#finds the set of continuations. Try:
#KidsP([{2,3},{3},{1,4},{1}],[1,2]);
KidsP:=proc(G,p) local e,j:
e:=G[p[nops(p)]]:
{seq([op(p),j], j in e)}:
end:



#coin1(p): Inputs a positive rational number p between 0 and 1 and outputs 1 with prob. p and 0 with prob. 1-p. Try:
#coin1(1/3);
coin1:=proc(p) local m,n,ra:
m:=numer(p):
n:=denom(p):
ra:=rand(1..n)():
if ra()<=m then
 RETURN(1):
else
 RETURN(0):
fi:
end:

#RG(n,p): A random directed graph with prob. p for an edge
RG:=proc(n,p) local L,i,j,L1:

L:=[]:

for i from 1 to n do
L1:={}:
 for j from 1 to n do
  if coin1(p)=1 then
   L1:=L1 union {j}:
  fi:
 od:
L:=[op(L),L1]:
od:
L:

end:



#Paths(G,i,k): inputs a directed graph G an edge i and a positive integer k, all the paths of length k that start at i. Try:
#G:=RG(10,1/2): Paths(G,1,5);
Paths:=proc(G,i,k) local gu,gu1:
option remember:
if k=0 then
 RETURN({[i]}):
fi:
gu:=Paths(G,i,k-1):
{seq(op(KidsP(G,gu1)), gu1 in gu)}:
end:




#KidsP1(G,p): Given a directed graph given as a list of sets where L[i] are the vertices reachable by an edge from and a path,
#finds the set of continuations. i, not using a previous vertex. Try:
#KidsP1([{2,3},{3},{1,4},{1}],[1,2]);
KidsP1:=proc(G,p) local e,j:
e:=G[p[nops(p)]] minus {op(p)}:
{seq([op(p),j], j in e)}:
end:



#PathsSA(G,i,k): inputs a directed graph G an edge i and a positive integer k, all the Self-avoiding paths of length k that start at i. Try:
#G:=RG(10,1/2): PathsSA(G,1,5);
PathsSA:=proc(G,i,k) local gu,gu1:
option remember:
if k=0 then
 RETURN({[i]}):
fi:
gu:=PathsSA(G,i,k-1):
{seq(op(KidsP1(G,gu1)), gu1 in gu)}:
end:


#AM(G): The adjacency matrix of the directed graph G
AM:=proc(G) local n,i,j,M,M1:
M:=[]:

n:=nops(G):

for i from 1 to n do
 M1:=[]:
 for j from 1 to n do
   if member(j,G[i]) then
    M1:=[op(M1),1] :
    else
  M1:=[op(M1),0] :
   fi:
 od:
  M:=[op(M),M1]:
od:

M:
end:




#AMs(G,a): The symbolic adjacency matrix of the directed graph G with a[i,j] meaning an edge from i to j
AMs:=proc(G,a) local n,i,j,M,M1:
M:=[]:

n:=nops(G):

for i from 1 to n do
 M1:=[]:
 for j from 1 to n do
   if member(j,G[i]) then
    M1:=[op(M1),a[i,j]] :
    else
  M1:=[op(M1),0] :
   fi:
 od:
  M:=[op(M),M1]:
od:

M:
end:

#NuPaths(G,i,j,k): Inputs a directed graph G vertices i and j (between 1 and nops(G)) and a positive inteter k, outputs the
#number of paths from i to j with k steps. Try:
#G:=RG(10,1/2); NuPaths(G,1,10,4);
NuPaths:=proc(G,i,j,k) local M:
M:=AM(G):
evalm(M&^k)[i,j]:
end:


#WtPaths(G,i,j,k,a): Inputs a directed graph G vertices i and j (between 1 and nops(G)) and a positive inteter k, and a symbol a, outputs the
#weight-enumerator of the number of paths from i to j with k steps, where the weight of an edge from i to j is a[i,j]. Try:
#G:=RG(10,1/2); WtPaths(G,1,10,4,a);
WtPaths:=proc(G,i,j,k,a) local M:
M:=AMs(G,a):
expand(evalm(M&^k)[i,j]):
end:





#AMs1(G,a,r): The symbolic adjacency matrix of the directed graph G with a[i,j,r], meaning an edge from i to j at step r. 
AMs1:=proc(G,a,r) local n,i,j,M,M1:
M:=[]:

n:=nops(G):

for i from 1 to n do
 M1:=[]:
 for j from 1 to n do
   if member(j,G[i]) then
    M1:=[op(M1),a[i,j,r]] :
    else
  M1:=[op(M1),0] :
   fi:
 od:
  M:=[op(M),M1]:
od:

M:
end:



#WtPaths1(G,i,j,k,a): Inputs a directed graph G vertices i and j (between 1 and nops(G)) and a positive inteter k, and a symbol a, outputs the
#weight-enumerator of the number of paths from i to j with k steps, where the weight of an edge from i to j at step r is a[i,j,r]. Try:
#G:=RG(10,1/2); WtPaths1(G,1,10,4,a);
WtPaths1:=proc(G,i,j,k,a) local M,r:
M:=AMs1(G,a,1):
for r from 2 to k do
M:=expand(evalm(M&*AMs1(G,a,r))):
od:

expand(M[i,j]):
end:


#IsMC(p,M,C,d): Is the pair p=[i,j] a legal vertex in a Missinaries and Cannibals puzzle?
IsMC:=proc(p,M,C,d) local i,j:
i:=p[1]:
j:=p[2]:

if not i>=0 and i<=M and j>=0 and j<=C then
 RETURN(false):
fi:

if i>0 and i-j<d then
 RETURN(false):
fi:

if M-i>0 and (M-i)-(C-j)<d then
 RETURN(false):
fi:
true:
end:



#CMg1(M,C,B,d): inputs posiitve integers M, C, B,d and finds all legal vertices and edges. Try:
#CMg1(3,3,2,0);
CMg1:=proc(M,C,B,d) local V,E,i,j,e1,e2:

V:={}:

for i from 0 to M do
 for j from 0 to C do
  if IsMC([i,j],M,C,d) then
    V:=V union {[i,j,0],[i,j,1]}:
  fi:
 od:
od:
E:={}:

for  e1 from 0 to B do
 for e2 from 0 to B-e1 do
  if (e1>0 or e2>0)  and not (e1>0 and e2>0 and e1-e2<d) then
   E:=E union {[e1,e2]}:
  fi:
 od:
od:

[V,E]:
end:



#MCgraph(M,C,B,d): The Missionaries and Cannibals graph with M in our data structure followed by the code. Try:
#The first vertex is the initial position [M,C,1] and the last one is the final position, [0,0,0]
#MCgraph(3,3,2,0);
MCgraph:=proc(M,C,B,d) local lu,V,V1,E,T,e,nei,ID1,i:
lu:=CMg1(M,C,B,d):
V:=convert(lu[1],set):
V1:=V minus  {[M,C,1],[0,0,0]}:
V:=[[M,C,1],op(V1),[0,0,0]]:
for i from 1 to nops(V) do
 ID1[V[i]]:=i:
od:
E:=lu[2]:

for i from 1 to nops(V) do
T[i]:={}:

if V[i][3]=1 then

 for e in E do
  nei:=[V[i][1]-e[1], V[i][2]-e[2],0]:
  if member(nei, {op(V)}) then
   T[i]:=T[i] union {ID1[nei]}:
  fi:
 od:

elif V[i][3]=0 then


 for e in E do
  nei:=[V[i][1]+e[1], V[i][2]+e[2],1]:
  if member(nei, {op(V)}) then
   T[i]:=T[i] union {ID1[nei]}:
  fi:
 od:
fi:
od:

[seq(T[i],i=1..nops(V))],V:

end:


#PathsG(G): Inputs a directed graph G and outputs all the shortest paths from vertex 1 to vertex nops(G).
PathsG:=proc(G) local gu,k,n,gu1,lu:
n:=nops(G):


for k from 1 to n do
 gu:=PathsSA(G,1,k):
 if gu={} then
   RETURN({}):
 fi:

if member( n,{seq(gu1[nops(gu1)],gu1 in gu)}) then
 lu:={}:
for gu1 in gu do
 if gu1[nops(gu1)]=n then
  lu:=lu union {gu1}:
 fi:
od:
RETURN(lu):
fi:

od:

{}:

end:


#SSpLA(G): Using linear algebra to find the length of the shortest path from 1 to nops(G) and how many there are. Try
#SSpLA(MCgraph(3,3,2,0)[1]):
SSpLA:=proc(G) local M,k,n,M1:
M:=AM(G):
n:=nops(G):
M1:=M:

for k from 1 to n do

 if M1[1,n]<>0 then
  RETURN([k, M1[1,n]]):
 fi:

 M1:=evalm(M&*M1):
od:

FAIL:
end:



#SSpLA(G): inputs a directed graph G, uses explicit construction to find the length of the shortest path from 1 to nops(G) and how many there are. Try
#SSp(MCgraph(3,3,2,0)[1]):
SSp:=proc(G) local gu:
gu:=PathsG(G):
if gu={} then
 RETURN(FAIL):
else
[nops(gu[1])-1,nops(gu)]:
fi:
end:


#SolveMC(M,C,B,d): Gives ONE solution of the Missionaries and Cannibals problem with M missionaries, C cannibals, boat capacity B, and safety margin d. Try:
#SolveMC(3,3,2,0);
SolveMC:=proc(M,C,B,d) local G,gu,L,i:

G:=MCgraph(M,C,B,d):
L:=G[2]:
G:=G[1]:
gu:=PathsG(G):
if gu={} then
 RETURN(FAIL):
fi:
gu:=gu[1]:
[seq(L[gu[i]],i=1..nops(gu))]:
end:


#SeqrBd(r,B,d,K): The sequence
#[SSp(MCgraph(i+r,i,B,d)[2],i=1..K)]:
#In other words the number of solutions with i+r missionaries, i cannibals, Boat size B, and safety margin d for i from 1 to K. Try:
#SeqrBd(5,3,1,5);
SeqrBd:=proc(r,B,d,K) local i,gu,lu:
gu:=[]:
for i from 1 to K do
lu:=SSp(MCgraph(i+r,i,B,d)[1]):
if lu=FAIL then
 lu:=0:
fi:
gu:=[op(gu),lu[2]]:
od:
 gu:
end: