# Ok to post homework
# Lucy Martinez, 03-05-2025, Assignment 12

with(combinat):

######################
#Problem 1: Write a procedure AveLuckyVertexFriends(n,K)
# that inputs a positive integer n, and for K random subsets of {1, ..., 2^n}
# with 2^(n-1)+1 elements gotten with combinat[randcomb](2^n,2^(n-1)+1)
# does the following
# 1. Finds the lucky vertex (in terms of the row number)
# 2. Finds the number of members of H that are joined by an edge to it
#    [Hint: The absolute value of An(n) is the adjacency matrix
#     of the n-dim unit cube so you have to sum the absolute value of
#     the entry in the row that correpond to the lucky vertex that belong
#     to columns of H]. 
#     Call that number FriendsOfLucky(n,H)
# 3. Make sure that it is always ≥ sqrt(n),
#    according to the sensitivity lemma. If you can find an H that violates it,
#    return FAIL.
# 4. Average it over all K random chices and output the average.

# What did you get for AveLuckyVertexFriends(9,100)?
# Answer: I get 9


AveLuckyVertexFriends:=proc(n,K) local co,i,N,H:
co:=0:
for i from 1 to K do
 H:=randcomb(2^n,2^(n-1)+1):
 N:=FriendsOfLucky(n,H):
 if N<evalf(sqrt(n)) then
  return(FAIL):
 else
  co:=co+N:
 fi:
od:
co/K:
end:


#FriendsOfLucky(n,H): Inputs an integer n and a set H.
# It first finds the lucky vertex (in terms of the row number) and
# then finds the number of members of H that are joined by an edge to it
# (to the lucky vertex)
# [Hint: The absolute value of An(n) is the adjacency matrix
#  of the n-dim unit cube so you have to sum the absolute value of
#  the entry in the row that correpond to the lucky vertex that belong
#  to columns of H]
FriendsOfLucky:=proc(n,H) local LV,A,i:
LV:=LuckyVertex(n,H):
A:=An(n):
add(abs(A[LV][i]),i=1..nops(A[1])):
end:

####################################previous classes
#C12.txt, March 03, 2025
#Continuing Hao Huang's proof (according to Knuth)
with(combinat):
with(linalg):


Help12:=proc():
print(`Verify1(n),Findy(n,H),LuckyVertex(n,H)`):
end:

#Verify1(n): checks that An(n)*Bn(n)=sqrt(n)*Bn(n)
Verify1:=proc(n):
evalb(Mul(An(n),Bn(n))=expand(sqrt(n)*Bn(n))):
end:


#Findy(n,H): inputs a positive integer n and a subset H with
# exactly 2^(n-1)+1 elements and outputs the vector y such that
# An(n)*y=sqrt(n)*y (y is an eigenvector corresponding to An(n))
Findy:=proc(n,H) local A,B,Hc,i,X,var,j,Xn,Y,var1,i1,v:
A:=An(n):
B:=Bn(n):
if nops(H)<>2^(n-1)+1 then
 print(`The set `, H, ` should have `, 2^(n-1)+1, `members`):
 return(FAIL):
fi:
Hc:={seq(i,i=1..2^n)} minus H:
#X is the vector we are looking for
X:=[seq(x[i],i=1..2^(n-1))]:
var:={seq(x[i],i=1..2^(n-1))}:
#the set of equations to determine what X is:
var:=solve({seq(add(B[i][j]*x[j],j=1..2^(n-1)),i in Hc)},var):
X:=subs(var,X):
#we need to seek for a specific solution since we might get infinitely many
var1:={}:
for v in var do
 if op(1,v)=op(2,v) then
  var1:=var1 union {op(1,v)}:
 fi:
od:
X:=subs({var1[1]=1,seq(var1[i1]=0,i1=2..nops(var1))},X):
Xn:=simplify(sqrt(add(X[i]^2,i=1..2^(n-1))),symbolic):
X:=expand(X/Xn):
Y:=[seq(add(B[i][j]*X[j],j=1..2^(n-1)),i=1..2^n)]:
if {seq(add(A[i][j]*Y[j],j=1..2^n)-sqrt(n)*Y[i],i=1..nops(A))}<>{0} then
 print(`Either Huang messed up or Knuth, or more likely we`):
 return(FAIL):
fi:
Y:
end:


#LuckyVertex(n,H): inputs a positive integer n and a subset H of size 2^(n-1)+1
# and outputs a member that is guaranteed to have at least sqrt(n) neighbors
# in H
LuckyVertex:=proc(n,H) local i,Y,Y1:
Y:=Findy(n,H):
Y1:=max[index]([seq(abs(Y[i]),i=1..nops(Y))]):
end:


##############################################
#C11.txt, Feb. 27, 2025
#The sensitivity conjecture of Hao Huang (according to Knuth)

Help11:=proc():
print(`NeiList(G),MaxNN(G,H),MinMaxNN(G,r),An(n),Bn(n),MinMaxNNpablo(G,r)`):
end:


#NeiList(G): Given a graph G=[n,E] outputs a list of length n-1
# whose i-th entry is the set of neighbors of vertex i
NeiList:=proc(G) local n,E,e,T,i:
n:=G[1]:
E:=G[2]:

for i from 1 to n do
 T[i]:={}:
od:

for e in E do
 T[e[1]]:=T[e[1]] union {e[2]}:
 T[e[2]]:=T[e[2]] union {e[1]}:
od:
[seq(T[i],i=1..n)]:
end:



#MaxNN(G,H): inputs a graph G=[n,E] and a subset H of {1,...,n}, and
# outputs the maximum number of neighbors a member that belong to H over all
# members 
MaxNN:=proc(G,H) local L,i:
L:=NeiList(G):
max(seq(nops(L[i] intersect H),i in H)):
end:


#MinMaxNN(G,r): the minimum of MaxNN(G,H) over all subsets H of {1,...,n}
# with cardinality r
# Note: the Hao Huang famous sensitivity conjecture that says that
# MinMaxNN(Ck(n),2^(n-1)+1)>=sqrt(n) 
MinMaxNN:=proc(G,r) local n,i,S,s:
n:=G[1]:
#all subsets of cardinality r
S:=choose({seq(i,i=1..n)},r):
min(seq(MaxNN(G,s),s in S)):
end:


#MinMaxNNpablo(G,r): minimum of MaxNN(G,H) over all subsets H of [n] with size r.
# more memory-efficient way suggested by Pablo Blanco 
# (after discussions with Aurora Hiveley)
MinMaxNNpablo:=proc(G,r) local n,i,c,mini,curr:
n:=G[1]:
c:=firstcomb(n,r):  
mini:=MaxNN(G,c): #check {1,..,r} first

while nextcomb(c,n)<>FAIL do:
#check all other r-subsets of [n] and update the minimum if they are better
 c:=nextcomb(c,n):
 curr:=MaxNN(G,c):
 if curr < mini then:
  mini:=curr:
 fi:
od:
mini:
end:



#An(n): the 2^n by 2^n matrix (given as a list of lists) in Knuth's write-up
# of Huang's proof
An:=proc(n) local A,i:
option remember:
if n=0 then
 return([[0]]):
fi:
A:=An(n-1):
[seq([op(A[i]),0$(i-1),1,0$(nops(A[i])-i)],i=1..nops(A)),
 seq([0$(i-1),1,0$(nops(A[i])-i), op(-A[i])],i=1..nops(A))
]:
end:


#Bn(n): the 2^(n-1) by 2^n matrix in the paper
Bn:=proc(n) local A,i:
A:=An(n-1):
[seq(A[i]+[0$(i-1),sqrt(n),0$(nops(A[i])-i)],i=1..nops(A)),
 seq([0$(i-1),1,0$(nops(A[i])-i)],i=1..nops(A))
]:
end:


#Mul(A,B): the product of matrix A and B (assuming that it exists)
Mul:=proc(A,B) local i,j,k,n:
[seq([seq(add(A[i][k]*B[k][j],k=1..nops(A[i])),j=1..nops(B[1]))],i=1..nops(A))]:
end:


######################################
#C2.txt: Jan. 27, 2025
Help2:=proc():
print(`LC(p), RG(n,p), Cliques(G,k)`):
end:


#LC(p): inputs a rational number between 0 and 1 and outputs true with prob. p
LC:=proc(p) local a,b,ra:
if not (type(p,fraction) and p>= 0 and p<=1) then
 return fail:
fi:

a:=numer(p):
b:=denom(p):

ra:=rand(1..b)():

if ra<=a then
 true:
else
 false:
fi:
end:


RG:=proc(n,p) local E,i,j:
E:={}:
for i from 1 to n do
 for j from i+1 to n do
  if LC(p) then
   E:=E union {{i,j}}:
  fi:
 od:
od:
[n,E]:
end:


#Cliques(G,k): inputs a graph G and a positive integer k, outputs
# the set of k-cliques
Cliques:=proc(G,k) local n, E,S,i,c,C:
n:=G[1]:
E:=G[2]:
S:={}:
C:=choose({seq(i,i=1..n)},k):
for c in C do
 if choose(c,2) minus E={} then
  S:= S union {c}:
 fi:
od:
S:
end:


###############################
#C1.txt: Jan. 23
Help1:=proc():
print(`Graphs(n),Tri(G),Comp(G),TotTri(G)`):
end:



#An undirected graph is a set of vertices V and a set of edges
#[V,E] and an edge is defined as e={i,j} where i and j belong to V
#Convention: our vertices are labeled {1,2,...,n}
#Our data structure is [n,E] where E is a set of edges
#[3,{{1,2},{1,3},{2,3}}];

#If there are n vertices, how many (undirected) graphs are there?
# answer: 2^(n choose 2) since you either include the vertex or not and
# then you choose which two to connect

#Graphs(n): inputs a non-neg. integer n and outputs the set of ALL
# graphs on {1,2,...,n}
Graphs:=proc(n) local i,j,S,E,s:
E:={seq(seq({i,j},j=i+1..n),i=1..n)}:
S:=powerset(E):
{seq([n,s],s in S)}:
end:


#Tri(G): inputs a graph [n,E] and outputs the set of all triples {i,j,k}
# such {{i,j},{i,k},{j,k}} is a subset of E
Tri:=proc(G) local n, S,E,i,j,k:
S:={}:
n:=G[1]:
E:=G[2]:
for i from 1 to n do
 for j from i+1 to n do
  for k from j+1 to n do
   #if member({i,j},E) and member({i,k},E) and member({j,k},E) then
   if {{i,j},{i,k},{j,k}} minus E={} then
	S:=S union {{i,j,k}}:
   fi:
  od:
 od:
od:
S:
end:

#Comp(G): the complement of G=[n,E]
Comp:=proc(G) local n,i,j,E:
n:=G[1]:
E:=G[2]:
[n,{seq(seq({i,j},j=i+1..n),i=1..n)} minus E]:
end:


#TotTri(G): the total number of love triangles and hate triangles
TotTri:=proc(G):
nops(Tri(G))+nops(Tri(Comp(G))):
end:


##############################
#C9.txt, Feb. 20, 2025
Help9:=proc():
print(`Gnd(n,d),AM(G),lam2(G),Vk(k),Nei(v),Ck(k)`):
end:

#Expanders


#Gnd(n,d): the graph whose vertices are {0,1,...,n-1}
# and edges {i,i+j mod n} where j=1..d. The vertices will have degree 2d
Gnd:=proc(n,d) local i,j,E:
E:={seq(seq({i,i+j mod n},j=1..d),i=0..n-1)}:
[n,subs({seq(i=i+1,i=0..n-1)},E)]:
end:


#AM(G): the adjancency matrix of the graph G (the stupid way)
AM:=proc(G) local n, E,e,i,j,A:
n:=G[1]:
E:=G[2]:
for i from 1 to n do
 for j from 1 to n do
  A[i,j]:=0:
 od:
od:

for e in E do
 A[e[1],e[2]]:=1:
 A[e[2],e[1]]:=1:
od:
[seq([seq(A[i,j],j=1..n)],i=1..n)]:
end:


#lam2(G): inputs a graph and returns FAIL if it is not regular, otherwise
# it returns [degree,Lambda_2] (Lambda_2 is the second largest eigenvalue 
# of the adjancencey matrix)
lam2:=proc(G) local A,x,d,S,i:
A:=AM(G):
S:={seq(add(A[i]),i=1..nops(A))}:
if nops(S)<>1 then
 return(FAIL):
fi:
d:=S[1]:
[d,fsolve(charpoly(A,x))[-2]]:
end:


#Vk(k): The list of all 0-1 vectors of length k in lex order 
# (starting at 00...0 and ends with 11..11
Vk:=proc(k) local S,i:
option remember:
if k=0 then
 return([[]]):
fi:

S:=Vk(k-1):
[seq([0,op(S[i])],i=1..nops(S)),seq([1,op(S[i])],i=1..nops(S))]:
end:


#Nei(v): inputs the 0-1 vector and outputs all the vectors where exactly
# one bit is changed (hamming distance is 1) from the vector v
Nei:=proc(v) local k,i:
k:=nops(v):
if {op(v)} minus {0,1} <> {} then
 return(FAIL):
fi:
{seq([op(1..i-1,v),1-v[i] ,op(i+1..k,v)],i=1..k)}
end:


#Ck(k): the k-dimensional unit cube as a graph in our format
Ck:=proc(k) local V,E,i,v:
V:=Vk(k):
E:={seq(seq({V[i],v},v in Nei(V[i])), i=1..nops(V))}:
E:=subs({seq(V[i]=i,i=1..nops(V))},E):
[nops(V),E]:
end: