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

with(combinat):

#2 
IntegerToVec:=proc(k,n) local i,bin:
bin:=[0$n]:
for i from 0 to n-1 do 
 if k mod 2^(i+1) >= 2^i  then 
   bin[n-i]:=1:
 fi:
od:
bin:
end:

VecToInteger:=proc(v,n) local i,k:
k:=0:
for i from 0 to n-1 do
k:=k+v[n-i]*2^i:
od:
k:
end:





#3
#outputs the alpha in the paper. (choose just one index alpha to give the max, as the paper)
#Important: alpha in the paper is in {0,1,...,2^n - 1}, which is incompatible with Findy, then there is a shift.

FindAlpha:=proc(H,n) max[index](abs~(Findy(n,H)))-1: end:





#4
#inputs n and a subset of {0,1,..,2^n - 1} with size 2^(n-1)+1
Haotest:=proc(n,S) local Vecs,i,alpha,center,nei,shift:
shift:=1+~S:
Vecs:={}:
for i from 1 to nops(S) do
Vecs:=Vecs union {IntegerToVec(S[i],n)}:
od:
alpha:=FindAlpha(shift,n):
center:=IntegerToVec(alpha,n):
nei:=Nei(Vecs,center):
print(`neighbor of alpha is`,nei,`sqrt of n is`,sqrt(n)):
if nei^2>=n then print(`correct`): fi:
end:

#Haotest(3,randcomb({seq(i, i = 0 .. 7)}, 5))
#Haotest(4,randcomb({seq(i, i = 0 .. 15)}, 9))
#Haotest(5,randcomb({seq(i, i = 0 .. 31)}, 17))
#Haotest(6,randcomb({seq(i, i = 0 .. 63)}, 33))
#Haotest(7,randcomb({seq(i, i = 0 .. 127)}, 65))
#Haotest(8,randcomb({seq(i, i = 0 .. 255)}, 129))
#they work perfectly




























#C21.txt: Lecture 21
Help21:=proc(): print(`In(n) , An(n), Check1(n), Bn(n), Check2(n), Findx(n,H), Findy(n,H) `):end:

with(linalg):
with(combinat):

#Glue(A): inputs a matrix of matrices given in block
#notation, "glues" them into one matrix
#For example
#Glue([ [[1,2],[1,3]], [[3,5],[1,13]]]]);
Glue:=proc(A) local i,j:
#TO BE CONTINUED IN HOMEWORK

end:

#In(n):n by n identity matrix
In:=proc(n) local L,i:

[seq([0$(i-1),1,0$(n-i)],i=1..n)]:

end:

#An(n): The matrix A_n in Knuth's note about Hunag's proof
An:=proc(n) local L1,L2:
if n=0 then
 RETURN([[0]]):
fi:

L1:=An(n-1):
L2:=In(2^(n-1)):

[
seq([op(L1[i]),op(L2[i])],i=1..nops(L1)),
seq([op(L2[i]),op(-L1[i])],i=1..nops(L1))
]:
end:

#Check1(n): checks that An(n)^2=n*In(2^n) (in Line 1 of the proof)
Check1:=proc(n):
evalb(convert(multiply(An(n),An(n)),list)=n*In(2^n)):
end:

#Bn(n): The 2^n by 2^(n-1) matrix B_n defined in the second paragraph of
#Knuth's note on Hao Huang's proof
Bn:=proc(n) local L1,L2:
L1:=An(n-1):
L2:=In(2^(n-1)):
[seq(L1[i]+sqrt(n)*~L2[i],i=1..nops(L1)), op(L2)]:
end:

#Check2(n):checks that A_n B_n = sqrt(n) B_n
Check2:=proc(n) local L,R,A,B:
A:=An(n):
B:=Bn(n):
L:=multiply(A,B):
R:=[seq(sqrt(n)*~B[i],i=1..nops(B))]:
convert(L,list),R:

end:

#Findx(n,H): inputs a pos. integer n and an arbitrary subset, H,
#of {1,...,2^n} of cardinality 2^(n-1)+1 finds the vector
#2^(n-1)-dimensional vector x such that all the rows of Bn(n) that correspond to those
#that do NOT belong to H in Bn(n)x are 0
#For example
#Findx(3,{1,3,5,6,8});
Findx:=proc(n,H) local B,x,i,var,eq,v,Nv:
if nops(H)<>2^(n-1)+1 then
  RETURN(FAIL):
fi:

B:=Bn(n):

var:={seq(x[i],i=1..2^(n-1))}:
eq:={}:

for i from 1 to 2^n do
 if not member(i,H) then
  eq:=eq union {add(B[i][j]*x[j],j=1..nops(B[i]))}:
 fi:
od:

var:=solve(eq,var):
v:=[seq(subs(var,x[i]),i=1..2^(n-1))]:
#get rid of the symbol
v:=subs({seq(x[i]=1,i=1..2^(n-1))},v):

#Nv is the length of v
Nv:=simplify(sqrt(add(v[i]^2,i=1..nops(v)))):
[seq(v[i]/Nv,i=1..nops(v))]:
end:


#Done after class
#Findy(n,H): The 2^n-dimensinal vector B_n x where x=Findx(n,H) (see above
Findy:=proc(n,H)  local x,y,i:
x:=Findx(n,H):
x:=[seq([x[i]],i=1..nops(x))]:

y:=convert(multiply(Bn(n),x),list):
[seq(y[i][1],i=1..nops(y))]:
end:

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

#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 minimimum degree
Hao:=proc(Dim,J) local B,S,cha,rec,cand,hope1,i:
B:=Box(Dim):
S:=choose(B,J):
cha:={}:
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: