> #Weiji Zheng, 10/25/2020, Assignment 14
;
> #OK TO POST HOMEWORK
;
> 
;
> #Code for understanding
;
> 
> Help14:=proc(): print(` LC(p) , LDie(A), Paths(m,n), GPaths(m,n), NuPaths(m,n), NuGPaths(m,n) , RPath(m,n)  , CycShi(L) , DvoMot(n), DrawP(P) `): end:
> 
> #LC(p): inputs a rational number p outputs 1 with probability p (and 0 otherwise)
> #LC is short for "Loaded Coin"
> LC:=proc(p) local a,b,r:
> a:=numer(p):
> b:=denom(p):
> r:=rand(1..b)():
> if  r<=a then
>   RETURN(1):
> else
>  RETURN(0):
> fi:
> end:
> 
> 
> #LDie(A): inputs a list of non-negative integers [a[1],..., a[k]] (say) outputs i with probability a[i]/(a[1]+...+a[k]). Try
> #LDie([1,2,3]);
> #LDie is short for "Loaded die"
> LDie:=proc(A) local N,i,j,r:
> 
> N:=add(A[i],i=1..nops(A)):
> r:=rand(1..N)():
> 
> for i from 1 to nops(A) while add(A[j],j=1..i-1)<r do od:
> i-1:
> 
> end:
> 
> #Paths(m,n): The set of lattice paths from [0,0] to [m,n] with atomic steps [1,0] (Called 1) , [0,1] (called -1)
> Paths:=proc(m,n) local S1,S2,s:
> option remember:
> 
> if m<0 or n<0 then
>  RETURN({}):
> fi:
> if m=0 and n=0 then
>  RETURN({[]}):
> fi:
> S1:=Paths(m-1,n):
> S2:= Paths(m,n-1):
> {seq([op(s),1], s in S1)  , seq([op(s),-1], s in S2)  }:
> end:
> 
> 
> 
> #GPaths(m,n): The set of lattice paths from [0,0] to [m,n] with atomic steps [1,0] (Called 1) , [0,1] (called -1)
> #that stay in m>=n
> GPaths:=proc(m,n) local S1,S2,s:
> option remember:
> 
> if ( m<0 or n<0 or m<n) then
>  RETURN({}):
> fi:
> if m=0 and n=0 then
>  RETURN({[]}):
> fi:
> S1:=GPaths(m-1,n):
> S2:= GPaths(m,n-1):
> {seq([op(s),1], s in S1)  , seq([op(s),-1], s in S2)  }:
> end:
> 
> 
> NuPaths:=proc(m,n) option remember:
> if m<0 or n<0 then 
> 0:
> elif m=0 and n=0 then
> 1:
> else
> NuPaths(m-1,n)+NuPaths(m,n-1):
> fi:
> end:
> 
> 
> NuGPaths:=proc(m,n) option remember:
> if (m<0 or n<0 or m<n)  then 
> 0:
> elif m=0 and n=0 then
> 1:
> else
> NuGPaths(m-1,n)+NuGPaths(m,n-1):
> fi:
> end:
> 
> 
> 
> #RPath(m,n): A (uniformaly random) path from [0,0] to [m,n]. Try:
> #RPath(6,5);
> RPath:=proc(m,n) local c:
> if m=0 and n=0 then
>  RETURN([]):
> fi:
> 
> c:=LDie([NuPaths(m-1,n),NuPaths(m,n-1)]):
> 
> if c=1 then
>  [op(RPath(m-1,n)),1]:
> else
>  [op(RPath(m,n-1)),2]:
> fi:
> end:
> 
> 
> #RGPath(m,n): A (uniformaly random) path from [0,0] to [m,n]. Try:
> #RGPath(6,5);
> RGPath:=proc(m,n) local c:
> if m=0 and n=0 then
>  RETURN([]):
> fi:
> 
> c:=LDie([NuGPaths(m-1,n),NuGPaths(m,n-1)]):
> 
> if c=1 then
>  [op(RGPath(m-1,n)),1]:
> else
>  [op(RGPath(m,n-1)),2]:
> fi:
> end:
> 
> 
> 
> #CycShi(L): The set of cyclic shifts of the list L
> CycShi:=proc(L) local i:
> {seq([op(i..nops(L),L),op(1..i-1,L)],i=1..nops(L))}:
> end:
> 
> 
> #DvoMot(n): illustration and empirical confirmation of the Dvoretzky-Motzkin proof that
> #the number of good paths (namely NuGPaths(n,n) ) is given by the Catalan number
> #(2*n+1)!/(n!*(n+1)!) by showing that the associated set of GPaths(n,n) when 1 is
> #appended to the end, all have distinct cyclic shifts (so 2*n+1 of them) and
> #their union is Paths(n+1,n) whose number of elements is binomial(2*n+1,n)=(2*n+1)!/(n!*(n+1)!)
> DvoMot:=proc(n) local S,s:
> 
> #S is the set of Good paths from [0,0] to [n,n]
> S:=GPaths(n,n):
> #Now we consider the associated set where we append 1 to each
> S:={seq([op(s),1], s in S)}:
> 
> #Now we check empirically that all the cyclic shifts are distinct
> 
> for s in S do
> 
> if nops(CycShi(s))<>2*n+1 then
>   print(`Something is wrong`):
>  RETURN(false):
> fi:
> od:
> 
> 
> #We now check that the union of the cyclic shifts of the members of S is Paths(n+1,n)
> 
> evalb({seq(op(CycShi(s)), s in S)}=Paths(n+1,n)):
> end:
> 
> #DrawP(P): Draws the path P
> DrawP:=proc(P) local P1,Pt,i:
> P1:=[[0,0]]:
> Pt:=[0,0]:
> for i from 1 to nops(P) do
> if P[i]=1 then
>  Pt:=Pt+[1,0]:
> else
>  Pt:=Pt+[0,1]:
> fi:
> P1:=[op(P1),Pt]:
> od:
> 
> plot(P1,axes=none):
> 
> end:

> #Q1
;
> #Use a procedure CountRuns(L,k)that inputs a list L of numbers (or for that matter, any objects) and outputs the number of "runs" of length k in L
> #[A run of length k is a string L is an occurence of an i such that L[i]=L[i+1]=...=L[i+k-1]
> #[Hint: a quick way to find whether location i is the beginning of a run of length k (where i is between 1 and nops(L)-k+1, is to see if nope({op(i..i+k-1,L)})=1]
;
> 
;
> 
;
> CountRuns := proc(L,k) local i, count:
> option remember:
> count := 0:
> for i from 1 to nops(L)-k+1 do
>     if nops({op(i..i+k-1,L)}) = 1 then
>        count := count + 1:
>     fi:
> od:
> 
> count:
> end:
> CountRuns([1,2,2,2,1,1,1,1,2,2,2],3)
4
;
> 
;
> #Q2
;
> #Use the above CountRuns(L,k) that you just wrote, combined with procedure RPath in M14.txt to write a procedure
> #AvePathRuns(m,n,k,N) that inputs positive inetgers m,n,k, and a LARGE positive integer N, generated N random paths, for each of them #records the number of runs, and finds the sample average and standard-deviation (the squre-root of the variance)
> #
> #Run
> #
> #AvePathRuns(200,200,5,3000)
> #
> #ten times and see what you get. Are they close to each other?
;
> 
;
> with(Statistics):
> AvePathRuns:=proc(m,n,k,N) local i,run,S,std,avg:
> 
> run := 0:
> S := []:
> for i from 1 to N do
>     run := CountRuns(RPath(m,n),k):
>     S := [op(S), run]:
> od:
> avg := Mean(S):
> std := StandardDeviation(S):
> print(avg): 
> print(std):
> 
> end:
> AvePathRuns(200,200,5,3000)
Warning,  computation interrupted
;
> #i takes too long
;
> 
;
> #Q3
;
> #use RGPath now
;
> AveGPathRuns:=proc(m,n,k,N) local i,run,S,std,avg:
> 
> run := 0:
> S := []:
> for i from 1 to N do
>     run := CountRuns(RGPath(m,n),k):
>     S := [op(S), run]:
> od:
> avg := Mean(S):
> std := StandardDeviation(S):
> print(avg): 
> print(std):
> 
> end:
> 
;
Warning,  computation interrupted
;
> #Q4
;
> #A maximal run is a run L[i]=L[i+1]=...L[i+k-1] AND L[i-1]<> L[i] and L[i+k]<>L[i]
> #[Hint: a quick way to find whether location i is the beginning of a MAXIMAL run of length k (where i is between 1 and nops(L)-k+1, is to #see if nope({op(i..i+k-1,L)})=1 and L[i-1]<>L[i] and L[i+k]<>L[i]]
> #
> #Write procedures
> #
> #CountMaximalRuns(L,k), AvePathMaximalRuns(m,n,k,N), and AveGPathMaximalRuns(m,n,k,N)
> #
> #and do the analogous thing for the above two problems, (that you did for Runs) for Maximal Runs.
;
> 
;
> CountMaximalRuns := proc(L,k) local n, i, S:
> option remember:
> S := []:
> n := nops(L)-k+1:
> 
> for i from 1 to n do
> if i=1 then
>    if nops({op(i..i+k-1, L)})=1 and L[i+k]<>L[i] then 
>    S := [op(S), i]:
>    fi:
> fi:
> if i=n then
>    if nops({op(i..i+k-1,L)})=1 and L[i]<>L[i-1] then
>    S := [op(S), i]:
>    fi:
> fi:
> if i<>1 and i<>n and nops({op(i..i+k-1,L)})=1 and L[i+k]<>L[i] and L[i]<>L[i-1] then
>    S:=[op(S), i]:
> fi:
> od:
> nops(S):
> end:
> 
;
> #use CountMaximalRunsuns now
;
> with(Statistics):
> AvePathRuns:=proc(m,n,k,N) local i,run,S,std,avg:
> 
> run := 0:
> S := []:
> for i from 1 to N do
>     run := CountMaximalRuns(RPath(m,n),k):
>     S := [op(S), run]:
> od:
> avg := Mean(S):
> std := StandardDeviation(S):
> print(avg): 
> print(std):
> 
> end:
> AveGPathRuns:=proc(m,n,k,N) local i,run,S,std,avg:
> 
> run := 0:
> S := []:
> for i from 1 to N do
>     run := CountMaximalRuns(RGPath(m,n),k):
>     S := [op(S), run]:
> od:
> avg := Mean(S):
> std := StandardDeviation(S):
> print(avg): 
> print(std):
> 
> end:
> 
;
> #Q5
;
> #Write a procedure GoodBrother(L)
> #input an arbitrary list of ODD length that adds up to 1 (so if the length is 2*n+1 it must have n+1 1's and n -1's), and outputs the UNIQUE cyclic shift that is a so-called Dyck path of length 2*n with 1 appended at the end.
> #For example, GoodBrother([[1,-1,-1,1,1,1,1,-1,-1,-1,1]=[1,1,1,-1,-1,-1,1,1,-1,-1,1]
;
> 
;
> GoodBrother := proc(L) local S,i,n:
> 
> S:=CycShi(L):
> end:
> 
;
> #Q6
;
> # if two cyclic shifts of such an L are identical, it means that there is some word w such that L is w repeated k times. The sum of L is k times the sum of w. Show that this can only happen if k=1
;