# hw4 Pablo Blanco
# OK to post

with(linalg):
#### implements Box 1.1 in the book
UM:=proc(alph,bet,gam,del) local A,B,C:
 A:=[[exp(-I*bet/2) ,0], [0,exp(I*bet/2)]]:
 B:=[[cos(gam/2),-sin(gam/2)], [sin(gam/2),cos(gam/2)]]:
 C:=[[exp(-I*del/2), 0],[0,exp(I*del/2)]]:
 #A,B,C:
 convert(multiply([[exp(I*alph),0],[0,exp(I*alph)]  ] ,A,B,C), list): 
end:

#### 



#################
# Testing:
# ra:=rand(0.0..100000.0):
#a:=ra();
#b:=ra();
#g:=ra();
#d:=ra();
#A:=UM(a,b,g,d);
#IsUM(A):
#
# Result: generally false.
#
#################


### Challenge:
## Claim: K_n has n^(n-2) labeled spanning trees. (prove using Matrix Tree Theorem).
## 
## Proof idea: Perform an LU decomposition on the n x n matrix M with n in the diagonal and -1 everywhere else. 
## This is a minor of the Laplacian of K_{n+1}. Then our goal is to prove the determinant of this matrix is (n+1)^(n-1).
## We decompose M=LU. Into L lower triangular matrix and U upper triangular.
## 
## L will have 1s in the diagonal, so it has det(L)=1. Note that det(M) = det(L)*det(U).
## Each diagonal entry d[k] of U is given by d[k] = (n+1)*(n-k+1)/(n-k+2). 
## This telescopes to give us that, the product of the diagonal entries of U, det(U)=(n+1)^(n-1).
## Thus, det(M)=det(U)=(n+1)^(n-1) and we are done.
##
## I explicitly construct M,L,U in the procedure below.

# outputs the LU decomposition of a MINOR of the Laplacian matrix of K_(n+1).
# In particular, outputs M,L,U.


with(LinearAlgebra):

# outputs the LU decomposition of a MINOR of the Laplacian matrix of K_(n+1).
# In particular, outputs M,L,U.
CayleyFrm:=proc(n) local M,L,U,i,j:
 M:=[[-1$n]$n]:
 L:=[[0$n]$n]:
 U:=[[0$n]$n]:
 
 # we use i to range over rows and j to range over columns.

 ### construct M ###
 for i from 1 to n do: 
  M[i][i]:=n:           # define all diagonal entries
 od:
 
 ### construct L ###
 for i from 1 to n do: 
  L[i][i]:=1:           # define all diagonal entries
 od: 
 # define lower triangular section. Note that the columns are identical.
 for j from 1 to n-1 do:
  for i from j+1 to n do:
   L[i][j]:= -1/(n-j+1):
  od: 
 od:  
 
 ### construct U ###
 for i from 1 to n do:
  U[i][i]:= (n+1)*(n-i+1)/(n-i+2):
 od:
 for i from 1 to n do:
  for j from i+1 to n do:
   U[i][j]:= -(n+1)/(n+2-i):
  od:
 od:

 M,L,U:
end:

## verifies that the previous procedure generates a correct LU decomposition. Assuming M is defined correctly.
## Also checks the determinant of U to make sure it matches Cayley's Formula.
## Outputs true or false
CheckProc:=proc(n) local S,M,L,U:
 S:=CayleyFrm(n):
 M:=Matrix(S[1]):
 L:=Matrix(S[2]):
 U:=Matrix(S[3]):
 Equal(M, Multiply(L,U)) and evalb(Determinant(U)=(n+1)^(n-1)):
end:





####################################################################################################################### End Pablo Code

#Feb. 3, 2025
Help4q:=proc(): print(` RM(n,K), CT(A), IsUM(A) ,PM()`): end:
with(combinat):
with(linalg):

#PM(): the three famous Pauli matrices sigma_x,sigma_y, sigma_z
PM:=proc():
[
[[0,1],[1,0]],
[[0,-I],[I,0]],
[[1,0],[0,-1]]
]:
end:

#RM(n,K): a random n by n matrix with complex entries from -K to K
RM:=proc(n,K) local ra,i,j:
ra:=rand(-K..K):
[seq([seq( ra()+I*ra(),j=1..n)],i=1..n)]:
end:

#CT(A): the conjugate transpose of A
CT:=proc(A): local n,i,j:
n:=nops(A):
[seq([seq(conjugate(A[j][i]),j=1..n)],i=1..n)]:
end:

#IsUM(A): Is the matrix A unitary
IsUM:=proc(A) local n,i,j, P:
n:=nops(A):
P:=multiply(A,CT(A)):
evalb([seq([seq(P[i,j],j=1..n)],i=1..n)]=[seq([0$(i-1),1,0$(n-i)],i=1..n)]):
end:

Help4:=proc(): print(` EstimateAverageClique(n,p,k,K) , STbf(G) `):end:

#STbf: Inputs a graph G=[n,E] and outputs the set of all spanning trees
#
STbf:=proc(G) local n,E,S,ST,g:
n:=G[1]:
E:=G[2]:
S:=choose(E,n-1):
ST:={}:
for g in S do
 if IsCo([n,g]) then
  ST:=ST union {[n,g]}:
 fi:
od:
ST:
end:

#NST(G): The number of spanning trees of G (using the Matrix Three Therorem
NST:=proc(G) local n,D1,A,i:
n:=G[1]:
A:=AM(G):
D1:=[seq(add(A[i]),i=1..n)]:
D1:=[seq([0$(i-1),D1[i],0$(n-i)], i=1..n)]:
A:=D1-A:
A:=[ seq([op(1..n-1,A[i])]   , i=1..n-1)]:
det(A):
end:

#code by Joseph K.
EstimateAverageClique := proc(n,p,k,K) local i:
  evalf([add(nops(Cliques(RG(n,p),k)), i=1..K) / K, binomial(n,k)*p^binomial(k,2)]):
end:




#old stuff
#C3.txt: Jan.30, 2025
Help3:=proc(): print(`AM(G), Neis(G), IsCo(G), NuCG(n), NuCGc(n) `): end:



#Neis(G): inputs a graph G=[n,E] and outputs a list of length n, N, such that
#N[i] is the set of neighbors of vertexi
Neis:=proc(G) local n,E,N,i,e:
n:=G[1]: 
E:=G[2]:
for i from 1 to n do
 N[i]:={}:
od:
for e in E do
N[e[1]]:=N[e[1]] union {e[2]}:
N[e[2]]:=N[e[2]] union {e[1]}:
od:
[seq(N[i],i=1..n)]:
end:

#C(G,i): The connected component of vertex i
CC:=proc(G,i) local n,C1,C2,C3,N,i1:
if not (type(G,list) and nops(G)=2 and type(G[1],integer) and G[1]>=0 and type(G[2],set)) then
  RETURN(FAIL):
fi:
n:=G[1]:
N:=Neis(G):

C1:={i}:
C2:=C1 union { seq(op(N[i1]), i1 in C1)}:
while C1<>C2 do
 C3:=C2 union { seq(op(N[i1]), i1 in C2)}:
C1:=C2:
C2:=C3:
od:
C2:
end:

#IsCo(G): Is the graph G connected?
IsCo:=proc(G): evalb(nops(CC(G,1))=G[1]):end:

#NuCG(n): the first n terms of the sequence enumerating CONNECTED labeled graphs on n vertices
#OEIS A001187  
NuCG:=proc(n) local i,g:
[seq(coeff(add(IsCo(g), g in Graphs(i)),true,1),i=1..n)]:
end:

NuCGc:=proc(n) local f,z,i:
#The number of labeled graphs on n vertices is 2^binomial(n,2)
f:=log(add(2^binomial(i,2)*z^i/i!,i=0..n)):
f:=taylor(f,z=0,n+1):
[seq(i!*coeff(f,z,i),i=1..n)]:
end:








## adjacency matrices Code by Aurora Hiveley
# AM(G): inputs a graph [n,E] and outputs the adjacency matrix, represented as a list of length n of lists of length n, 
# such that M[i][j]=1 if {i,j} belongs to E and 0 otherwise. 
# For example AM([2,{{1,2}}]); should output [[0,1],[1,0]] .

AM := proc(G) local n,E,e,M:
n := G[1]:
E := G[2]:

M := [[0$n]$n]:     # initialize n x n matrix of all 0's

for e in E do
    M[e[1]][e[2]]++:
    M[e[2]][e[1]]++:
od:

M:
end:



#C2.txt: Jan. 27, 2025
Help2:=proc(): print(`LC(p), RG(n,p), Cliques(G,k) `):end:
with(combinat):
#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 pos. 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:

###From C1
#C1.txt: Jan. 23, 2025 Exp Math (Dr. Z.)
Help1:=proc(): print(`Graphs(n), Tri(G) , TotTri(G) `):   end:



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

#If there are n vertices how many (undirected) graphs there

#Graphs(n): inputs a non-neg. integer and outputs the set of ALL
#graphs on {1,...,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:
n:=G[1]:
E:=G[2]:
#S is the set of love triangles
S:={}:
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:
#Tot(G): the total number of love triangles and hate triangles
TotTri:=proc(G)
nops(Tri(G))+nops(Tri(Comp(G))):
end: