# Please do not post homework
# Lucy Martinez, 02-07-2025, Assignment 5

with(combinat):

######################

# Problem 1: Write a procedure RandomState(n,K) that inputs
#  a pos. integer n and a pos. integer K and outputs a random state vector
#  with n complex components whose real and imaginary parts are from -K to K,
#  then normalize it.
RandomState:=proc(n,K) local ra,v,i,k:
ra:=rand(-K..K):
v:=[seq(ra()+I*ra(),i=1..n)]:
k:=sqrt(add(abs(v[i])^2,i=1..n)):
[seq(v[i]/k,i=1..n)]:
end:


# Problem 2: Experiment with three different random random Hermitian
#  matrices using RHM(n,K) with n=3 and K=10 and check that
#  the eigenvalues are real and that the eigenvectors (pure states)
#  are all orhogonal to each other.



# Problem 3: Experiment with five different random random Hermitian
#  matrices using RHM(n,K) with n=3 and K=10 and corresponding
#  five random state vector and apply ProbDist(L,A) to them.
#  For each case check that the proba. distribution consists
#  of non-negative numbers that add up to 1.

# Here is what I get:

#H11:=RHM(3,10):
#h11:=RandomState(3,10):
#P11:=ProbDist(H11,h11)[1]:
#evalf(add(P11[i],i=1..nops(P11))); outputs 0.9999999980

#H12:=RHM(3,10):
#h12:=RandomState(3,10):
#P12:=ProbDist(H12,h12)[1]:
#evalf(add(P12[i],i=1..nops(P12))); outputs 1.000000000

#H13:=RHM(3,10):
#h13:=RandomState(3,10):
#P13:=ProbDist(H13,h13)[1]:
#evalf(add(P13[i],i=1..nops(P13))); outputs 1.000000000


# Problem 4: Write a procedure GeneralPauli(n)
#  that inputs a NORMALIZED vector n=[nx,ny,nz]
#  that constructs the 2 by 2 matrix. at the bottom of p. 349
#  of the Susskind-Frieman Quantum Mechanics book.
#  Find its eigenvalues and eigenvectors.
GeneralPauli:=proc(n):
[[n[3],n[1]-I*n[2]],[n[1]+I*n[2],-n[3]]]:
end:

# Answer: The eigenvalues are the following
#  {sqrt(nx^2 + ny^2 + nz^2), -sqrt(nx^2 + ny^2 + nz^2)}
# The eigenvectors are the following
#  {[(sqrt(nx^2 + ny^2 + nz^2) + nz)/((nx + ny*I)*sqrt(1 + abs((sqrt(nx^2 + ny^2 + nz^2) + nz)/(nx + ny*I))^2)), 1/sqrt(1 + abs((sqrt(nx^2 + ny^2 + nz^2) + nz)/(nx + ny*I))^2)],
#    [-(sqrt(nx^2 + ny^2 + nz^2) - nz)/((nx + ny*I)*sqrt(1 + abs((-sqrt(nx^2 + ny^2 + nz^2) + nz)/(nx + ny*I))^2)), 1/sqrt(1 + abs((-sqrt(nx^2 + ny^2 + nz^2) + nz)/(nx + ny*I))^2)] }



######################################Previous classes:
#C5.txt, Feb. 6, 2025 Experimental Mathematics (Rutgers, Dr. Z.)
Help:=proc(): print(` PM(), RHM(n,K) , EV(L), EVC1(L,g), EVC(L), IP(A,B), ProbDist(L,A) `): end:

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:

#RHM(n,K): a random n by n Hermitian matrix with entriries from -K to K
RHM:=proc(n,K) local ra,i,j,L:
ra:=rand(-K..K):
for i from 1 to n do
 for j from 1 to i-1 do
  L[i,j]:=ra()+I*ra():
  L[j,i]:=conjugate(L[i,j]):
 od:
L[i,i]:=ra():
od:
[seq([seq(L[i,j],j=1..n)],i=1..n)]:
end:

#EV(L): inputs a Hermitian matrix (say n by n) and outputs the LIST of eigenvalues
EV:=proc(L) local  x,L1,i,n:
n:=nops(L):
L1:=[seq(L[i]-[0$(i-1),x,0$(n-i)],i=1..n)]:
[solve(det(L1),x)]:
end:

#EVC1(L,g): inputs a matrix L and an eigenvalue g finds the NORMALIZED eigenvector
EVC1:=proc(L,g) local v,i,j,var,eq,V,n,L1,V1,va,k:
n:=nops(L):
L1:=[seq(L[i]-[0$(i-1),g,0$(n-i)],i=1..n)]:
if simplify(det(L1))<>0 then
  RETURN(FAIL):
fi:

V:=[seq(v[i],i=1..n)]:
var:={seq(v[i],i=1..n)}:
eq:={seq(add(L1[i,j]*V[j],j=1..n),i=1..n)}:
var:=solve(eq,var):
#V1 the set of unknowns that is arbitrary
V1:={}:
for va in var do
  if op(1,va)=op(2,va) then
  V1:=V1 union {op(1,va)}:
  fi:
od:
if nops(V1)<>1 then
  RETURN(FAIL):
fi:
V:=subs(var,V):
V:=subs(V1[1]=1,V):
k:=sqrt(add(abs(V[i])^2,i=1..n)):
[seq(V[i]/k,i=1..n)]:
end:

#EVC(L): inputs a Hermitian matrix and outputs the pair ListOfEigenvalues, ListOfEigenvectors
EVC:=proc(L) local n,ev,evc,i:
n:=nops(L):
ev:=EV(L):
evc:=[seq(EVC1(L,ev[i]),i=1..n)]:
[ev,evc]:
end:

#IP(A,B): the inner product of state-vectors A and B
IP:=proc(A,B) local n,i:
n:=nops(A):
if nops(B)<>n then
  RETURN(FAIL):
fi:
simplify(add(A[i]*conjugate(B[i]),i=1..n)):
end:

#ProbDist(L,A): Given an observable represented by the Hermitian matrix L and a state vector A
#outputs first all the eigenvectors (pure states) and the prob. of winding in the respective pure states
#followed by the expected value of A (i.e. if you make the observation L and the state is A the average in the long run
ProbDist:=proc(L,A) local evc,i,A1,k,PD,n:
n:=nops(L):
k:=sqrt(add(abs(A[i])^2,i=1..n)):
A1:=[seq(A[i]/k,i=1..n)]:
evc:=EVC(L):
PD:=[seq(abs(IP(A1,evc[2][i]))^2,i=1..n)]:
[PD,add(PD[i]*evc[1][i],i=1..n)]:
end:



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


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


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


# IsUM(A): is the matrix A unitary? i.e. A^T*A=I (id matrix)
IsUM:=proc(A) local n, P,i,j:
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:


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


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


with(linalg):


# EstimateAverageClique(n,p,k,K): estimates the average number of cliques
#  of size k in a random graph given by RG(n,p)
#  by picking K of them and averaging
EstimateAverageClique:=proc(n,p,k,K) local co,G,i:
co:=0:
for i from 1 to K do
 G:=RG(n,p):
 co:=co+nops(Cliques(G,k)):
od:
[evalf(co/K),evalf(binomial(n,k)*p^binomial(k,2))]:
end:


# STbf(G): 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 Tree Thm)
NST:=proc(G) local D1,A,n,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:


########################################################

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

# Neis(G): inputs a graph G=[n,E] and outputs a list N of lenght n such that
#  N[i] is the set of neighbors of vertex i
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:



#CC(G,i): The connected component of vertex i
CC:=proc(G,i) local n,N,C1,C2,C3,i1:

if not (type(G,list) and nops(G)=2 and type(G[1],posint) and type(G[2],set)) then
 return(FAIL):
fi:

n:=G[1]:
N:=Neis(G):
#C1 is the current friend
C1:={i}:
#C2 is the current friends
C2:=C1 union {seq(op(N[i1]),i1 in C1)}:

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



#IsCo(G): is the graph G=[n,E] connected?
IsCo:=proc(G) local n,i,V:
n:=G[1]:
V:={seq(i,i=1..n)}:
evalb(CC(G,1)=V):
end:



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


# NuCGc(n): the first n terms of the sequence enumerating
#  CONNECTED labeled graphs on n vertices, in other words
#  sequence OEIS A001187 , done cleverly,
#  using Exponential Generting Function, the fact that
#  Sum(CC[i]/i!*z^i,i=0..infinity)= log(Sum(2^binomial(i,2)/i!*z^i,i=0..infinity)
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:


# AM(G): inputs a graph [n,E] and outputs the adjacency matrix,
#  represented as a list of length of lists
#  of length n, such that M[i][j]=1 if {i,j} belongs to E and 0 otherwise
AM:=proc(G) local L,e,i:
L:=[]:
for i from 1 to G[1] do
 L:=[op(L),[seq(0,i=1..G[1])]]:
od:

#iterate through every edge
for e in G[2] do
  L[e[1],e[2]]:=1:
  L[e[2],e[1]]:=1:
od:
L:
end:


########################################################

#C2.txt: Jan. 27, 2025
Help2:=proc():
print(`LC(p), RG(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: