<jupyter_start><jupyter_text>Efficient-simulation-of-non-markovian-dynamics-using-Max-Ent-states One of the fundamental problems in Quantum Information Theory and Quantum Computing is the accurate representation of quantum composite systems, in particular: their's states and dynamics. Said composite systems are univoquely represented by mathematical objects, the density operators, which containt the information of all possible n-body correlations present. A notable exception are the Gaussian dynamics. In these dynamics, achievable for bosonic systems, the dynamics is closed over the set of Gaussian states. Said set of Gaussian states are parameterizable in terms of pairwise correlations, thus forming a $n$ or $n^2$-dimensional Riemannian differentiable manifold; with a metric given by Hilbert-Schmidt inner product of operators. This has motivated to search for generalizations of the Gaussian states to still bosonic systems.In this work, a generalization based on the Max-Ent property of Gaussian states is proposed, in which we will consider families of states that maximize the entropy of the system (Max-Ent principle), under the restriction of fixing the mean values of a certain set of independent observables. Strategies to build approximations within this family, which represent arbitrary states, will then be discussed. As an application case, we will study the relative entropy between the states that results from the dynamics in the Dicke model with its corresponding estimates as MaxEnt states defined by their local values and two-body correlations. * We'll compare rho(t) with its max-ent state associated with a base of observables,* compare rho(t) with its projected state, using the sc corresponding to the initial state, associated to a base of observables, * and compare rho(t) with its projected state, using the sc corresponding to the instantaneous state, associated to a base of observables.<jupyter_code>import qutip
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt
import pickle

dim=20

def prod_basis(b1, b2):
    return [qutip.tensor(b,s) for b in b1 for s in b2]

def scalar_prod(op1,op2,rho0=None):
    if op1.dims[0][0]!=op1.dims[0][0]:
        return None
    if rho0 is None:
        rho0 = qutip.qeye(op1.dims[0])/op1.dims[0][0]
    return ((op1.dag()*op2+op2.dag()*op1)*rho0).tr()
        

def base_orto(ops,rho0):
    dim = ops[0].dims[0][0]
    base = []
    # hacer gramm schmidt
    for op in ops:
        coeff = [scalar_prod(op2,op, rho0) for op2 in base]
        op_mod = op - sum([ c*op2 for c, op2 in zip(coeff, base)])
        op_mod = op_mod/np.sqrt(scalar_prod(op_mod,op_mod, rho0))
        base.append(op_mod)
    return base


def proj_op(K,base,rho0):
    return sum([   scalar_prod(b, K,rho0) * b for b in base])

def logM(rho):
    vals, vecs = rho.eigenstates()
    return sum([np.log(val)*vec*vec.dag() for val,vec in zip(vals, vecs) if val>0])

def sqrtM(rho):
    vals, vecs = rho.eigenstates()
    return sum([ (abs(val)**.5)*vec*vec.dag() for val,vec in zip(vals, vecs)])

def rel_entropy(rho, sigma):
    val = (rho*(logM(rho)-logM(sigma))).tr()
    if abs(val.imag)>1.e-6:
        print("rho or sigma not positive")
        print(rho.eigenstates())
        print(sigma.eigenstates())
    return val.real


def bures(rho, sigma):
    val = abs((sqrtM(rho)*sqrtM(sigma)).tr())
    val = max(min(val,1.),-1.)
    return np.arccos(val)/np.pi
        
def maxent_rho(rho, basis):   
    def test(x, rho, basis):
        k = sum([-u*b for u,b in zip(x, basis)])        
        sigma = (.5*(k+k.dag())).expm()
        sigma = sigma/sigma.tr()
        return rel_entropy(rho, sigma)    
    res = opt.minimize(test,np.zeros(len(basis)),args=(rho,basis))
    k = sum([-u*b for u,b in zip(res.x, basis)])        
    sigma = (.5*(k+k.dag())).expm()
    sigma = sigma/sigma.tr()
    return sigma
    
    
def error_maxent_state(rho, basis, distance=bures):
    try:
        sigma = maxent_rho(rho, basis)
        return distance(rho,sigma)
    except:
        print("fail")
        return None
    
    
def error_proj_state(rho, rho0, basis, distance=bures):
    try:
        basis = base_orto(basis, rho0)
        sigma = proj_op(logM(rho), basis, rho0).expm()
        return distance(rho, sigma)
    except:
        print("fail")
        return None
    
class Result(object):
    def __init__(self, ts=None, states=None):
        self.ts = ts
        self.states = states
        self.max_ent_app = None
        self.projrho0_app = None
        self.projrho_inst_app = None 
        
        

def simul(omega_bos=3, omega_s=3, temp=1, gaussian=False, deltat=10., tmax=500., distance=bures):
    basis_bos = [qutip.qeye(dim), qutip.create(dim),qutip.create(dim).dag(),qutip.num(dim)]
    H_bos = qutip.tensor(qutip.num(dim), qutip.qeye(2))
    H_i = qutip.tensor(.5*(qutip.create(dim)+qutip.destroy(dim)), qutip.sigmax())
    # Estado inicial
    rho0 = qutip.tensor((-qutip.num(dim)/temp).expm(), qutip.qeye(2)/2.)
    rho0 = rho0/rho0.tr()
    # Base
    if gaussian:
        basis_spin = [qutip.qeye(2), qutip.sigmax()]
        H_s = qutip.tensor(qutip.qeye(dim),qutip.sigmax())
    else:
        basis_spin = [qutip.qeye(2), qutip.sigmax(),qutip.sigmay(),qutip.sigmaz()]
        H_s = qutip.tensor(qutip.qeye(dim),qutip.sigmaz())
        
    basis = base_orto(prod_basis(basis_bos, basis_spin), rho0)
    # Hamiltoniano
    H = omega_bos * H_bos + omega_s * H_s + H_i

    sampling = int(10*max(1,omega_bos, omega_s)*deltat)
    
    states = [rho0]
    rho = rho0    
    ts= [0]
    for i in range(int(tmax/deltat)):
        result = qutip.mesolve(H, states[-1], np.linspace(0,deltat, sampling))
        states.append(result.states[-1])
        ts.append(deltat*i)
    
    result = Result(ts, states)
    result.times = ts
    result.states = states
    result.max_ent_app = np.array([error_maxent_state(rho, basis, distance) for rho in states])
    result.projrho0_app = np.array([error_proj_state(rho, rho0, basis,distance) for rho in states])
    result.projrho_inst_app = np.array([error_proj_state(rho, qutip.tensor(rho.ptrace([0]),rho.ptrace([1])), 
                                                         basis, distance) for rho in states])
    
    if gaussian:
        title = distance.__name__ + f" - Dinámica gaussiana    dim={dim} wb={omega_bos} dw={abs(omega_s-omega_bos)} "
    else:
        title = distance.__name__ + f" - Dinámica no gaussiana dim={dim} wb={omega_bos} dw={abs(omega_s-omega_bos)} " 

    with open(title+".pkl","wb") as f:
        pickle.dump(result, f)
    return result, title
## Dinámica Gaussiana, resonante


result, title = simul(omega_bos=3., omega_s=3., temp=1, 
                      gaussian=True, 
                      deltat=10., tmax=500., 
                      distance=bures)

plt.plot(result.times, result.max_ent_app, label="max-ent")
plt.plot(result.times, result.projrho0_app, label="proj rho0")
plt.plot(result.times, result.projrho_inst_app, label="proj rho(t)")
plt.xlabel("t")
plt.ylabel("Arccos(F)")

plt.legend()
plt.title(title + f" dim={dim}")
plt.savefig(title + f" dim={dim}.svg")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12371277483996661, 'fraction_numerical': 0.013080990815474534, 'mean_word_length': 3.245126993502658, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '500895', 'license': 'CNRI-Python', 'n_tokens_mistral': 2536, 'n_tokens_neox': 2419, 'n_words': 660, 'path': 'thesis_notebook.ipynb', 'repo_name': 'licTomasPerez/-Code-Thesis-Non-markovian-Dynamics', 'size': '7186'}