#OK to post homework
#Blair Seidler, 2021-02-07, Assignment 4

with(combinat):

Help:=proc():
print(` CheckNuSnk() `,` RnkFast(n,k) `): 
print(`EstimateAveMax(n,k,K)`,`CalcAveMax(n,k)`,`CalcAveMaxBF(n,k)`):
print(`EstimateAveSum(n,k,K)`,`CalcAveSum(n,k)`,`CalcAveSumBF(n,k)`):
print(`Wnk(a,b)`,`NuWnk(a,b)`,`RWnk(a,b)`):
print(`W3nk(a,b,c)`,`NuW3nk(a,b,c)`,`RW3nk(a,b,c)`):
print(`WORDS(L)`,`NuWORDS(L)`,`RandWORD(L)`):
end:

# 1.

#CheckNuSnk(): Prove that NuSnk(n,k)=n!/(k!*(n-k)!). You can do it by hand, but you are
#welcome to use Maple to verify that G(n,k):=n!/(k!*(n-k)!) satisfies the equation
#G(n,k)=G(n-1,k)+G(n-1,k-1)
CheckNuSnk:=proc() local G,lhs,rhs:

G:=(n,k)->n!/(k!*(n-k)!):
lhs:=G(n,k):
rhs:=G(n-1,k)+G(n-1,k-1):

print(lhs=rhs):
print(lhs=simplify(rhs)):
evalb(lhs=simplify(rhs)):

end:
#### This procedure returns true, so G satisfies the same recurrence as NuSnk ####

# 2.

#RnkFast(n,k): Inputs NON-NEG. integer n and integer k and outputs
#the a (Uniformly at Random) k-element subset of {1,...,n}
#RnkFast(3,2) should output each of {1,2},{1,3},{2,3} with prob. 1/3
RnkFast:=proc(n,k) local c:
if not (type(n,integer) and type(k,integer) and n>=0) then
 RETURN(FAIL):
fi:
if k<0 or k>n then
 RETURN(0):
fi:

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

c:=LD([n-k,k]):
if c=1 then
 RETURN(RnkFast(n-1,k)):
else
 RETURN(RnkFast(n-1,k-1) union {n}):
fi:

end:

#### Why is this valid? ####
#Using the notation from question 1, G(n-1,k)=(1/k)*G(n-1,k-1) and
#G(n-1,k-1)=(1/(n-k))*G(n-1,k-1), so the ratio G(n-1,k)/G(n-1,k-1)=(n-k)/k
#This means that an object chosen uniformly at random from the union of
#Snk(n-1,k) and Snk(n-1,k-1) has a probability of (n-k)/n of being in the first
#set and k/n of being from the second.

#### Is RnkFast(5000,2000) faster than Rnk(5000,2000)? ####
#Certainly true for the first run, where RnkFast takes from 0.3 to 5.3 seconds and
#Rnk takes between 4 and 6 minutes. Because I have "option remember" in NuSnk,
#after the first run Rnk and RnkFast run in about the same time.

# 3.

#EstimateAveMax(n,k,K): uses RnkFast(n,k) K times, and estimates the average of the quantity
#(max(S)) over all k-elemenet subsets of {1,...,n}. Experiment it with various n and k, and K=1000.
EstimateAveMax:=proc(n,k,K) local i:

add(seq(max(RnkFast(n,k)),i=1..K))/K:

end:

#### I don't know if this counts as finding the exact value as an expression in n and k ####
#### but it does seem to calculate the exact value. It also runs an awful lot faster    ####
#### than the brute force calculation in the procedure below it.                        ####
#CalcAveMax(n,k): exactly calculates the average of the quantity (max(S)) over all
#k-elemenet subsets of {1,...,n}.
CalcAveMax:=proc(n,k) local i:

sum(binomial(i-1,k-1)*i,i=k..n)/binomial(n,k):

end:

#CalcAveMaxBF(n,k): exactly calculates the average of the quantity (max(S)) over all
#k-elemenet subsets of {1,...,n} by brute force.
CalcAveMaxBF:=proc(n,k) local S,s:

S:=Snk(n,k):
add(seq(max(s),s in op(S)))/nops(S):

end:

# 4.

#EstimateAveSum(n,k,K): uses RnkFast(n,k) K times, and estimates the average of the quantity
#(sum(S)) over all k-elemenet subsets of {1,...,n}. Experiment it with various n and k, and K=1000.
EstimateAveSum:=proc(n,k,K) local i:

add(seq(add(op(RnkFast(n,k))),i=1..K))/K:

end:

#### There are k elements in each subset, and the average value of an item across all ####
#### subsets is just (n+1)/2. The brute force routine below gives the same results,   ####
#### just much more slowly.                                                           ####
#CalcAveSum(n,k): exactly calculates the average of the quantity (sum(S)) over all
#k-elemenet subsets of {1,...,n}.
CalcAveSum:=proc(n,k):

(k*(n+1))/2:

end:

#CalcAveSumBF(n,k): exactly calculates the average of the quantity (max(S)) over all
#k-elemenet subsets of {1,...,n} by brute force.
CalcAveSumBF:=proc(n,k) local S,s:

S:=Snk(n,k):
add(seq(add(op(s)),s in op(S)))/nops(S):

end:

# 5.

#### I decided to do this one with the "words" being integers rather than a list of
#### 1's and 2's. If I want the list form, I can use WORDS([a,b]), etc.

#Wnk(a,b): input non-negative integers a,b, and outputs, the set of all words in {1,2}
#with a 1s, b 2s
#For example
#Wnk(1,2) should be {122,212,221}
Wnk:=proc(a,b) local S,A,B,i:
option remember:
if not (type(a,integer) and type(a,integer) and a>=0 and b>=0) then
 RETURN(FAIL):
fi:

if a+b=1 then
 RETURN({1+b}):
fi:

S:={}:
if a>0 then
  A:=Wnk(a-1,b):
  S:={seq(10^(a+b-1)+A[i],i=1..nops(A))}:
fi:
if b>0 then
  B:=Wnk(a,b-1):
  S:=S union {seq(2*10^(a+b-1)+B[i],i=1..nops(B))}:
fi:
S:

end:

#NuWnk(a,b): input non-negative integers a,b, and outputs, the number of words in {1,2}
#with a 1s, b 2s
#NuWnk(1,2) should be 3
NuWnk:=proc(a,b) local s:
option remember:
if not (type(a,integer) and type(b,integer) and a>=0 and b>=0) then
 RETURN(FAIL):
fi:

if a+b=1 then
 RETURN(1):
fi:

s:=0:
if a>0 then
  s:=NuWnk(a-1,b):
fi:
if b>0 then
  s:=s+NuWnk(a,b-1):
fi:
s:

end:

#RWnk(n,k): input non-negative integers a,b, and outputs a (Uniformly at Random) word in {1,2}
#with a 1s, b 2s
#RWnk(1,2) should output each of 122,212,221 with prob. 1/3
RWnk:=proc(a,b) local i,c:
if not (type(a,integer) and type(b,integer) and a>=0 and b>=0) then
  RETURN(FAIL):
fi:

if a+b=0 then
  RETURN(0):
fi:

if a=0 then
  RETURN(2*10^(b-1)+RWnk(0,b-1)):
fi:
if b=0 then
 RETURN(10^(a-1)+RWnk(a-1,0)):
fi:

c:=LD([NuWnk(a-1,b),NuWnk(a,b-1)]):
if c=1 then
 RETURN(10^(a+b-1)+RWnk(a-1,b)):
else
 RETURN(2*10^(a+b-1)+RWnk(a,b-1)):
fi:

end:

# 6.
#Once I decided to do number 7, I didn't really need to code this one.
W3nk:=proc(a,b,c): RETURN(WORDS([a,b,c])): end:
NuW3nk:=proc(a,b,c): RETURN(NuWORDS([a,b,c])): end:
RW3nk:=proc(a,b,c): RETURN(RandWORD([a,b,c])): end:

# 7.

#### Procedures: WORDS(L), NuWORDS(L), RandWORD(L)
# input a list L, of length k, say, (i.e. nops(L)=k) of non-negative numbers, and
# output, respectively, the set of all words in the alphabet {1,2,...k} with 
# exactly L[1] 1s, L[2] 2s, ..., L[k] k-s, their number, and a random such word.

#WORDS(L)
WORDS:=proc(L) local S,T,i,w:

if add(op(L))=0 then
  RETURN([[]]):
fi:

S:=([]):
for i from 1 to nops(L) do
  if L[i]>0 then
    T:=L:
    T[i]:=T[i]-1:
    S:=[op(S),seq([i,op(w)],w in WORDS(T))]:
  fi:
od:
S:
end:

#NuWORDS(L)
NuWORDS:=proc(L) local s,T,i:

if add(op(L))=1 then
  RETURN(1):
fi:

s:=0:
for i from 1 to nops(L) do
  if L[i]>0 then
    T:=L:
    T[i]:=T[i]-1:
    s:=s+NuWORDS(T):
  fi:
od:
s:
end:

#RandWORD(L)
RandWORD:=proc(L) local i,c,t,T,P:

if add(op(L))=0 then
  RETURN([]):
fi:

#Note: LD(L) works just fine if some of the entries are 0
P:=[]:
for i from 1 to nops(L) do
  if L[i]>0 then
    t:=i: #Saving this in case L consists of a single 1 with the rest 0's
    T:=L:
    T[i]:=T[i]-1:
    P:=[op(P),NuWORDS(T)]:
  else
    P:=[op(P),0]:
  fi:
od:

if add(op(P))=0 then RETURN([t]): fi:
c:=LD(P):
T:=L:
T[c]:=T[c]-1:

[c,op(RandWORD(T))]:

end:


#### Code included from C4.txt ####
#Snk(n,k): Inputs NON-NEG. integer n and integer k and outputs
#the set of k-element subsets of {1,...,n}
#For example
#Snk(3,2) should be {{1,2},{1,3},{2,3}}
Snk:=proc(n,k) local S1,S2:
option remember:
if not (type(n,integer) and type(k,integer) and n>=0) then
 RETURN(FAIL):
fi:
if k<0 or k>n then
 RETURN({}):
fi:

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

S1:=Snk(n-1,k):
S2:=Snk(n-1,k-1):
S1 union {seq(s union {n}, s in S2)}:

end:

#NuSnk(n,k): Inputs NON-NEG. integer n and integer k and outputs
#the number of k-element subsets of {1,...,n}
#For example
#Snk(3,2) should be 3
NuSnk:=proc(n,k) option remember:
if not (type(n,integer) and type(k,integer) and n>=0) then
 RETURN(FAIL):
fi:
if k<0 or k>n then
 RETURN(0):
fi:

if n=0 then
 RETURN(1):
fi:

NuSnk(n-1,k)+NuSnk(n-1,k-1):

end:

#LD(L): Inputs a list of positive integers L (of n:-nops(L) members)
#outputs an integer i from 1 to n with the probability of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6, 2 with prob. 1/3, 3 with prob. 1/2
LD:=proc(L) local n,i,r,su:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:
for i from 1 to n do
  su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:

end:


#Rnk(n,k): Inputs NON-NEG. integer n and integer k and outputs
#the a (Uniformly at Random) k-element subset of {1,...,n}
#Rnk(3,2) should output each of {1,2},{1,3},{2,3} with prob. 1/3
Rnk:=proc(n,k) local c:
if not (type(n,integer) and type(k,integer) and n>=0) then
 RETURN(FAIL):
fi:
if k<0 or k>n then
 RETURN(0):
fi:

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

c:=LD([NuSnk(n-1,k),NuSnk(n-1,k-1)]):
if c=1 then
 RETURN(Rnk(n-1,k)):
else
 RETURN(Rnk(n-1,k-1) union {n}):
fi:

end:
#### End of code included from C4.txt ####