Merge branch 'Noise_Characterization' of https://github.com/cda-tum/mqt-yaqs into Noise_Characterization

Author: MaxFroehlich1410

src/yaqs/noise_char/optimization.py

@@ -5,10 +5,17 @@
 
 def trapezoidal(y, x):
 
-    integral = np.zeros(len(y))
+    len_x = len(x)
+    len_y = len(y)
+
+    integral = np.zeros(len_y)
 
     integral[0] = 0
 
+    if len_x != len_y:
+        raise ValueError("Mismatch in the number of elements between x and y. len(x) = {len_x} and len(y) = {len_y}")
+
+
     for i in range(1,len(y)):
         integral[i] = integral[i-1] + 0.5*(x[i] - x[i-1])*(y[i] + y[i-1])
 

src/yaqs/noise_char/propagation.py

@@ -1,6 +1,8 @@
 import numpy as np
 import qutip as qt
 
+import scipy as sp
+
 from yaqs.core.data_structures.networks import MPO, MPS
 from yaqs.core.data_structures.noise_model import NoiseModel
 from yaqs.core.data_structures.simulation_parameters import Observable, PhysicsSimParams
@@ -9,6 +11,10 @@
 
 from yaqs.noise_char.optimization import trapezoidal
 
+import scikit_tt.tensor_train as tt
+from scikit_tt.tensor_train import TT
+import scikit_tt.solvers.ode as ode
+import scikit_tt
 
 
 
@@ -23,6 +29,10 @@ class SimulationParameters:
     gamma_deph: float = 0.1
     observables = ['x','y','z']
 
+    # For scikit_tt
+    N:int = 100
+    rank: int= 8
+
 
 
 # def qutip_traj(sim_params_class: SimulationParameters):
@@ -363,3 +373,151 @@ def tjm(sim_params_class: SimulationParameters, N=1000):
 
     return t, tjm_exp_vals
 
+
+
+
+def construct_Ank(O_list, L_list):
+    """Construct observable array (same for all sites)"""
+    # define necessary tensors for compact construction
+    L = np.asarray(L_list)
+    L_dag = np.conjugate(L).transpose([0,2,1])
+    L_prod = np.einsum('ijk,ikl->ijl', L_dag, L)
+    O = np.asarray(O_list)
+    I = np.eye(2)
+    # define A_nk
+    A_nk = np.einsum('ijk,ilm->ijklm', L_dag, L) - 0.5*np.einsum('ij,klm->kijlm',I,L_prod) - 0.5*np.einsum('ijk,lm->ijklm',L_prod,I)
+    A_nk = np.einsum('ijklm,nkl->ijnm', A_nk, O)
+    return A_nk
+
+
+def evaluate_Ank(A_nk, state):
+    """Compute expected values of A_nk for each site (assuming that state is right-orthonormal)"""
+    n_jumps = A_nk.shape[0]
+    n_obs = A_nk.shape[2]
+    n_sites = state.order
+    A = np.zeros([n_sites, n_jumps, n_sites, n_obs], dtype=complex)
+
+    A[0,:,0,:] = np.einsum('ijk, mjnl, ilk -> mn', np.conjugate(state.cores[0][:,:,0,:]), A_nk, state.cores[0][:,:,0,:])
+    
+    for i in range(state.order-1):
+        # apply SVD
+        [u, s, v] = sp.linalg.svd(state.cores[i].reshape(state.ranks[i]*state.row_dims[i]*state.col_dims[i],state.ranks[i + 1]), full_matrices=False, overwrite_a=True, check_finite=False)
+        # define updated rank and core
+        state.ranks[i + 1] = u.shape[1]
+        state.cores[i] = u.reshape(state.ranks[i], state.row_dims[i], state.col_dims[i], state.ranks[i + 1])
+        # shift non-orthonormal part to next core
+        state.cores[i + 1] = np.tensordot(np.diag(s).dot(v), state.cores[i + 1], axes=(1, 0))
+
+        A[i+1,:,i+1,:] = np.einsum('ijk, mjnl, ilk -> mn', np.conjugate(state.cores[i+1][:,:,0,:]), A_nk, state.cores[i+1][:,:,0,:])
+
+    return A.transpose([1,0,3,2]).reshape([n_sites*n_jumps, n_sites*n_obs])
+
+
+
+
+
+
+def scikit_tt_traj(sim_params_class: SimulationParameters):
+
+
+    T = sim_params_class.T
+    dt = sim_params_class.dt
+    L = sim_params_class.L
+    J = sim_params_class.J
+    g = sim_params_class.g
+    gamma_rel = sim_params_class.gamma_rel
+    gamma_deph = sim_params_class.gamma_deph
+
+    rank = sim_params_class.rank
+    N = sim_params_class.N
+
+
+    t = np.arange(0, T + dt, dt) 
+    timesteps=len(t)-1
+
+
+    # Parameters
+    X = np.array([[0, 1], [1, 0]])
+    Y = np.array([[0, -1j], [1j, 0]])
+    Z = np.array([[1, 0], [0, -1]])
+    I = np.eye(2)
+    L_1 = np.array([[0, 1], [0, 0]])
+    L_2 = np.array([[1, 0], [0, -1]])
+    O_list = [X, Y, Z]
+    L_list = [L_1, L_2]
+
+    # A_nk = construct_Ank(O_list, L_list)
+
+    
+    cores = [None] * L
+    cores[0] = tt.build_core([[-g * X, - J * Z, I]])
+    for i in range(1, L - 1):
+        cores[i] = tt.build_core([[I, 0, 0], [Z, 0, 0], [-g * X, - J * Z, I]])
+    cores[-1] = tt.build_core([I, Z, -g*X])
+    hamiltonian = TT(cores)# jump operators and parameters
+
+    jump_operator_list = [[L_1, L_2] for _ in range(L)]
+    jump_parameter_list = [[np.sqrt(gamma_rel), np.sqrt(gamma_deph)] for _ in range(L)]
+
+
+    obs_list=[]
+
+    for pauli in O_list:
+        for j in range(L):
+            obs= tt.eye(dims=[2]*L)
+            obs.cores[j]=np.zeros([1,2,2,1], dtype=complex)
+            obs.cores[j][0,:,:,0]=pauli
+            obs_list.append(obs)
+
+    
+
+    n_obs= len(obs_list)
+    n_jump= 2*L
+
+
+
+    exp_vals = np.zeros([n_obs,timesteps+1])
+
+    A_nk=construct_Ank(O_list, L_list)
+
+
+    A_kn_numpy=np.zeros([n_jump, n_obs, timesteps+1],dtype=complex)
+
+    for k in range(N):
+        initial_state = tt.unit([2] * L, [0] * L)
+        for i in range(rank - 1):
+            initial_state += tt.unit([2] * L, [0] * L)
+        initial_state = initial_state.ortho()
+        initial_state = (1 / initial_state.norm()) * initial_state
+
+        
+        for j in range(n_obs):
+            exp_vals[j,0] += initial_state.transpose(conjugate=True)@obs_list[j]@initial_state
+        
+        A_kn_numpy[:,:,0]+= evaluate_Ank(A_nk, initial_state)
+        
+        
+        
+        for i in range(timesteps):
+            initial_state = ode.tjm(hamiltonian, jump_operator_list, jump_parameter_list, initial_state, dt, 1)[-1]
+
+            for j in range(n_obs):                
+                exp_vals[j,i+1] += initial_state.transpose(conjugate=True)@obs_list[j]@initial_state
+
+            A_kn_numpy[:,:,i+1] += evaluate_Ank(A_nk, initial_state)
+            
+
+
+    exp_vals = (1/N)*exp_vals
+    
+    A_kn_numpy = (1/N)*A_kn_numpy
+    
+
+
+    ## The .real part is added as a workaround 
+    d_On_d_gk = [ [trapezoidal(A_kn_numpy[i,j].real,t)  for j in range(n_obs)] for i in range(n_jump) ]
+
+
+
+
+    return t, exp_vals, d_On_d_gk