# Set up the cluster for parallel computation
from IPython.parallel import Client
client = Client(profile='lsf')
len(client.ids)
34
dview = client[:]
dview.use_dill()
lview = client.load_balanced_view();
%%px --local
%matplotlib inline
import qutip
import numpy as np
import scipy
import matplotlib as mpl
import matplotlib.pyplot as plt
from IPython import display
import sys
import time
from itertools import product
from types import FunctionType, BuiltinFunctionType
from functools import partial
def expand_lindblad_operator_to_qutip_std_form(*summands):
"""Expands the superoperator L(A+B+...) to sum of independent dissipator operators.
Warning: only real coefficients are supported.
The input `[matA, coefA], [matB, coefB], ...` represents the Lindblad dissipator
`L(matA*coefA + matB*coefB + ...)`. `mat*` is an operator, `coef*` is a time
dependent coefficient (in any of the qutip supported formats).
If a summand is constant just input it on its own,
e.g. `[matA, coefA], matConst, [matB, coefB]`.
The output is a list of independent dissipator superoperators in form that can be
used directly with `mesolve`'s `c_ops`."""
matrix_coeff_pairs = [_ if isinstance(_, list) else [_, None] for _ in summands]
if all(isinstance(_, np.ndarray) or _ is None for __, _ in matrix_coeff_pairs):
coef_type = 'array'
elif all(isinstance(_, str) or _ is None for __, _ in matrix_coeff_pairs):
coef_type = 'str'
elif all(isinstance(_, (FunctionType, BuiltinFunctionType, partial)) or _ is None for __, _ in matrix_coeff_pairs):
coef_type = 'function'
else:
raise TypeError("Incorrect format for the coefficients")
c_ops = sum((properly_conjugated_lindblad(*_, coef_type=coef_type) for _ in product(matrix_coeff_pairs, repeat=2)),[])
return sorted(c_ops, key=lambda _:isinstance(_, list), reverse=True)
def properly_conjugated_lindblad(left, right, coef_type):
'''Return the two components of the Lindblad superoperator with properly conjugated coeffs.
For input `(A, ca), (B, cb)` return
`[[pre(A)*post(Bdag), ca*conj(cb)],[-0.5*(pre(Adag*B)+post(Adag*B)), conj(ca)*cb]]`.
It supports constant operators (just set `ca` or `cb` to `None`).
'''
if coef_type == 'array':
conj = lambda _: np.conj(_)
prod = lambda a,b: a*b
elif coef_type == 'str':
conj = lambda _: 'conj(%s)'%_
prod = lambda a,b: '(%s)*(%s)'%(a,b)
elif coef_type == 'function':
raise NotImplementedError
m_left, c_left = left
m_right, c_right = right
A = qutip.spre(m_left) * qutip.spost(m_right.dag())
ld_r = m_left.dag()*m_right
B = - 0.5 * qutip.spre(ld_r) - 0.5 * qutip.spost(ld_r)
if c_left is None and c_right is None:
return [A+B]
elif c_left is None:
return [[A, conj(c_right)],
[B, c_right]]
elif c_right is None:
return [[A, c_left],
[B, conj(c_left)]]
else:
return [[A, prod( c_left, conj(c_right))],
[B, prod(conj(c_left), c_right)]]
%%px --local
def do_all_for_given_circle_size(radius):
N = 40
samples = 1600
T = 200
title = 'two_blobs_two_circles_%.3f'%radius
unit_interval = np.linspace(0,1,samples)
centered_interval = np.linspace(-1,1,samples)
one = np.ones(samples)
zero = np.zeros(samples)
gaussian = np.exp(-centered_interval**2)
a_op = qutip.destroy(N)
adag_op = qutip.create(N)
id_op = qutip.identity(N)
zero_op = qutip.zero_oper(N)
num_op = qutip.num(N)
alpha0 = 2 + 0j
beta0 = -2 + 0j
init_state = (qutip.coherent(N,alpha0) + qutip.coherent(N,beta0)).unit()
alphas = alpha0 + radius*(1-np.cos( unit_interval*2*np.pi)) - radius*1j*np.sin(unit_interval*2*np.pi)
betas = beta0 - radius*(1-np.cos(-unit_interval*2*np.pi)) + radius*1j*np.sin(-unit_interval*2*np.pi)
c_ops = expand_lindblad_operator_to_qutip_std_form(
a_op**2,
[-a_op, alphas+betas],
[id_op, alphas*betas]
)
ret = qutip.mesolve(zero_op,init_state,unit_interval*T,
c_ops=c_ops,
e_ops=[],
progress_bar=qutip.ui.progressbar.TextProgressBar())
def animate_states(states, paths=[], every=5, extent=[[-4,4],[-4,4]], title=''):
plt.close('all')
for i, s in enumerate(states[::every]):
f, a = qutip.WignerDistribution(s, extent=np.array(extent)*np.sqrt(2)).visualize()
for p in paths:
pp = p*np.sqrt(2)
a.plot(np.real(pp),np.imag(pp),'gx-')
a.plot(np.real(pp[i*every]),np.imag(pp[i*every]),'ro')
f.savefig(title+'_ani_%06d.png'%(i*every))
print('%d/%d'%(i+1,len(states[::every])))
sys.stdout.flush()
plt.close('all')
def plot_xyz(ret, alphas, betas, title=''):
plt.close('all')
proj_pop_coh = [(
qutip.expect(s,qutip.coherent_dm(N,a)-qutip.coherent_dm(N,b)),
qutip.expect(s,qutip.coherent(N,a)*qutip.coherent(N,b).dag())
)
for s, a, b in zip(ret.states, alphas, betas)]
proj_pop_coh = np.array(proj_pop_coh)
z = proj_pop_coh[:,0]
x = 2*np.real(proj_pop_coh[:,1])
y = 2*np.imag(proj_pop_coh[:,1])
f = plt.figure()
plt.plot(ret.times,np.real(z), label='z')
plt.plot(ret.times,x, label='x')
plt.plot(ret.times,y, label='y')
plt.plot(ret.times,np.sqrt(np.abs(z)**2+x**2+y**2), label='norm(x,y,z)')
plt.xlabel('Time')
plt.title('x y z components of the logical qubit')
plt.legend()
f.savefig(title+'_xyz.png')
return f
def plot_purity(ret, title=''):
plt.close('all')
f = plt.figure()
plt.plot(ret.times,[(s**2).tr() for s in ret.states])
plt.title('Purity')
plt.xlabel('Time')
f.savefig(title+'_purity.png')
return f
animate_states(ret.states, paths=[alphas,betas], every=10, extent=[[-3,3],[-2,2]], title=title)
plot_xyz(ret, alphas, betas, title=title)
plot_purity(ret, title=title)
return ret
radii = np.linspace(-0.5,1.,num=40)
results = lview.map(do_all_for_given_circle_size, radii)
results.wait_interactive()
40/40 tasks finished after 1729 s done
result_vals = results.result
last_states = [r.states[-1] for r in result_vals]
last_states_pop_coh = [(
qutip.expect(s,qutip.coherent_dm(40,2)-qutip.coherent_dm(40,-2)),
qutip.expect(s,qutip.coherent(40,2)*qutip.coherent(40,-2).dag())
)
for s in last_states]
last_states_pop_coh = np.array(last_states_pop_coh)
z = last_states_pop_coh[:,0]
x = 2*np.real(last_states_pop_coh[:,1])
y = 2*np.imag(last_states_pop_coh[:,1])
plt.plot(radii, x, '.-', label='x')
plt.plot(radii, y, '.-', label='y')
plt.plot(radii, z, label='z')
plt.legend()
plt.title('x y z components')
plt.xlabel('radius')
plt.savefig('xyz_radius.png')
/gpfs/home/fas/jiang/sk943/virtenv3.4/lib/python3.4/site-packages/numpy/core/numeric.py:460: ComplexWarning: Casting complex values to real discards the imaginary part return array(a, dtype, copy=False, order=order)
area = np.pi * radii**2 * np.sign(radii)
plt.plot(area, x, '.-', label='x')
plt.plot(area, y, '.-', label='y')
plt.plot(area, z, label='z')
plt.legend()
plt.title('x y z components')
plt.xlabel(r'$\pi r^2$')
plt.savefig('xyz_area.png')
argument = np.angle(x+1j*y)
plt.plot(area, argument, '-o')
plt.legend()
plt.title(r'$\arg{(x+iy)}$')
plt.xlabel(r'$\pi r^2$')
plt.savefig('argxy_area.png')
/gpfs/home/fas/jiang/sk943/virtenv3.4/lib/python3.4/site-packages/matplotlib/axes.py:4749: UserWarning: No labeled objects found. Use label='...' kwarg on individual plots. warnings.warn("No labeled objects found. "
last_purity = [(r.states[-1]**2).tr() for r in result_vals]
plt.plot(radii, last_purity)
plt.title('Purity')
plt.xlabel('radius')
plt.savefig('purity_radius.png')