Skip to content

Commit d7001d3

Browse files
aleramos119aleramos119
committed
qutip_traj computes only the A_kn different than zero
d_On_d_gk has shape (n_jump_site, n_obs_site, L, n_t) loss_function changed accordingly
1 parent 6a44563 commit d7001d3

File tree

2 files changed

+43
-213
lines changed

2 files changed

+43
-213
lines changed

‎src/yaqs/noise_char/optimization.py‎

Lines changed: 15 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -52,37 +52,33 @@ def loss_function(sim_params, ref_traj, traj_der):
5252
tjm_time = end_time - start_time
5353
# print(f"TJM time -> {tjm_time:.4f}")
5454

55-
# Initialize loss
56-
loss = 0.0
55+
5756

5857
# Ensure both lists have the same structure
59-
if len(ref_traj) != len(exp_vals_traj):
58+
if np.shape(ref_traj) != np.shape(exp_vals_traj):
6059
raise ValueError("Mismatch in the number of sites between qt_exp_vals and tjm_exp_vals.")
6160

62-
# Compute squared distance for each site
63-
for ref_vals, tjm_vals in zip(ref_traj, exp_vals_traj):
64-
loss += np.sum((np.array(ref_vals) - np.array(tjm_vals)) ** 2)
65-
6661

67-
n_jump = len(d_On_d_gk)
68-
n_obs = len(d_On_d_gk[0])
69-
n_t = len(d_On_d_gk[0][0])
62+
n_jump_site, n_obs_site, L, nt = np.shape(d_On_d_gk)
7063

71-
n_gr = n_jump//2
7264

65+
# Initialize loss
66+
loss = 0.0
7367

7468
dJ_d_gr = 0
7569
dJ_d_gd = 0
7670

7771

78-
for i in range(n_obs):
79-
for j in range(n_t):
80-
# I have to add all the derivatives with respect to the same gamma_relaxation and gamma_dephasing
81-
for k in range(n_gr):
82-
# The initial half of the jump operators are relaxation operators
83-
dJ_d_gr += 2*(exp_vals_traj[i][j] - ref_traj[i][j]) * d_On_d_gk[k][i][j]
84-
# The second half of the jump operators are dephasing operators
85-
dJ_d_gd += 2*(exp_vals_traj[i][j] - ref_traj[i][j]) * d_On_d_gk[n_gr + k][i][j]
72+
for i in range(n_obs_site):
73+
for j in range(L):
74+
for k in range(nt):
75+
76+
loss += (exp_vals_traj[i,j,k] - ref_traj[i,j,k])**2
77+
78+
# I have to add all the derivatives with respect to the same gamma_relaxation and gamma_dephasing
79+
dJ_d_gr += 2*(exp_vals_traj[i,j,k] - ref_traj[i,j,k]) * d_On_d_gk[0,i,j,k]
80+
81+
dJ_d_gd += 2*(exp_vals_traj[i,j,k] - ref_traj[i,j,k]) * d_On_d_gk[1,i,j,k]
8682

8783

8884

‎src/yaqs/noise_char/propagation.py‎

Lines changed: 28 additions & 194 deletions
Original file line numberDiff line numberDiff line change
@@ -35,150 +35,6 @@ class SimulationParameters:
3535

3636

3737

38-
# def qutip_traj(sim_params_class: SimulationParameters):
39-
40-
# T = sim_params_class.T
41-
# dt = sim_params_class.dt
42-
# L = sim_params_class.L
43-
# J = sim_params_class.J
44-
# g = sim_params_class.g
45-
# gamma_rel = sim_params_class.gamma_rel
46-
# gamma_deph = sim_params_class.gamma_deph
47-
48-
49-
# t = np.arange(0, T + dt, dt)
50-
51-
# '''QUTIP Initialization + Simulation'''
52-
53-
# # Define Pauli matrices
54-
# sx = qt.sigmax()
55-
# sy = qt.sigmay()
56-
# sz = qt.sigmaz()
57-
58-
# # Construct the Ising Hamiltonian
59-
# H = 0
60-
# for i in range(L-1):
61-
# H += -J * qt.tensor([sz if n==i or n==i+1 else qt.qeye(2) for n in range(L)])
62-
# for i in range(L):
63-
# H += -g * qt.tensor([sx if n==i else qt.qeye(2) for n in range(L)])
64-
65-
66-
67-
# # Construct collapse operators
68-
# c_ops = []
69-
# gammas = []
70-
71-
# # Relaxation operators
72-
# for i in range(L):
73-
# c_ops.append(np.sqrt(gamma_rel) * qt.tensor([qt.destroy(2) if n==i else qt.qeye(2) for n in range(L)]))
74-
# gammas.append(gamma_rel)
75-
76-
# # Dephasing operators
77-
# for i in range(L):
78-
# c_ops.append(np.sqrt(gamma_deph) * qt.tensor([sz if n==i else qt.qeye(2) for n in range(L)]))
79-
# gammas.append(gamma_deph)
80-
81-
# #c_ops = [rel0, rel1, rel2,... rel(L-1), deph0, deph1,..., deph(L-1)]
82-
83-
# # Initial state
84-
# psi0 = qt.tensor([qt.basis(2, 0) for _ in range(L)])
85-
86-
# # Define measurement operators
87-
# sx_list = [qt.tensor([sx if n == i else qt.qeye(2) for n in range(L)]) for i in range(L)]
88-
# sy_list = [qt.tensor([sy if n == i else qt.qeye(2) for n in range(L)]) for i in range(L)]
89-
# sz_list = [qt.tensor([sz if n == i else qt.qeye(2) for n in range(L)]) for i in range(L)]
90-
91-
92-
# # obs_list = [x0, x1, x2,..., x(L-1), y0, y1,..., y(L-1), z0, z1,..., z(L-1)]
93-
# obs_list = sx_list + sy_list + sz_list
94-
95-
96-
97-
98-
# A_kn_list = []
99-
# # number of sites
100-
101-
# for site in range(L):
102-
# # For each site, get the two collapse operators and their corresponding gamma values.
103-
# c_op_rel = c_ops[site] # relaxation collapse operator for this site
104-
# gamma_rel = gammas[site]
105-
# c_op_deph = c_ops[site + L] # dephasing collapse operator for this site
106-
# gamma_deph = gammas[site + L]
107-
108-
# # For each observable type (x, y, z), find the corresponding operator from obs_list.
109-
# # We assume obs_list is ordered: [sx_0, sx_1, ..., sx_{L-1}, sy_0, ..., sy_{L-1}, sz_0, ..., sz_{L-1}]
110-
# obs_x = obs_list[site]
111-
# obs_y = obs_list[site + L]
112-
# obs_z = obs_list[site + 2*L]
113-
114-
# # Compute A_kn for relaxation on this site for each observable type.
115-
# for obs in (obs_x, obs_y, obs_z):
116-
# A_kn = (1 / gamma_rel) * (c_op_rel.dag() * obs * c_op_rel -
117-
# 0.5 * obs * c_op_rel.dag() * c_op_rel -
118-
# 0.5 * c_op_rel.dag() * c_op_rel * obs)
119-
# A_kn_list.append(A_kn)
120-
121-
# # Compute A_kn for dephasing on this site for each observable type.
122-
# for obs in (obs_x, obs_y, obs_z):
123-
# A_kn = (1 / gamma_deph) * (c_op_deph.dag() * obs * c_op_deph -
124-
# 0.5 * obs * c_op_deph.dag() * c_op_deph -
125-
# 0.5 * c_op_deph.dag() * c_op_deph * obs)
126-
# A_kn_list.append(A_kn)
127-
128-
129-
# # A_kn_list = [x0rel0,y0rel0,z0rel0,x0deph0,y0deph0,z0deph0,x1rel1,y1rel1,...,z(L-1)deph(L-1)]
130-
131-
132-
# new_obs_list = obs_list + A_kn_list
133-
134-
135-
136-
137-
# n_obs= len(obs_list)
138-
# n_jump= len(c_ops)
139-
140-
141-
# # Exact Lindblad solution
142-
# result_lindblad = qt.mesolve(H, psi0, t, c_ops, new_obs_list, progress_bar=True)
143-
144-
# exp_vals = []
145-
# for i in range(len(new_obs_list)):
146-
# exp_vals.append(result_lindblad.expect[i])
147-
148-
149-
# # Separate original and new expectation values
150-
# original_exp_vals = exp_vals[:n_obs]
151-
# new_exp_vals = exp_vals[n_obs:]
152-
153-
154-
# n_types = len(obs_list) // L # number of observable types, e.g. 3 for x,y,z
155-
# n_noise = 2 # since you have relaxation and dephasing
156-
# n_Akn_per_site = n_noise * n_types # e.g., 2*3 = 6
157-
158-
# # new_exp_vals (for the A_kn part) should have length = L * n_Akn_per_site.
159-
# # Reshape it into a list of L blocks, each containing n_Akn_per_site arrays.
160-
# A_kn_exp_vals = [ new_exp_vals[site * n_Akn_per_site : (site + 1) * n_Akn_per_site]
161-
# for site in range(L) ]
162-
163-
# # Compute the derivative for each expectation value using trapezoidal integration.
164-
# d_On_d_gk = [ [ trapezoidal(A_kn_exp_vals[site][j], t)
165-
# for j in range(n_Akn_per_site) ]
166-
# for site in range(L) ]
167-
168-
169-
170-
# n_obs = len(obs_list) # still 3L
171-
# n_sites = sim_params_class.L # number of sites
172-
# # new_exp_vals now has length = 6 * L (since A_kn_list now has 6 elements xireli,yireli,zireli,xidephi,yidephi,zidephi for site i)
173-
174-
# # Reshape new_exp_vals into a list of L lists, each containing 6 entries.
175-
# A_kn_exp_vals = [ new_exp_vals[site * 6 : (site + 1) * 6] for site in range(n_sites) ]
176-
177-
# # Then compute the derivative for each of these 6 expectation values per site:
178-
# d_On_d_gk = [ [ trapezoidal(A_kn_exp_vals[site][j], t) for j in range(6) ] for site in range(n_sites) ]
179-
180-
181-
# return t, original_exp_vals, d_On_d_gk, A_kn_exp_vals
18238

18339
def qutip_traj(sim_params_class: SimulationParameters):
18440

@@ -193,6 +49,8 @@ def qutip_traj(sim_params_class: SimulationParameters):
19349

19450
t = np.arange(0, T + dt, dt)
19551

52+
n_t = len(t)
53+
19654
'''QUTIP Initialization + Simulation'''
19755

19856
# Define Pauli matrices
@@ -244,45 +102,29 @@ def qutip_traj(sim_params_class: SimulationParameters):
244102
obs_list.extend([qt.tensor([sz if n == i else qt.qeye(2) for n in range(L)]) for i in range(L)])
245103

246104

247-
# Create new set of observables by multiplying every operator in obs_list with every operator in c_ops
248-
A_kn_list= []
249-
for i,c_op in enumerate(c_ops):
250-
for obs in obs_list:
251-
A_kn_list.append( (1/gammas[i]) * (c_op.dag()*obs*c_op - 0.5*obs*c_op.dag()*c_op - 0.5*c_op.dag()*c_op*obs) )
252105

253106

254107

255-
# A_kn_list = []
256-
# n_types = len(sim_params_class.observables) # number of observable types
257-
# for site in range(L):
258-
# # For each site, get the two collapse operators and their corresponding gamma values.
259-
# c_op_rel = c_ops[site] # relaxation collapse operator for this site
260-
# gamma_rel = gammas[site]
261-
# c_op_deph = c_ops[site + L] # dephasing collapse operator for this site
262-
# gamma_deph = gammas[site + L]
108+
jump_site_list = [ qt.destroy(2) , sz]
109+
110+
obs_site_list = [sx, sy, sz]
111+
112+
113+
A_kn_site_list = []
263114

264-
# # For each observable type, get the corresponding operator from obs_list.
265-
# # The operator for the k-th observable type at this site is:
266-
# # obs_list[site + k*L]
267-
# for k in range(n_types):
268-
# obs_current = obs_list[site + k * L]
269-
# # Compute A_kn for relaxation on this site for this observable type.
270-
# A_kn = (1 / gamma_rel) * (c_op_rel.dag() * obs_current * c_op_rel -
271-
# 0.5 * obs_current * c_op_rel.dag() * c_op_rel -
272-
# 0.5 * c_op_rel.dag() * c_op_rel * obs_current)
273-
# A_kn_list.append(A_kn)
274115

275-
# # And now for dephasing on this site for each observable type.
276-
# for k in range(n_types):
277-
# obs_current = obs_list[site + k * L]
278-
# A_kn = (1 / gamma_deph) * (c_op_deph.dag() * obs_current * c_op_deph -
279-
# 0.5 * obs_current * c_op_deph.dag() * c_op_deph -
280-
# 0.5 * c_op_deph.dag() * c_op_deph * obs_current)
281-
# A_kn_list.append(A_kn)
282-
# # A_kn_list = [x0rel0,y0rel0,z0rel0,x0deph0,y0deph0,z0deph0,x1rel1,y1rel1,...,z(L-1)deph(L-1)]
116+
n_jump_site = len(jump_site_list)
117+
n_obs_site = len(obs_site_list)
283118

284119

285-
new_obs_list = obs_list + A_kn_list
120+
for lk in jump_site_list:
121+
for on in obs_site_list:
122+
for k in range(L):
123+
A_kn_site_list.append( qt.tensor([ lk.dag()*on*lk - 0.5*on*lk.dag()*lk - 0.5*lk.dag()*lk*on if n == k else qt.qeye(2) for n in range(L)]) )
124+
125+
126+
127+
new_obs_list = obs_list + A_kn_site_list
286128

287129
n_obs= len(obs_list)
288130
n_jump= len(c_ops)
@@ -301,33 +143,25 @@ def qutip_traj(sim_params_class: SimulationParameters):
301143
original_exp_vals = exp_vals[:n_obs]
302144
new_exp_vals = exp_vals[n_obs:] # these correspond to the A_kn operators
303145

304-
# # Determine parameters:
305-
# n_types = len(sim_params_class.observables) # e.g., 3 for ['x','y','z']
306-
# n_noise = 2 # since you have relaxation and dephasing
307-
# n_Akn_per_site = n_noise * n_types # e.g., 2*3 = 6
308-
309-
# # new_exp_vals should have a total length of L * n_Akn_per_site.
310-
# # Reshape it into a list of L lists, each containing n_Akn_per_site arrays.
311-
# A_kn_exp_vals = [new_exp_vals[site * n_Akn_per_site : (site + 1) * n_Akn_per_site]
312-
# for site in range(sim_params_class.L)]
313146

314-
# # Compute the derivative for each A_kn expectation value using trapezoidal integration.
315-
# d_On_d_gk = [
316-
# [trapezoidal(A_kn_exp_vals[site][j], t) for j in range(n_Akn_per_site)]
317-
# for site in range(sim_params_class.L)
318-
# ]
319-
320-
# Reshape new_exp_vals to be a list of lists with dimensions n_jump times n_obs
321-
A_kn_exp_vals = [new_exp_vals[i * n_obs:(i + 1) * n_obs] for i in range(n_jump)]
322147

323148
# Compute the integral of the new expectation values to obtain the derivatives
324-
d_On_d_gk = [ [trapezoidal(A_kn_exp_vals[i][j],t) for j in range(n_obs)] for i in range(n_jump) ]
149+
d_On_d_gk = [ trapezoidal(new_exp_vals[i],t) for i in range(len(A_kn_site_list)) ]
150+
151+
d_On_d_gk = np.array(d_On_d_gk).reshape(n_jump_site, n_obs_site, L, n_t)
152+
original_exp_vals = np.array(original_exp_vals).reshape(n_obs_site, L, n_t)
325153

326154
# return t, original_exp_vals, d_On_d_gk, A_kn_exp_vals
327155
return t, original_exp_vals, d_On_d_gk
328156

329157

330158

159+
160+
161+
162+
163+
164+
331165
def tjm(sim_params_class: SimulationParameters, N=1000):
332166

333167
T = sim_params_class.T

0 commit comments

Comments
 (0)