###Generalized Gambler’s Ruin Problem###
#Victoria Chayes, Chun Lau, Yuxuan Yang, Yukun Yao and Doron Zeilberger


##########################################################################
## GamblersRuin.txt Save this file as #GamblersRuin.txt                 ##
## to use it, stay in the same directory, get into Maple                ##
## (by typing: maple <Enter> ) and then type:                           ##
## read GamblersRuin.txt <Enter>                                        ##
## Then follow the instructions given there                             ##
##                                                                      ##
## Written by Victoria Chayes, Chun Lau, Yuxuan Yang, Yukun Yao         ##
## and Doron Zeilberger, Rutgers University                             ## 
## vc362@math.rutgers.edu                                               ##
## larryl@cs.rutgers.edu                                                ##
## yy458@math.rutgers.edu                                               ##
## yao@math.rutgers.edu                                                 ##
## doronzeil@gmail.com                                                  ##
##########################################################################

with(combinat):
with(plots):
with(Statistics):

print(`First Written: May. 2020`):
print(`Version 2020:`):
print():
print(`This is GamblersRuin.txt, A Maple package`):
print(`accompanying by Victoria Chayes, Chun Lau, Yuxuan Yang, Yukun Yao and Doron Zeilberger’s article`):   
print(`”Generalized Gambler’s Ruin Problem” `):

print():
print(`The most current version is available on WWW at:`):
print(`https://sites.math.rutgers.edu/~zeilberg/EM20/projs.html .`):
print(`Please report all bugs to: yao at math dot rutgers dot edu .`):
print():
print(`For an overview of the MAIN functions type “Help()”`):


Help:=proc():
print(`  StatAnal(P,t,K) , ProbWinning(N,p), ExpDuration(N,p), PGF(N,p,t), `): 
print(`  ProbWinningCF(N,n,p)  , ExpDurationCF(N,n,p), PGFcf(N,n,p,t) `):
print(` StatAnalPure(P,t,K), StatAnalPureL(P,t,K), GP1(L,n,d), GP(L,n), MomsGR1(n,N,K) , MomsGR(n,N,K,N1,N2) `):
print(`  CheckStra(L),TiS(N), BoS(N), RaS(N), GRsim1(L,p,n) , GRsim(L,p,n,M)   `): 
print(` ProbWinningStr(L,p), ExpDurationStr(L,p), PGFStr(L,p,t), RW(d,K), PolCons(d,K,M) `):
print(` PlotGR1(M,N,p,q,S,iter),PlotGR2(M,N,p,q,g,S,iter),PWNGR1(M,N,p,q,n),PWNGR2(M,N,p,q,n),EDNGR1(M,N,p,q),EDNGR2(M,N,p,q,g) `):
print(` TestEDdiff(M,N,p,q),PGFDNGR1(M,N,p,q,t),PGFDNGR2(M,N,p,q,g,t),TestPGdiff(M,N,p,q),ProbWinning2D(M,N,S,p,q,node)`):
print(` ProbWinningGraph(A,SetOfLosingSites,SetOfWinningSites), ExpDurationGraph(A,SetOfLosingSites,SetOfWinningSites), PGFGraph(A,SetOfLosingSites,SetOfWinningSites,t) `):
print(` LineGraph(N, p), ExpDurationR(N,n,p1,p2,p3,p4), Torus(m,n,p,q,r,s), ProbWinningTR(M,N,p,q,r,s,SetOfLosingSites,SetOfWinningSites), ExpDurationTR(M,N,p,q,r,s,SetOfLosingSites,SetOfWinningSites), PGFTR(M,N,p,q,r,s,SetOfLosingSites,SetOfWinningSites,t) `):
print(` GWinProbOne(N1,N2,N3,p,q,r), GExpDurOne(N1,N2,N3,p,q,r), GPGFExpDurOne(N1,N2,N3,p,q,r) `):
print(` GWinProbAll(N1,N2,N3,p,q,r), GExpDurAll(N1,N2,N3,p,q,r), GPGFExpDurAll(N1,N2,N3,p,q,r) `):

end:



######
#Exact expressions for the probability of winning  expected duration and higher moments of duration
######


#modified to get around "division by 0" by taking the limit as t goes to 1
#StatAnal(P,t,K): Inputs a polynomial (or even a function whose Taylor expansion has positive coefficients)
#describing the (not-necessarily normalized) probability generating function of some discrete random variable
#outputs  a list of length K whose
#(i) First entry is the expectation (alias mean, alias average)
#(ii) Second entry is the variance (aka as second moment about the mean)
#(iii) third-through K entries are the scaled moments about the mean
#in particular the third entry is the skewness and the fourth moment is the kurtosis. For example, try:
#StatAnal((1+t)^1000,t,6);
StatAnal:=proc(P,t,K) local P1,L,k,mu:

#First let's normalize it (just in case)
if limit(P,t=1)=0 then
 RETURN(FAIL):
else
 P1:=P/limit(P,t=1):
fi:

#mu is the average
mu:=limit(diff(P1,t),t=1):
L:=[mu]:

#Now let's CENTRALIZE, get the prob. gen. function for the "deviation from the mean"
P1:=P1/t^mu:

#We will now apply the operation td/dt repeatedly and plug-in t=1
P1:=t*diff(P1,t):

#The expectation of the "deviation from the mean" should be 0, let's check it
if normal(limit(P1,t=1))<>0 then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

for k from 2 to K do
P1:=t*diff(P1,t):
#We append the current moment-about the mean
L:=[op(L),limit(P1,t=1)]:
od:

#Finally we adjust the third-through the K-th moment by dividing the k-th moment by
#the L[2]^(k/2) (the k-th power of the standard deviation)

if L[2]=0 then
 print(`The variance is 0, the unscaled moments are`):
 print(L):
 RETURN(FAIL):
else
[L[1],L[2],seq(L[k]/L[2]^(k/2),k=3..K)]:
fi:

end:


#ProbWinning(N,p): inputs a positive integer N and a prob. p (either numeric or symbolic),
#outputs a list of length N-1 whose i-th entry is the probability of exiting a winner in
#Gambler's ruin game where you start with a capital of i dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ProbWinning(10,1/2);
#ProbWinning(10,1/3);
#ProbWinning(10,2/3);
#ProbWinning(10,p);
ProbWinning:=proc(N,p) local eq,var,x,i:

#We call the prob. of winning if you are currently with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently in possession of i dollars, then with probability p next round you would have i+1 dollars
#and with prob. 1-p next time you would have i-1 dollars. By CONDITIONING on the event that you
#currently have i dollars, we get the equation
#x[i]=p*x[i+1]+(1-p)*x[i-1]. 
#Also if right now you have 0 dollars your prob. of winning is 0
#and if you currently have N dollars than your prob. of winning the game is 1 (THESE ARE THE BOUNDARY CONDITIONS)

eq:={seq(x[i]=p*x[i+1]+(1-p)*x[i-1],i=1..N-1),x[0]=0,x[N]=1}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

[seq(subs(var,x[i]),i=1..N-1)]:

end:




#ExpDuration(N,p): inputs a positive integer N and a prob. p (either numeric or symbolic),
#outputs a list of length N-1 whose i-th entry is the expected duration of the a
#Gambler's ruin game where you start with a capital of i dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ExpDuration(10,1/2);
#ExpDuration(10,1/3);
#ExpDuration(10,2/3);
#ExpDuration(10,p);
ExpDuration:=proc(N,p) local eq,var,x,i:

#We call the expected remaining duration if you currently have  with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently in possession of i dollars, then with probability p next round you would have i+1 dollars
#and with prob. 1-p next time you would have i-1 dollars. By CONDITIONAL EXPECTATION on the event that you
#currently have i dollars, we get the equation
#x[i]=p*x[i+1]+(1-p)*x[i-1]+1
#(Note that the +1 accounts that the duration today is one more than the duration tomorrow)
#Also if right now you have 0 dollars or N dollars your duration (and hence expected duration) is 0
eq:={seq(x[i]=p*x[i+1]+(1-p)*x[i-1]+1,i=1..N-1),x[0]=0,x[N]=0}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

[seq(subs(var,x[i]),i=1..N-1)]:

end:




#PGF(N,p,t): inputs a positive integer N and a prob. p (either numeric or symbolic),
#and a variable name t,
#outputs a list of length N-1 whose i-th entry is the PROBABILITY GENERATING FUNCTION FOR DURATION
#(the rational function in t whose coefficient of t^i is the probability that the game will end exactly after i round)
#In a Gambler's ruin game where you start with a capital of i dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#PGF(10,1/2,t);
#PGF(10,1/3,t);
#PGF(10,2/3,t);
#PGF(10,p,t);
PGF:=proc(N,p,t) local eq,var,x,i:

#We call the prob. generating functions for the remaining duration if you currently have  with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently have i dollars, then with probability p next round you would have i+1 dollars
#and with prob. 1-p next time you would have i-1 dollars. By CONDITIONAL EXPECTATION on the event that you
#currently have i dollars, we get the equation
#x[i]=t*(p*x[i+1]+(1-p)*x[i-1])
#(Note that the multiplication by t  accounts that the duration today is one more than the duration tomorrow)
#Also if right now you have 0 dollars or N dollars your duration the prob. generating function is t^0=1 
#(I will definitely take you 0 rounds to finish the game)
eq:={seq(x[i]=t*(p*x[i+1]+(1-p)*x[i-1]),i=1..N-1),x[0]=1,x[N]=1}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

normal([seq(subs(var,x[i]),i=1..N-1)]):

end:



#ProbWinningCF(N,n,p): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic),
#outputs the CLOSED-FORM EXPRESSION for the prob. of winning in
#Gambler's ruin game where you start with a capital of n dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ProbWinningCF(N,n,p);
ProbWinningCF:=proc(N,n,p) local x,F,A,A0:
#We let Maple solve the recurrence equation x(n)=(1-p)*x(n-1)+p*x(n+1), but it can only handle
#INITIAL conditions, and so far we have no clue what x(1) is, so we call it A

F:=rsolve({x(n)=(1-p)*x(n-1)+p*x(n+1),x(0)=0, x(1)=A},x(n)):

#Now we use the fact that x(N)=1, and solve for A
A0:=solve(subs(n=N,F)=1,A):

#Now we plug that expression for A0 into F

normal(simplify(subs(A=A0,F))):

end:


#ExpDurationCF(N,n,p): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic),
#outputs the CLOSED-FORM EXPRESSION for the Expected duration
#Gambler's ruin game where you start with a capital of n dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ExpDurationCF(N,n,p);
ExpDurationCF:=proc(N,n,p) local x,F,A,A0:
#We let Maple solve the recurrence equation x(n)=(1-p)*x(n-1)+p*x(n+1)+1, but it can only handle
#INITIAL conditions, and so far we have no clue what x(1) is, so we call it A

F:=rsolve({x(n)=(1-p)*x(n-1)+p*x(n+1)+1,x(0)=0, x(1)=A},x(n)):

#Now we use the fact that x(N)=0, and solve for A
A0:=solve(subs(n=N,F)=0,A):

#Now we plug that expression for A0 into F

normal(simplify(subs(A=A0,F))):

end:


#PGCcf(N,n,p,t): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic), and a variable name (symbol) t;
#outputs the CLOSED-FORM EXPRESSION for the Prob. Generating function in t, of the duration
#Gambler's ruin game where you start with a capital of n dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#PGFcf(N,n,p,t);
PGFcf:=proc(N,n,p,t) local x,F,A,A0:
#We let Maple solve the recurrence equation x(n)=t*((1-p)*x(n-1)+p*x(n+1)), but it can only handle
#INITIAL conditions, and so far we have no clue what x(1) is, so we call it A


F:=rsolve({x(n)=t*((1-p)*x(n-1)+p*x(n+1)),x(0)=1, x(1)=A},x(n)):

#Now we use the fact that x(N)=1, and solve for A
A0:=solve(subs(n=N,F)=1,A):

#Now we plug that expression for A0 into F

normal(simplify(subs(A=A0,F))):

end:



#StatAnalPureL(P,t,K): The expecation, variance, and the UNSCALDED moments about the mean
#the first entry is the expectation, the second the variance, and then the third, ..., up to the
#K-th moment about the mean. 
#It uses limit(...,t=1) instead of subs(t=1, ...) MUCH SLOWER
#Try:
#StatAnalPureL((1+t)^n,t,10);
StatAnalPureL:=proc(P,t,K) local P1,L,k,mu:

#First let's normalize it (just in case)
if limit(P,t=1)=0 then
 RETURN(FAIL):
else
 P1:=P/limit(P,t=1):
fi:

#mu is the average
mu:=limit(diff(P1,t),t=1):
L:=[mu]:

#Now let's CENTRALIZE, get the prob. gen. function for the "deviation from the mean"
P1:=P1/t^mu:

#We will now apply the operation td/dt repeatedly and plug-in t=1
P1:=t*diff(P1,t):

#The expectation of the "deviation from the mean" should be 0, let's check it
if normal(limit(P1,t=1))<>0 then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

for k from 2 to K do
P1:=t*diff(P1,t):
#We append the current moment-about the mean
L:=[op(L),limit(subs(P1),t=1)]:
od:
L:

end:


#StatAnalPure(P,t,K): The expecation, variance, and the UNSCALDED moments about the mean
#the first entry is the expectation, the second the variance, and then the third, ..., up to the
#K-th moment about the mean. Try:
#StatAnalPure((1+t)^n,t,10);
StatAnalPure:=proc(P,t,K) local P1,L,k,mu:

#First let's normalize it (just in case)
if subs(t=1,P)=0 then
 RETURN(FAIL):
else
 P1:=P/subs(t=1,P):
fi:

#mu is the average
mu:=subs(t=1,diff(P1,t)):
L:=[mu]:

#Now let's CENTRALIZE, get the prob. gen. function for the "deviation from the mean"
P1:=P1/t^mu:

#We will now apply the operation td/dt repeatedly and plug-in t=1
P1:=t*diff(P1,t):

#The expectation of the "deviation from the mean" should be 0, let's check it
if normal(limit(P1,t=1))<>0 then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

for k from 2 to K do
P1:=t*diff(P1,t):
#We append the current moment-about the mean
L:=[op(L),subs(t=1,P1)]:
od:
L:

end:


######
#Guessing
######

#GP1(L,n,d): inputs a list L , a variable name n, and a positive integer d, guesses a polynomial of degree <=d in n,
#let's call it P(n) such that P(i)=L[i]. The length of L must be at least d+3. If it fails it returns FAIL.
#Try:
#GP1([seq(i^3,i=1..10)],n,3);
GP1:=proc(L,n,d) local a,i,var,eq,P:

#Any d+1 data points can fit a polynomial of degree d. We want at least 2 "confirmation" points
if nops(L)<=d+2 then
 print(`The list must be of length at least`, d+3):
 RETURN(FAIL):
fi:

#We define a generic polynomial P of n, of degree d, with UNDETERMINED coefficients. Our goal is to determine them

P:=add(a[i]*n^i,i=0..d):

#This is the set of unknowns with d+1 unknowns
var:={seq(a[i],i=0..d)}:

#We set-up a system with d+3 equations fitting the polynomial to the data
eq:={seq(subs(n=i,P)-L[i],i=1..d+3)}:

#Maple tries to solve the system, we rename var the solution
var:=solve(eq,var):

#If var=NULL then there is no solution, in other words, no polynomial of degree d fits the data

if var=NULL then
 RETURN(FAIL):
fi:

#Now we plug the solution var into the generic polynomial P, making it an actual polynomial (that fits the data)
P:=subs(var,P):

#If nops(L) is larger than d+3, we check that the polynomial fits the remaining data
for i from d+4 to nops(L) do
 if expand(subs(n=i,P))-L[i]<>0 then
  print(P, `does not agree with n=`, i):
  RETURN(FAIL):
 fi:
od:

#If it agrees we return the polynomial
P:

end:


  

#GP(L,n): inputs a list L , a variable name n, guesses a polynomial of degree <=nops(L)-3 in n,
#let's call it P(n) such that P(i)=L[i] for all i. If it fails it returns FAIL.
#Try:
#GP([seq(i^3,i=1..10)],n);

GP:=proc(L,n) local d,P:

for d from 0 to nops(L)-3 do
#We start with d=0 and keep trying with higher and higher degree
 P:=GP1(L,n,d):
 if P<>FAIL then
  RETURN(P):
 fi:
od:
#If there is no such polynomial of degree <=nopd(L)-3 it returns FAIL
FAIL:
end:


#MomsGR1(n,N,K): inputs a variable name n, and a positive integer N (sufficiently large) and a positive integer K
#tries to find polynomial expression in n for the first K moments of the r.v. duration of a Gambler's
#Ruin game with a fair coin with exit capital N and starting capital n. If it can only partially do it,
#it returns what it can. Otherwise it returns all the K moments
#Try:
#MomsGR1(n,20,6);

MomsGR1:=proc(n,N,K) local t,L,i,M,j,P:
#We use the semi-numerical, semi-symbolic procedure (that uses simple linear algebra)
#to find the list of N-1 prob. generating functions for the duration stating with n dollars
L:=PGF(N,1/2,t):
#We rename L to be the list of lists each of length K, describing the first K moments for each of n=1..N-1

L:=[seq(StatAnalPure(L[i],t,K),i=1..nops(L))]:

M:=[]:

#We try to guess polynomial expressions for the first K moments one by one, each is a polynomial in n

for i from 1 to K do
 P:=GP([seq(L[j][i],j=1..nops(L))],n):
 #If we failed to get a polynomial expression for the i-th moment, we return what we got 

if P=FAIL then
   RETURN(M):
 fi:

#otherwise we append the new moment
 M:=[op(M),P]:
od:

#If all goes well we got polynomial expressions in n for the first K moments of the Gambler's ruin with exit capital N (N numeric)
M:
end:



#MomsGR(n,N,K,N1,N2): tries to conjecture polynomial expressions in n and N for the first K moments,
#N=N1 to N=N2 is  where we collect data
MomsGR:=proc(n,N,K,N1,N2) local L1,P,N3,M,i:

#L is the list of lists of length K from N1 to N2, we first make sure that with N=N1 we get the full K moments

L1:=MomsGR1(n,N1,K):

if nops(L1)<K then
 print(` Make `, N1, `larger `):
fi:


M:=[]:

for i from 1 to K do

#P is the  polynomial in N (and n) for the i-th moment

P:=GP([seq(MomsGR1(n,N3,K)[i],N3=N1..N2)],N):

#GP starts at 1, since here the data starts at N1 we have to change the variable
P:=expand(subs(N=N-N1+1,P)):

if P=FAIL then
  print(` Make `, N2, `bigger `):
  RETURN(M):
else
M:=[op(M),P]:
fi:
od:

M:
end:



with(plots):

#plotDist(f,x,K): Given a prob. gen. function f(x) that has a 
#Taylor series
#for a discrete r.v.
#plots its normalized version (X-mu)/sig between mu-K1*sig and
#mu+K2*sig
#For example, try:
#plotDist((1+x)^40,x,5,5);
plotDist:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1:
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
gu:=f1/x^mu:
sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))):

lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10):

lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]:

lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]):

plot(lu,scaling=constrained):

end:


######
#Strategic Gambling
######


#TiS(N): The timid strategy with maximum capital N
TiS:=proc(N) local i:
 [1$(N-1)]:
end:

#BoS(N): The bold strategy with maximum capital N
BoS:=proc(N) local i:
[seq(min(i,N-i),i=1..N-1)]:
end:

