#OK to post
#Yuxuan Yang, April 10th, Assignment 20

with(combinat):

#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,m:
if S={} then
 RETURN(FAIL):
fi:
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 i,sum:
sum:=0:
for i from 1 to K do
sum:=sum+AvHD(randcomb(Box(Dim),J)):
od:
sum/K:
end:
#seq(evalf(AvAv([2,2,2,2,2,2,2],65,10000)),i=1..10) ;
#3.529073846, 3.526649231, 3.528236923, 3.528753846, 3.528723077, 3.525726154, 3.525292308, 3.527375385, 3.525821538, 3.525396923


#3
#n*2^(n-1)              *   (binomial(2^n-2,J-2))/(binomial(2^n,J))                  *   2                           /    J
#number of edges in cube*    probability: incident points to be in the set           *   contribution to the deg sum /   number of points
Realav:=proc(n,J) local i:
n*2^(n-1)* (binomial(2^n-2,J-2))/(binomial(2^n,J)) * 2 / J:
end:


#Realav(7,65) gives 448/127 = 3.527559055...



#SimuHao(Dim,J,K):Inputs a list Dim, and an integer J and a large positive integer K, 
#uses combina[randcomb](Box(Dim),J) to generate a random subset of size J, K times, 
#and outputs the "champion" and the "record" among the K sampled sets. 
SimuHao:=proc(Dim,J,K) local i,cha,cand,rec,hope1:

cha:={}:
cand:=randcomb(Box(Dim),J):
cha:={cand}: rec:=MaxHD(cand):

for i from 1 to K do
cand:=randcomb(Box(Dim),J):
hope1:=MaxHD(cand):
if hope1<rec then
 rec:=hope1:
 cha:={cand}:
 elif hope1=rec then
  cha:=cha union {cand}:
fi:
od:
[rec,cha]:
end:

#SimuHao([2$10],513,10000)?
#My Maple does NOT accept this. Maybe there are too many champions to print.










#C20.txt, April 5, 2021
Help20:=proc(): print(`Box(Dim), HamD(u,v) , Nei(S,v), MaxHD(S), Hao(Dim,J)`):end:

#Box(Dim): inputs a list of pos. integers DIM
#and outputs the set of vectors
#{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 vetor 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 the 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(w,v)=1 then
  co:=co+1:
fi:
od:
co:
end:

#MaxHD(S): inputs a set of vectors S and outputs the
#max 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:

with(combinat):