#1. 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]

#For example

#CountRuns([1,2,2,2,1,1,1,1,2,2,2],3)

#is 4 for the runs starting at i=2, i=5, i=6, i=9. Note that it can be part of a longer run.

#[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 isRun, n, i, j:
n:=0:
for i from 1 to nops(L) do:
 isRun:=true:
 for j from 1 to k-1 do:
  if i+j>nops(L) then
   isRun:=false:
  elif i+j<nops(L) and L[i+j]<>L[i] then
   isRun:=false:
  fi:
 od:
 if isRun=true then
  n:=n+1:
 fi:
 isRun:=true:
od:
n:
end:

#2. 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):

WtEn:=proc(L,x) local i:
add(x^L[i],i=1..nops(L)):
end:

AveClever:=proc(L) local i,f,x: 
f:=WtEn(L,x):
subs(x=1,diff(f,x))/subs(x=1,f):
end:

AvePathRuns:=proc(m,n,k,N) local i, random, runs, ave, stdev: 
runs:=[]:
for i from 1 to N do:
 random:=RPath(m,n):
 runs:=[op(runs),CountRuns(random,k)]:
od:
ave:=AveClever(runs):
stdev:=StandardDeviation(Array(runs)):
[ave,stdev]:
end:

DoTenTimes:=proc(m,n,k,N) local i, l:
l:=[]:
for i from 1 to 10 do:
 l:=[op(l),AvePathRuns(m,n,k,N)]:
od:
l:
end:

#evalf(DoTenTimes(200,200,5,3000))

#[[24.06900000, 7.390925580],
#[24.31766667, 7.441093861],
#[24.22833333, 7.354424501],
#[24.16066667, 7.635506087],
#[24.27100000, 7.420551228],
#[24.17133333, 7.508713060],
#[24.05866667, 7.504976112],
#[24.20566667, 7.624483292],
#[24.11966667, 7.292809809],
#[24.15100000, 7.598867016]]

#Yeah they are pretty close to each other.

#3. Do the analogous thing for Good paths, writing
#AveGPathRuns(m,n,k,N)

#[Hint: you only have to change one word in the code]

#Run

#AvePathGRuns(200,200,5,3000)

#ten times. How does it compare with the previous output for all paths?

AveGPathRuns:=proc(m,n,k,N) local i, random, runs, ave, stdev:
runs:=[]:
for i from 1 to N do:
 random:=RGPath(m,n):
 runs:=[op(runs),CountRuns(random,k)]:
od:
ave:=AveClever(runs):
stdev:=StandardDeviation(Array(runs)):
[ave,stdev]:
end:

DoTenTimesG:=proc(m,n,k,N) local i, l:
l:=[]:
for i from 1 to 10 do:
 l:=[op(l),AveGPathRuns(m,n,k,N)]:
od:
l:
end:

#evalf(DoTenTimesG(200,200,5,3000))

#[[24.30066667, 7.58982661494345],
#[24.27500000, 7.57811658331519],
#[23.96033333, 7.36864023018062],
#[24.11666667, 7.36626156826195],
#[24.00066667, 7.28164406818044],
#[24.15533333, 7.44314035438004],
#[24.20000000, 7.48264636780030],
#[24.19066667, 7.46491043827761],
#[24.01900000, 7.39681530292205],
#[24.32233333, 7.41656499787321]]

#It is also very close to the previous output for all the paths.

#4. 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 isMaxRun, n, i, j:
n:=0:
for i from 1 to nops(L) do:
 isMaxRun:=true:
 if i <> 1 and L[i-1]=L[i] then
  isMaxRun:=false
 elif i+k<nops(L) and L[i+k]=L[i] then
  isMaxRun:=false
 else
  for j from 1 to k-1 do:
   if i+j>nops(L) then
    isMaxRun:=false:
   elif i+j<nops(L) and L[i+j]<>L[i] then
    isMaxRun:=false:
   fi:
  od:
 fi:
 if isMaxRun=true then
  n:=n+1:
 fi:
 isMaxRun:=true:
od:
n:
end:

AvePathMaximalRuns:=proc(m,n,k,N) local i, random, runs, ave, stdev:
runs:=[]:
for i from 1 to N do:
 random:=RPath(m,n):
 runs:=[op(runs),CountMaximalRuns(random,k)]:
od:
ave:=AveClever(runs):
stdev:=StandardDeviation(Array(runs)):
[ave,stdev]:
end:

AveGPathMaximalRuns:=proc(m,n,k,N) local i, random, runs, ave, stdev:
runs:=[]:
for i from 1 to N do:
 random:=RGPath(m,n):
 runs:=[op(runs),CountMaximalRuns(random,k)]:
od:
ave:=AveClever(runs):
stdev:=StandardDeviation(Array(runs)):
[ave,stdev]:
end:

DoTenTimesMax:=proc(m,n,k,N) local i, l:
l:=[]:
for i from 1 to 10 do:
 l:=[op(l),AvePathMaximalRuns(m,n,k,N)]:
od:
l:
end:

DoTenTimesMaxG:=proc(m,n,k,N) local i, l:
l:=[]:
for i from 1 to 10 do:
 l:=[op(l),AveGPathMaximalRuns(m,n,k,N)]:
od:
l:
end:

#evalf(DoTenTimesMax(200,200,5,3000)):

#[[6.181666667, 2.30718549836451],
#[6.205000000, 2.33925893766790],
#[6.310666667, 2.40362476818470],
#[6.277333333, 2.38823242313736],
#[6.298666667, 2.36276514037933],
#[6.279333333, 2.37567972875416],
#[6.268000000, 2.40820316435864],
#[6.255000000, 2.39928333470088],
#[6.232000000, 2.38594640483581],
#[6.223000000, 2.33432803231554]]

#All test data points are close to each other, as before. 

#evalf(DoTenTimesMaxG(200,200,5,3000)):

#[[6.234333333, 2.36740844873530],
#[6.296333333, 2.34542538337189],
#[6.262666667, 2.36238576677461],
#[6.304000000, 2.33910405980599],
#[6.264000000, 2.34139517747981],
#[6.288333333, 2.35861221330393],
#[6.250666667, 2.34854666677875],
#[6.246333333, 2.40060565789770],
#[6.324000000, 2.36365951343077],
#[6.164000000, 2.34327427149217]]

#Correspondingly, these points are very similar to the those tested for all the paths, like before.

#5. Write a procedure
#GoodBrother(L)

#that inputs 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]

GoodShift:=proc(L) local i, k:
k:=0:
for i in L do:
 k:=k+i:
 if k<0 then
   return false:
 fi:
od:
true:
end:

GoodBrother:=proc(L) local i, shifts:
shifts:=CycShi(L):
for i in shifts do:
 if i[nops(L)]=1 and GoodShift(i) then
  return i:
 fi:
od:
end:

#Test
#GoodBrother([1, -1, -1, 1, 1, 1, 1, -1, -1, -1, 1]);
#            [1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1] 

#6. [Human exercise] Prove that the all the cyclic shifts of a word in {1,-1} that adds up to 1 are distinct.
#[Hint: 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.

#Proof by contradiction
#Suppose that there are some cyclic shifts of L which are not distinct.
#That means L is one of its cyclic shifts w repeated k times where k is not equal to 1 (otherwise l=w and then there are no repeated cyclic shifts)
#Then L= w_1 w_2 ... w_k
#Sum(L)=Sum(w)*k	
#1=Sum(w)*k		L adds up to 1 so Sum(L)=1
#Since Sum(w) and k must both be integers, the only way this can occur is if Sum(w) and k are both 1.
#This is a contradiction since we assumed that k is not equal to 1.
#Thus proved that all the cyclic shifts of a word in {1,-1} that adds up to 1 are distinct.

#7. [Optional Challenge] Use the Wilf methodology discussed in the lecture to write a procedure
#RandSetPartition(n)

#That inputs a positive integer n and outputs a set-partion of {1, ..., n} where each of them has the same probability of being picked.

#Hint: One way of doing it is first write a prelininary procedure

#RandSetPartition1(n,k)

#that inputs positive integers n and k (k between 1 and n) and outputs a random set-partition with exactly k parts, using Snk(n,k) of M11.txt doing 
#a LDie([Snk(n-1,k-1), k*Snk(n-1,k)]) and if the die lands on Snk(n-1,k-1) recursively generate a set partions of {1, ..., n-1} with k-1 members 
#[calling RandSetPartition1(n-1,k-1) ] and adjoin {n} to it, and in the latter case generate a random set partion of {1, ..., n-1} with k members, 
#[calling RandSetPartition1(n-1,k) ] and then rand(1..k) to decide which of the k sets should invite n to join.

#Having done that RandSetPartition(n) would use

#LDie([seq(Snk(n,k),k=1..n)])

#To pick a k, then use the above procedure RandSetPartition1(n,k) with that k