#RaS(N): a random legal strategy  maximum capital N
RaS:=proc(N) local i,L:
L:=BoS(N):
[seq(rand(1..L[i])(),i=1..N-1)]:
end:


#CheckStra(L): inputs a list L (of length nops(L)) snf lry N:=nops(L)+1, of positive integers indicating a gambling strategy
#such that that if you currently have i dollars how much you would stake (it has to be less than min(i,N-i)). Try:
#CheckStra([1,2,3,2,1]);
#CheckStra([1,2,4,2,1]);

CheckStra:=proc(L) local N,i:

if not type(L,list) then
 RETURN(false):
fi:


N:=nops(L)+1:

for i from 1 to nops(L) do
 if L[i]>min(L[i],N-L[i]) or L[i]<1 then
  RETURN(false):
 fi:
od:
true:
end:


#GRsim1(L,p,n): Simulating ONE Gambler's game with initial capital of n dollars,  and exit capital N
#and gambling strategy L, where nops(L)=N-1 and L[i] is the amount you would stake if you had i dollars
#At each turn the gambler's win the staked amount with probability p and loses it with probability 1-p
#the game ends as soon as it reaches N dollars or 0 dollars (when the gambler's is ruined).
#The output is a pair [0,m] or [1,m], according to whether the gambler's left the casino broke or
#rich, and m is the number of rounds. Try:
#GRsim1([1$9],1/2,5);
#GRsim1(10,10/21,5);
GRsim1:=proc(L,p,n) local N,x,U,coin,count:

N:=nops(L)+1:

if not (type(p,numeric) and p>=0 and p<=1 and type(n, integer) and n>=1 and n<=N) then
 print(`Bad input`):
 RETURN(FAIL):
fi:



if not CheckStra(L) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


U:=RandomVariable(Bernoulli(p)):

x:=n:
count:=0:
while x>0 and x<N do

count:=count+1:
coin:=Sample(U,1)[1]:

 if coin=0 then
  x:=x-L[x]:
 else
  x:=x+L[x]:
 fi:
od:

if x=0 then
 RETURN([0,count]):
elif x=N then
 RETURN([1,count]):
else
 print(`Something is wrong`):
 RETURN(FAIL):
fi:


end:



#GRsim(L,p,n,M): Inputs a gambling strategy L, prob. of winning n, runs GRsim1(L,p,n) M times
#outputs the ratio of winning and average duration.
#Try:
#GRsim(TiS(30),1/2,20,100);
GRsim:=proc(L,p,n,M) local i:

evalf(add(GRsim1(L,p,n),i=1..M)/M):
end:


#ProbWinningStr(L,p): inputs a  winning strategy L,  a prob. p (either numeric or symbolic),
#outputs a list of length nops(L) whose i-th entry is the probability of exiting a winner in
#a game following strategy L, where you start with a capital of i dollars, and then at each round
#you win  your stake, L[i] with prob. p and lose it  with prob. 1-p. Try:
#ProbWinningStr(TiS(10),1/2);
#ProbWinningStr(BoS(10),1/3);
ProbWinningStr:=proc(L,p) local N,eq,var,x,i:

if not CheckStra(L) then
 RETURN(FAIL):
fi:

N:=nops(L)+1:


#We call the prob. of winning if you are currently with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently in possession of i dollars, then with probability p next round you would have i+L[i] dollars
#and with prob. 1-p next time you would have i-L[i] dollars. By CONDITIONING on the event that you
#currently have i dollars, we get the equation
#x[i]=p*x[i+L[i]]+(1-p)*x[i-L[i]]. 
#Also if right now you have 0 dollars your prob. of winning is 0
#and if you currently have N dollars than your prob. of winning the game is 1 (THESE ARE THE BOUNDARY CONDITIONS)

eq:={seq(x[i]=p*x[i+L[i]]+(1-p)*x[i-L[i]],i=1..N-1),x[0]=0,x[N]=1}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

[seq(subs(var,x[i]),i=1..N-1)]:

end:




#ExpDurationStr(L,p): inputs a gambling strategy, L, (Let N:=nops(L)+1), and a prob. p (either numeric or symbolic),
#outputs a list of length N-1 whose i-th entry is the expected duration of the a
#Gambler's ruin game where you start with a capital of i dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ExpDurationStr(TiS(10),1/2);
#ExpDurationStr(BoS(10),1/3);
ExpDurationStr:=proc(L,p) local eq,var,x,i,N:

if not CheckStra(L) then
 RETURN(FAIL):
fi:

N:=nops(L)+1:

#We call the expected remaining duration if you currently have  with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently in possession of i dollars, then with probability p next round you would have i+L[i] dollars
#and with prob. 1-p next time you would have i-L[i] dollars. By CONDITIONAL EXPECTATION on the event that you
#currently have i dollars, we get the equation
#x[i]=p*x[i+L[i]]+(1-p)*x[i-L[i]]+1
#(Note that the +1 accounts that the duration today is one more than the duration tomorrow)
#Also if right now you have 0 dollars or N dollars your duration (and hence expected duration) is 0
eq:={seq(x[i]=p*x[i+L[i]]+(1-p)*x[i-L[i]]+1,i=1..N-1),x[0]=0,x[N]=0}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

[seq(subs(var,x[i]),i=1..N-1)]:

end:




#PGFStr(L,p,t): inputs a strategy L and a prob. p (either numeric or symbolic),
#and a variable name t,
#outputs a list of length N-1 whose i-th entry is the PROBABILITY GENERATING FUNCTION FOR DURATION
#(the rational function in t whose coefficient of t^i is the probability that the game will end exactly after i round)
#In a  game with strategy L, where you start with a capital of i dollars, and then at each round
#you win L[i] dollars with prob. p and lose them with prob. 1-p. Try:
#PGFStr(TiS(10),1/2,t);
#PGFStr(BoS(10),1/3,t);
PGFStr:=proc(L,p,t) local eq,var,x,i,N:

if not CheckStra(L) then
 RETURN(FAIL):
fi:

N:=nops(L)+1:
#We call the prob. generating functions for the remaining duration if you currently have  with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently have i dollars, then with probability p next round you would have i+L[i] dollars
#and with prob. 1-p next time you would have i-L[i] dollars. By CONDITIONAL EXPECTATION on the event that you
#currently have i dollars, we get the equation
#x[i]=t*(p*x[i+L[i]]+(1-p)*x[i-L[i]])
#(Note that the multiplication by t  accounts that the duration today is one more than the duration tomorrow)
#Also if right now you have 0 dollars or N dollars your duration the prob. generating function is t^0=1 
#(I will definitely take you 0 rounds to finish the game)
eq:={seq(x[i]=t*(p*x[i+L[i]]+(1-p)*x[i-L[i]]),i=1..N-1),x[0]=1,x[N]=1}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

normal([seq(subs(var,x[i]),i=1..N-1)]):

end:





######
#Random Walk
######

#RW(d,K): the number of steps for a  simple random walk in d-dimensional space until it returns to the origin
#it outputs the number of steps or FAIL
RW:=proc(d,K) local pt,co,i,j,ra1,ra2:
pt:=[0$d]:
ra1:=rand(1..d):
ra2:=rand(0..1):

for co from 1 to K do
i:=ra1():
j:=ra2():

if j=0 then
 pt:=[op(1..i-1,pt),pt[i]-1,op(i+1..d,pt)]:
else
 pt:=[op(1..i-1,pt),pt[i]+1,op(i+1..d,pt)]:
fi:

if pt=[0$d] then
 RETURN(co):
fi:

od:

FAIL:

end:

#PolCons(d,K,M): estimating Polya's random walk constant in d dimension, pretending that K is infinity
#it also estimates the expected time of return
#doing it M times
#Try:
#PolCons(3,10000,1000);
PolCons:=proc(d,K,M) local co1,co2,i,W:
co1:=0:
co2:=0:
for i from 1 to M do
  W:=RW(d,K):
   if W<>FAIL then
     co1:=co1+1:
      co2:=co2+W:
   fi:
od:

evalf([co1/M,co2/M]):

end:



################################
# VISUALISING 2D GABLER'S RUIN #
################################

#PlotGR1(M,N,p,q,S,iter): computes and plots successive iterations of the probability density of 2D Gambler's Ruin Game Type 1 (player simultaneously bets vertically and horizontally every round)
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# a list S=[s1, s2] of coordinates for the starting vertical and horizontal amounts
# an integer iter of number of rounds the game should run
# outputs for each round a matrix of probabilities of the amount of money the gambler has, and if p and q are numeric, a density plot

PlotGR1:=proc(M,N,p,q,S,iter) local A,i,B,B1,r,v,c,n,j,Bnew,f:
A:=Matrix(M,N,0):
A[S[1],S[2]]:=1:
print(A):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->A[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
B:=Matrix(M,N,0):
for i from 1 to iter do
  Bnew:=Matrix(M,N,0):
  if i=1 then
  r,c,v:=SearchArray(A):
  else
  r,c,v:=SearchArray(B):
  fi:
  n:=numelems(r):
  for j from 1 to n do
  B1:=Matrix(M,N,0):
    if r[j]=1 then
      if c[j]=1 then
         B1[1,1]:=v[j]:
      elif c[j]=N then
         B1[1,N]:=v[j]:
      else
         B1[1,c[j]+1]:=v[j]*q:
         B1[1,c[j]-1]:=v[j]*(1-q):
      fi:
    elif r[j]=M then
      if c[j]=1 then
         B1[M,1]:=v[j]:
      elif c[j]=N then
         B1[M,N]:=v[j]:
      else
         B1[M,c[j]+1]:=v[j]*q:
         B1[M,c[j]-1]:=v[j]*(1-q):
      fi:
    elif c[j]=1 then
         B1[r[j]+1,1]:=v[j]*p:
         B1[r[j]-1,1]:=v[j]*(1-p):
    elif c[j]=N then
         B1[r[j]+1,N]:=v[j]*p:
         B1[r[j]-1,N]:=v[j]*(1-p):
    else
         B1[r[j]+1,c[j]+1]:=v[j]*p*q:
         B1[r[j]-1,c[j]+1]:=v[j]*q*(1-p): 
         B1[r[j]+1,c[j]-1]:=v[j]*(1-q)*p:
         B1[r[j]-1,c[j]-1]:=v[j]*(1-p)*(1-q):
    fi:
   Bnew:=Bnew+B1:
od:
B:=Bnew:
print(B):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->B[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
od:
end:

#PlotGR2(M,N,p,q,g,S,iter): computes and plots successive iterations of the probability density of 2D Gambler's Ruin Game Type 2 (player bets horizontally with probability g, or vertically with probability (1-g), each round)
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# a numeric or symbolic probability g playing the vertical game each round 
# a list S=[s1, s2] of coordinates for the starting vertical and horizontal amounts
# an integer iter of number of rounds the game should run
# outputs for each round a matrix of probabilities of the amount of money the gambler has, and if p, q, and g are symbolic, a density plot of the matrix

PlotGR2:=proc(M,N,p,q,g,S,iter) local A,i,B,B1,r,v,c,n,j,Bnew,f,ra:
A:=Matrix(M,N,0):
A[S[1],S[2]]:=1:
print(A):
if type(p, numeric) and type(q, numeric) and type(g, numeric) then
f:=(x,y)->A[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
B:=Matrix(M,N,0):
for i from 1 to iter do
  Bnew:=Matrix(M,N,0):
  if i=1 then
  r,c,v:=SearchArray(A):
  else
  r,c,v:=SearchArray(B):
  fi:
  n:=numelems(r):
  for j from 1 to n do
  B1:=Matrix(M,N,0):
    if r[j]=1 then
      if c[j]=1 then
         B1[1,1]:=v[j]:
      elif c[j]=N then
         B1[1,N]:=v[j]:
      else
         B1[1,c[j]+1]:=v[j]*q:
         B1[1,c[j]-1]:=v[j]*(1-q):
      fi:
    elif r[j]=M then
      if c[j]=1 then
         B1[M,1]:=v[j]:
      elif c[j]=N then
         B1[M,N]:=v[j]:
      else
         B1[M,c[j]+1]:=v[j]*q:
         B1[M,c[j]-1]:=v[j]*(1-q):
      fi:
    elif c[j]=1 then
         B1[r[j]+1,1]:=v[j]*p:
         B1[r[j]-1,1]:=v[j]*(1-p):
    elif c[j]=N then
         B1[r[j]+1,N]:=v[j]*p:
         B1[r[j]-1,N]:=v[j]*(1-p):
    else
         B1[r[j]+1,c[j]]:=v[j]*g*p:
         B1[r[j]-1,c[j]]:=v[j]*g*(1-p): 
         B1[r[j],c[j]+1]:=v[j]*(1-g)*q:
         B1[r[j],c[j]-1]:=v[j]*(1-g)*(1-q):
    fi:
   Bnew:=Bnew+B1:
od:
B:=Bnew:
print(B):
if type(p, numeric) and type(q, numeric) and type(g, numeric) then
f:=(x,y)->B[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
od:
end:

#############################
# NUMERIC 2D GAMBLER'S RUIN #
#############################

#PWNGR1(M,N,p,q,n): computes and plots the probability of ending a game of 2D Gambler's Ruin Game Type 1 in each possible node from each potential starting position
#inputs: an integer M>1 for the win condition of $M-1 in the vertical game
# an integer N>1 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# an integer 1-5 indicating the output is for: 1: the node [1,1]; 2: the node [M,1]; 3: the node [1,N]; 4: the node [M,N]; or 5: all four nodes.
#outputs: a matrix with each entry indicating the probability of ending each game in the given node from that starting position, and if p and q are numeric, a density plot of said matrix. 

PWNGR1:=proc(M,N,p,q,n) local eq,var,x,i,j,sol,f:
if n=1 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,1]=1, seq(x[1,i]=q*x[1,i+1]+(1-q)*x[1,i-1],i=2..N-1),seq(x[i,1]=p*x[i+1,1]+(1-p)*x[i-1,1],i=2..M-1),seq(x[i,N]=0,i=1..M),seq(x[M,i]=0,i=1..N-1),seq(seq(x[i,j]=x[i+1,j+1]*p*q+x[i+1,j-1]*p*(1-q)+x[i-1,j+1]*(1-p)*q+x[i-1,j-1]*(1-p)*(1-q),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)', x=0.01..N, y=0.01..M, style=patchnogrid)):
fi:
fi:

if n=2 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[M,1]=1, seq(x[M,i]=q*x[M,i+1]+(1-q)*x[M,i-1],i=2..N-1),seq(x[i,1]=p*x[i+1,1]+(1-p)*x[i-1,1],i=2..M-1),seq(x[i,N]=0,i=1..M),seq(x[1,i]=0,i=1..N-1),seq(seq(x[i,j]=x[i+1,j+1]*p*q+x[i+1,j-1]*p*(1-q)+x[i-1,j+1]*(1-p)*q+x[i-1,j-1]*(1-p)*(1-q),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)', x=0.01..N, y=0.01..M, style=patchnogrid)):
fi:
fi:

if n=3 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,N]=1, seq(x[1,i]=q*x[1,i+1]+(1-q)*x[1,i-1],i=2..N-1),seq(x[i,N]=p*x[i+1,N]+(1-p)*x[i-1,N],i=2..M-1),x[1,1]=0,seq(x[i,1]=0,i=1..M), seq(x[M,i]=0,i=1..N),seq(seq(x[i,j]=x[i+1,j+1]*p*q+x[i+1,j-1]*p*(1-q)+x[i-1,j+1]*(1-p)*q+x[i-1,j-1]*(1-p)*(1-q),i=2..M-1),j=2..N-1) }:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
fi:

if n=4 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[M,N]=1, seq(x[M,i]=q*x[M,i+1]+(1-q)*x[M,i-1],i=2..N-1),seq(x[i,N]=p*x[i+1,N]+(1-p)*x[i-1,N],i=2..M-1),seq(x[i,1]=0,i=1..M),seq(x[1,i]=0,i=2..N),seq(x[1,i]=0,i=1..N-1), seq(seq(x[i,j]=x[i+1,j+1]*p*q+x[i+1,j-1]*p*(1-q)+x[i-1,j+1]*(1-p)*q+x[i-1,j-1]*(1-p)*(1-q),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
fi:
end:


#PWNGR2(M,N,p,q,n): computes and plots the probability of ending a game of 2D Gambler's Ruin Game Type 2 in each possible node from each potential starting position
#inputs: an integer M>1 for the win condition of $M-1 in the vertical game
# an integer N>1 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# a numeric or symbolic probability g playing the vertical game each round 
# an integer 1-5 indicating the output is for: 1: the node [1,1]; 2: the node [M,1]; 3: the node [1,N]; 4: the node [M,N]; or 5: all four nodes.
#outputs: a matrix with each entry indicating the probability of ending each game in the given node from that starting position, and if p and q are numeric, a density plot of said matrix. 

PWNGR2:=proc(M,N,p,q,g,n) local eq,var,x,i,j,sol,f:
if n=1 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,1]=1, seq(x[1,i]=q*x[1,i+1]+(1-q)*x[1,i-1],i=2..N-1),seq(x[i,1]=p*x[i+1,1]+(1-p)*x[i-1,1],i=2..M-1),seq(x[i,N]=0,i=1..M),seq(x[M,i]=0,i=1..N-1),seq(seq(x[i,j]=x[i+1,j]*p*g+x[i-1,j]*g*(1-p)+x[i,j+1]*(1-g)*q+x[i,j-1]*(1-g)*(1-q),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)', x=0.01..N, y=0.01..M, style=patchnogrid)):
fi:
fi:

if n=2 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[M,1]=1, seq(x[M,i]=q*x[M,i+1]+(1-q)*x[M,i-1],i=2..N-1),seq(x[i,1]=p*x[i+1,1]+(1-p)*x[i-1,1],i=2..M-1),seq(x[i,N]=0,i=1..M),seq(x[1,i]=0,i=1..N-1),seq(seq(x[i,j]=x[i+1,j]*p*g+x[i-1,j]*g*(1-p)+x[i,j+1]*(1-g)*q+x[i,j-1]*(1-g)*(1-q),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)', x=0.01..N, y=0.01..M, style=patchnogrid)):
fi:
fi:

if n=3 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,N]=1, seq(x[1,i]=q*x[1,i+1]+(1-q)*x[1,i-1],i=2..N-1),seq(x[i,N]=p*x[i+1,N]+(1-p)*x[i-1,N],i=2..M-1),x[1,1]=0,seq(x[i,1]=0,i=1..M), seq(x[M,i]=0,i=1..N),seq(seq(x[i,j]=x[i+1,j]*p*g+x[i-1,j]*g*(1-p)+x[i,j+1]*(1-g)*q+x[i,j-1]*(1-g)*(1-q),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
fi:

if n=4 or n=5 then 
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[M,N]=1, seq(x[M,i]=q*x[M,i+1]+(1-q)*x[M,i-1],i=2..N-1),seq(x[i,N]=p*x[i+1,N]+(1-p)*x[i-1,N],i=2..M-1),seq(x[i,1]=0,i=1..M),seq(x[1,i]=0,i=2..N),seq(x[1,i]=0,i=1..N-1), seq(seq(x[i,j]=x[i+1,j]*p*g+x[i-1,j]*g*(1-p)+x[i,j+1]*(1-g)*q+x[i,j-1]*(1-g)*(1-q),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
print(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
fi:
end:


#EDNGR1(M,N,p,q): computes and plots the expected duration of a game of 2D Gambler's Ruin Type 1 from each potential starting position
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
#outputs: a matrix with each entry indicating the expected duration of the game from that starting position, and if p and q are numeric, a density plot of said matrix. 

EDNGR1:=proc(M,N,p,q) local eq,var,sol,x,i,j,f:
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,1]=0, x[1,N]=0, x[M,1]=0, x[M,N]=0, seq(x[1,i]=q*x[1,i+1]+(1-q)*x[1,i-1]+1,i=2..N-1), seq(x[M,i]=q*x[M,i+1]+(1-q)*x[M,i-1]+1,i=2..N-1), seq(x[i,1]=p*x[i+1,1]+(1-p)*x[i-1,1]+1,i=2..M-1),seq(x[i,N]=p*x[i+1,N]+(1-p)*x[i-1,N]+1,i=2..M-1), seq(seq(x[i,j]=x[i+1,j+1]*p*q+x[i+1,j-1]*p*(1-q)+x[i-1,j+1]*(1-p)*q+x[i-1,j-1]*(1-p)*(1-q)+1,i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
sol:
end:


#EDNGR2(M,N,p,q,g): computes and plots the expected duration of a game of 2D Gambler's Ruin Type 2 from each potential starting position
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# a numeric or symbolic probability g playing the vertical game each round 
#outputs: a matrix with each entry indicating the expected duration of the game from that starting position, and if p and q are numeric, a density plot of said matrix. 

EDNGR2:=proc(M,N,p,q,g) local eq,var,sol,x,i,j,f:
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,1]=0, x[1,N]=0, x[M,1]=0, x[M,N]=0, seq(x[1,i]=q*x[1,i+1]+(1-q)*x[1,i-1]+1,i=2..N-1), seq(x[M,i]=q*x[M,i+1]+(1-q)*x[M,i-1]+1,i=2..N-1), seq(x[i,1]=p*x[i+1,1]+(1-p)*x[i-1,1]+1,i=2..M-1),seq(x[i,N]=p*x[i+1,N]+(1-p)*x[i-1,N]+1,i=2..M-1),seq(seq(x[i,j]=x[i+1,j]*p*g+x[i-1,j]*g*(1-p)+x[i,j+1]*(1-g)*q+x[i,j-1]*(1-g)*(1-q)+1,i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=[seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]:
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->sol[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
sol:
end:


#TestEDdiff(M,N,p,q): computes and plots the difference in expected duration between a game of 2D Gambler's Ruin Type 1 and Type 2 from each potential starting position
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
#outputs: a matrix with each entry indicating the difference in expected durations of the games from that starting position, and if p and q are numeric, a density plot of said matrix. 

TestEDdiff:=proc(M,N,p,q) local A,B,f:
A:=EDNGR1(M,N,p,q):
B:=EDNGR2(M,N,p,q,g):
A:=B-A:
if type(p, numeric) and type(q, numeric) then
f:=(x,y)->A[ceil(y),ceil(x)]:
print(densityplot('f(x,y)',x=0.01..N,y=0.01..M,style=patchnogrid)):
fi:
A:
end:


#PGFDNGR1(M,N,p,q,t): computes the probability generating functions for the duration of a 2D Gambler's Ruin Game Type 1 from each starting position
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# a variable symbol t
#outputs: a matrix with each entry as the probability generating function for duration (the rational function in t whose coefficient of t^i is the probability that the game will end exactly after i round) from that starting position. 

PGFDNGR1:=proc(M,N,p,q,t) local eq,var,sol,x,i,j,f:
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,1]=1, x[1,N]=1, x[M,1]=1, x[M,N]=1, seq(x[1,i]=t*(q*x[1,i+1]+(1-q)*x[1,i-1]),i=2..N-1), seq(x[M,i]=t*(q*x[M,i+1]+(1-q)*x[M,i-1]),i=2..N-1), seq(x[i,1]=t*(p*x[i+1,1]+(1-p)*x[i-1,1]),i=2..M-1),seq(x[i,N]=t*(p*x[i+1,N]+(1-p)*x[i-1,N]),i=2..M-1), seq(seq(x[i,j]=t*(x[i+1,j+1]*p*q+x[i+1,j-1]*p*(1-q)+x[i-1,j+1]*(1-p)*q+x[i-1,j-1]*(1-p)*(1-q)),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=normal([seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]):
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
end:


#PGFDNGR2(M,N,p,q,g,t): computes the probability generating functions for the duration of a 2D Gambler's Ruin Game Type 2 from each starting position
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# a numeric or symbolic probability g playing the vertical game each round 
# a variable symbol t
#outputs: a matrix with each entry as the probability generating function for duration (the rational function in t whose coefficient of t^i is the probability that the game will end exactly after i round) from that starting position. 

PGFDNGR2:=proc(M,N,p,q,g,t) local eq,var,sol,x,i,j,f:
var:={seq(seq(x[i,j],i=1..M),j=1..N)}:
eq:={x[1,1]=1, x[1,N]=1, x[M,1]=1, x[M,N]=1, seq(x[1,i]=t*(q*x[1,i+1]+(1-q)*x[1,i-1]),i=2..N-1), seq(x[M,i]=t*(q*x[M,i+1]+(1-q)*x[M,i-1]),i=2..N-1), seq(x[i,1]=t*(p*x[i+1,1]+(1-p)*x[i-1,1]),i=2..M-1),seq(x[i,N]=t*(p*x[i+1,N]+(1-p)*x[i-1,N]),i=2..M-1), seq(seq(x[i,j]=t*(x[i+1,j]*p*g+x[i-1,j]*g*(1-p)+x[i,j+1]*(1-g)*q+x[i,j-1]*(1-g)*(1-q)),i=2..M-1),j=2..N-1)}:
var:=solve(eq,var):
sol:=normal([seq([seq(subs(var,x[i,j]),i=1..M)],j=1..N)]):
sol:=convert(sol, Matrix):
sol:=Transpose(sol):
end:


#TestPGdiff(M,N,p,q): computes and plots the difference in expected duration between a game of 2D Gambler's Ruin Type 1 and Type 2 from each potential starting position
#inputs: an integer M>0 for the win condition of $M-1 in the vertical game
# an integer N>0 for the win condition of $N-1 in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
#outputs: a matrix with each entry indicating the difference in probability generating function of durations of the games from that starting position

TestPGdiff:=proc(M,N,p,q,t) local A,B,f:
A:=PGFDNGR1(M,N,p,q,t):
B:=PGFDNGR2(M,N,p,q,g,t):
simplify(B-A):
end:


##############################
# SYMBOLIC 2D GAMBLER'S RUIN #
##############################

#ProbWinning2D(M,N,S,p,q,node): computes the probability of ending a 2D Gambler's Ruin game in each node with all symbolic location input, and numeric or symbolic probability. 
#inputs: a symbol M for the win condition of $M in the vertical game
# a symbol N for the win condition of $N in the horizontal game 
# a list of symbols S=[m,n] for starting with $m in the vertical game and $n in the horizontal game 
# a numeric or symbolic probability p for gaining $1 in the vertical game 
# a numeric or symbolic probability q for gaining $1 in the horizontal game 
# an integer 1-5 indicating the output is for: 1: the node [1,1]; 2: the node [M,1]; 3: the node [1,N]; 4: the node [M,N]; or 5: all four nodes.
#outputs: a matrix with each entry indicating the probability of ending each game in the given node from that starting position, and if p and q are numeric, a density plot of said matrix. 

ProbWinning2D:=proc(M,N,S,p,q,node) local H,V:
V:=ProbWinningCF(M,S[1],p):
H:=ProbWinningCF(N,S[2],q):
if node=1 or node=5 then
print(simplify((1-V)*(1-H))):
fi:
if node=2 or node=5 then
print(simplify((V)*(1-H))):
fi:
if node=3 or node=5 then
print(simplify((1-V)*(H))):
fi:
if node=4 or node=5 then
print(simplify((V)*(H))):
fi:
end:



######
#Ring Graph and Torus
######

#inputs a matrix A of numbers (the transition matrix) given as a list of lists, such that A[i][j] is the probability that if today are at i, 
#then tomorrow you would be at j (of course it has to be a stochastic matrix, each row should have non-negative entries, and they should add-up to 1). 
#The walk is on that graph with vertices {1, ..., N=nops(A)}. StartingPoint is location from 1 to N, 
#and SetOfLosingSites and SetOfWinningSites are disjoint subsets of {1,..,N}.
#The rows correpsponding to members of SetOfLosingSites and SetOfWinningSites should be unit vectors that definitely point to the same place. 
#The output should be a vector of length N that tells you the probability of winning at each starting point. 
#Of course the locations occupied by SetOfLosingSites is automatically 0 and the locations occupied by SetOfWinningSites is automatically 1 

ProbWinningGraph:=proc(A,SetOfLosingSites,SetOfWinningSites) local x,n,i,j,k,l,freex,A0:
n:=nops(A):
freex:={seq(i,i=1..n)} minus {op(SetOfLosingSites)} minus {op(SetOfWinningSites)}:
A0:=solve({seq(x[j]=0,j in SetOfLosingSites),seq(x[k]=1,k in SetOfWinningSites),seq(add(A[l][m]*x[m],m=1..n)=x[l],l in freex)},{seq(x[i],i=1..n)}):
end:

ExpDurationGraph:=proc(A,SetOfLosingSites,SetOfWinningSites) local x,n,i,j,k,l,freex,A0:
n:=nops(A):
freex:={seq(i,i=1..n)} minus {op(SetOfLosingSites)} minus {op(SetOfWinningSites)}:
A0:=solve({seq(x[j]=0,j in SetOfLosingSites),seq(x[k]=0,k in SetOfWinningSites),seq(add(A[l][m]*x[m],m=1..n)+1=x[l],l in freex)},{seq(x[i],i=1..n)}):
A0:
end:

PGFGraph:=proc(A,SetOfLosingSites,SetOfWinningSites,t)  local x,n,i,j,k,l,freex,A0:
n:=nops(A):
freex:={seq(i,i=1..n)} minus {op(SetOfLosingSites)} minus {op(SetOfWinningSites)}:
A0:=solve({seq(x[j]=1,j in SetOfLosingSites),seq(x[k]=1,k in SetOfWinningSites),seq(add(A[l][m]*x[m],m=1..n)*t=x[l],l in freex)},{seq(x[i],i=1..n)}):
end:


LineGraph:=proc(N,p) local i,j,A:
A:=[seq([seq(0,i=1..N+1)],j=1..N+1)]:
A[1][1]:=1:
A[N+1][N+1]:=1:
for i from 2 to N do
A[i][i-1]:=1-p:
A[i][i+1]:=p:
A:
od:
end:




#ExpDurationR(N,n,p1,p2,p3,p4): inputs SYMBOLS n, N and prob. p1,p2,p3,p4 (either numeric or symbolic),
#outputs the CLOSED-FORM EXPRESSION for the expectation of duration in the 
#game on Zn(start at n, aim at 0,+2 with prob p1 ,+1 with p2, -1 with p3. -2 with p4)
ExpDurationR:=proc(N,n,p1,p2,p3,p4) local x,F,A,B,C,A0:
F:=rsolve({x(n)=p4*x(n-2)+p3*x(n-1)+p2*x(n+1)+p1*x(n+2)+1,x(0)=0, x(1)=A, x(2)=B, x(-1)=C},x(n)):
#Now we use the fact that x(N)=0, X(N+1)=x(1)=A, X(N+2)=x(2)=B, X(N-1)=x(-1)=C, and solve for A,B,C
A0:=solve({subs(n=N,F)=0,subs(n=N+1,F)=A,subs(n=N-1,F)=C},{A,B,C}):
#Now we plug that expression for A0 into F
normal(simplify(subs(A0,F))):
end:
#Experiment with this.
#F := [seq(evalf[200](subs({N = 40, n = i}, ExpDurationR(N, n, 19/40, 1/40, 1/40, 19/40))), i = 0 .. 40)]



#Torus(m,n,p,q,r,s): the transition matrix of a torus.
Torus:=proc(m,n,p,q,r,s) local i,j,k,l,A:
A:=Array(1..m*n,1..m*n):

for k from 1 to m-1 do 
 for l from 1 to n-1 do
  A[(k-1)*n+l,(k-1)*n+l+1]:=p:
  A[(k-1)*n+l,(k)*n+l]:=q:
  A[(k-1)*n+l+1,(k-1)*n+l]:=r:
  A[(k)*n+l,(k-1)*n+l]:=s:
 od:
od:
for k from 1 to m-1 do 
  A[(k)*n,(k-1)*n+1]:=p:
  A[(k)*n,(k+1)*n]:=q:
  A[(k-1)*n+1,(k)*n]:=r:
  A[(k+1)*n,(k)*n]:=s:
od:

for l from 1 to n-1 do
A[(m-1)*n+l,(m-1)*n+l+1]:=p:
A[(m-1)*n+l,l]:=q:
A[(m-1)*n+l+1,(m-1)*n+l]:=r:
A[l,(m-1)*n+l]:=s:
od:
 
A[m*n,(m-1)*n+1]:=p:
A[m*n,n]:=q:
A[(m-1)*n+1,m*n]:=r:
A[n,m*n]:=s:
convert(A,listlist):
end:


#On a torus lattice, with C_M vertically times C_N horizontally, for each step it has prob. p to go right, prob. q to go up, prob. r to go left, prob. s to go down.
#The following functions give the prob. of winning, expectation of duration and prob. generating function of duration.
ProbWinningTR:=proc(M,N,p,q,r,s,SetOfLosingSites,SetOfWinningSites) local i,j,A,lose,win:
lose:=[seq((i[1]-1)*N+i[2],i in SetOfLosingSites)]:
win:=[seq((i[1]-1)*N+i[2],i in SetOfWinningSites)]:
A:=ProbWinningGraph(Torus(M,N,p,q,r,s),lose,win):
Matrix([seq([seq(rhs(A[(j-1)*N+i]),j=1..M)],i=1..N)]):
end:

ExpDurationTR:=proc(M,N,p,q,r,s,SetOfLosingSites,SetOfWinningSites) local i,j,A,lose,win:
lose:=[seq((i[1]-1)*N+i[2],i in SetOfLosingSites)]:
win:=[seq((i[1]-1)*N+i[2],i in SetOfWinningSites)]:
A:=ExpDurationGraph(Torus(M,N,p,q,r,s),lose,win):
Matrix([seq([seq(rhs(A[(j-1)*N+i]),j=1..M)],i=1..N)]):
end:

PGFTR:=proc(M,N,p,q,r,s,SetOfLosingSites,SetOfWinningSites,t) local i,j,A,lose,win:
lose:=[seq((i[1]-1)*N+i[2],i in SetOfLosingSites)]:
win:=[seq((i[1]-1)*N+i[2],i in SetOfWinningSites)]:
A:=PGFGraph(Torus(M,N,p,q,r,s),lose,win,t):
Matrix([seq([seq(rhs(A[(j-1)*N+i]),j=1..M)],i=1..N)]):
end:

#e.g.ExpDurationTR(20, 20, 1/4, 1/4, 1/4, 1/4, [[1, 1]], [[10, 10]])
#ProbWinningTR(20, 20, 1/4, 1/4, 1/4, 1/4, [[1, 1]], [[10, 10]])
#listplot3d can plot the data above.


######
#High Dimensional Game
######


#GWinProbOne(N1,N2,N3,p,q,r): inputs positive integers N1, N2, and N3 and a probabilities p,q and r (either numeric or symbolic),
#outputs a list of length (N1+1)*(N2+1)*(N3+1) where each entry is the probability of exiting a winner in
#Generalized Gambler's Ruin Game where a player starts with a capital of (i,j,k) dollars
#At each round, a player randomly chooses a subgame to play, earning 1 dollar with probability p_i
#and losing 1 dollar with probability (1-p_i) 
#A player wins if the player wins at least one subgame.

GWinProbOne:=proc(N1,N2,N3,p,q,r) local x,i,j,k,var,var1,eq:
var:=[seq(seq(seq(x[i,j,k],i=0..N1),j=0..N2),k=0..N3)]:
eq:={seq(seq(seq(x[i,j,k]=(1/3)*(p*x[i+1,j,k]+(1-p)*x[i-1,j,k])
                         +(1/3)*(q*x[i,j+1,k]+(1-q)*x[i,j-1,k])
                         +(1/3)*(r*x[i,j,k+1]+(1-r)*x[i,j,k-1]),i=1..N1-1),j=1..N2-1),k=1..N3-1),
     x[0,0,0]=0,
     x[N1,0,0]=1,
     x[0,N2,0]=1,
     x[0,0,N3]=1,
     x[N1,N2,0]=1,
     x[0,N2,N3]=1,
     x[N1,0,N3]=1,
     x[N1,N2,N3]=1}:
for i from 1 to N1-1 do
  eq:=eq union {x[i,0,0]=p*x[i+1,0,0]+(1-p)*x[i-1,0,0]}:
  eq:=eq union {x[i,N2,0]=1}:
  eq:=eq union {x[i,0,N3]=1}:
  eq:=eq union {x[i,N2,N3]=1}:
od:
for j from 1 to N2-1 do
  eq:=eq union {x[0,j,0]=q*x[0,j+1,0]+(1-q)*x[0,j-1,0]}:
  eq:=eq union {x[N1,j,0]=1}:
  eq:=eq union {x[0,j,N3]=1}:
  eq:=eq union {x[N1,j,N3]=1}:
od:
for k from 1 to N3-1 do
  eq:=eq union {x[0,0,k]=r*x[0,0,k+1]+(1-r)*x[0,0,k-1]}:
  eq:=eq union {x[N1,0,k]=1}:
  eq:=eq union {x[0,N2,k]=1}:
  eq:=eq union {x[N1,N2,k]=1}:
od:
eq:=eq union {seq(seq(x[i,j,0]=.5*(p*x[i+1,j,0]+(1-p)*x[i-1,j,0]) 
                              +.5*(q*x[i,j+1,0]+(1-q)*x[i,j-1,0]),i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,j,N3]=1,i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,0,k]=.5*(p*x[i+1,0,k]+(1-p)*x[i-1,0,k]) 
                              +.5*(r*x[i,0,k+1]+(1-r)*x[i,0,k-1]),i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[i,N2,k]=1,i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[0,j,k]=.5*(q*x[0,j+1,k]+(1-q)*x[0,j-1,k]) 
                              +.5*(r*x[0,j,k+1]+(1-r)*x[0,j,k-1]),j=1..N2-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[N1,j,k]=1,j=1..N2-1),k=1..N3-1)}:
var:=solve(eq,var):
var[1]:
end:


#GExpDurOne(N1,N2,N3,p,q,r): inputs positive integers N1, N2, and N3 and a probabilities p,q and r (either numeric or symbolic),
#outputs a list of length (N1+1)*(N2+1)*(N3+1) where each entry is the expected duration of a
#Generalized Gambler's Ruin Game where a player starts with a capital of (i,j,k) dollars
#At each round, a player randomly chooses a subgame to play, earning 1 dollar with probability p_i
#and losing 1 dollar with probability (1-p_i) 
#A player wins if the player wins at least one subgame.

GExpDurOne:=proc(N1,N2,N3,p,q,r) local x,i,j,k,var,eq:
var:=[seq(seq(seq(x[i,j,k],i=0..N1),j=0..N2),k=0..N3)]:
eq:={seq(seq(seq(x[i,j,k]=(1/3)*(p*x[i+1,j,k]+(1-p)*x[i-1,j,k])
                         +(1/3)*(q*x[i,j+1,k]+(1-q)*x[i,j-1,k])
                         +(1/3)*(r*x[i,j,k+1]+(1-r)*x[i,j,k-1])+1,i=1..N1-1),j=1..N2-1),k=1..N3-1),
     x[0,0,0]=0,
     x[N1,0,0]=0,
     x[0,N2,0]=0,
     x[0,0,N3]=0,
     x[N1,N2,0]=0,
     x[0,N2,N3]=0,
     x[N1,0,N3]=0,
     x[N1,N2,N3]=0}:
for i from 1 to N1-1 do
  eq:=eq union {x[i,0,0]=p*x[i+1,0,0]+(1-p)*x[i-1,0,0]+1}:
  eq:=eq union {x[i,N2,0]=0}:
  eq:=eq union {x[i,0,N3]=0}:
  eq:=eq union {x[i,N2,N3]=0}:
od:
for j from 1 to N2-1 do
  eq:=eq union {x[0,j,0]=q*x[0,j+1,0]+(1-q)*x[0,j-1,0]+1}:
  eq:=eq union {x[N1,j,0]=0}:
  eq:=eq union {x[0,j,N3]=0}:
  eq:=eq union {x[N1,j,N3]=0}:
od:
for k from 1 to N3-1 do
  eq:=eq union {x[0,0,k]=r*x[0,0,k+1]+(1-r)*x[0,0,k-1]+1}:
  eq:=eq union {x[N1,0,k]=0}:
  eq:=eq union {x[0,N2,k]=0}:
  eq:=eq union {x[N1,N2,k]=0}:
od:
eq:=eq union {seq(seq(x[i,j,0]=.5*(p*x[i+1,j,0]+(1-p)*x[i-1,j,0]) 
                              +.5*(q*x[i,j+1,0]+(1-q)*x[i,j-1,0])+1,i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,j,N3]=0,i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,0,k]=.5*(p*x[i+1,0,k]+(1-p)*x[i-1,0,k]) 
                              +.5*(r*x[i,0,k+1]+(1-r)*x[i,0,k-1])+1,i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[i,N2,k]=0,i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[0,j,k]=.5*(q*x[0,j+1,k]+(1-q)*x[0,j-1,k]) 
                              +.5*(r*x[0,j,k+1]+(1-r)*x[0,j,k-1])+1,j=1..N2-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[N1,j,k]=0,j=1..N2-1),k=1..N3-1)}:
var:=solve(eq,var):
var[1]:
end:


#GPGFExpDurOne(N1,N2,N3,p,q,r): inputs positive integers N1, N2, and N3 and a probabilities p,q and r (either numeric or symbolic),
#outputs a list of length (N1+1)*(N2+1)*(N3+1) where each entry is the probability generating function
#of the expected duration of a Generalized Gambler's Ruin Game
#where a player starts with a capital of (i,j,k) dollars
#At each round, a player randomly chooses a subgame to play, earning 1 dollar with probability p_i
#and losing 1 dollar with probability (1-p_i) 
#A player wins if the player wins at least one subgame.

GPGFExpDurOne:=proc(N1,N2,N3,p,q,r,t) local x,i,j,k,var,var1,eq:
var:=[seq(seq(seq(x[i,j,k],i=0..N1),j=0..N2),k=0..N3)]:
eq:={seq(seq(seq(x[i,j,k]=t*((1/3)*(p*x[i+1,j,k]+(1-p)*x[i-1,j,k])
                         +(1/3)*(q*x[i,j+1,k]+(1-q)*x[i,j-1,k])
                         +(1/3)*(r*x[i,j,k+1]+(1-r)*x[i,j,k-1])),i=1..N1-1),j=1..N2-1),k=1..N3-1),
     x[0,0,0]=1,
     x[N1,0,0]=1,
     x[0,N2,0]=1,
     x[0,0,N3]=1,
     x[N1,N2,0]=1,
     x[0,N2,N3]=1,
     x[N1,0,N3]=1,
     x[N1,N2,N3]=1}:
for i from 1 to N1-1 do
  eq:=eq union {x[i,0,0]=t*(p*x[i+1,0,0]+(1-p)*x[i-1,0,0])}:
  eq:=eq union {x[i,N2,0]=1}:
  eq:=eq union {x[i,0,N3]=1}:
  eq:=eq union {x[i,N2,N3]=1}:
od:
for j from 1 to N2-1 do
  eq:=eq union {x[0,j,0]=t*(q*x[0,j+1,0]+(1-q)*x[0,j-1,0])}:
  eq:=eq union {x[N1,j,0]=1}:
  eq:=eq union {x[0,j,N3]=1}:
  eq:=eq union {x[N1,j,N3]=1}:
od:
for k from 1 to N3-1 do
  eq:=eq union {x[0,0,k]=t*(r*x[0,0,k+1]+(1-r)*x[0,0,k-1])}:
  eq:=eq union {x[N1,0,k]=1}:
  eq:=eq union {x[0,N2,k]=1}:
  eq:=eq union {x[N1,N2,k]=1}:
od:
eq:=eq union {seq(seq(x[i,j,0]=t*(.5*(p*x[i+1,j,0]+(1-p)*x[i-1,j,0]) 
                              +.5*(q*x[i,j+1,0]+(1-q)*x[i,j-1,0])),i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,j,N3]=1,i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,0,k]=t*(.5*(p*x[i+1,0,k]+(1-p)*x[i-1,0,k]) 
                              +.5*(r*x[i,0,k+1]+(1-r)*x[i,0,k-1])),i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[i,N2,k]=1,i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[0,j,k]=t*(.5*(q*x[0,j+1,k]+(1-q)*x[0,j-1,k]) 
                              +.5*(r*x[0,j,k+1]+(1-r)*x[0,j,k-1])),j=1..N2-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[N1,j,k]=1,j=1..N2-1),k=1..N3-1)}:
var:=solve(eq,var):
var[1]:
end:



#GWinProbAll(N1,N2,N3,p,q,r): inputs positive integers N1, N2, and N3 and a probabilities p,q and r (either numeric or symbolic),
#outputs a list of length (N1+1)*(N2+1)*(N3+1) where each entry is the probability of exiting a winner in
#Generalized Gambler's Ruin Game where a player starts with a capital of (i,j,k) dollars
#At each round, a player randomly chooses a subgame to play, earning 1 dollar with probability p_i
#and losing 1 dollar with probability (1-p_i) 
#A player wins only if the player wins every subgame.

GWinProbAll:=proc(N1,N2,N3,p,q,r) local x,i,j,k,var,var1,eq:
var:=[seq(seq(seq(x[i,j,k],i=0..N1),j=0..N2),k=0..N3)]:
eq:={seq(seq(seq(x[i,j,k]=(1/3)*(p*x[i+1,j,k]+(1-p)*x[i-1,j,k])
                         +(1/3)*(q*x[i,j+1,k]+(1-q)*x[i,j-1,k])
                         +(1/3)*(r*x[i,j,k+1]+(1-r)*x[i,j,k-1]),i=1..N1-1),j=1..N2-1),k=1..N3-1),
     x[0,0,0]=0,
     x[N1,0,0]=0,
     x[0,N2,0]=0,
     x[0,0,N3]=0,
     x[N1,N2,0]=0,
     x[0,N2,N3]=0,
     x[N1,0,N3]=0,
     x[N1,N2,N3]=1}:
for i from 1 to N1-1 do
  eq:=eq union {x[i,0,0]=p*x[i+1,0,0]+(1-p)*x[i-1,0,0]}:
  eq:=eq union {x[i,N2,0]=p*x[i+1,N2,0]+(1-p)*x[i-1,N2,0]}:
  eq:=eq union {x[i,0,N3]=p*x[i+1,0,N3]+(1-p)*x[i-1,0,N3]}:
  eq:=eq union {x[i,N2,N3]=p*x[i+1,N2,N3]+(1-p)*x[i-1,N2,N3]}:
od:
for j from 1 to N2-1 do
  eq:=eq union {x[0,j,0]=q*x[0,j+1,0]+(1-q)*x[0,j-1,0]}:
  eq:=eq union {x[N1,j,0]=q*x[N1,j+1,0]+(1-q)*x[N1,j-1,0]}:
  eq:=eq union {x[0,j,N3]=q*x[0,j+1,N3]+(1-q)*x[0,j-1,N3]}:
  eq:=eq union {x[N1,j,N3]=q*x[N1,j+1,N3]+(1-q)*x[N1,j-1,N3]}:
od:
for k from 1 to N3-1 do
  eq:=eq union {x[0,0,k]=r*x[0,0,k+1]+(1-r)*x[0,0,k-1]}:
  eq:=eq union {x[N1,0,k]=r*x[N1,0,k+1]+(1-r)*x[N1,0,k-1]}:
  eq:=eq union {x[0,N2,k]=r*x[0,N2,k+1]+(1-r)*x[0,N2,k-1]}:
  eq:=eq union {x[N1,N2,k]=r*x[N1,N2,k+1]+(1-r)*x[N1,N2,k-1]}:
od:
eq:=eq union {seq(seq(x[i,j,0]=.5*(p*x[i+1,j,0]+(1-p)*x[i-1,j,0]) 
                              +.5*(q*x[i,j+1,0]+(1-q)*x[i,j-1,0]),i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,j,N3]=.5*(p*x[i+1,j,N3]+(1-p)*x[i-1,j,N3]) 
                              +.5*(q*x[i,j+1,N3]+(1-q)*x[i,j-1,N3]),i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,0,k]=.5*(p*x[i+1,0,k]+(1-p)*x[i-1,0,k]) 
                              +.5*(r*x[i,0,k+1]+(1-r)*x[i,0,k-1]),i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[i,N2,k]=.5*(p*x[i+1,N2,k]+(1-p)*x[i-1,N2,k]) 
                              +.5*(r*x[i,N2,k+1]+(1-r)*x[i,N2,k-1]),i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[0,j,k]=.5*(q*x[0,j+1,k]+(1-q)*x[0,j-1,k]) 
                              +.5*(r*x[0,j,k+1]+(1-r)*x[0,j,k-1]),j=1..N2-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[N1,j,k]=.5*(q*x[N1,j+1,k]+(1-q)*x[N1,j-1,k]) 
                              +.5*(r*x[N1,j,k+1]+(1-r)*x[N1,j,k-1]),j=1..N2-1),k=1..N3-1)}:
var:=solve(eq,var):
var[1]:
end:

#GExpDurAll(N1,N2,N3,p,q,r): inputs positive integers N1, N2, and N3 and a probabilities p,q and r (either numeric or symbolic),
#outputs a list of length (N1+1)*(N2+1)*(N3+1) where each entry is the expected duration of a
#Generalized Gambler's Ruin Game where a player starts with a capital of (i,j,k) dollars
#At each round, a player randomly chooses a subgame to play, earning 1 dollar with probability p_i
#and losing 1 dollar with probability (1-p_i) 
#A player wins only if the player wins every subgame.

GExpDurAll:=proc(N1,N2,N3,p,q,r) local x,i,j,k,var,eq:
var:=[seq(seq(seq(x[i,j,k],i=0..N1),j=0..N2),k=0..N3)]:
eq:={seq(seq(seq(x[i,j,k]=(1/3)*(p*x[i+1,j,k]+(1-p)*x[i-1,j,k])
                         +(1/3)*(q*x[i,j+1,k]+(1-q)*x[i,j-1,k])
                         +(1/3)*(r*x[i,j,k+1]+(1-r)*x[i,j,k-1])+1,i=1..N1-1),j=1..N2-1),k=1..N3-1),
     x[0,0,0]=0,
     x[N1,0,0]=0,
     x[0,N2,0]=0,
     x[0,0,N3]=0,
     x[N1,N2,0]=0,
     x[0,N2,N3]=0,
     x[N1,0,N3]=0,
     x[N1,N2,N3]=0}:
for i from 1 to N1-1 do
  eq:=eq union {x[i,0,0]=p*x[i+1,0,0]+(1-p)*x[i-1,0,0]+1}:
  eq:=eq union {x[i,N2,0]=p*x[i+1,N2,0]+(1-p)*x[i-1,N2,0]+1}:
  eq:=eq union {x[i,0,N3]=p*x[i+1,0,N3]+(1-p)*x[i-1,0,N3]+1}:
  eq:=eq union {x[i,N2,N3]=p*x[i+1,N2,N3]+(1-p)*x[i-1,N2,N3]+1}:
od:
for j from 1 to N2-1 do
  eq:=eq union {x[0,j,0]=q*x[0,j+1,0]+(1-q)*x[0,j-1,0]+1}:
  eq:=eq union {x[N1,j,0]=q*x[N1,j+1,0]+(1-q)*x[N1,j-1,0]+1}:
  eq:=eq union {x[0,j,N3]=q*x[0,j+1,N3]+(1-q)*x[0,j-1,N3]+1}:
  eq:=eq union {x[N1,j,N3]=q*x[N1,j+1,N3]+(1-q)*x[N1,j-1,N3]+1}:
od:
for k from 1 to N3-1 do
  eq:=eq union {x[0,0,k]=r*x[0,0,k+1]+(1-r)*x[0,0,k-1]+1}:
  eq:=eq union {x[N1,0,k]=r*x[N1,0,k+1]+(1-r)*x[N1,0,k-1]+1}:
  eq:=eq union {x[0,N2,k]=r*x[0,N2,k+1]+(1-r)*x[0,N2,k-1]+1}:
  eq:=eq union {x[N1,N2,k]=r*x[N1,N2,k+1]+(1-r)*x[N1,N2,k-1]+1}:
od:
eq:=eq union {seq(seq(x[i,j,0]=.5*(p*x[i+1,j,0]+(1-p)*x[i-1,j,0]) 
                              +.5*(q*x[i,j+1,0]+(1-q)*x[i,j-1,0])+1,i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,j,N3]=.5*(p*x[i+1,j,N3]+(1-p)*x[i-1,j,N3]) 
                              +.5*(q*x[i,j+1,N3]+(1-q)*x[i,j-1,N3])+1,i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,0,k]=.5*(p*x[i+1,0,k]+(1-p)*x[i-1,0,k]) 
                              +.5*(r*x[i,0,k+1]+(1-r)*x[i,0,k-1])+1,i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[i,N2,k]=.5*(p*x[i+1,N2,k]+(1-p)*x[i-1,N2,k]) 
                              +.5*(r*x[i,N2,k+1]+(1-r)*x[i,N2,k-1])+1,i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[0,j,k]=.5*(q*x[0,j+1,k]+(1-q)*x[0,j-1,k]) 
                              +.5*(r*x[0,j,k+1]+(1-r)*x[0,j,k-1])+1,j=1..N2-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[N1,j,k]=.5*(q*x[N1,j+1,k]+(1-q)*x[N1,j-1,k]) 
                              +.5*(r*x[N1,j,k+1]+(1-r)*x[N1,j,k-1])+1,j=1..N2-1),k=1..N3-1)}:
var:=solve(eq,var):
var[1]:
end:

#GPGFExpDurAll(N1,N2,N3,p,q,r): inputs positive integers N1, N2, and N3 and a probabilities p,q and r (either numeric or symbolic),
#outputs a list of length (N1+1)*(N2+1)*(N3+1) where each entry is the probability generating function
#of the expected duration of a Generalized Gambler's Ruin Game 
#where a player starts with a capital of (i,j,k) dollars
#At each round, a player randomly chooses a subgame to play, earning 1 dollar with probability p_i
#and losing 1 dollar with probability (1-p_i) 
#A player wins only if the player wins every subgame.

GPGFExpDurAll:=proc(N1,N2,N3,p,q,r,t) local x,i,j,k,var,var1,eq:
var:=[seq(seq(seq(x[i,j,k],i=0..N1),j=0..N2),k=0..N3)]:
eq:={seq(seq(seq(x[i,j,k]=t*((1/3)*(p*x[i+1,j,k]+(1-p)*x[i-1,j,k])
                         +(1/3)*(q*x[i,j+1,k]+(1-q)*x[i,j-1,k])
                         +(1/3)*(r*x[i,j,k+1]+(1-r)*x[i,j,k-1])),i=1..N1-1),j=1..N2-1),k=1..N3-1),
     x[0,0,0]=1,
     x[N1,0,0]=1,
     x[0,N2,0]=1,
     x[0,0,N3]=1,
     x[N1,N2,0]=1,
     x[0,N2,N3]=1,
     x[N1,0,N3]=1,
     x[N1,N2,N3]=1}:
for i from 1 to N1-1 do
  eq:=eq union {x[i,0,0]=t*(p*x[i+1,0,0]+(1-p)*x[i-1,0,0])}:
  eq:=eq union {x[i,N2,0]=t*(p*x[i+1,N2,0]+(1-p)*x[i-1,N2,0])}:
  eq:=eq union {x[i,0,N3]=t*(p*x[i+1,0,N3]+(1-p)*x[i-1,0,N3])}:
  eq:=eq union {x[i,N2,N3]=t*(p*x[i+1,N2,N3]+(1-p)*x[i-1,N2,N3])}:
od:
for j from 1 to N2-1 do
  eq:=eq union {x[0,j,0]=t*(q*x[0,j+1,0]+(1-q)*x[0,j-1,0])}:
  eq:=eq union {x[N1,j,0]=t*(q*x[N1,j+1,0]+(1-q)*x[N1,j-1,0])}:
  eq:=eq union {x[0,j,N3]=t*(q*x[0,j+1,N3]+(1-q)*x[0,j-1,N3])}:
  eq:=eq union {x[N1,j,N3]=t*(q*x[N1,j+1,N3]+(1-q)*x[N1,j-1,N3])}:
od:
for k from 1 to N3-1 do
  eq:=eq union {x[0,0,k]=t*(r*x[0,0,k+1]+(1-r)*x[0,0,k-1])}:
  eq:=eq union {x[N1,0,k]=t*(r*x[N1,0,k+1]+(1-r)*x[N1,0,k-1])}:
  eq:=eq union {x[0,N2,k]=t*(r*x[0,N2,k+1]+(1-r)*x[0,N2,k-1])}:
  eq:=eq union {x[N1,N2,k]=t*(r*x[N1,N2,k+1]+(1-r)*x[N1,N2,k-1])}:
od:
eq:=eq union {seq(seq(x[i,j,0]=t*(.5*(p*x[i+1,j,0]+(1-p)*x[i-1,j,0]) 
                              +.5*(q*x[i,j+1,0]+(1-q)*x[i,j-1,0])),i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,j,N3]=t*(.5*(p*x[i+1,j,N3]+(1-p)*x[i-1,j,N3]) 
                              +.5*(q*x[i,j+1,N3]+(1-q)*x[i,j-1,N3])),i=1..N1-1),j=1..N2-1)}:
eq:=eq union {seq(seq(x[i,0,k]=t*(.5*(p*x[i+1,0,k]+(1-p)*x[i-1,0,k]) 
                              +.5*(r*x[i,0,k+1]+(1-r)*x[i,0,k-1])),i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[i,N2,k]=t*(.5*(p*x[i+1,N2,k]+(1-p)*x[i-1,N2,k]) 
                              +.5*(r*x[i,N2,k+1]+(1-r)*x[i,N2,k-1])),i=1..N1-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[0,j,k]=t*(.5*(q*x[0,j+1,k]+(1-q)*x[0,j-1,k]) 
                              +.5*(r*x[0,j,k+1]+(1-r)*x[0,j,k-1])),j=1..N2-1),k=1..N3-1)}:
eq:=eq union {seq(seq(x[N1,j,k]=t*(.5*(q*x[N1,j+1,k]+(1-q)*x[N1,j-1,k]) 
                              +.5*(r*x[N1,j,k+1]+(1-r)*x[N1,j,k-1])),j=1..N2-1),k=1..N3-1)}:
var:=solve(eq,var):
var[1]:
end: