Notebook

In [1]:

# Set up the cluster for parallel computation
from IPython.parallel import Client
client = Client(profile='lsf')

In [2]:

len(client.ids)

Out[2]:

In [3]:

dview = client[:]
dview.use_dill()
lview = client.load_balanced_view();

In [1]:

#%%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)]]

In [19]:

%%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

In [2]:

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))

In [21]:

%%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

In [ ]:

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)

In [23]:

results_two_circles_both_move_standard.wait_interactive()

 304/304 tasks finished after 7214 s
done

In [24]:

import pickle
with open("results_two_circles_both_move_standard", "wb") as f:
    pickle.dump(results_two_circles_both_move_standard.result, f)

In [25]:

results_two_circles_one_moves_sqrt2_standard.wait_interactive()

 304/304 tasks finished after 11433 s
done

In [26]:

with open("results_two_circles_one_moves_sqrt2_standard", "wb") as f:
    pickle.dump(results_two_circles_one_moves_sqrt2_standard.result, f)

In [27]:

results_two_circles_one_moves_standard.wait_interactive()

 304/304 tasks finished after 14926 s
done

In [28]:

with open("results_two_circles_one_moves_standard", "wb") as f:
    pickle.dump(results_two_circles_one_moves_standard.result, f)

In [3]:

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)

In [4]:

len([_ for _ in results_two_circles_both_move_standard if _[1]])

Out[4]:

In [55]:

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')

In [56]:

f = plot_stuff(results_two_circles_one_moves_standard, 'Two Blobs: Only One Moves')