Add files via upload

Fig2a_plot.py

@@ -0,0 +1,389 @@
+def onorm(M):
+    MM=np.conj(M).T@M
+    eigenvalues,eigenvectors=np.linalg.eig(MM)
+
+    return np.sqrt(np.max(eigenvalues))
+
+
+
+
+
+
+
+import math as math
+import numpy as np
+import sympy as sp
+from mpmath import coth
+import scipy.linalg as scln
+from scipy.integrate import trapz
+import itertools
+from sympy.tensor.array.expressions import ArraySymbol
+from sympy.abc import i, j, k
+import h5py
+import matplotlib.pyplot as plt
+from engine.inverse.inverse import *
+plt.rc("font", family="serif")
+plt.rc("xtick")
+
+
+plt.rc("ytick")
+fig = plt.figure(  )  # width=8 inches, height=6 inches
+
+ax = fig.add_subplot(1, 1, 1)
+N=500
+wantdelt=0.1
+t=np.linspace(0,wantdelt*(N-1),num=N)
+
+Idval=np.array([[1,0],[0,1]],dtype = 'complex')
+Hs=np.array([[0,1],[1,0]],dtype = 'complex')
+delt=wantdelt
+G=scln.expm(np.kron(1j*Hs*delt/2,Idval)+np.kron(Idval,(-1j*Hs*delt/2)))
+filename = f'Data/Fig 2 Data/jobgqme0.11.0/dyn0.1beta5delt0.1ohmicity1.0.h5'
+
+h5 = h5py.File(filename,'r')
+kmax=14
+#MK=np.zeros((4,4,kmax+1),dtype='complex')
+M = h5["M"]  # VSTOXX futures data
+A = np.array(h5["10AT"])
+print(list(A))
+UU=np.zeros((4,4,N),dtype = 'complex')
+UU[:,:,1]=M[:,:,1]
+Mnorm=np.zeros((kmax+1),dtype = 'complex')
+MMM=np.zeros((4,4,kmax+1),dtype = 'complex')
+
+for k in range(2,kmax+1):
+    tempM=np.zeros((4,4),dtype = 'complex')
+    for j in range(1,k):
+        tempM += np.matmul(M[:, :, j],UU[:, :, k - j])
+    (UU[:,:,(k)])=M[:,:,(k)]+tempM
+for i in range(1,kmax+1):
+    UU[:,:,i]=np.matmul(G,np.matmul(UU[:,:,i],G))
+    MMM[:,:,i]=M[:,:,i]
+MK=KfromU(UU,kmax)
+#ax.plot(np.arange(0,14)*wantdelt,MK[1,2,:],color="k", ls="solid",label="Dyck/INFPI $\mathcal{K}_{12}^N$")
+#ax.plot(np.arange(0,14)*wantdelt,MK[0,0,:],color=(0.6,0.6,0.6), ls="solid",label="Dyck/INFPI $\mathcal{K}_{00}^N$")
+#ax.plot(np.arange(0,14)*wantdelt,MK[0,3,:],color=(0.2,0.2,0.2), ls="solid",label="Dyck/INFPI $\mathcal{K}_{03}^N$")
+for i in range(kmax+1):
+    Mnorm[i]=onorm(M[:,:,i]/(delt**2))
+#ax.plot(t,A[0,:],color="k")
+ax.semilogy(np.arange(2,kmax+1)*wantdelt,Mnorm[2:],color="k", ls="none",marker='.',markersize=7)
+
+h5.close()
+
+ax.set_facecolor((0.9,1,1))
+plt.xlabel("Time",fontsize=15)
+#ax.set_ylabel(r"$\sigma_z$")
+#ax.set_ylim([-0.2,1])
+ax.set_ylim([0.0001,15])
+ax.set_xlim([0,4])
+
+#plt.legend()
+#plt.savefig("plots/K051.png",dpi=1000)
+#plt.show()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+import math as math
+import numpy as np
+import sympy as sp
+from mpmath import coth
+import scipy.linalg as scln
+from scipy.integrate import trapz
+import itertools
+from sympy.tensor.array.expressions import ArraySymbol
+from sympy.abc import i, j, k
+import h5py
+import matplotlib.pyplot as plt
+from engine.inverse.inverse import *
+filename = f'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial1.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A1 = np.transpose(h5["rho"])
+
+filename = f'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial2.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A2 = np.transpose(h5["rho"])
+
+filename = f'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial3.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A3 = np.transpose(h5["rho"])
+
+filename = 'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial4.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A4 = np.transpose(h5["rho"])
+
+Uex=UfromA(A1,A2,A3,A4)
+A5=A1.copy()
+for i in range(1,len(A1)):  
+    A5[:,i]=np.matmul(Uex[:,:,i],[1,0,0,0])
+mmax=150
+Mex=KfromU(Uex,mmax)
+
+t=np.arange(0,0.1*(mmax+1),0.1)
+ddelt=0.1
+Mnormex=np.zeros((mmax+1),dtype = 'complex')
+for i in range(mmax+1):
+    Mnormex[i]=onorm(Mex[:,:,i]/(ddelt**2))
+ax.plot(t[2:],Mnormex[2:],color=(0,0,0),alpha=1, ls="solid",label="HEOM $\mathcal{K}_{12}^N$")
+plt.xticks(fontsize=15)
+
+plt.yticks(fontsize=15)
+#plt.savefig("plots/K051HEOM.png",dpi=1000)
+#plt.plot(t[0:500],A1[0,:500])
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+import math as math
+import numpy as np
+import sympy as sp
+import scipy as sc
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import scipy.linalg as scln
+from scipy.integrate import trapz, quad
+import itertools
+from sympy.tensor.array.expressions import ArraySymbol
+from sympy.abc import i, j, k
+mpl.rcParams.update(mpl.rcParamsDefault)
+beta=5
+kmax=40
+N=500
+wantdelt=0.1
+kondo=0.1
+ga=kondo* np.pi / 2
+No=20000
+hbar=1
+om=np.linspace(-100,100,No)
+yy1=np.zeros(No,dtype = "complex")
+yy2=np.zeros(No,dtype = "complex")
+tarr = np.linspace(0,wantdelt*(N-1),num=N)
+delt=(tarr[1]-tarr[0])
+Hs=np.array([[0,1],[1,0]],dtype = "complex")
+wcut = 7.5
+JJ=om.copy()
+ohmicity=1
+wc=wcut
+def J(xs):  # Spectral density (Ohmic bath with exponential cutoff)
+    output = np.zeros_like(xs)
+    output[xs > 0] = ga * (xs**ohmicity)[xs > 0] * np.exp(-xs[xs > 0] / wc)/(wc**(ohmicity-1))
+    output[xs <= 0] = - ga * ((-xs)**ohmicity)[xs <= 0] * np.exp(xs[xs <= 0] / wc)/(wc**(ohmicity-1))
+    #output=np.heaviside(xs,0.5)ga*(xs**ohmicity)/(wc**(1-ohmicity))*np.exp(-(xs)/wc)
+    return output
+for i in enumerate(om):
+    JJ[i[0]]=J(i[1])
+
+etakkd=np.zeros(kmax+1,dtype="complex")
+
+
+delt2=delt*1
+for iii in enumerate(om):   
+    yy1[iii[0]]=1/(2*np.pi)*J(om[iii[0]])/(om[iii[0]]**2)*np.exp(beta*hbar*om[iii[0]]/2)/(math.sinh(beta*hbar*om[iii[0]]/2))*(1-np.exp(-1j*om[iii[0]]*delt2))
+etanul=trapz(yy1,om)
+
+
+onetwo=[1,-1]
+def yy2f(x,diffk): 
+    return 2/(1*np.pi)*J(x)/(x**2)*np.exp(beta*hbar*x/2)/(math.sinh(beta*hbar*x/2))*((math.sin(x*delt2/2))**2)*np.exp(-1j*x*delt2*(diffk))
+for difk in np.linspace(1,kmax,kmax,dtype="int"):
+    for i in enumerate(om):
+        #print(om[i[0]])
+        yy2[i[0]]=2/(1*np.pi)*J(om[i[0]])/(om[i[0]]**2)*np.exp(beta*hbar*om[i[0]]/2)/(math.sinh(beta*hbar*om[i[0]]/2))*((math.sin(om[i[0]]*delt2/2))**2)*np.exp(-1j*om[i[0]]*delt2*(difk))
+    #aaa=quad(yy2f,-np.inf,np.inf,args=difk)
+    etakkd[difk]=trapz(yy2,om)
+def Influencediff(dx1,dx2,dx3,dx4,diff):
+    x1=onetwo[dx1]
+    x2=onetwo[dx2]
+    x3=onetwo[dx3]
+    x4=onetwo[dx4]
+    Sum=-1/hbar*(x3-x4)*(etakkd[diff]*x1-np.conjugate(etakkd[diff])*x2)     # eq 12 line 1 Nancy quapi I           
+    return np.exp(Sum)
+
+def Influencenull(dx1,dx2): # eq 12 line 2 Nancy quapi I    
+    x1=onetwo[dx1]
+
+    x2=onetwo[dx2]
+    Sum=-1/hbar*(x1-x2)*((etanul)*x1-np.conjugate(etanul)*x2)                
+    return np.exp(Sum)
+
+
+I0_val=np.array([Influencenull(0,0),Influencenull(0,1),Influencenull(1,0),Influencenull(1,1)], dtype = "complex")
+
+TI_val=np.zeros((4,4,kmax+1),dtype = "complex")  
+I_val=np.zeros((4,4,kmax+1),dtype = "complex")  
+A=np.zeros((4,len(tarr)),dtype = "complex")
+A[0,0]=1
+P_val=np.zeros((4,4),dtype = "complex")
+
+
+for i in np.arange(0,kmax+1):
+    I_val[0,0,i]=Influencediff(0,0,0,0,i)
+    I_val[0,1,i]=Influencediff(1,0,0,0,i)
+    I_val[0,2,i]=Influencediff(0,1,0,0,i)
+    I_val[0,3,i]=Influencediff(1,1,0,0,i)
+
+    I_val[1,0,i]=Influencediff(0,0,1,0,i)
+    I_val[1,1,i]=Influencediff(1,0,1,0,i)
+    I_val[1,2,i]=Influencediff(0,1,1,0,i)
+    I_val[1,3,i]=Influencediff(1,1,1,0,i)
+
+    I_val[2,0,i]=Influencediff(0,0,0,1,i)
+    I_val[2,1,i]=Influencediff(1,0,0,1,i)
+    I_val[2,2,i]=Influencediff(0,1,0,1,i)
+    I_val[2,3,i]=Influencediff(1,1,0,1,i)
+
+    I_val[3,0,i]=Influencediff(0,0,1,1,i)
+    I_val[3,1,i]=Influencediff(1,0,1,1,i)
+    I_val[3,2,i]=Influencediff(0,1,1,1,i)
+    I_val[3,3,i]=Influencediff(1,1,1,1,i)
+for i in range(kmax+1):
+    TI_val[:,:,i]=I_val[:,:,i]-np.ones(4)
+#plt.rc('text', usetex=True)
+
+#plt.rc("text", usetex=True)
+
+
+
+I0_valm=np.zeros((4,4),dtype="complex")
+for v in range(4):
+    I0_valm[v,v]=I0_val[v]  
+#fnorm[0]=np.linalg.norm(I_val)
+print(np.linalg.norm(I0_valm))
+x=np.log(np.arange(1,kmax+1))
+for i in range(kmax):
+    x[i]=onorm(TI_val[:,:,i])
+#print(sc.stats.linregress(x, y=np.log(fnorm[1:]), alternative="two-sided"))
+#ax.semilogy(range(1,kmax+1)*delt,np.abs(I_val[1,2,1:]-1),color="k", ls="solid",label="Re${I}_{12}^{\Delta k}$")
+#ax.semilogy(range(1,kmax+1)*delt,np.abs(I_val[1,3,1:]-1),color=(0.6,0.6,0.6), ls="solid",label="Re${I}_{13}^{\Delta k}$")
+#ax.semilogy(range(1,kmax+1)*delt,np.abs(I_val[0,0,1:]-1),color=(0.2,0.2,0.2), ls="solid",label="Re${I}_{00}^{\Delta k}$")
+
+#ax.semilogy(range(1,14)*delt,np.abs(I_val[1,2,1:14]-1),color="k", ls="none",marker='x')
+#ax.semilogy(range(1,14)*delt,np.abs(I_val[1,3,1:14]-1),color=(0.6,0.6,0.6), ls="none",marker='x')
+#ax.semilogy(range(1,14)*delt,np.abs(I_val[0,0,1:14]-1),color=(0.2,0.2,0.2), ls="none",marker='x')
+
+#ax.set_xlabel("Time")
+#ax.set_xlim([delt,delt*kmax])
+#ax.set_ylim([0.99,1.01])
+#plt.plot(np.log(abs(etakkd[1:])))
+#plt.plot(2.5*1/x**0.5
+#ax.set_ylabel(r'$\|\mathcal{K}_{N \Delta t}\|_O$', color='k',fontsize=15)
+#ax2 = ax.twinx()
+ax.semilogy(range(2,kmax+1)*delt,x[1:], color='#5E091A')
+ax.semilogy(range(2,14+1)*delt,x[1:14], color='#5E091A', ls="none",marker='x',markersize=7)
+#ax.set_ylabel(r'$\|\tilde{I}\|_O$', color='#5E091A',fontsize=15)
+ax.tick_params('y', color='#5E091A')
+#fnorm=np.zeros((kmax+1),dtype="complex")
+#matrixone=np.ones((4),dtype="complex")
+#plt.legend()
+
+plt.xticks(fontsize=15)
+plt.yticks(fontsize=15)
+filename="Data/Fig 2 Data/crestcrest6/dyn0.15707963267948966ohmicity1.0.h5"
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+sett = h5["sett"]  # VSTOXX futures data
+nnn = h5["norm"]
+nnnn=np.zeros(17,dtype='complex')
+aaa=np.zeros(17,dtype='complex')
+bbb=np.zeros((4,4,18),dtype='complex')
+#plt.plot(range(2,18)*delt,list(sett[1:]),alpha=0.9,color='black', ls="dashed",label=r'$\|\mathcal{K}_{N,crest}\|$')
+for i in range(17):
+
+    bbb[:,:,i]=sett[:,:,i]
+    #aaa[i]=onorm(bbb[:,:,i]-MMM[:,:,i+1]/0.01)
+    nnnn[i]=nnn[i]
+h5.close()
+print(Mex[2,1,1:8]/0.01)
+print(bbb[2,1,:7])
+
+plt.plot(range(2,18)*delt,(nnnn[1:]),alpha=0.9,color='black', ls="dashed",label=r'$\|\mathcal{K}_{N,crest}-\mathcal{K}_{N}\|$')
+#plt.plot(range(2,16)*delt,(aaa[0:-3]),alpha=0.9,color='black', ls="dashed",label=r'$\|\mathcal{K}_{N,crest}-\mathcal{K}_{N}\|$')
+plt.savefig(f"Fig1a.eps",format='eps',bbox_inches="tight")
+
+plt.show()