# Set up the cluster for parallel computation
from IPython.parallel import Client
client = Client(profile='lsf')
len(client.ids)
48
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 mpl_toolkits.mplot3d import Axes3D
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 two_circles_both_move_radius_one_half(alpha0=2, T=200, samples_factor=1):
assert np.abs(alpha0) < 5.1, 'alpha > 6 at any point in time goes out of the cutoff for the hilbert space'
N = 80
radius = 0.5
samples = int(15000*samples_factor)
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)
beta0 = -alpha0
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]
)
import signal
def signal_handler(signum, frame):
raise Exception("Timed out!")
signal.signal(signal.SIGALRM, signal_handler)
signal.alarm(2*60*60)
try:
ret = qutip.mesolve(zero_op,init_state,unit_interval*T,
c_ops=c_ops,
e_ops=[],
progress_bar=qutip.ui.progressbar.TextProgressBar())
return ret.states[-1], len(ret.states)
except Exception:
signal.alarm(0)
return init_state, 1
alpha0s = np.linspace(2,5,num=16)
Ts = np.concatenate([np.linspace(20,200,num=10)[:-1],np.linspace(200,2000,num=10)])
args = list(product(alpha0s, Ts))
%%px --local
def two_circles_one_moves_radius_one_half_sqrt2(alpha0=2, T=200, samples_factor=1):
assert np.abs(alpha0) < 5.1, 'alpha > 6 at any point in time goes out of the cutoff for the hilbert space'
N = 80
radius = 0.5*(2**0.5)
samples = int(15000*samples_factor)
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)
beta0 = -alpha0
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 * np.ones_like(alphas)
c_ops = expand_lindblad_operator_to_qutip_std_form(
a_op**2,
[-a_op, alphas+betas],
[id_op, alphas*betas]
)
import signal
def signal_handler(signum, frame):
raise Exception("Timed out!")
signal.signal(signal.SIGALRM, signal_handler)
signal.alarm(2*60*60)
try:
ret = qutip.mesolve(zero_op,init_state,unit_interval*T,
c_ops=c_ops,
e_ops=[],
progress_bar=qutip.ui.progressbar.TextProgressBar())
return ret.states[-1], len(ret.states)
except Exception:
signal.alarm(0)
return init_state, 1
def two_circles_one_moves_radius_one_half(alpha0=2, T=200, samples_factor=1):
assert np.abs(alpha0) < 5.1, 'alpha > 6 at any point in time goes out of the cutoff for the hilbert space'
N = 80
radius = 0.5
samples = int(15000*samples_factor)
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)
beta0 = -alpha0
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 * np.ones_like(alphas)
c_ops = expand_lindblad_operator_to_qutip_std_form(
a_op**2,
[-a_op, alphas+betas],
[id_op, alphas*betas]
)
import signal
def signal_handler(signum, frame):
raise Exception("Timed out!")
signal.signal(signal.SIGALRM, signal_handler)
signal.alarm(2*60*60)
try:
ret = qutip.mesolve(zero_op,init_state,unit_interval*T,
c_ops=c_ops,
e_ops=[],
progress_bar=qutip.ui.progressbar.TextProgressBar())
return ret.states[-1], len(ret.states)
except Exception:
signal.alarm(0)
return init_state, 1
results_two_circles_both_move_standard = lview.map(lambda a: two_circles_both_move_radius_one_half(*a), args)
results_two_circles_one_moves_sqrt2_standard = lview.map(lambda a: two_circles_one_moves_radius_one_half_sqrt2(*a), args)
results_two_circles_one_moves_standard = lview.map(lambda a: two_circles_one_moves_radius_one_half(*a), args)
results_two_circles_both_move_standard.wait_interactive()
304/304 tasks finished after 7214 s done
import pickle
with open("results_two_circles_both_move_standard", "wb") as f:
pickle.dump(results_two_circles_both_move_standard.result, f)
results_two_circles_one_moves_sqrt2_standard.wait_interactive()
304/304 tasks finished after 11433 s done
with open("results_two_circles_one_moves_sqrt2_standard", "wb") as f:
pickle.dump(results_two_circles_one_moves_sqrt2_standard.result, f)
results_two_circles_one_moves_standard.wait_interactive()
304/304 tasks finished after 14926 s done
with open("results_two_circles_one_moves_standard", "wb") as f:
pickle.dump(results_two_circles_one_moves_standard.result, f)
import pickle
with open("results_two_circles_both_move_standard", "rb") as f:
results_two_circles_both_move_standard = pickle.load(f)
with open("results_two_circles_one_moves_sqrt2_standard", "rb") as f:
results_two_circles_one_moves_sqrt2_standard = pickle.load(f)
with open("results_two_circles_one_moves_standard", "rb") as f:
results_two_circles_one_moves_standard = pickle.load(f)
len([_ for _ in results_two_circles_both_move_standard if _[1]])
304
def plot_stuff(results, title=''):
steps_purity = np.array([(_[0]**2).tr() if _[1]==15000 else 0 for _ in results])
steps_exchange = np.array([(qutip.coherent(80,-a).dag()*_[0]*qutip.coherent(80,a)).tr() if _[1]==15000 else 0
for _, (a, t) in zip(results,args)])
mask = np.array([_[1] for _ in results])
steps_purity.shape = len(alpha0s), len(Ts)
steps_exchange.shape = len(alpha0s), len(Ts)
mask.shape = len(alpha0s), len(Ts)
#plt.imshow(mask, interpolation='nearest', origin='lower', cmap=plt.get_cmap('RdBu'))
#mask = np.triu(mask[:,::-1],k=5)[:,::-1]
#plt.imshow(mask, interpolation='nearest', origin='lower', cmap=plt.get_cmap('RdBu'))
steps_purity[mask!=15000] = np.nan
steps_purity = np.ma.MaskedArray(steps_purity, mask!=15000)
steps_exchange[mask!=15000] = np.nan
steps_exchange = np.ma.MaskedArray(steps_exchange, mask!=15000)
f = plt.figure(figsize=(16,16))
pvst = f.add_subplot(2,2,1)
pvst.plot(Ts, steps_purity.T[:,:], '.-')
pvst.set_xlabel(r'$T$')
pvst.set_title(r'Purity-vs-$T$ for various $\alpha_0$')
pvst.legend(list(map(str, alpha0s)), loc=4, title=r'$\alpha_0$')
pvst.set_ylim(0.9,1)
pvsa = f.add_subplot(2,2,2)
pvsa.plot(alpha0s, steps_purity[:,:], '.-')
pvsa.set_xlabel(r'$\alpha_0$')
pvsa.set_title(r'Purity-vs-$\alpha_0$ for various $T$')
pvsa.legend(list(map(str, Ts[:10])), loc=4, title='$T$')
pvsa.set_ylim(0.9,1)
pvsboth = f.add_subplot(2,2,3, projection='3d')
am,tm = np.meshgrid(alpha0s,Ts)
pvsboth.plot_wireframe(am,tm,steps_purity.T)
pvsboth.set_ylim(20,200)
pvsboth.set_xlabel(r'$\alpha_0$')
pvsboth.set_ylabel(r'$T$')
pvsboth.set_title(r'Purity vs $\alpha_0$ and $T$')
pvsboth.set_zlim(0.9,1)
phasevsboth = f.add_subplot(2,2,4, projection='3d')
phasevsboth.plot_wireframe(am,tm,np.angle(steps_exchange).T)
phasevsboth.set_ylim(20,200)
phasevsboth.set_xlabel(r'$\alpha_0$')
phasevsboth.set_ylabel(r'$T$')
phasevsboth.set_title(r'$arg(\langle -\alpha_0 | \rho_{final} | \alpha_0 \rangle)$ vs $\alpha_0$ and $T$')
f.suptitle(title)
return f
f = plot_stuff(results_two_circles_both_move_standard, 'Two Blobs: Both Move on Opposite Circles')
f = plot_stuff(results_two_circles_one_moves_standard, 'Two Blobs: Only One Moves')