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

with(linalg):
with(combinat):

Help:=proc(): print(`IntegerToVec(k,n), VecToInteger(v,n),FindAlpha(H,n),CheckKnuth(H,n)`): end:

# 2.
#IntegerToVec(k,n): inputs an integer k in [1,2^n] (and n) and outputs the vector of
#length n that corresponds to the binary represenation of k-1 (written as a list)
#(with zeros padded at the front to make it of length n, if necessary)
IntegerToVec:=proc(k,n) local i,rem,L:
rem:=k-1:
L:=[0$n]:
for i from n to 1 by -1 do
  L[i]:=rem mod 2:
  rem:=trunc(rem/2):
od:
L:
end:

#VecToInteger(v,n): inverse of IntegerToVec(k,n)
VecToInteger:=proc(v,n) local i:
sum(2^(n-i)*v[i],i=1..n)+1:
end:

# 3.
#FindAlpha(H,n): inputs a subset of Box([2$n]) and computes the member(s) of H predicted by taking
#the index (indices) with the largest absolute value of Findy(n,H).
FindAlpha:=proc(H,n) local S,L,i,y,m:
S:={seq(VecToInteger(H[i],n),i=1..nops(H))}:
y:=Findy(n,S):
m:=max(seq(abs(y[i]),i=1..nops(y))):
L:=[]:
for i from 1 to nops(y) do
  if abs(y[i])=m then
    L:=[op(L),IntegerToVec(i,n)]:
  fi:
od:
L:
end:

#Check that its number of neighbors is indeed the number of neighbors in H (number of other vectors
#in H whose Hamming distance with it is 1). In other words, the procedure checks the last step in
#Knuth's version of Hao Huang's proof.
CheckKnuth:=proc(H,n) local L,c,nb,i:
L:=FindAlpha(H,n):
c:=ceil(sqrt(n)):
#print(`ceil(sqrt(n)) is `,c):
for i from 1 to nops(L) do
  nb:=Nei(H,L[i]):
#  print(L[i],` has `, nb ,` neighbors`):
  if nb<c then RETURN(false): fi:
od:
return(true):
end:


# 4. I ran my CheckKnuth with a bunch of random subsets of [2$8] and everything seems to work

# Note that I wrote a new version of Box called Box0(Dim) which indexes each dimension of the
# box starting at 0 instead of 1. Wow I miss programming in C where my array indices always
# started at 0 like they should. Yes, I believe that 0 is a natural number.




#### From C20.txt ####
#Box0(Dim): inputs a list of pos. integers DIM
#and outputs the set of vectors
#{0,1,2,3.., Dim[1]} x....x {0,1,2,... Dim[nops(Dim)-1]}
Box0:=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:=Box0(Dim1):

{seq(seq([op(v) , i  ], i=0..Dim[d]-1), 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:



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


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