Quantum Random Walk

qrw.py

# import NumPy
from numpy import *
from matplotlib.pyplot import *

# add a line to print out floating point numbers to 3 decimal places

set_printoptions(precision=3) # float precision for printing

N = 100 # define number of random steps
P = 2*N+1 # define the number of positions

# define a coin state c0=|0> and c1=|1> using the matric representation

c0 = array([1, 0]) # |0>
c1 = array([0, 1]) # |1>

Coin00 = outer(c0, c0) # define the Hademard coin operator |0><0|

#print Coin00
#[[1 0]
# [0 0]]

# define Coin11 to be |1><1|

Coin11 = outer(c1, c1)

#print Coin11
#[[0 0]
# [0 1]]


Coin10 = outer(c1, c0)
Coin01 = outer(c0, c1)

# define coin operator C_hat that "toosses/rotates" a quantum coin
C_hat = ( 1/sqrt(2) ) * ( Coin00 + Coin01 + Coin10 - Coin11 )

#print C_hat
#[[ 0.707  0.707]
# [ 0.707 -0.707]]

# shift (step) operator

#print eye(P)
#print roll(eye(P),N)

# define the shift which shifts down the arguments in the matrix by a certain position 
shift_d = roll(eye(P), -1, axis=1)
shift_u = roll(eye(P), 1, axis=1)

#print shift

shift_down = kron( shift_d, Coin00 )
shift_up = kron( shift_u, Coin11 )

#print shift_up

# define full shift operator

S_hat = shift_down + shift_up

#print S_hat

# define the combined walk operator

combo_walk = kron( eye(P), C_hat )

#print combo_walk

#define full walk operator

U_hat = dot( S_hat, combo_walk )

#print U_hat

# define initial state of the coin and walker position as a tensor product

posn0 = zeros(P) 
posn0[N] = 1

tensor = ( 1/(sqrt(2) ) * kron( posn0, ( c0 + ( c1 * 1j ) ) ) )

#print tensor

psi = dot( linalg.matrix_power(U_hat, N)  , tensor )

#print final_state

# apply measurement operator M 
M_hat = zeros( (2*P,2*P,2*P) ) 
for p in range(P): 
 # M[p] = Id2 (x) |p><p| 
 posn = zeros(P) 
 posn[p] = 1 
 M_hat[p] = kron( outer(posn,posn), eye(2)) 
# project, then square 
prob = empty(P) 
for p in range(P): 
 proj = M_hat[p].dot(psi) 
 prob[p] = proj.dot(proj.conjugate()).real 
 # prob is real, by construction

#print prob

fig = figure() 
ax = fig.add_subplot(111) 
 
plot(arange(P), prob) 
#Location of ticks 
loc = range (0, P, P / 10) 
xticks(loc) 
xlim(0, P) 
ax.set_xticklabels(range (-N, N+1, P / 10)) 
 
show()