#OK to post homework
#Blair Seidler, 2021-04-11 Assignment 20

with(combinat):

Help:=proc(): print(`AvHD(S), AvAv(Dim,J,K), SimuHao(Dim,J,K)`): end:

# 1.

#AvHD(S): inputs a set of vectors S and outputs the average number of neighbors a
#vertex has, rather than the maximum
AvHD:=proc(S) local v:
if S={} then RETURN(FAIL): fi:
evalf(add(seq(Nei(S,v),v in S))/nops(S)):
end:

# 2.
#AvAv(Dim,J,K): inputs Dim, J (like in Hao(Dim,J)) and a LARGE integer K, and by using 
#combinat[randcomb](Box(Dim),J) picks K random subsets, for each computes AvHD(S), and
#averages over all of them.
AvAv:=proc(Dim,J,K) local B:
B:=Box(Dim):
evalf(add(seq(AvHD(randcomb(B,J)),i=1..K))/K):
end:

#What are your estimates for (run it serval times and see whether you get close answers
#to each other) AvAv([2,2,2,2,2,2,2],65,10000)

#I ran three trials, and these are close enough that I'm sold.
#                          3.527384539
#                          3.529172197
#                          3.528547577

# 3.

#I'm assuming that in the notation we were using, you meant AvAv([2$n],J), because [1$n] only has
#one element and is therefore boring.

#I think this is just a linearity of expectation problem. If I am some vertex in a J-subset
#(and we are all indistinguishable by symmetry), let X_i be the indicator function for
#"v_i is my neighbor". I have exactly n neighbors, and there are 2^n-1 vertices which are not me,
#so the expected value of X_i is just n/(2^n-1). There are J-1 elements in the subset which
#aren't me, so my expected number of neighbors is just (J-1)*n/(2^n-1).

#For AvAv([2,2,2,2,2,2,2],65,10000) from problem 2, this would give (65-1)*7/(128-1) which is
#about 3.52756. This is extremely close to the values for the three trials I ran for Problem 2.

# 4.

#SimuHao(Dim,J,K): Inputs a list Dim and an integer J between
#1 and nops(Box(Dim)) and outputs
#The champion and the record in the competition
#which of K randomly generated subsets of Box(Dim) with J elements has the minimum degree
SimuHao:=proc(Dim,J,K) local B,S,cha,rec,i,cand,hope1:
B:=Box(Dim):
cha:={randcomb(B,J)}: rec:=MaxHD(op(cha)):
for i from 2 to K do
  cand:=randcomb(B,J): hope1:=MaxHD(cand):
  if hope1<rec then
    rec:=hope1:
    cha:={cand}:
#  elif hope1=rec then
### The output was too long for my computer if I tried to save all of the "champions"
#    cha:=cha union {cand}:
  fi:
od:
[rec,cha]:
end:

# I ran SimuHao([2$10],513,10000) twice, and it took several hours to run each time.
# The "champion" on both runs was only an 8, rather than the theoretical value of 4.

# It is not terribly surprising that this experimental result is off. The total number
# of subsets is around 2*10^299 and even on two runs I only examined 2*10^4 of them.

#### From C20.txt ####

Help20:=proc(): print(`Box(Dim), HamD(u,v)`): end:

#Box(Dim): inputs a list of positive integers DIM
#{1,2,3,...,Dim[1]}x...x{1,2,...,Dim[nops(Dim)]}
Box:=proc(Dim) local d, Dim1, S, v, i:
option remember:
d:=nops(Dim):
if d=0 then RETURN({[]}): fi:
Dim1:=[op(1..d-1,Dim)]:
S:=Box(Dim1):
{seq(seq([op(v),i],i=1..Dim[d]),v in S)}
end:

#HamD(u,v): The Hamming distance between vector u and v
HamD:=proc(u,v) local i,co:
co:=0:
for i from 1 to nops(u) do
  if u[i]<>v[i] then co:=co+1: fi:
od:
co:
end:

#Nei(S,v): inputs a set of vectors S and a vector v
#outputs number of members of S whose Hamming distance to v is 1
Nei:=proc(S,v) local co,w:
co:=0:
for w in S do
  if HamD(v,w)=1 then co:=co+1: fi:
od:
co:
end:

#MinHD(S): inputs a set of vectors S and outputs the
#minimal number of neighbors a member of S can have
MinHD:=proc(S) local v,m:
if S={} then RETURN(FAIL): fi:
min(seq(Nei(S,v),v in S)):
end:

#MaxHD(S): inputs a set of vectors S and outputs the
#maximal number of neighbors a member of S can have
MaxHD:=proc(S) local v,m:
if S={} then RETURN(FAIL): fi:
max(seq(Nei(S,v),v in S)):
end:


#Hao(Dim,J): Inputs a list Dim and an integer J between
#1 and nops(Box(Dim)) and outputs
#The champion and the record in the competition
#which subset of Box(Dim) with J elements has the minimum degree
Hao:=proc(Dim,J) local B,S,cha,rec,i,cand,hope1:
B:=Box(Dim):
S:=choose(B,J):
cha:={S[1]}: rec:=MaxHD(S[1]):
for i from 2 to nops(S) do
  cand:=S[i]: hope1:=MaxHD(cand):
  if hope1<rec then
    rec:=hope1:
    cha:={cand}:
  elif hope1=rec then
    cha:=cha union {cand}:
  fi:
od:
[rec,cha]:
end: