######################################################################
##UmbralMarkov.txt: Save this file as  UmbralMarkov.txt              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  UmbralMarkov.txt<Enter>                            #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: July 2019
 
print(`Created:  July , 2019`):
print(` This is UmbralMarkov.txt `):
print(`It is one of the packages that accompany the article `):
print(`Using Symbolic Computation to analyze some Childish Board Games`):
print(`by Shalosh B. Ekhad  and Doron Zeilberger`):
print(`published in the Personal J. of Ekhad and Zeilberger and arxiv.org`):
print(``):
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 procedures epecific to the Counnt Your Chickens! game type ezraCh();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

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


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):
Digits:=20:




ezraR:=proc()

if args=NULL then
 print(` The random procedures are:  `):
 print(` RandDG, RandSet1, RandWMCe  `):

else
ezra(args):
fi:

end:

ezraCh:=proc()

if args=NULL then
 print(` The Count Your Chickens! procedures are:  CCb, CCMC, FindNext,  NeiState `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  CCb, FindNext,  HarogKatan, HarogGadol, NeiState, Pashet, PreUmb, PreUmbT, Umb, UmbT `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Alpha, Info, InfoT, WMCSer, WMCSerT `):

 print(` `):





elif nops([args])=1 and op(1,[args])=Alpha then
print(`Alpha(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond)`):
print(`returns a list, of length N, whose `):
print(`(i) First entry is the average`):
print(` (ii): Second entry is the variance `):
print(` for i=3...N, the i-th entry is the so-called alpha-coefficients `):
print(` that is the i-th moment about the mean divided by the `):
print(` variance to the power i/2 (For example, i=3 is the skewness `):
print(` and i=4 is the Kurtosis) `):
print(`If f(1) is not 1, than it first normalizes it by dividing`):
print(`by f(1) (if it is not 0) .`):
print(` For example, try: `):
print(` Alpha(((1+x)/2)^100,x,4); `):

elif nops([args])=1 and op(1,[args])=CCb then
print(`CCb(C,D,P,S,T): The "Count Your Chickens board". It is a triple , the labels , followed by the board, followed by the set of blue squares`):
print(`An emptry square is indicated by 0. It is assumed that the last square (nops(gu[2])+1) has all of them. Try:`):
print(`CCb(C,D,P,S,T);`):



elif nops([args])=1 and op(1,[args])=CCMC then
print(`CCMC(CB,SP): inputs a general Count Your Chickens! Style game, with the probability distribution on the set of animals`):
print(`outputs the weighted absorbing Markov Chain. The format of a weighted absorbing Markov Chain is`):
print(`a list L of length N, say, where there are N states labeled 1, ...,N, and the i-th entry is`):
print(` a set of triples [j,Prob(i->j),gain], and the convention is that the absorbing states are labelled N+1, .., N+r. Try: `):
print(` CCMC([[S,C],[0,0,S,C,0,C,0,S],{3,6}],[1/3,1/3]);`):
print(`CCMC(CCb1(),[(1/6)$5]);`):


elif nops([args])=1 and op(1,[args])=FindNext then
print(`FindNext(CB,i,j): inputs a "Count your Chicken board, and a location i, and an outcome of the spinner  from 1 to 6, outputs`):
print(` the next location of the chicken. Try: `):
print(` FindNext(CCb(C,D,P,S,T),4,4); `):



elif nops([args])=1 and op(1,[args])=HarogGadol then
print(`HarogGadol(f,t,N): replaces each power larger than N by t^N. Try:`):
print(`HarogGadol(a*t+b*t^3+d*t^4,t,2);`):

elif nops([args])=1 and op(1,[args])=HarogKatan then
print(`HarogKatan(f,t,N): replaces each power smaller than N by t^N. Try:`):
print(`HarogKatan(a*t+b*t^3+d*t^4,t,2);`):


elif nops([args])=1 and op(1,[args])=Info then
print(`Info(WMC,K,st): Inputs a weighted Markov Chain with absorbing states WMC (according to our convention)`):
print(`a positive integer K, and an intial state st (from 1 to nops(WMC)), outputs a list of`):
print(` assuming that it terminated after K steps `):
print(` (i) The probabilities of lending in each of the absorbing states `):
print(` (ii) The list of conditional expectation of the gain at each state `):
print(` (iii) The list of conditional expectation of the number of rounds at each state `):
print(` (iv) The probability of finishing after K steps. Here we do NOT use truncation, i.e. there is no limit to your income and you can be in debt`):
print(`Try: `):
print(`W:=RandWMCe(13,4,2,4,0,20); Info(W,50,1);`):
print(`Info(CCMC(CCb1(),[(1/6)$5]),50,1);`):


elif nops([args])=1 and op(1,[args])=InfoT then
print(`InfoT(WMC,K,st,Kat,Gad): Inputs a weighted Markov Chain with absorbing states WMC (according to our convention)`):
print(`a positive integer K, and an intial state st (from 1 to nops(WMC)), outputs a list of`):
print(` assuming that it terminated after K steps `):
print(` (i) The probabilities of lending in each of the absorbing states `):
print(` (ii) The list of conditional expectation of the gain at each state `):
print(` (iii) The list of conditional expectation of the number of rounds at each state `):
print(` (iv) The probability of finishing after K steps.`):
print(` Here we DO use truncation, i.e. minimum fortune is Kat and maximal fortune is Gad`):
print(`Try: `):
print(`W:=RandWMCe(13,4,2,4,0,20); InfoT(W,50,1,0,100);`):
print(` InfoT(CCMC(CCb1(),[(1/6)$5]),50,1,0,40);`):

elif nops([args])=1 and op(1,[args])=NeiState then
print(`NeiState(CB,SP,st): Given a Count Your Chicken Board, CB, a prob. dist. on the set of Animals`):
print(`and a state a, where a is the location 1<=a<=nops(CB[2])+1, and b is`):
print(` outputs the  states that come from it together with the gain of chicks.`):
print(`Try:`):
print(`NeiState([[S,C],[0,0,S,C,0,C,0,S],{3,6}],[1/3,1/3],1);`):
print(`NeiState(CCb1(),[(1/6)$5],1);`):

elif nops([args])=1 and op(1,[args])=Pashet then
print(`Pashet(f,x,t,eps): inputs a polynomial f in the variables x and t and a cut-op eps, discards`):
print(`terms whose absolute value is less than eps. Try:`):
print(`Pashet(x*t+1/10^11*x-1/10^11,x,t,1/10^10);`):


elif nops([args])=1 and op(1,[args])=PreUmb then
print(`PreUmb(f,s,t,WMC): Given a monomial f in s and t applies the pre-umbra operation implied by the Weighted Markov Chain WMC, w/o Truncation`):
print(`Try:`):
print(`PreUmb(s^5*t^11,s,t,CCMC(CCb1(),[(1/6)$5])); `):

elif nops([args])=1 and op(1,[args])=PreUmbT then
print(`PreUmbT(f,s,t,WMC,K,G): Given a monomial f in s and t applies the pre-umbra operation implied by the Weighted Markov Chain WMC, with Truncation`):
print(`where powers less than t^K are replaced by t^K and powers larger than t^G are replaced by t^G.`):
print(`Try:`):
print(`PreUmbT(s^40*t^11,s,t,CCMC(CCb1(),[(1/6)$5]),0,40); `):


elif nops([args])=1 and op(1,[args])=RandDG then
print(`RandDG(N,s,deg1,deg2): a random directed graph with N vertices labelled 1, ..., N that have positive in-degree and s vertices`):
print(`labelled N+1, .., N+s,the output is a list of length N whose i-th entry is the the set of out-neighbors of i`):
print(`with each vertex with a random outdegree between deg1 deg2. Try:`):
print(`RandDG(7,2,2,3);`):

elif nops([args])=1 and op(1,[args])=RandSet1 then
print(`RandSet1(S,k): a random subset with k elements of the set S, try:`):
print(`RandSet1({1,2,3,4},2);`):

