# OK to post homework
# Aurora Hiveley, 2/6/25, Assignment 5

Help := proc(): print(`RandomState(n,K), GeneralPauli(n)`): end:

with(linalg):
with(combinat):

## random state vector
# RandomState(n,K):  inputs a pos. integer n and a pos. integer K 
# outputs a random state vector with n complex components whose real and imaginary parts 
# are from -K to K, which is then normalized

RandomState:=proc(n,K) local ra,i,V,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:


## experimenting with eigenvalues & eigenvectors

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

# for i from 1 to 3 do 
#   A := RHM(3,10):
#   eigen := EVC(A):
  
#   for v in eigen[1] do
#     v := simplify(v):
#     if v<>conjugate(v) then 
##      print(v);
##      RETURN(FAIL):   # if an eigenvalue isn't real, returns FAIL
#    fi:
#  od:

#  pairs := choose(eigen[2],2): 
#  for p in pairs do 
#    if IP(p[1],p[2])<>0 then 
#      RETURN(FAIL):    # if a pair of eigenvectors isn't orthogonal, returns FAIL
#    fi:
#  od:
# od:

# some eigenvalues appear to come out imaginary, still. i suspect this is a result of maple's difficulty with simplification? 
# however, orthogonality seems to behave as expected!



## experimenting with prob dist 

# Experiment with five different 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 prob dist consists of non-negative numbers that add up to 1.

# for i from 1 to 5 do 
#   M := RHM(3,10):
#   V := RandomState(3,10):

#   P := ProbDist(M,V):
#   for p in P do 
#     if not p>=0 then 
#       RETURN(FAIL):
#     fi:
#   od:

#   if simplify(add(p,p in P))<>1 then 
#     RETURN(FAIL):
#   fi:
# od:

# there's also some issues here with imaginary numbers, perhaps because of the simplify command? i'm not certain



## general pauli matrices
# GeneralPauli(n): inputs a NORMALIZED vector n=[nx,ny,nz], then constructs the 2x2 matrix
# at the bottom of p. 349 of Susskind-Frieman 

GeneralPauli := proc(n) local i : 
if nops(n)<>3 or add(abs(n[i])^2,i=1..3)<>1 then    # verifies that input vector is normalized and of length 3
  RETURN(FAIL):
fi:

[[n[3], n[1]-I*n[2]], [n[1]+I*n[2], -n[3]]]:
end:


## Find its eigenvalues and eigenvectors.
# example: 
# M := GeneralPauli(RandomState(3,10)):
# EVC(N);








### copied from C5 
#C5.txt, Feb. 6, 2025 Experimental Mathematics (Rutgers, Dr. Z.)
Help5:=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 obervalble 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: