#Homework by Mike Murr
#OK to post
#hw20
read(`C20.txt`);

#Box([2,2,2,2]);

#HamD([1,2,1,2],[1,2,2,2]);

#S:=Box([2,2,2,2]);

#Nei(S,[1,2,2,2]);

#MaxHD(S);

#AvHD(S): inputs a set of vectors S and outputs the
#average number of neighbors a vertex has.
#
# We are assuming that each member of S represents a vertex
# of a hypercube.
#
AvHD:=proc(S) local v, dimHypercube:
if S={} then
 RETURN(FAIL):
fi:
dimHypercube:=nops(S[1]);
#print(dimHypercube);
add(seq(Nei(S,v), v in S))/(2^dimHypercube);
end;

S:=Box([2,2,2,2]);

AvHD(S);

#Hao([2,2,2,2],3);

#Hao([2,2,2,2],4);

#Hao([2,2,2,2],7);

#Hao([2,2,2,2],8);

#Hao([2,2,2,2],9);

with(combinat);

#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, S, i, oneAvHD, sumAvHD;

B:=Box(Dim);

sumAvHD:=0;

for i from 1 to K do

  S:=randcomb(B,J);
#  print(S);
  oneAvHD:=AvHD(S);
#  if oneAvHD > 1 then
#    print(S, oneAvHD);
#  fi;
  sumAvHD:=sumAvHD + oneAvHD;

od;

sumAvHD/K;

end proc;

#AvAv([2,2,2],5,10);

#evalf(AvAv([2,2],3,10000));

#evalf(AvAv([2,2,2],5,10000));

#evalf(AvAv([2,2,2,2],9,10000));

#evalf(AvAv([2,2,2,2,2],17,10000));

#evalf(AvAv([2,2,2,2,2,2],33,10000));

#evalf(AvAv([2,2,2,2,2,2,2],65,10000));

# Result: 1.7918…
# Result: 1.792150000
# Result: 1.7913078…

#SimuHao(Dim,J,K): inputs a list Dim, and an integer J and a large
# positive integer K, uses combinat[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 B,S,cha,rec,cand,hope1,i:
B:=Box(Dim):
#S:=choose(B,J):
cha:={}:
cand:=randcomb(B,J);
cha:={cand}: 
rec:=MaxHD(cand):

for i from 2 to K do
cand:=randcomb(B,J);
hope1:=MaxHD(cand):
if hope1<rec then
 rec:=hope1:
 cha:={cand}:
# The next two lines can be restored if desired.
# elif hope1=rec then
#  cha:=cha union {cand}:
fi:
od:
[rec,cha]:
end proc;

#SimuHao([2,2,2],5,10);
# Result = 2

#SimuHao([2,2,2,2],9,10);
# Result = 3, repeated 6 times. (Expected result = 2)

#SimuHao([2,2,2,2],9,1000);
# Result = 2, repeated 2 times.

#SimuHao([2$5],17,1000);
# Result = 3, repeated many times.

#SimuHao([2$6],33,1000);
# Result = 4, repeated many times.

#SimuHao([2$6],33,10000);
# Result = 4, repeated many, many times. (Ceiling(sqrt(6)) = 3)

#SimuHao([2$7],65,10000);
# Result = 5, repeated many, many times. (Ceiling(sqrt(7)) = 3)

#SimuHao([2$8],129,10000);
# Result = 6 (Ceiling(sqrt(8)) = 3)

SimuHao([2$9],257,10000);
# Result = 7 (Celing(sqrt(9)) = 3)