#Nathan Fox
#Homework 21
#I give permission for this file to be posted online

##Read old files
read(`C21.txt`):

#Help procedure
Help:=proc():
print(` FixRatioXY(R,S,T,P,p,Ratio) , RandWalk(M,n) `):
print(` RandUnif() , RandStochMat(m) , TestRandWalk(m,n,k) `):
end:

##Problem 1
#FixRatioXY(R,S,T,P,p,Ratio): inputs numbers (or symbols!)
#R,S,T,P, a symbol p, and a number Ratio, and outputs a general
#vector p=[p[1],p[2],p[3],p[4]] such that
#sX(R,S,T,P,p,q)=Ratio*sY(R,S,T,P,p,q):
#For ALL q.
FixRatioXY:=proc(R,S,T,P,p,Ratio) local eq,i,freeman,dyson,q,tent:

eq:={seq(sY(R,S,T,P,p,[0$(i-1),1,0$(4-i)])*Ratio=
    sX(R,S,T,P,p,[0$(i-1),1,0$(4-i)]),i=1..4),
    sY(R,S,T,P,p,[1/3,1/4,3/5,2/5])*Ratio=
    sX(R,S,T,P,p,[1/3,1/4,3/5,2/5])}:

freeman:={solve(eq,{seq(p[i],i=1..4)})}:

dyson:={}:

for i from 1 to nops(freeman) do
 tent:=subs(freeman[i],[p[1],p[2],p[3],p[4]]):
 if normal(sY(R,S,T,P,tent,q)*Ratio-sX(R,S,T,P,tent,q))=0 then
  dyson:=dyson union {tent}:
 fi:
od:

dyson:

end:

##Problem 2
#RandUnif(): random floating-point number in the interval [0,1)
RandUnif:=proc()
rand()/1000000000000:
end:

#RandWalk(M,n): inputs a stochastic matrix M and a positive integer
#n. Simulate a random walk on a graph given by transition matrix M,
#starting at a random vertex.  Outputs a vector of the fraction of
#time spent at each vertex
RandWalk:=proc(M,n) local ret, vtx, i, m, rval, sm, j:

m:=nops(M):
vtx:=rand(1..m)():
ret:=[seq(0, i=1..m)]:
ret[vtx]:=ret[vtx] + 1:

for i from 1 to n-1 do
    rval:=RandUnif():
    sm:=0:
    for j from 1 to m do
        sm:=sm + M[vtx][j]:
        if sm > rval then
            vtx:=j:
            break:
        fi:
    od:
    ret[vtx]:=ret[vtx] + 1:
od:

[seq(ret[i]/n,i=1..m)]:

end:

#RandStochMat(m): a random m by m stochastic matrix
RandStochMat:=proc(m) local M,i,j,sm,rval:

M:=[seq([seq(0,i=1..m)],j=1..m)]:

for i from 1 to m do
    sm:=0:
    for j from 1 to m-1 do:
        rval:=RandUnif()*(1-sm):
        M[i][j]:=rval:
        sm:=sm + rval:
    od:
    M[i][m]:=1 - sm:
od:

M:

end:

#TestRandWalk(m,n,k): inputs positive integers m, n, and k. Outputs
#the average (arithmetic mean) l2 norm of k vectors
#SS(M)-RandWalk(M,n), where M is a random m by m stochastic matrix.
TestRandWalk:=proc(m,n,k) local sm, ss, rw, M, i, j, l2, v:

sm:=0:

for i from 1 to k do
    M:=RandStochMat(m):
    ss:=SS(M):
    rw:=RandWalk(M,n):
    v:=ss - rw:
    l2:=0:
    for j from 1 to m do
        l2:=l2 + v[j]^2:
    od:
    sm:=sm + sqrt(l2):
od:

sm/k:

end:

#evalf(TestRandWalk(5,100,100)); returned 0.06684220851
#evalf(TestRandWalk(10,100,100)); returned 0.06765910220
#evalf(TestRandWalk(2,300,300)); returned 0.02953965825
#These numbers are all close to 0, as they should be