elif nops([args])=1 and op(1,[args])=RandWMCe then
print(`RandWMCe(N,s,deg1,deg2,T1,T2): a random Weighted Markov Chain where the weights are integers drawn randomly between T1 and T2`):
print(`and the underlying directed graph with s sinks is RandDG(N,s,deg1,deg2) (q.v.). Try:`):
print(`RandWMCe(7,2,2,3,0,10); `):

elif nops([args])=1 and op(1,[args])=Umb then
print(`Umb(f,s,t,WMC): Given a polynomial f in s and t applies the umbral operation implied by the weighted Markov Chain w/o transaction.`):
print(` Try: `):
print(`Umb(s^5*t^11/2+s^8*t^16/2,s,t, CCMC(CCb1(),[(1/6)$5]));`):

elif nops([args])=1 and op(1,[args])=UmbT then
print(`UmbT(f,s,t,WMC,K,G): Given a monomial f in s and t applies the umbral operation implied by the Weighted Markov Chain WMC, with Truncation`):
print(`where powers less than t^K are replaced by t^K and powers larger than t^G are replaced by t^G.`):
print(`Try:`):
print(`UmbT(s^40*t^11,s,t,CCMC(CCb1(),[(1/6)$5]),0,40); `):


elif nops([args])=1 and op(1,[args])=WMCSer then
print(`WMCSer(WMC,t,X,K,st): The inputs are (i) WMC: Weighted Markov chain according to our convention`):
print(` (ii) variables t and X `):
print(` (iii) a  positive integer K `):
print(` (iv): the starting state, usually 1, in application, but any positive integer beteen 1 and nops(WMC) `):
print(` Let N be the number of non-absorbing states in WMC. `):
print(` By convention that absorbing states are N+1,N+2, .., N+r where is is the number of absorbing states (found from the WMC thatmust adhere to our convention)`):
print(`The output is a list of length r, where the j-th entry is the`):
print(`polynomial of degree K in X where the coeff. of X^i is the`):
print(`prob. generating function, according to gain, if it ended after i rounds and winded up at state N+j, starting at state st.`):
print(` Try:`):
print(` WMCSer(CCMC([[S,C],[0,0,S,C,0,C,0,S],{3,6}],[1/3,1/3]),t,X,30,1);`):
print(`WMCSer(CCMC(CCb1(),[(1/6)$5]),t,X,30,1);`):

elif nops([args])=1 and op(1,[args])=WMCSerT then
print(`WMCSerT(WMC,t,X,K,st, Kat,Gad)`):
print(`Like WMCSer(WMC,t,X,K,st) (q.v.) but at each step the polynomials in t are adjusted so as to replace powers lower than t^Kat by t^Kat`):
print(`and powers higher than t^Gad by t^Gad. `):
print(` Try:`):
print(` WMCSerT(CCMC([[S,C],[0,0,S,C,0,C,0,S],{3,6}],[1/3,1/3]),t,X,30,1,0,8);`):
print(`WMCSerT(CCMC(CCb1(),[(1/6)$5]),t,X,30,1,0,40);`):

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



##########From HISTABRUT

#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:

#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:

##########End From HISTABRUT




#Tsamek1(L): Given a  of list of destinations with their probability [j,Pr[i->j]]
#ony mentions those destination with non-zero prob. Try:
#Tsamek1([[2,1/2],[3,1/2],[4,0],[5,0]]));
Tsamek1:=proc(L) local gu, i:
gu:=[]:

for i from 1 to nops(L) do
 if L[i][2]<>0 then
  gu:=[op(gu), L[i]]:
 fi:
od:
gu:
end:



#Tsamek(L): Given a list of lists L where L[i] is the list of destinations with their probability [j,Pr[i->j]]
#shrinks it as to only mention  destination with non-zero probabilities. Try:
#Tsamek([[[2,1/3],[3,2/3],[4,0]]]);
Tsamek:=proc(L) local i:
[seq(Tsamek1(L[i]),i=1..nops(L))]:
end:


CCb1:=proc():
[[1,2,3,4,5],
[0, 0, 4, 3, 5, 1, 2, 3, 1, 2, 4, 5, 0, 1, 3, 0, 0, 0, 5, 0, 5, 2, 4, 1, 2, 3, 5, 0, 4, 1, 0, 0, 5, 3, 4, 2, 0, 4, 1, 3], 
{5, 9, 23, 36, 40}]
end:


#CCb(C,D,P,S,T): The "Count Your Chickens board". It is a triple , the labels , followed by the board, followed by the set of blue squares
#An emptry square is indicated by 0. It is assumed that the last square (nops(gu[2])+1) has all of them. Try:
#CCb();
CCb:=proc(C,D,P,S,T) local gu1,gu2,gu3:
gu1:=[C,D,P,S,T]:
gu2:=[0,0,S,P,T,C,D,P,C,D,
      S,T,0,C,P,0,0,0,T,0,
      T,D,S,C,D,P,T,0,S,C,
      0,0,T,P,S,D,0,S,C,P]:
gu3:={5,9,23,36,40}:
[gu1,gu2,gu3]:
end:

#FindNext(CB,i,j): inputs a "Count your Chicken board, and a location i, and an outcome of the spinner  from 1 to 6, outputs
#the next location of the chicken. Try:
#FindNext(CCb(C,D,P,S,T),4,4);

FindNext:=proc(CB,i,j) local mu,i1:


if j=nops(CB[1])+1 then
 RETURN(i):
fi:

mu:=CB[1][j]:

for i1 from i+1 to nops(CB[2]) while CB[2][i1]<>mu do od:
i1:
end:


#NeiState(CB,SP,a): Given a Count Your Chicken Board, CB, with a prob. distribution on the animals, and a state a the location 1<=a<=nops(CB[2])+1,
#outputs the six  states that come from it, together with the proba. and the  gain of chicks. 
#The absorbing Losing state is indicated by 0 and winning by 1
#Try:
#NeiState(CCb1(),[(1/6)$5],1);
NeiState:=proc(CB,SP,a) local j,lu,a1,Pfox:

if a>=nops(CB[2])+1 then
 RETURN([a]):
fi:

Pfox:=1-convert(SP,`+`):

lu:=[]:

for j from 1 to nops(CB[1]) do
 a1:=FindNext(CB,a,j):
 
 if member(a1,CB[3]) then
  lu:=[op(lu),[a1,SP[j],a1-a+1] ]:
else
  lu:=[op(lu),[a1,SP[j],a1-a]]:
 fi:

od:

 lu:=[op(lu),[a,Pfox,-1]]:

lu:

end:


#PreUmb(f,s,t,WMC): Given a monomial f in s and t applies the pre-umbra operation implied by the Weighted Markov Chain w/o Truncation
#Try:
#PreUmb(s^5*t^11,s,t,CCMC(CCb1(),[(1/6)$5]););
PreUmb:=proc(f,s,t,CB) local d1,c,i,lu:

d1:=degree(f,s):
 c:=normal(f/(s^d1)):
if (normal(diff(c,s))<>0) then
 RETURN(FAIL):
fi:

if not (1<=d1 and d1<=nops(CB)) then
 RETURN(FAIL):
fi:


lu:=c*add(s^CB[d1][i][1]*CB[d1][i][2]*t^CB[d1][i][3],i=1..nops(CB[d1])):

lu:=expand(lu):

#lu:=c*add(s^CB[i][1]*CB[i][2]*t^CB[i][3],i=1..nops(CB)):

#lu:=expand(KillGadol(KillNeg(lu,t),t,nops(CB[2]))):

end:



#Umb(f,s,t,WMC): Given a polynomial f in s and t applies the umbral operation implied by the weighted Markov chain WMC
#Try:
#Umb(s^5*t^11/2+ s^7*t^12/2,s,t,CCMC(CCb1(),[(1/6)$5]) );
Umb:=proc(f,s,t,WMC) local i:

if not type(f,`+`) then
 RETURN(PreUmb(f,s,t,WMC)):
fi:

add(PreUmb(op(i,f),s,t,WMC),i=1..nops(f)):
end:


#CCMC(CB,SP): inputs a general Count Your Chickens! Style game, with the probability distribution on the set of animals
#outputs the weighted absorbing Markov Chain. The format of a weighted absorbing Markov Chain is
#a list L of length N, say, where there are N states labeled 1, ...,N, and the i-th entry is
# a set of triples [j,Prob(i->j),gain], and the convention is that the absorbing states are labelled N+1, .., N+r. Try: 
#CCMC(CCb1(),[(1/6)$5]);
#CCMC(CCb1(),[(1/6)$5]);

CCMC:=proc(CB,SP) local gu,i:

if nops(CB[1])<>nops(SP) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=[]:

for i from 1 to nops(CB[2]) do
 gu:=[op(gu),NeiState(CB,SP,i)]:
od:

gu:

end:



#HarogGadol(f,t,N): replaces each power larger than N by t^N. Try:
#HarogGadol(a*t+b*t^3+d*t^4,t,2);
HarogGadol:=proc(f,t,N) local i:
if degree(f,t)<=N then
  RETURN(f):
else
add(coeff(f,t,i)*t^i,i=ldegree(f,t)..N)+add(coeff(f,t,i),i=N+1..degree(f,t))*t^N:
fi:
end:


#HarogKatan(f,t,N): replaces each power less than N by t^N. Try:
#HarogKatan(a*t+b*t^3+d*t^4,t,2);
HarogKatan:=proc(f,t,N) local i:
if ldegree(f,t)>=N then
  RETURN(f):
else
add(coeff(f,t,i)*t^N,i=ldegree(f,t)..N-1)+add(coeff(f,t,i)*t^i,i=N..degree(f,t)):
fi:
end:





#PreUmbT(f,s,t,WMC,K,G): Given a monomial f in s and t applies the pre-umbra operation implied by the Weighted Markov Chain with Truncation
#where powers smaller than K are replaced by t^K and powers larger than G are replaced by t^G. Try:
#Try:
#PreUmbT(s^5*t^11,s,t,CCMC(CCb1(),[(1/6)$5],0,40););
PreUmbT:=proc(f,s,t,CB,K,G) local d1,c,i,lu:

d1:=degree(f,s):
 c:=normal(f/(s^d1)):
if (normal(diff(c,s))<>0) then
 RETURN(FAIL):
fi:

if not (1<=d1 and d1<=nops(CB)) then
 RETURN(FAIL):
fi:


lu:=c*add(s^CB[d1][i][1]*CB[d1][i][2]*t^CB[d1][i][3],i=1..nops(CB[d1])):

lu:=expand(lu):

expand(HarogGadol(HarogKatan(lu,t,K),t,G)):


end:

#UmbT(f,s,t,WMC,K,G): Given a polynomial f in s and t applies the umbral operation implied by the weighted Markov chain WMC
#with low-truncation K and high-truncation G
#Try:
#UmbT(s^5*t^11/2+ s^7*t^12/2,s,t,CCMC(CCb1(),[(1/6)$5]),0,40 );
UmbT:=proc(f,s,t,WMC,K,G) local i:

if not type(f,`+`) then
 RETURN(PreUmbT(f,s,t,WMC,K,G)):
fi:

add(PreUmbT(op(i,f),s,t,WMC,K,G),i=1..nops(f)):
end:



#Pashet(f,x,t,eps): inputs a polynomial f in the variables x and t and a cut-op eps, discards
#terms whose absolute value is less than eps. Try:
#Pashet(x*t+1/10^11*x-1/10^11,x,t,1/10^10);
Pashet:=proc(f,x,t,eps) local i,j,gu,f1:
gu:=0:
for i from ldegree(f,x) to degree(f,x) do
 f1:=coeff(f,x,i):
 for j from ldegree(f1,t) to degree(f1,t) do
  if abs(coeff(f1,t,j))>eps then
    gu:=gu+coeff(f1,t,j)*x^i*t^j:
  fi:
 od:
od:
gu:
end:



#WMCSer(WMC,t,X,K,st): The inputs are (i) WMC: Weighted Markov chain according to our convention
#(ii) variables t and X
#(iii) a  positive integer K
#(iv) the starting state (normally 1)
#Let N be the number of non-absorbing states in WMC.
#By convention that absorbing states are N+1,N+2, .., N+r where is is the number of absorbing states (found from the WMC thatmust adhere to our convention)
#The output is a list of length r, where the j-th entry is the
#polynomial of degree K in X where the coeff. of X^i is the
#prob. generating function, according to gain, if it ended after i rounds and winded up at state N+j.
#Try:
#WMCSer(CCMC(CCb1(),[(1/6)$5]),t,X,30,1);
WMCSer:=proc(WMC,t,X,K,st) local N, ABS,gu,mu,i,j,s,r,sofi:
N:=nops(WMC):

if not (st>=1 and st<=N) then
 RETURN(FAIL):
fi:

ABS:= {seq(seq(WMC[i][j][1],j=1..nops(WMC[i])),i=1..nops(WMC))} minus {seq(i,i=1..N)}:

r:=nops(ABS):

if ABS<>{seq(N+i,i=1..r)} then
 RETURN(FAIL):
fi:




gu:=0:

mu:=s^st:


for i from 1 to K do


 mu:=Umb(mu,s,t,WMC):


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

sofi:=add(coeff(mu,s,N+j)*s^(N+j),j=1..r):


 gu:=expand(gu+sofi*X^i):



 mu:=expand(mu-sofi):
od:


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



gu:=expand(normal(gu/subs({X=1,s=1,t=1},gu))):


[[seq(coeff(gu,s,N+j),j=1..r)],evalf(subs({s=1,t=1},mu))]:

end:



#WMCSerT(WMC,t,X,K,st,Kat,Gad): Like WMCSet(WMC,t,X,K,t) (q.v.) but with the additional inputs K and G
#indicating that powers of t lower than t^Kat are idenitified with t^Kat, and
#powers of t higher than t^Gad are idenitified with t^Gad, 
#Try:
#WMCSerT(CCMC(CCb1(),[(1/6)$5]),t,X,30,1,0,40);
WMCSerT:=proc(WMC,t,X,K,st,Kat,Gad) local N, ABS,gu,mu,i,j,s,r,sofi:
N:=nops(WMC):

if not (st>=1 and st<=N) then
 RETURN(FAIL):
fi:

ABS:= {seq(seq(WMC[i][j][1],j=1..nops(WMC[i])),i=1..nops(WMC))} minus {seq(i,i=1..N)}:

r:=nops(ABS):

if ABS<>{seq(N+i,i=1..r)} then
 RETURN(FAIL):
fi:




gu:=0:

mu:=s^st:


for i from 1 to K do


 mu:=UmbT(mu,s,t,WMC,Kat,Gad):


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

sofi:=add(coeff(mu,s,N+j)*s^(N+j),j=1..r):


 gu:=expand(gu+sofi*X^i):



 mu:=expand(mu-sofi):
od:


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



gu:=expand(normal(gu/subs({X=1,s=1,t=1},gu))):


[[seq(coeff(gu,s,N+j),j=1..r)],evalf(subs({s=1,t=1},mu))]:

end:




#RandSet1(S,k): a random subset with k elements of the set S, try:
#RandSet1({1,2,3,4},2);
RandSet1:=proc(S,k) local gu,S1,i,a:
if k>nops(S) then
 RETURN(FAIL):
fi:

if k=0 then
 RETURN({}):
fi:

gu:={}:
S1:=S:

for i from 1 to k do
 a:=S1[rand(1..nops(S1))()]:
 gu:=gu union {a}:
 S1:=S1 minus {a}:
od:
gu:
end:


#RandDG(N,s,deg1,deg2): a random directed graph with N vertices labelled 1, ..., N that have positive in-degree and s vertices
#labelled N+1, .., N+s,the output is a list of length N whose i-th entry is the the set of out-neighbors of i
#with each vertex with a random outdegree between deg1 deg2. Try:
#RandDG(7,2,2,3);
RandDG:=proc(N,s,deg1,deg2) local i,S1,S2,d,T,i1,j,ra:

if deg2>N-2 then
 RETURN(FAIL):
fi:

ra:=rand(deg1..deg2):
S1:={seq(i,i=1..N)}:
S2:={seq(i,i=N+1..N+s)}:


for i from 1 to N do
 d:=ra():
T[i]:=RandSet1(S1 union S2,d):
od:

for j from 1 to s do
i1:=rand(1..N):
T[i1]:=T[i1] union {N+j}:
od:

[seq(T[i1],i1=1..N)]:

end:

#RandWMCe(N,s,deg1,deg2,T1,T2): a random Weighted Markov Chain where the weights are integers drawn randomly between T1 and T2
#and the underlying directed graph with s sinks is RandDG(N,s,deg1,deg2) (q.v.). Try:
#RandWMCe(7,2,2,3,0,10); 
RandWMCe:=proc(N,s,deg1,deg2,T1,T2) local G,gu,j,i,ra:
G:=RandDG(N,s,deg1,deg2):

ra:=rand(T1..T2):

gu:=[]:

for i from 1 to nops(G) do
gu:=[op(gu),{seq([G[i][j],1/nops(G[i]),ra()],j=1..nops(G[i]))}]:
od:
gu:

end:



#Info(WMC,K,st): Inputs a weighted Markov Chain with absorbing states WMC (according to our convention)
#a positive integer K, and an intial state st (from 1 to nops(WMC)), outputs a list of
#assuming that it terminated after K steps
#(i) The probabilities of lending in each of the absorbing states
#(ii) The list of conditional expectation of the gain at each state
#(iii) The list of conditional expectation of the number of rounds at each state
#(iv) The probability of finishing after K steps. Here we do NOT use truncation, i.e. there is no limit to your income and you can be in debt
#Try
#Info(CCMC(CCb1(),[(1/6)$5]),50,1);
Info:=proc(WMC,K,st) local X,t,gu,eps,mu,mu1,mu2,j:

gu:=WMCSer(WMC,t,X,K,st):

eps:=gu[2]:

gu:=gu[1]:



mu:=evalf(subs({X=1,t=1},gu)):



mu1:=[]:

for j from 1 to nops(gu) do

if subs({X=1,t=1},gu[j])=0 then
   mu1:=[op(mu1),0]:
else
  mu1:=evalf([op(mu1),subs({X=1,t=1},diff(gu[j],t))/subs({X=1,t=1},gu[j])]):
fi:
od:

mu2:=[]:
for j from 1 to nops(gu) do

if subs({X=1,t=1},gu[j])=0 then
   mu2:=[op(mu2),0]:
else
  mu2:=evalf([op(mu2),subs({X=1,t=1},diff(gu[j],X))/subs({X=1,t=1},gu[j])]):
fi:
od:

[mu,mu1,mu2,eps]:

end:




#InfoT(WMC,K,st,Kat,Gad): Inputs a weighted Markov Chain with absorbing states WMC (according to our convention)
#a positive integer K, and an intial state st (from 1 to nops(WMC)), outputs a list of
#assuming that it terminated after K steps
#(i) The probabilities of lending in each of the absorbing states
#(ii) The list of conditional expectation of the gain at each state
#(iii) The list of conditional expectation of the number of rounds at each state
#(iv) The probability of finishing after K steps. 
#Here we DO use truncation where fortune less than Kat becomes Kat and fortune larger than Gad becomes Gad
#Try
#InfoT(CCMC(CCb1(),[(1/6)$5]),50,1,0,40);
InfoT:=proc(WMC,K,st,Kat,Gad) local X,t,gu,eps,mu,mu1,mu2,j:

gu:=WMCSerT(WMC,t,X,K,st,Kat,Gad):

eps:=gu[2]:

gu:=gu[1]:



mu:=evalf(subs({X=1,t=1},gu)):



mu1:=[]:

for j from 1 to nops(gu) do

if subs({X=1,t=1},gu[j])=0 then
   mu1:=[op(mu1),0]:
else
  mu1:=evalf([op(mu1),subs({X=1,t=1},diff(gu[j],t))/subs({X=1,t=1},gu[j])]):
fi:
od:

mu2:=[]:
for j from 1 to nops(gu) do

if subs({X=1,t=1},gu[j])=0 then
   mu2:=[op(mu2),0]:
else
  mu2:=evalf([op(mu2),subs({X=1,t=1},diff(gu[j],X))/subs({X=1,t=1},gu[j])]):
fi:
od:

[mu,mu1,mu2,eps]:

end: