Skip to content

Conversation

@aleramos119
Copy link
Contributor

Description

These pull request adds the necessary functionality to do noise characterization for an Ising model, such as propagation and computation of loss and gradients.

Checklist:

  • [x ] The pull request only contains commits that are focused and relevant to this change.
  • [ x] I have added appropriate tests that cover the new/changed functionality.
  • [ x] I have updated the documentation to reflect these changes.
  • [ x] I have added entries to the changelog for any noteworthy additions, changes, fixes, or removals.
  • [ x] I have added migration instructions to the upgrade guide (if needed).
  • [ x] The changes follow the project's style guidelines and introduce no new warnings.
  • [x ] The changes are fully tested and pass the CI checks.
  • [ x] I have reviewed my own code changes.

MaxFroehlich1410 and others added 30 commits February 19, 2025 15:38
force the parameters to be positives
force the parameters to be positive
d_On_d_gk has shape (n_jump_site, n_obs_site, L, n_t)
loss_function changed accordingly
@aaronleesander aaronleesander self-requested a review August 13, 2025 11:15
@aaronleesander aaronleesander added feature New feature or request major Major version update labels Aug 13, 2025
@aaronleesander
Copy link
Member

@aleramos119 Send me a ping when you are done and ready for me to review this. Looking forward to adding this, nice work!

@MaxFroehlich1410 MaxFroehlich1410 self-requested a review August 13, 2025 11:45
@codecov
Copy link

codecov bot commented Aug 15, 2025

Codecov Report

❌ Patch coverage is 93.50181% with 18 lines in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
src/mqt/yaqs/noise_char/optimization.py 92.1% 16 Missing ⚠️
src/mqt/yaqs/noise_char/propagation.py 97.2% 2 Missing ⚠️

📢 Thoughts on this report? Let us know!

Copy link
Member

@aaronleesander aaronleesander left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@aleramos119 Hey I added some comments throughout. I would generally prefer to create something that uses a characterizer..py similar to the simulator.py so that the user can use all these parts without directly getting into the details. Can you make the changes I requested if possible, then maybe think about the usability? I may also be able to edit it directly if you want (but not sure if I have permissions right now).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would remove this and leave it local for you.

Comment on lines +85 to +110
"""A base LossClass to track optimization history and compute averages."""

n_avg = 20

t: np.ndarray
exp_vals_traj: np.ndarray
work_dir: Path
history_file_name: Path
history_avg_file_name: Path

d: int

def __init__(self, *, print_to_file: bool = False) -> None:
"""Initializes the LossClass with default values.

Args:
print_to_file (bool, optional): If True, enables printing output to a file. Defaults to False.
"""
self.n_eval = 0
self.x_history: list[np.ndarray] = []
self.f_history: list[float] = []
self.x_avg_history: list[np.ndarray] = []
self.diff_avg_history: list[float] = []
self.grad_history: list[np.ndarray] = []
self.print_to_file: bool = print_to_file
self.work_dir: Path = Path()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this area, can you change the variable names to be more descriptive instead of single letter variables? I know n, t, etc. probably makes sense in terms of optimization, so if you want to leave it, can you at least give a good description in the docstring?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah, I can do that

Comment on lines +327 to +588
class LossClass2(LossClass):
"""LossClass2 represents the loss for a Ising model with the same noise parameters for each site.

This class encapsulates the objective function and its gradient computation for optimizing noise parameters
(relaxation and dephasing rates) in quantum system simulations. It compares simulated trajectories
to a reference trajectory and provides the sum of squared differences as the loss, along with its gradient
with respect to the noise parameters. It is designed for the case of the same noise parameters for each site,
in total 2 parameters.

Attributes:
print_to_file : bool
Flag indicating whether to print output to a file.
ref_traj : np.ndarray
Copy of the reference trajectory data.
traj_der : Callable[[SimulationParameters], tuple]
Function to compute trajectory derivatives given simulation parameters.
sim_params : SimulationParameters
Deep copy of the simulation parameters, including noise rates.
n_gamma_rel : int
Number of relaxation rates in the simulation parameters.
n_gamma_deph : int
Number of dephasing rates in the simulation parameters.
d : int
Total number of noise parameters (relaxation + dephasing).

Methods:
__init__(sim_params, ref_traj, traj_der, print_to_file=False)
Initializes the loss_class_2 instance with simulation parameters, reference trajectory,
and trajectory derivative function.
__call__(x: np.ndarray) -> tuple
Evaluates the objective function and its gradient for the given noise parameters.
Updates the simulation parameters, runs the trajectory simulation and its derivatives,
computes the difference between simulated and reference trajectories, and calculates:
- The objective function value (sum of squared differences)
- The gradient with respect to the noise parameters
- The simulation time
- The average of the minimum and maximum trajectory times
Array containing the noise parameters to be optimized. The first n_gamma_rel elements
correspond to relaxation rates, and the remaining elements correspond to dephasing rates.
Value of the objective function (sum of squared differences).
Gradient of the objective function with respect to the noise parameters.
Time taken to run the simulation, in seconds.
Average of the minimum and maximum trajectory times from the simulation.
"""

def __init__(
self,
sim_params: SimulationParameters,
ref_traj: np.ndarray,
traj_der: Callable[[SimulationParameters], tuple[np.ndarray, np.ndarray, np.ndarray, list[None]]],
*,
print_to_file: bool = False,
) -> None:
"""Initializes the optimization class for noise characterization.

Args:
sim_params (SimulationParameters): The simulation parameters to be used.
ref_traj (np.ndarray): Reference trajectory as a NumPy array.
traj_der (Callable[[SimulationParameters], tuple]): A callable that computes the trajectory derivative
given simulation parameters.
print_to_file (bool, optional): If True, output will be printed to a file. Defaults to False.

Attributes:
d (int): Dimensionality, set to 2.
print_to_file (bool): Indicates whether to print output to a file.
ref_traj (np.ndarray): Copy of the reference trajectory.
traj_der (Callable): Function to compute trajectory derivative.
sim_params (SimulationParameters): Deep copy of the simulation parameters.
"""
super().__init__(print_to_file=print_to_file)

self.d = 2

self.ref_traj = ref_traj.copy()
self.traj_der = traj_der
self.sim_params = copy.deepcopy(sim_params)

def __call__(self, x: np.ndarray) -> tuple[float, np.ndarray, float, list[None] | list[float]]:
"""Evaluates the objective function and its gradient for the given parameters.

This method updates the simulation parameters with the provided gamma values,
runs the trajectory simulation and its derivative, computes the loss (sum of squared
differences between the simulated and reference trajectories), and calculates the gradient
of the loss with respect to the gamma parameters. It also measures the simulation time
and retrieves the average minimum and maximum trajectory times.

Args:
x (np.ndarray): Array of gamma parameters to be set in the simulation.

Returns:
tuple:
- f (float): The value of the objective function (sum of squared differences).
- grad (np.ndarray): The gradient of the objective function with respect to gamma parameters.
- sim_time (float): The time taken to run the simulation (in seconds).
- avg_min_max_traj_time (Any): Average, minimum and maximum trajectory running times.
"""
self.sim_params.set_gammas(x[0], x[1])

start_time = time.time()

self.t, self.exp_vals_traj, self.d_on_d_gk, avg_min_max_traj_time = self.traj_der(self.sim_params)

end_time = time.time()

_n_jump_site, n_obs_site, sites, nt = np.shape(self.d_on_d_gk)

diff = self.exp_vals_traj - self.ref_traj

f: float = np.sum(diff**2)

# I reshape diff so it has a shape compatible with d_on_d_gk (n_jump_site, n_obs_site, sites, nt)
# to do elemtwise multiplication. Then I sum over the n_obs_site, sites and nt dimensions to
# get the gradient for each gamma, returning a vector of shape (n_jump_site)
grad = np.sum(2 * diff.reshape(1, n_obs_site, sites, nt) * self.d_on_d_gk, axis=(1, 2, 3))

self.post_process(x.copy(), f, grad.copy())

sim_time = end_time - start_time # Simulation time

return f, grad, sim_time, avg_min_max_traj_time


class LossClass2L(LossClass):
"""LossClass2L represents the loss for a Ising model with site-independent noise parameters.

This class encapsulates the objective function and its gradient computation for optimizing noise parameters
(relaxation and dephasing rates) in quantum system simulations. It compares simulated trajectories
to a reference trajectory and provides the sum of squared differences as the loss, along with its gradient
with respect to the noise parameters. It is designed for the case of independent noise parameters for each site,
in total 2*sites parameters.

Attributes:
print_to_file : bool
Flag indicating whether to print output to a file.
ref_traj : np.ndarray
Copy of the reference trajectory data.
traj_der : Callable[[SimulationParameters], tuple]
Function to compute trajectory derivatives given simulation parameters.
sim_params : SimulationParameters
Deep copy of the simulation parameters, including noise rates.
n_gamma_rel : int
Number of relaxation rates in the simulation parameters.
n_gamma_deph : int
Number of dephasing rates in the simulation parameters.
d : int
Total number of noise parameters (relaxation + dephasing).

Methods:
__init__(sim_params, ref_traj, traj_der, print_to_file=False)
Initializes the loss_class_2L instance with simulation parameters,
reference trajectory, and trajectory derivative function.
__call__(x: np.ndarray) -> tuple
Evaluates the objective function and its gradient for the given noise parameters.
Updates the simulation parameters, runs the trajectory simulation and its derivatives,
computes the difference between simulated and reference trajectories, and calculates:
- The objective function value (sum of squared differences)
- The gradient with respect to the noise parameters
- The simulation time
- The average of the minimum and maximum trajectory times
Array containing the noise parameters to be optimized. The first n_gamma_rel elements
correspond to relaxation rates, and the remaining elements correspond to dephasing rates.
Value of the objective function (sum of squared differences).
Gradient of the objective function with respect to the noise parameters.
Time taken to run the simulation, in seconds.
Average of the minimum and maximum trajectory times from the simulation.
"""

def __init__(
self,
sim_params: SimulationParameters,
ref_traj: np.ndarray,
traj_der: Callable[[SimulationParameters], tuple[np.ndarray, np.ndarray, np.ndarray, list[None]]],
*,
print_to_file: bool = False,
) -> None:
"""Initializes the optimization class for noise characterization.

Args:
sim_params (SimulationParameters): Simulation parameters containing relaxation and dephasing rates.
ref_traj (np.ndarray): Reference trajectory data as a NumPy array.
traj_der (Callable[[SimulationParameters], tuple]): Callable that computes the trajectory derivative
given simulation parameters.
print_to_file (bool, optional): If True, enables printing output to a file. Defaults to False.

Attributes:
print_to_file (bool): Flag indicating whether to print output to a file.
ref_traj (np.ndarray): Copy of the reference trajectory.
traj_der (Callable): Function to compute trajectory derivatives.
sim_params (SimulationParameters): Deep copy of the simulation parameters.
n_gamma_rel (int): Number of relaxation rates in the simulation parameters.
n_gamma_deph (int): Number of dephasing rates in the simulation parameters.
d (int): Total number of noise parameters (relaxation + dephasing).
"""
super().__init__(print_to_file=print_to_file)

self.ref_traj = ref_traj.copy()
self.traj_der = traj_der
self.sim_params = copy.deepcopy(sim_params)

self.n_gamma_rel = len(self.sim_params.gamma_rel)
self.n_gamma_deph = len(self.sim_params.gamma_deph)

self.d = self.n_gamma_rel + self.n_gamma_deph

def __call__(self, x: np.ndarray) -> tuple[float, np.ndarray, float, list[None] | list[float]]:
"""Evaluates the objective function and its gradient for the given parameters.

This method updates the simulation parameters with the provided gamma values,
runs the trajectory simulation and its derivatives, computes the difference
between the simulated and reference trajectories, and calculates the objective
function (sum of squared differences) and its gradient with respect to the
gamma parameters. It also records the simulation time and average trajectory
time statistics.
Parameters.
----------
x : np.ndarray
Array containing the gamma parameters to be optimized. The first
`self.n_gamma_rel` elements correspond to one set of gammas, and the
remaining elements correspond to another set.

Returns:
-------
f : float
The value of the objective function (sum of squared differences between
simulated and reference trajectories).
grad : np.ndarray
The gradient of the objective function with respect to the gamma
parameters, flattened into a 1D array.
sim_time : float
The time taken to run the simulation, in seconds.
avg_min_max_traj_time : float
The average of the minimum and maximum trajectory times, as returned by
the trajectory simulation.
"""
self.sim_params.set_gammas(x[: self.n_gamma_rel], x[self.n_gamma_rel :])

start_time = time.time()

t, exp_vals_traj, d_on_d_gk, avg_min_max_traj_time = self.traj_der(self.sim_params)

end_time = time.time()

self.t = t.copy()
self.exp_vals_traj = exp_vals_traj.copy()

_n_jump_site, n_obs_site, sites, nt = np.shape(d_on_d_gk)

diff = exp_vals_traj - self.ref_traj

f: float = np.sum(diff**2)

# I reshape diff so it has a shape compatible with d_on_d_gk (n_jump_site, n_obs_site, sites, nt)
# to do elemtwise multiplication. Then I sum over the n_obs_site and nt dimensions to get
# the gradient for each gamma for each site, returning a matrix of shape (n_jump_site, sites)
# which I then flatten obtaining a vector of shape (n_jump_site*sites)
grad = np.sum(2 * diff.reshape(1, n_obs_site, sites, nt) * d_on_d_gk, axis=(1, 3)).flatten()

self.post_process(x.copy(), f, grad.copy())

sim_time = end_time - start_time # Simulation time

return f, grad, sim_time, avg_min_max_traj_time
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm confused what these are for exactly. For YAQS, it would be best to put the general components in the library, so maybe the first LossClass is fine, then you can add all the Ising stuff as an exxample for someone to build on. I would like to be able to input my Hamiltonian + Noise model and then everything else is taken care of.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I mean normally you want to give the user the possibility to define a different loss function. That's why I did like this. They inherit common functionality from LossClass, but still gives you the flexibility to define your own. But I can make it that somehow you input the Hamiltonian and Noise model.

Comment on lines +22 to +122
class SimulationParameters:
"""A class to encapsulate simulation parameters for open quantum system simulations.

Attributes:
T (float): Total simulation time. Default is 5.
dt (float): Time step for the simulation. Default is 0.1.
J (float): Coupling constant for the Ising Hamiltonian. Default is 1.
g (float): Transverse field strength. Default is 0.5.
observables (list): List of observables to measure (e.g., ['x', 'y', 'z']).
threshold (float): Threshold for truncation in tensor network simulations. Default is 1e-6.
max_bond_dim (int): Maximum bond dimension for tensor network simulations. Default is 4.
order (int): Order of the integration method. Default is 2.
N (int): Number of trajectories or samples for stochastic simulations. Default is 100.
rank (int): Rank for tensor train decompositions. Default is 8.
req_cpus (int): Number of CPUs requested for parallel simulations. Default is 1.
scikit_tt_solver (dict): Dictionary specifying the solver and method for scikit_tt simulations.

Args:
L (int): Number of sites in the system.
gamma_rel (list or float): Relaxation rates for each site, or a single float for all sites.
gamma_deph (list or float): Dephasing rates for each site, or a single float for all sites.

Methods:
set_gammas(gamma_rel, gamma_deph):
Sets the relaxation and dephasing rates for the system. Accepts either a list of length L or a single float.
set_solver(solver='tdvp1', local_solver='krylov_5'):
Sets the solver and local solver method for scikit_tt simulations.

Args:
solver (str): 'tdvp1' or 'tdvp2'.
local_solver (str): 'krylov_<number>' or 'exact'.
"""

T: float = 5
dt: float = 0.1
J: float = 1
g: float = 0.5

threshold: float = 1e-6
max_bond_dim: int = 4
order: int = 2

# For scikit_tt
N: int = 100
rank: int = 8

def __init__(self, sites: int, gamma_rel: list[float] | float, gamma_deph: list[float] | float) -> None:
"""Defines the system with the number of sites and noise parameters.

Parameters
----------
sites : int
The number of sites in the system
gamma_rel : list[float] | float | np.ndarray
Relaxation rates. If a float is provided, the same value is used for all sites (length L).
If a list or array is provided, it must have length L.
gamma_deph : list[float] | float | np.ndarray
Dephasing rates. If a float is provided, the same value is used for all sites (length L).
If a list or array is provided, it must have length L.
"""
self.L = sites

self.set_gammas(gamma_rel, gamma_deph)

def set_gammas(
self, gamma_rel: np.ndarray | list[float] | float, gamma_deph: np.ndarray | list[float] | float
) -> None:
"""Set the relaxation (gamma_rel) and dephasing (gamma_deph) rates for the system.

Args:
gamma_rel (list[float] | float | np.ndarray): Relaxation rates. If a float is
provided, the same value is used for all sites (length L). If a list or array
is provided, it must have length L.
gamma_deph (list[float] | float | np.ndarray): Dephasing rates. If a float is
provided, the same value is used for all sites (length L). If a list or array
is provided, it must have length L.

Raises:
ValueError: If ``gamma_rel`` is a list or array and its length does not match ``L``.
ValueError: If ``gamma_deph`` is a list or array and its length does not match ``L``.

Notes:
This method sets the attributes ``gamma_rel`` and ``gamma_deph`` as lists of
length L.
"""
if (isinstance(gamma_rel, (list, np.ndarray))) and len(gamma_rel) != self.L:
msg = "gamma_rel must be a list of length L."
raise ValueError(msg)
if (isinstance(gamma_deph, (list, np.ndarray))) and len(gamma_deph) != self.L:
msg = "gamma_deph must be a list of length L."
raise ValueError(msg)

if isinstance(gamma_rel, float):
self.gamma_rel = [gamma_rel] * self.L
else:
self.gamma_rel = list(gamma_rel)

if isinstance(gamma_deph, float):
self.gamma_deph = [gamma_deph] * self.L
else:
self.gamma_deph = list(gamma_deph)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would prefer not to have a separate SimulationParameters in the characterization. Could you either

  1. Find a way to get all of this into the AnalogSimParams directly
  2. Create a new CharSimParams in the src/core/data_structures/sim_params.py file?

Comment on lines +125 to +245
def tjm_traj(sim_params_class: SimulationParameters) -> tuple[np.ndarray, np.ndarray, np.ndarray, list[None]]:
"""Simulates the time evolution of an open quantum system using the Lindblad master equation with TJM.

This function constructs the system Hamiltonian and collapse operators for a spin chain with
relaxation and dephasing noise, initializes the system state, and computes the expectation values
of specified observables and their derivatives with respect to the noise parameters over time.
Parameters.
----------
sim_params_class : SimulationParameters
An instance of SimulationParameters containing simulation parameters:
- T (float): Total simulation time.
- dt (float): Time step.
- L (int): Number of sites (spins) in the chain.
- J (float): Ising coupling strength.
- g (float): Transverse field strength.
- gamma_rel (array-like): Relaxation rates for each site.
- gamma_deph (array-like): Dephasing rates for each site.
- observables (list of str): List of observables to measure ('x', 'y', 'z').

Returns:
-------
t : numpy.ndarray
Array of time points at which the system was evolved.
original_exp_vals : numpy.ndarray
Expectation values of the specified observables at each site and time, shape (n_obs_site, L, n_t).
d_on_d_gk : numpy.ndarray
Derivatives of the observables with respect to the noise parameters, shape (n_jump_site, n_obs_site, L, n_t).
avg_min_max_traj_time : list
Placeholder list [None, None, None] for compatibility with other interfaces.

Notes:
-----
- The function uses QuTiP for quantum object and solver operations.
- The system is initialized in the ground state |0>^{⊗L}.
- The Hamiltonian is an Ising model with a transverse field.
- Collapse operators are constructed for both relaxation and dephasing noise.
- The function computes both the expectation values of observables
and their derivatives with respect to noise parameters
using the Lindblad master equation.
"""
sim_time = sim_params_class.T
dt = sim_params_class.dt
sites = sim_params_class.L
coupl = sim_params_class.J
g = sim_params_class.g
gamma_rel = sim_params_class.gamma_rel
gamma_deph = sim_params_class.gamma_deph
ntraj = sim_params_class.N

threshold = sim_params_class.threshold
max_bond_dim = sim_params_class.max_bond_dim
order = sim_params_class.order

t = np.arange(0, sim_time + dt, dt)
n_t = len(t)

# Define the system Hamiltonian
h_0 = MPO()

h_0.init_ising(sites, coupl, g)
# Define the initial state
state = MPS(sites, state="zeros")

noise_model = NoiseModel(
[{"name": "lowering", "sites": [i], "strength": gamma_rel[i]} for i in range(sites)]
+ [{"name": "pauli_z", "sites": [i], "strength": gamma_deph[i]} for i in range(sites)]
)

obs_list = (
[Observable(X(), site) for site in range(sites)]
+ [Observable(Y(), site) for site in range(sites)]
+ [Observable(Z(), site) for site in range(sites)]
)

jump_site_list = [Destroy(), Z()]

obs_site_list = [X(), Y(), Z()]

a_kn_site_list: list[Observable] = []

n_jump_site = len(jump_site_list)
n_obs_site = len(obs_site_list)

for lk in jump_site_list:
for on in obs_site_list:
a_kn_site_list.extend(
Observable(lk.dag() * on * lk - 0.5 * on * lk.dag() * lk - 0.5 * lk.dag() * lk * on, k)
for k in range(sites)
)

new_obs_list = obs_list + a_kn_site_list

sim_params = AnalogSimParams(
observables=new_obs_list,
elapsed_time=sim_time,
dt=dt,
num_traj=ntraj,
max_bond_dim=max_bond_dim,
threshold=threshold,
order=order,
sample_timesteps=True,
)
simulator.run(state, h_0, sim_params, noise_model)

exp_vals = [observable.results for observable in sim_params.observables]

# Separate original and new expectation values from result_lindblad.
n_obs = len(obs_list) # number of measurement operators (should be sites * n_types)
original_exp_vals = exp_vals[:n_obs]
new_exp_vals = exp_vals[n_obs:] # these correspond to the A_kn operators
assert all(v is not None for v in new_exp_vals)

# Compute the integral of the new expectation values to obtain the derivatives
d_on_d_gk = [trapezoidal(new_exp_vals[i], t) for i in range(len(a_kn_site_list))]

d_on_d_gk = np.array(d_on_d_gk).reshape(n_jump_site, n_obs_site, sites, n_t)
original_exp_vals = np.array(original_exp_vals).reshape(n_obs_site, sites, n_t)

avg_min_max_traj_time = [None, None, None] # Placeholder for average, min, and max trajectory time

return t, original_exp_vals, d_on_d_gk, avg_min_max_traj_time
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would also be better as an example in the docs like I mentioned above.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

feature New feature or request major Major version update

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants