# Copyright © 2019-2023 HQS Quantum Simulations GmbH. All Rights Reserved.
"""Bath fitting routines."""
import numpy as np
import scipy.linalg as la
from struqture_py import mixed_systems, spins, bosons, fermions
from hqs_noise_app import ( # type: ignore
coupling_to_spectral_function,
SpinBRNoiseOperator,
calculate_excitation_spectral_function,
FermionBRNoiseOperator,
integrate_over_frequency_intervals,
)
from typing import Optional, List, Dict, Tuple, Any, Union
from lmfit import Parameters, minimize
import warnings
from copy import deepcopy
from scipy.interpolate import PchipInterpolator as pchip
[docs]
class BathFitter:
"""Utility for fitting effective open quantum systems to original open quantum systems."""
[docs]
def __init__(
self,
number_boson_modes: int,
spins_per_bosonic_mode: int = 1,
broadening_constraint: Optional[List[float]] = None,
background_broadening_ratio: float = 0.0,
minimum_eigenfrequencies: Optional[float] = None,
maximum_eigenfrequencies: Optional[float] = None,
fitting_window: Optional[Tuple[float, float, int]] = None,
coupling_types: Optional[Union[Dict[Tuple[int, int], List[str]], List[str]]] = None,
coupling_indices: Optional[List[Tuple[int, int]]] = None,
max_fitting_iterations: int = 5,
max_fitting_error: float = 0.05,
) -> None:
"""Initialized Fitting utility.
Args:
number_boson_modes (int): Number of bosonic modes used for the fit of bosonic baths.
spins_per_bosonic_mode (int): Number of spin modes used to represent one bosonic mode.
broadening_constraint (Optional[List[float]]): The optional broadening constraints.
When None broadenings of bosonic modes
are fitted freely.
When given the relative broadening of all
modes is fixed and only a prefactor is
fitted. The prefactor corresponds to
the Trotter timestep in a quantum circuit.
background_broadening_ratio (float): Adds a constant background offset
to the diagonal spectral functions when fitting.
Given in terms of (as a ratio of) the average broadening.
minimum_eigenfrequencies (Optional[float]): Minimal value allowed
for bath eigenfrequencies.
maximum_eigenfrequencies (Optional[float]): Maximum value allowed
for bath eigenfrequencies.
fitting_window (Optional[Tuple[float,float, int]]): The frequency window used for the
fitting (start, end , steps).
If no values are provided the functions
uses the whole frequency range to determine
the fit.
coupling_types (Optional[Union[Dict[Tuple[int, int], List[str]], List[str]]]): A list
of the couplings to include. If None, all
the couplings are used: X, Y, Z.
coupling_indices (Optional[List[Tuple[int, int]]]): A list of allowed fermionic \
hopping operators of the form `c^dagger_j c_k` that are allowed to couple\
to bosonic modes.\
For example [(0,0), (0,1)] only allows coupling operators\
`c^dagger_0 c_0` and `c^dagger_0 c_1` \
If None, all the couplings are allowed. The default is None.
"""
self.number_boson_modes = number_boson_modes
self.spins_per_bosonic_mode = spins_per_bosonic_mode
self.broadening_constraint = broadening_constraint
self.background_broadening_ratio = background_broadening_ratio
self.minimum_eigenfrequencies = minimum_eigenfrequencies
self.maximum_eigenfrequencies = maximum_eigenfrequencies
self.fitting_window = fitting_window
self.coupling_indices = coupling_indices
self._max_fitting_iterations = max_fitting_iterations
self._max_fitting_error = max_fitting_error
# The cached last fitted broadenings.
# When using None broadening constraints the broadenings are fitted but not automatically
# returned. This member stores the last fitted broadenings and enables to use
# the last_fitted_broadenings property to retrieve them.
self.__cached_broadenings: Optional[List[float]] = None
if coupling_types is None:
self.coupling_types: Union[Dict[Tuple[int, int], List[str]], List[str]] = [
"X",
"Y",
"Z",
]
else:
self.coupling_types = coupling_types
@property
def max_fitting_iterations(self) -> int:
"""The number of retries allowed when fitting the spectrum.
The bath fitter uses a simple metric for the quality of the fit:
Let A be the sum of squares of the difference between the
fitter and target, and B be the sum of squares of the fitted spectral function.
The quality of the fit is defined as the ratio A/B where a small ratio
corresponds to a good fit. By default, a deviation of 5% is allowed.
This value can be changed using the `max_fitting_error` setter.
If the criterion is not met, the fitting is retried and if the number of retries
exceeds the maximum, the fit fails.
Returns:
(int): The number of retries allowed when fitting the spectrum.
"""
return self._max_fitting_iterations
@max_fitting_iterations.setter
def max_fitting_iterations(self, value: int) -> None:
self._max_fitting_iterations = value
@property
def last_fitted_broadenings(self) -> Optional[List[float]]:
"""The broadenings from the last fit of the spectrum.
Returns:
Optional[List[float]]: The broadenings if a fit has been performed.
"""
return self.__cached_broadenings
def _set_cached_broadenings(
self,
params: Parameters,
) -> None:
"""Set the cached broadenings to the values of the parameters.
Args:
params: the parameters to extract the broadenings from
"""
new_broadenings = []
for m in range(self.number_boson_modes):
new_broadenings.append(float(params[f"gam_{m}"]))
self.__cached_broadenings = new_broadenings
@property
def max_fitting_error(self) -> float:
"""The maximum allowed fitting error in fitting the spectrum.
The bath fitter uses a simple metric for the quality of the fit:
Let A be the sum of squares of the difference between the
fitter and target and B be the sum of squares of the fitted spectral function.
The quality of the fit is defined as the ratio A/B where a small ratio
corresponds to a good fit. By default, a deviation of 5% is allowed.
If the criterion is not met, the fitting is retried and if the number of retries
exceeds the maximum the fit fails. By default, the number of retries is 5.
This value can be changed using the `max_fitting_iterations` setter.
Returns:
(float): The maximum allowed fitting error as a ratio between spectrum deviation\
and the spectrum.
"""
return self._max_fitting_error
@max_fitting_error.setter
def max_fitting_error(self, value: float) -> None:
self._max_fitting_error = value
def _get_coupling_types(
self, number_particles: int, is_spin: bool
) -> Dict[Tuple[int, int], List[str]]:
if isinstance(self.coupling_types, list):
if not is_spin:
warnings.warn(
"Fermionic systems have only one coupling type."
+ " Ignoring coupling types setting."
)
else:
tmp_couplings = deepcopy(self.coupling_types)
coupling_types = {}
for particle in range(number_particles):
for mode in range(self.number_boson_modes):
coupling_types[(particle, mode)] = tmp_couplings
return coupling_types
else:
return self.coupling_types
[docs]
def to_json(self) -> Dict[str, Any]:
"""Returns a JSON representation of the object.
Returns:
(Dict[str, Any]): JSON representation of the object.
"""
return {
"number_boson_modes": self.number_boson_modes,
"spins_per_bosonic_mode": self.spins_per_bosonic_mode,
"broadening_constraint": self.broadening_constraint,
"background_broadening_ratio": self.background_broadening_ratio,
"minimum_eigenfrequencies": self.minimum_eigenfrequencies,
"maximum_eigenfrequencies": self.maximum_eigenfrequencies,
"fitting_window": self.fitting_window,
"coupling_types": self.coupling_types,
"coupling_indices": self.coupling_indices,
"max_fitting_iterations": self.max_fitting_iterations,
"max_fitting_error": self.max_fitting_error,
}
[docs]
@classmethod
def from_json(cls, json_dict: Dict[str, Any]) -> "BathFitter":
"""Initializes the object from a JSON representation.
Args:
json_dict (Dict[str, Any]): JSON representation of the object.
"""
self = cls(
number_boson_modes=json_dict["number_boson_modes"],
spins_per_bosonic_mode=json_dict["spins_per_bosonic_mode"],
broadening_constraint=json_dict["broadening_constraint"],
background_broadening_ratio=json_dict["background_broadening_ratio"],
minimum_eigenfrequencies=json_dict["minimum_eigenfrequencies"],
maximum_eigenfrequencies=json_dict["maximum_eigenfrequencies"],
fitting_window=json_dict["fitting_window"],
coupling_types=json_dict["coupling_types"],
coupling_indices=json_dict["coupling_indices"],
max_fitting_iterations=json_dict["max_fitting_iterations"],
max_fitting_error=json_dict["max_fitting_error"],
)
return self
[docs]
def fit_boson_bath_to_boson_bath(
self,
original_system: mixed_systems.MixedLindbladOpenSystem,
frequencies: np.ndarray,
) -> Tuple[mixed_systems.MixedLindbladOpenSystem, Optional[float]]:
"""Fits a model Spin-Boson system to an original Spin-Boson system.
Args:
original_system (MixedLindbladOpenSystem): The spin-boson system bath model
the effective model is fitted to
frequencies (ndarray): The frequencies for the spectral functions
that is calculated from the Bosonic bath
Returns:
(MixedHamiltonianSystem, Optional[float]): The Spin-Boson system and
for constrained broadenings
the fitted prefactor.
"""
spectrum = coupling_to_spectral_function(original_system, frequencies)
new_system = spins.SpinHamiltonianSystem(original_system.number_spins()[0])
original_system = original_system.system()
for index in original_system.keys():
# Filter out all indices where all boson indices in Mixed
# indices are empty, in other words the system terms
if all(str(x) == "I" for x in index.bosons()):
# new index that contains the spin part of the old index as
new_system.set(index.spins()[0], original_system.get(index))
return self.fit_boson_bath_to_spectral_function(
new_system,
spectral_function=spectrum,
)
[docs]
def fit_fermion_boson_system_to_fermion_boson_system(
self,
original_system: mixed_systems.MixedLindbladOpenSystem,
frequencies: np.ndarray,
) -> Tuple[mixed_systems.MixedLindbladOpenSystem, Optional[float]]:
"""Fits a model Fermion-Boson system to an original Fermion-Boson system.
Args:
original_system (MixedLindbladOpenSystem): The fermion-boson system bath model
the effective model is fitted to
frequencies (ndarray): The frequencies for the spectral functions
that is calculated from the Bosonic bath
Returns:
(MixedHamiltonianSystem, Optional[float]): The Fermion-Boson system and
for constrained broadenings
the fitted prefactor.
"""
spectrum = coupling_to_spectral_function(original_system, frequencies)
new_system = fermions.FermionHamiltonianSystem(original_system.number_fermionic_modes()[0])
original_system = original_system.system()
for index in original_system.keys():
# Filter out all indices where all boson indices in Mixed
# indices are empty, in other words the system terms
if all(str(x) == "I" for x in index.bosons()):
# new index that contains the spin part of the old index as
new_system.set(index.fermions()[0], original_system.get(index))
return self.fit_fermion_boson_system_to_spectral_function(
new_system,
spectral_function=spectrum,
)
[docs]
def fit_spin_bath_to_boson_bath(
self,
noise_app: Any,
original_system: mixed_systems.MixedLindbladOpenSystem,
frequencies: np.ndarray,
device: Any,
) -> Tuple[mixed_systems.MixedHamiltonianSystem, Optional[float]]:
"""Fits a model Spin-Spin system to an original Spin-Boson system.
Args:
noise_app (HqsNoiseApp): The noise app used to construct simulation programs
from the result of the fitting
original_system (MixedLindbladOpenSystem): The spin-boson system bath model
the effective model is fitted to
frequencies (ndarray): The frequencies for the spectral functions
is calculated from the Bosonic bath
device (Device): The qoqo device on which the spin-spin system bath model
should be simulated. Used to calculate the circuit depth
necessary to determine the Trotter timestep
Returns:
(MixedHamiltonianSystem, Optional[float]): The Spin-Spin system and
for constrained broadenings
the trotter timestep that
fixes the correct broadening
prefactor.
"""
spectrum = coupling_to_spectral_function(original_system, frequencies)
new_system = spins.SpinHamiltonianSystem(original_system.number_spins()[0])
original_system = original_system.system()
for index in original_system.keys():
# Filter out all indices where all boson indices in Mixed
# indices are empty, in other words the system terms
if all(str(x) == "I" for x in index.bosons()):
new_system.set(index.spins()[0], original_system.get(index))
return self.fit_spin_bath_to_spectral_function(
noise_app, new_system, spectral_function=spectrum, device=device
)
[docs]
def fit_fermion_spin_system_to_fermion_boson_system(
self,
noise_app: Any,
original_system: mixed_systems.MixedLindbladOpenSystem,
frequencies: np.ndarray,
device: Any,
) -> Tuple[mixed_systems.MixedHamiltonianSystem, Optional[float]]:
"""Fits a model Fermion-Spin system to an original Fermion-Boson system.
Args:
noise_app (HqsNoiseApp): The noise app used to construct simulation programs
from the result of the fitting
original_system (MixedLindbladOpenSystem): The fermion-boson system bath model
the effective model is fitted to
frequencies (ndarray): The frequencies for the spectral functions
is calculated from the Bosonic bath
device (Device): The qoqo device on which the fermion-spin system bath model
should be simulated. Used to calculate the circuit depth
necessary to determine the Trotter timestep
Returns:
(MixedHamiltonianSystem, Optional[float]): The Fermion-Spin system and
for constrained broadenings
the trotter timestep that
fixes the correct broadening
prefactor.
"""
spectrum = coupling_to_spectral_function(original_system, frequencies)
new_system = fermions.FermionHamiltonianSystem(original_system.number_fermionic_modes()[0])
original_system = original_system.system()
for index in original_system.keys():
# Filter out all indices where all boson indices in Mixed
# indices are empty, in other words the system terms
if all(str(x) == "I" for x in index.bosons()):
new_system.set(index.fermions()[0], original_system.get(index))
(
fermion_boson_model,
prefactor,
) = self.fit_fermion_boson_system_to_spectral_function(
fermionic_system=new_system,
spectral_function=spectrum,
)
system = self._create_fermion_system_from_fitting(fermion_boson_system=fermion_boson_model)
if prefactor is not None:
execution_time = noise_app.system_bath_execution_time(system, 1.0, device)
# The effective broadening scales with execution time.
# It increases the effective broadening and therefore rescales the prefactor
trotter_timestep = 1 / (prefactor / execution_time)
else:
trotter_timestep = None
return (system, trotter_timestep)
[docs]
def fit_boson_bath_to_fermion_bath(
self,
original_system: mixed_systems.MixedHamiltonianSystem,
temperature: float,
number_frequency_points: int,
) -> Tuple[Tuple[mixed_systems.MixedLindbladOpenSystem, Optional[float]], List[float]]:
"""Fits a model Spin-Boson system to an original Spin-Fermion system.
Args:
original_system (MixedHamiltonianSystem): The spin-fermion system bath model
the effective model is fitted to
temperature (float): The fermionic bath temperature assumed to obtain the
spectral function
number_frequency_points (int): The number of frequency points used to discretize
the spectral function obtained from the Fermion bath
Returns:
((MixedLindbladOpenSystem, Optional[float]), List[float]): The Spin-Boson system,
the fitted prefactor (for constrained
broadenings) and the corresponding frequencies.
"""
(spectrum, frequencies) = calculate_excitation_spectral_function(
original_system, temperature, number_frequency_points
)
new_system = spins.SpinHamiltonianSystem(original_system.number_spins()[0])
for index in original_system.keys():
# Filter out all indices where all boson indices in Mixed
# indices are empty, in other words the system terms
if all(str(x) == "I" for x in index.fermions()):
# new index that contains the spin part of the old index as
# first subsystem index.
new_system.set(index.spins()[0], original_system.get(index))
return (
self.fit_boson_bath_to_spectral_function(
new_system,
spectral_function=spectrum,
),
frequencies,
)
[docs]
def fit_spin_bath_to_fermion_bath(
self,
noise_app: Any,
original_system: mixed_systems.MixedHamiltonianSystem,
temperature: float,
number_frequency_points: int,
device: Any,
) -> Tuple[Tuple[mixed_systems.MixedHamiltonianSystem, Optional[float]], List[float]]:
"""Fits a model Spin-Spin system to an original Spin-Fermion system.
Args:
noise_app (HqsNoiseApp): The noise app used to construct simulation programs
from the result of the fitting
original_system (MixedHamiltonianSystem): The spin-fermion system bath model
the effective model is fitted to
temperature (float): The fermionic bath temperature assumed to obtain the
spectral function
number_frequency_points (int): The number of frequency points used to discretize
the spectral function opained from the Fermion bath
device (Device): The qoqo device on which the spin-spin system bath model
should be simulated. Used to calculate the circuit depth
necessary to determine the Trotter timestep
Returns:
((MixedHamiltonianSystem, Optional[float]), List[float]): The Spin-Spin system,
the trotter timestep that
fixes the correct broadening
prefactor for constrained broadenings,
and the corresponding frequencies.
"""
(spectrum, frequencies) = calculate_excitation_spectral_function(
original_system, temperature, number_frequency_points
)
new_system = spins.SpinHamiltonianSystem(original_system.number_spins()[0])
for index in original_system.keys():
# Filter out all indices where all boson indices in Mixed
# indices are empty, in other words the system terms
if all(str(x) == "I" for x in index.fermions()):
new_system.set(index.spins()[0], original_system.get(index))
return (
self.fit_spin_bath_to_spectral_function(
noise_app, new_system, spectral_function=spectrum, device=device
),
frequencies,
)
[docs]
def fit_boson_bath_to_spectral_function(
self,
spin_system: spins.SpinHamiltonianSystem,
spectral_function: SpinBRNoiseOperator,
) -> Tuple[mixed_systems.MixedLindbladOpenSystem, Optional[float]]:
"""Fits a model Spin-Boson system to a Spin spectral function.
Args:
spin_system (SpinHamiltonianSystem): The coherent spin system for which the
spectral function determines the noise
spectral_function (SpinBRNoiseOperator): The Bloch-Redfield type spectral function
that determines the decoherence in the
open system the effective system is fitted to.
Returns:
(MixedLindbladOpenSystem, Optional[float]): The Spin-Boson system and
for constrained broadenings
the fitted prefactor.
"""
if self.fitting_window is not None:
new_frequencies = np.linspace(
self.fitting_window[0], self.fitting_window[1], self.fitting_window[2]
)
spectral_function = spectral_function.resample(new_frequencies)
if self.broadening_constraint is not None:
fitter_initial_guess = BathFitter(
number_boson_modes=self.number_boson_modes,
spins_per_bosonic_mode=self.spins_per_bosonic_mode,
broadening_constraint=None,
background_broadening_ratio=self.background_broadening_ratio,
minimum_eigenfrequencies=self.minimum_eigenfrequencies,
maximum_eigenfrequencies=self.maximum_eigenfrequencies,
fitting_window=self.fitting_window,
coupling_types=self.coupling_types,
coupling_indices=self.coupling_indices,
max_fitting_iterations=self.max_fitting_iterations,
max_fitting_error=self.max_fitting_error,
)
(_, _, initial_guess) = fitter_initial_guess._internal_fitting(
spectrum=spectral_function,
coherent_system=spin_system,
is_spin=True,
initial_parameter_guess=None,
)
fitted_broadening = np.zeros(self.number_boson_modes)
for i in range(self.number_boson_modes):
fitted_broadening[i] = initial_guess.get(f"gam_{i}")
fitted_couplings = np.zeros(
(self.number_boson_modes, 3 * spin_system.number_spins()), dtype=complex
)
coupling_types = self._get_coupling_types(
number_particles=spin_system.number_spins(), is_spin=True
)
for j in range(spin_system.number_spins()):
for m in range(self.number_boson_modes):
for added, coupling in enumerate(["X", "Y", "Z"]):
if coupling in coupling_types.get((j, m), []):
fitted_couplings[m, 3 * j + added] = initial_guess.get(
f"g_{m}_{j}{coupling}_re", 0.0
) + 1j * initial_guess.get(f"g_{m}_{j}{coupling}_im", 0.0)
fitted_frequencies = np.zeros(self.number_boson_modes)
for i in range(self.number_boson_modes):
fitted_frequencies[i] = initial_guess.get(f"x0_{i}")
# Sort both and remember sorting for both, then line them up
broadening_index_array_i = np.argsort(self.broadening_constraint)
broadening_index_array_j = np.argsort(fitted_broadening)
new_initial_broadenings = np.zeros(self.number_boson_modes)
new_initial_couplings = np.zeros(fitted_couplings.shape, dtype=complex)
new_initial_frequencies = np.zeros(self.number_boson_modes)
for new, old in zip(broadening_index_array_i, broadening_index_array_j):
new_initial_broadenings[new] = fitted_broadening[old]
new_initial_couplings[new, :] = fitted_couplings[old, :]
new_initial_frequencies[new] = fitted_frequencies[old]
initial_guess_dict = {
"couplings": new_initial_couplings.tolist(),
"bosonic_frequencies": new_initial_frequencies.tolist(),
"initial_broadenings": new_initial_broadenings.tolist(),
}
else:
initial_guess_dict = None
(spin_boson_model, prefactor, _) = self._internal_fitting(
spectrum=spectral_function,
coherent_system=spin_system,
is_spin=True,
initial_parameter_guess=initial_guess_dict,
)
return (spin_boson_model, prefactor)
[docs]
def fit_fermion_boson_system_to_spectral_function(
self,
fermionic_system: fermions.FermionHamiltonianSystem,
spectral_function: FermionBRNoiseOperator,
) -> Tuple[mixed_systems.MixedLindbladOpenSystem, Optional[float]]:
"""Fits a model Fermion-Boson system to a Fermion spectral function.
Args:
fermion_system (FermionHamiltonianSystem): The coherent fermion system for which the
spectral function determines the noise
spectral_function (FermionBRNoiseOperator): The Bloch-Redfield type spectral function
that determines the decoherence in the
open system the effective system is fitted to.
Returns:
(MixedLindbladOpenSystem, Optional[float]): The Fermion-Boson system and
for constrained broadenings
the fitted prefactor.
"""
if self.fitting_window is not None:
new_frequencies = np.linspace(
self.fitting_window[0], self.fitting_window[1], self.fitting_window[2]
)
spectral_function = spectral_function.resample(new_frequencies)
if self.broadening_constraint is not None:
fitter_initial_guess = BathFitter(
number_boson_modes=self.number_boson_modes,
spins_per_bosonic_mode=self.spins_per_bosonic_mode,
broadening_constraint=None,
background_broadening_ratio=self.background_broadening_ratio,
minimum_eigenfrequencies=self.minimum_eigenfrequencies,
maximum_eigenfrequencies=self.maximum_eigenfrequencies,
fitting_window=self.fitting_window,
coupling_types=self.coupling_types,
coupling_indices=self.coupling_indices,
max_fitting_iterations=self.max_fitting_iterations,
max_fitting_error=self.max_fitting_error,
)
(_, _, initial_guess) = fitter_initial_guess._internal_fitting(
spectrum=spectral_function,
coherent_system=fermionic_system,
is_spin=False,
initial_parameter_guess=None,
)
fitted_broadening = np.zeros(self.number_boson_modes)
for i in range(self.number_boson_modes):
fitted_broadening[i] = initial_guess.get(f"gam_{i}")
fitted_couplings = np.zeros(
(self.number_boson_modes, fermionic_system.number_modes() ** 2),
dtype=complex,
)
for j in range(fermionic_system.number_modes()):
for m in range(self.number_boson_modes):
for k in range(j, fermionic_system.number_modes()):
if self.coupling_indices is None or (j, k) in self.coupling_indices:
fitted_couplings[m, j * fermionic_system.number_modes() + k] = (
initial_guess.get(f"g_{m}_{j}_{k}_re", 0.0)
+ 1j * initial_guess.get(f"g_{m}_{j}_{k}_im", 0.0)
)
fitted_frequencies = np.zeros(self.number_boson_modes)
for i in range(self.number_boson_modes):
fitted_frequencies[i] = initial_guess.get(f"x0_{i}")
# Sort both and remember sorting for both, then line them up
broadening_index_array_i = np.argsort(self.broadening_constraint)
broadening_index_array_j = np.argsort(fitted_broadening)
new_initial_broadenings = np.zeros(self.number_boson_modes)
new_initial_couplings = np.zeros(fitted_couplings.shape, dtype=complex)
new_initial_frequencies = np.zeros(self.number_boson_modes)
for new, old in zip(broadening_index_array_i, broadening_index_array_j):
new_initial_broadenings[new] = fitted_broadening[old]
new_initial_couplings[new, :] = fitted_couplings[old, :]
new_initial_frequencies[new] = fitted_frequencies[old]
initial_guess_dict = {
"couplings": new_initial_couplings.tolist(),
"bosonic_frequencies": new_initial_frequencies.tolist(),
"initial_broadenings": new_initial_broadenings.tolist(),
}
else:
initial_guess_dict = None
(fermion_boson_model, prefactor, _) = self._internal_fitting(
spectrum=spectral_function,
coherent_system=fermionic_system,
is_spin=False,
initial_parameter_guess=initial_guess_dict,
)
return (fermion_boson_model, prefactor)
[docs]
def fit_spin_bath_to_spectral_function(
self,
noise_app: Any,
spin_system: spins.SpinHamiltonianSystem,
spectral_function: SpinBRNoiseOperator,
device: Any,
) -> Tuple[mixed_systems.MixedHamiltonianSystem, Optional[float]]:
"""Fits a model Spin-Spin system to a Spin spectral function.
Args:
noise_app (HqsNoiseApp): The noise app used to construct simulation programs
from the result of the fitting
spin_system (SpinHamiltonianSystem): The coherent spin system for which the
spectral function determines the noise
spectral_function (SpinBRNoiseOperator): The Bloch-Redfield type spectral function
that determines the decoherence in the
open system the effective system is fitted to.
device (Device): The qoqo device on which the spin-spin system bath model
should be simulated. Used to calculate the circuit depth
necessary to determine the Trotter timestep
Returns:
(MixedHamiltonianSystem, Optional[float]): The Spin-Spin system and
for constrained broadenings
the trotter timestep that
fixes the correct broadening
prefactor.
"""
(spin_boson_model, prefactor) = self.fit_boson_bath_to_spectral_function(
spin_system=spin_system,
spectral_function=spectral_function,
)
system, trotter_timestep = self.spin_bath_trotterstep_from_boson_bath(
noise_app=noise_app,
spin_boson_system=spin_boson_model,
prefactor=prefactor,
device=device,
)
return (system, trotter_timestep)
[docs]
def spin_bath_trotterstep_from_boson_bath(
self,
noise_app: Any,
spin_boson_system: mixed_systems.MixedLindbladOpenSystem,
prefactor: Optional[float],
device: Any,
) -> Tuple[mixed_systems.MixedHamiltonianSystem, Optional[float]]:
"""Creates Spin-Spin Hamiltonian and Trotterstep from Spin-Boson system.
Converts a Spin-Boson open system and a fitting prefactor
for the Bosonic mode broadening
to a Spin-Spin system and a Trotter timestep
that can be used to construct a quantum circuit
to simulate the Spin part of the Spin-Boson circuit with an effective bath.
Args:
noise_app (HqsNoiseApp): The noise app used to construct simulation programs
from the result of the fitting
spin_boson_system (MixedLindbladOpenSystem): The Spin-Boson open system
that is transformed to a Spin-Spin system
prefactor (Optional[float]): The fitted prefactor of the bosonic mode broadenings
Is used to derive the Trotter timestep
device (Device): The qoqo device on which the spin-spin system bath model
should be simulated. Used to calculate the circuit depth
necessary to determine the Trotter timestep
Returns:
(MixedHamiltonianSystem, Optional[float]): The Spin-Spin system and
for constrained broadenings
the trotter timestep that
fixes the correct broadening
prefactor.
"""
system = self._create_spin_system_from_fitting(spin_boson_system=spin_boson_system)
if prefactor is not None:
execution_time = noise_app.system_bath_execution_time(system, 1.0, device)
# The effective broadening scales with execution time.
# It increases the effective broadening and therefore rescales the prefactor
trotter_timestep = 1 / (prefactor / execution_time)
else:
trotter_timestep = None
return (system, trotter_timestep)
def _create_spin_system_from_fitting(
self,
spin_boson_system: mixed_systems.MixedLindbladOpenSystem,
) -> mixed_systems.MixedHamiltonianSystem:
"""Create a spin-spin MixedHamiltonianSystem from a fitted spin-boson system-bath.
Args:
spin_boson_system (MixedLindbladOpenSystem): The coupling system.
Returns:
MixedHamiltonianSystem: The MixedHamiltonianSystem with the fitted bath.
Raises:
ValueError: Bosonic coupling operator must be given by single annihilator.
ValueError: Bath coupling of the given form is not supported.
"""
spin_boson_system = spin_boson_system.system()
new_system = mixed_systems.MixedHamiltonianSystem(
[
spin_boson_system.number_spins()[0],
spin_boson_system.number_bosonic_modes()[0] * self.spins_per_bosonic_mode,
],
[],
[],
)
for index in spin_boson_system.keys():
value = spin_boson_system.get(index)
boson_op = index.bosons()[0]
spin_op = index.spins()[0]
# Adding pure system terms
if boson_op == bosons.BosonProduct([], []):
new_system.set(
mixed_systems.HermitianMixedProduct([spin_op, spins.PauliProduct()], [], []),
value,
)
# Adding pure bath terms
elif spin_op == spins.PauliProduct():
try:
boson_index = boson_op.annihilators()[0]
except IndexError:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
try:
boson_index_2 = boson_op.creators()[0]
except IndexError:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
if boson_index != boson_index_2:
raise ValueError(
"Bath coupling of form c_{}a_{} is not supported".format(
boson_index_2, boson_index
)
) from None
for spin_offset in range(self.spins_per_bosonic_mode):
spin_bath_index = boson_index * self.spins_per_bosonic_mode + spin_offset
spin_index = spins.PauliProduct().z(spin_bath_index)
set_index = mixed_systems.HermitianMixedProduct([spin_op, spin_index], [], [])
new_system.set(set_index, value * 0.5)
# Adding system-bath terms
else:
try:
boson_index = boson_op.annihilators()[0]
except IndexError:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
if len(boson_op.creators()) > 0:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
for spin_offset in range(self.spins_per_bosonic_mode):
spin_bath_index = boson_index * self.spins_per_bosonic_mode + spin_offset
spin_index = spins.PauliProduct().x(spin_bath_index)
set_index = mixed_systems.HermitianMixedProduct([spin_op, spin_index], [], [])
new_system.set(set_index, value / np.sqrt(self.spins_per_bosonic_mode))
return new_system
def _create_fermion_system_from_fitting(
self,
fermion_boson_system: mixed_systems.MixedLindbladOpenSystem,
) -> mixed_systems.MixedHamiltonianSystem:
"""Create a fermion-spin MixedHamiltonianSystem from a fitted fermion-boson system-bath.
Args:
fermion_boson_system (MixedLindbladOpenSystem): The coupling system.
Returns:
MixedHamiltonianSystem: The MixedHamiltonianSystem with the fitted bath.
Raises:
ValueError: Bosonic coupling operator must be given by single annihilator.
ValueError: Bath coupling of the given form is not supported.
"""
fermion_boson_system = fermion_boson_system.system()
new_system = mixed_systems.MixedHamiltonianSystem(
[
fermion_boson_system.number_bosonic_modes()[0] * self.spins_per_bosonic_mode,
],
[],
[
fermion_boson_system.number_fermionic_modes()[0],
],
)
for index in fermion_boson_system.keys():
value = fermion_boson_system.get(index)
boson_op = index.bosons()[0]
fermion_op = index.fermions()[0]
# Adding pure system terms
if boson_op == bosons.BosonProduct([], []):
new_system.set(
mixed_systems.HermitianMixedProduct([spins.PauliProduct()], [], [fermion_op]),
value,
)
# Adding pure bath terms
elif fermion_op == fermions.FermionProduct([], []):
try:
boson_index = boson_op.annihilators()[0]
except IndexError:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
try:
boson_index_2 = boson_op.creators()[0]
except IndexError:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
if boson_index != boson_index_2:
raise ValueError(
"Bath coupling of form c_{}a_{} is not supported".format(
boson_index_2, boson_index
)
) from None
for spin_offset in range(self.spins_per_bosonic_mode):
spin_bath_index = boson_index * self.spins_per_bosonic_mode + spin_offset
spin_index = spins.PauliProduct().z(spin_bath_index)
set_index = mixed_systems.HermitianMixedProduct([spin_index], [], [fermion_op])
new_system.set(set_index, value * 0.5)
# Adding system-bath terms
else:
try:
boson_index = boson_op.annihilators()[0]
except IndexError:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
if len(boson_op.creators()) > 0:
raise ValueError(
"Bosonic coupling operator must be given by single annihilator"
) from None
for spin_offset in range(self.spins_per_bosonic_mode):
spin_bath_index = boson_index * self.spins_per_bosonic_mode + spin_offset
spin_index = spins.PauliProduct().x(spin_bath_index)
set_index = mixed_systems.HermitianMixedProduct([spin_index], [], [fermion_op])
new_system.set(set_index, value / np.sqrt(self.spins_per_bosonic_mode))
return new_system
def _internal_fitting(
self,
spectrum: Union[SpinBRNoiseOperator, FermionBRNoiseOperator],
coherent_system: Union[fermions.FermionHamiltonianSystem, spins.SpinHamiltonianSystem],
is_spin: bool,
initial_parameter_guess: Optional[Dict[str, List[float]]],
) -> Tuple[mixed_systems.MixedLindbladOpenSystem, Optional[float], Parameters]:
"""Fits a system bath model to a spectral function of Spin-Boson model.
Args:
spectrum (Union[SpinBRNoiseOperator, FermionBRNoiseOperator]): The spectrum that the
System-Bath model is
fitted to.
coherent_system (Union[FermionHamiltonianSystem, SpinHamiltonianSystem]): The spin
system or
fermionic
system.
initial_parameter_guess (Optional[Dict[str, List[float]]]): An initial guess
for the couplings.
Returns:
Tuple[mixed_systems.MixedLindbladOpenSystem, Optional[float], Parameters]
Raises:
ValueError: The number of broadenings is
not commensurate with the number of boson modes.
ValueError: Broadening input not recognised.
"""
number_particles = (
coherent_system.number_spins() if is_spin else coherent_system.number_modes()
)
minimum_frequency_diff = np.max(spectrum.frequencies())
for i in range(len(spectrum.frequencies()) - 1):
if spectrum.frequencies()[i + 1] - spectrum.frequencies()[i] < minimum_frequency_diff:
minimum_frequency_diff = spectrum.frequencies()[i + 1] - spectrum.frequencies()[i]
if 5 * self.number_boson_modes > len(spectrum.frequencies()):
warnings.warn(
"There are less than 5 points per boson mode in the spectrum frequencies. This may"
+ " not be enough points and could lead to fitting issues. Please consider either "
+ "providing more data, or using the `resample` method of the BRNoiseOperator."
)
# Prepare input guess:
if initial_parameter_guess is None:
(
couplings,
bosonic_frequencies,
initial_broadenings_tmp,
) = self._prepare_initial_guess_bugfix(
spectrum=spectrum,
number_particles=number_particles,
is_spin=is_spin,
)
initial_broadenings = [initial_broadenings_tmp for _ in range(self.number_boson_modes)]
else:
(couplings, bosonic_frequencies, initial_broadenings) = (
np.array(initial_parameter_guess["couplings"]),
initial_parameter_guess["bosonic_frequencies"],
initial_parameter_guess["initial_broadenings"],
)
# Construct input to lmfit
fit_params = Parameters()
# If required, constrain the broadening widths to be equal.
if self.broadening_constraint is not None:
average_constraint = np.mean(self.broadening_constraint)
average_width_of_initial_guess = np.mean(initial_broadenings) / 2
prefactor_guess = average_width_of_initial_guess / average_constraint
fit_params.add("prefactor", value=prefactor_guess, min=1e-10, vary=True)
else:
pass
for m in range(self.number_boson_modes):
# Specify the centres (constrained to the available frequencies)
fit_params.add(
f"x0_{m}",
value=bosonic_frequencies[m],
min=self.minimum_eigenfrequencies,
max=self.maximum_eigenfrequencies,
)
# Specify the broadenings, can be default/None, a float for all or a list
if self.broadening_constraint is None:
if initial_broadenings[m] < minimum_frequency_diff:
warnings.warn(
"The smallest frequency difference in the spectrum is larger than the "
+ f"broadening for the boson mode {m}. This may lead to fitting issues. "
+ "Please providing more data, or using the `resample` method of the "
+ "BRNoiseOperator."
)
fit_params.add(
f"gam_{m}",
value=initial_broadenings[m] / 2,
min=1e-8,
vary=True,
)
# No maximum, must be positive
elif isinstance(self.broadening_constraint, list) or isinstance(
self.broadening_constraint, np.ndarray
):
if len(self.broadening_constraint) == self.number_boson_modes:
if self.broadening_constraint[m] * prefactor_guess < minimum_frequency_diff:
warnings.warn(
"The smallest frequency difference in the spectrum is larger than the "
+ f"broadening constraint for the boson mode {m}. This may lead to "
+ "fitting issues. Please providing more data, or using the `resample`"
+ " method of the BRNoiseOperator."
)
fit_params.add(
f"gam_{m}",
expr=f"{self.broadening_constraint[m]} * prefactor",
)
# No maximum, must be positive
else:
raise ValueError(
f"The number of broadenings {len(self.broadening_constraint)} is not "
+ "commensurate with the number of boson modes"
+ f" {self.number_boson_modes}."
)
else:
raise ValueError("Broadening input not recognised")
# Set all the initial couplings as derived from our guess
# Can be negative, zero, positive
# if is_spin there are 3 * number_particles free couplings
# 3 coupling operators for each spin mode
if is_spin:
coupling_operators = ["X", "Y", "Z"]
coupling_types = self._get_coupling_types(
number_particles=number_particles,
is_spin=is_spin,
)
for j in range(number_particles):
for added, coup in enumerate(coupling_operators):
if coup in coupling_types.get((j, m), []):
fit_params.add(
f"g_{m}_{j}{coup}_re",
value=couplings[m, 3 * j + added].real,
vary=True,
)
if not np.isclose(
couplings[m, 3 * j + added].imag,
0.0,
rtol=1e-09,
atol=1e-09,
):
fit_params.add(
f"g_{m}_{j}{coup}_im",
value=couplings[m, 3 * j + added].imag,
vary=True,
)
# if it is not a spin system there are number_particles ** 2
# free couplings (for each combination j `c` k `a`)
else:
for j in range(number_particles):
for k in range(j, number_particles):
if self.coupling_indices is None or (j, k) in self.coupling_indices:
fit_params.add(
f"g_{m}_{j}_{k}_re",
value=couplings[m, j * number_particles + k].real,
vary=True,
)
if not np.isclose(
couplings[m, j * number_particles + k].imag,
0.0,
rtol=1e-09,
atol=1e-09,
):
fit_params.add(
f"g_{m}_{j}_{k}_im",
value=couplings[m, j * number_particles + k].imag,
vary=True,
)
fitting_found = False
for opt_counter in range(self._max_fitting_iterations):
if not fitting_found:
out = minimize(
self._objective,
fit_params,
kws={
"spectrum": spectrum,
"number_particles": number_particles,
"number_boson_modes": self.number_boson_modes,
"background_broadening_ratio": self.background_broadening_ratio,
"is_spin": is_spin,
},
)
if is_spin:
empty_spectrum = SpinBRNoiseOperator(spectrum.frequencies())
else:
empty_spectrum = FermionBRNoiseOperator(spectrum.frequencies())
helper_resid = self._objective(
out.params,
spectrum=empty_spectrum,
number_particles=number_particles,
number_boson_modes=self.number_boson_modes,
background_broadening_ratio=self.background_broadening_ratio,
is_spin=is_spin,
)
ratio_residuals_values = np.sum(np.array(out.residual) ** 2) / np.sum(
np.array(helper_resid) ** 2
)
if out.success and ratio_residuals_values < self._max_fitting_error:
fitting_found = True
break
else:
warnings.warn(
(
"Bad fit found. Ratio of deviation in spectrum squared to spectrum"
+ f" squared: {ratio_residuals_values}."
+ f" Trying again. Try {opt_counter}"
+ f" of {self._max_fitting_iterations}"
)
)
if not fitting_found:
raise RuntimeError(
(
"Could not find a good fit to the data."
+ " Ratio of deviation in spectrum squared to spectrum"
+ " squared always larger than"
+ f" {self._max_fitting_error}."
)
)
if is_spin:
result_model = self._create_spin_model(out.params, coherent_system)
else:
result_model = self._create_fermionic_model(out.params, coherent_system)
# Caching the broadenings that where fitted.
self._set_cached_broadenings(out.params)
return (result_model, out.params.get("prefactor", None), out.params)
def _objective(
self,
params: Parameters,
spectrum: Union[SpinBRNoiseOperator, FermionBRNoiseOperator],
number_particles: int,
number_boson_modes: int,
background_broadening_ratio: float,
is_spin: bool,
) -> np.ndarray:
"""Calculate total residual for fits of Gaussians to several data sets.
Args:
params (Parameters): The parameters of the model.
spectrum (Union[SpinBRNoiseOperator, FermionBRNoiseOperator]): Spin/Fermion spectrum
of the problem.
number_particles (int): The number of spins or fermions in the model.
number_boson_modes (int): The number of bosons in the model.
background_broadening_ratio (float): Adds a constant background offset
to the diagonal spectral functions when fitting.
Given as a ratio of the average broadening.
Returns:
numpy.ndarray: residual calculated for the fits.
"""
broadenings: List[float] = []
for m in range(number_boson_modes):
broadenings.append(params[f"gam_{m}"])
background = background_broadening_ratio * np.mean(broadenings)
model = self._params_to_spectrum(
params,
spectrum.frequencies(),
number_particles,
float(background),
is_spin,
)
resid = spectrum._optimization_difference(model, number_particles)
resid = np.array(resid)
return resid.view(float)
def _params_to_spectrum(
self,
params: Parameters,
frequencies: np.ndarray,
number_particles: int,
background: float,
is_spin: bool,
) -> Union[SpinBRNoiseOperator, FermionBRNoiseOperator]:
"""Construct model and extract spectrum.
Args:
params (Parameters): The parameters of the model.
frequencies (numpy.ndarray): The frequencies.
number_particles (int): The number of spins in the model.
number_boson_modes (int): The number of bosons in the model.
background (float): The constant background offset
added to diagonal spectral functions.
Returns:
Union[SpinBRNoiseOperator, FermionBRNoiseOperator]: constructed from the provided parameters.
"""
if is_spin:
spin_system = spins.SpinHamiltonianSystem(number_particles)
model = self._create_spin_model(params, spin_system)
else:
fermion_system = fermions.FermionHamiltonianSystem(number_particles)
model = self._create_fermionic_model(params, fermion_system)
spec_function = coupling_to_spectral_function(
model, frequencies, background, allow_empty_coupling=True
)
return spec_function
def _prepare_initial_guess_bugfix(
self,
spectrum: Union[SpinBRNoiseOperator, FermionBRNoiseOperator],
number_particles: int,
is_spin: bool,
) -> Tuple[np.ndarray, List[float], float]:
"""Fits a system bath model to a spectral function of Fermion/Spin-Boson model.
Args:
spectrum (SpinBRNoiseOperator): The spectrum that the System-Bath model is fitted to.
number_particles (int): The number of the spins in the system.
Returns:
Tuple[numpy.ndarray, List[float], float]: The couplings of the resampled spectrum,
the frequencies of the input spectrum and their width.
"""
# Prepare input guess:
start_spectrum = (
spectrum.frequencies()[0]
if self.minimum_eigenfrequencies is None
else self.minimum_eigenfrequencies
)
end_spectrum = (
spectrum.frequencies()[-1]
# if self.maximum_eigenfrequencies is None
# else self.maximum_eigenfrequencies
)
end_points = np.linspace(
start_spectrum,
end_spectrum,
self.number_boson_modes + 1,
)
if len(end_points) > 1:
widths = end_points[1] - end_points[0]
else:
widths = 0.0
# Resample spectrum with new points
avg_spectra, mid_points = self._bath_fitter_resampling(
spectrum, number_particles, end_points, is_spin
)
# Get couplings initial guess
# For spin spectra there are 3*number_particles couplings (three coupling types per spin mode)
if is_spin:
couplings = np.zeros((self.number_boson_modes, 3 * number_particles), dtype=complex)
else:
couplings = np.zeros((self.number_boson_modes, number_particles**2), dtype=complex)
tol = 1e-8
for ii in range(self.number_boson_modes):
D, V = la.eigh(avg_spectra[:, :, ii])
D[D < tol] = tol
max_index = np.argmax(D)
couplings[ii, :] = np.sqrt(D[max_index]) * V[:, max_index].T.conj()
return (couplings, mid_points, widths)
def _bath_fitter_resampling(self, spectrum, number_particles, end_points, is_spin):
if is_spin:
coupling_operators = ["X", "Y", "Z"]
spectrum_array = np.zeros(
(
3 * number_particles,
3 * number_particles,
len(spectrum.frequencies()),
),
dtype=complex,
)
for ii in range(number_particles):
for added_left, coup_left in enumerate(coupling_operators):
for jj in range(number_particles):
for added_right, coup_right in enumerate(coupling_operators):
left = f"{coup_left}{ii}"
right = f"{coup_right}{jj}"
spectrum_array[3 * ii + added_left, 3 * jj + added_right, :] = (
spectrum.get((left, right))
)
spectrum_interp = pchip(
spectrum.frequencies(), spectrum_array, axis=2, extrapolate=None
)
avg_spectra = np.zeros(
(3 * number_particles, 3 * number_particles, self.number_boson_modes),
dtype=complex,
)
mid_points = []
for ii in range(self.number_boson_modes):
mid_points.append(0.5 * (end_points[ii] + end_points[ii + 1]))
avg_spectra[:, :, ii] = spectrum_interp.integrate(
end_points[ii], end_points[ii + 1], extrapolate=None
)
else:
spectrum_array = np.zeros(
(
number_particles**2,
number_particles**2,
len(spectrum.frequencies()),
),
dtype=complex,
)
for ii in range(number_particles):
for jj in range(ii, number_particles):
for kk in range(number_particles):
for ll in range(kk, number_particles):
if self.coupling_indices is None or (
(ii, jj) in self.coupling_indices
and (kk, ll) in self.coupling_indices
):
left = f"c{ii}a{jj}"
right = f"c{kk}a{ll}"
spectrum_array[
ii * number_particles + jj,
kk * number_particles + ll,
:,
] = spectrum.get((left, right))
spectrum_interp = pchip(
spectrum.frequencies(), spectrum_array, axis=2, extrapolate=None
)
avg_spectra = np.zeros(
(number_particles**2, number_particles**2, self.number_boson_modes),
dtype=complex,
)
mid_points = []
for ii in range(self.number_boson_modes):
mid_points.append(0.5 * (end_points[ii] + end_points[ii + 1]))
avg_spectra[:, :, ii] = spectrum_interp.integrate(
end_points[ii], end_points[ii + 1], extrapolate=None
)
return avg_spectra, mid_points
def _create_spin_model(
self,
params: Parameters,
spin_system: spins.SpinHamiltonianSystem,
) -> mixed_systems.MixedLindbladOpenSystem:
"""Construct model and extract spectrum.
Args:
params (Parameters): The parameters of the model.
spin_system (SpinHamiltonianSystem): The coherent part of the system.
Returns:
MixedLindbladOpenSystem: created based on the model parameters.
Raises:
ValueError: Coupling type not supported
"""
model = mixed_systems.MixedLindbladOpenSystem(
[spin_system.number_spins()], [self.number_boson_modes], []
)
for m in range(self.number_boson_modes):
index = mixed_systems.HermitianMixedProduct(
# Identity spin operator
[spins.PauliProduct()],
# Create a Boson occupation operator
[bosons.BosonProduct([m], [m])],
[],
)
model.system_set(index, params[f"x0_{m}"])
index = mixed_systems.MixedDecoherenceProduct(
# Identity spin operator
[spins.PauliProduct()],
# Create a Boson occupation operator
[bosons.BosonProduct([], [m])],
[],
)
model.noise_set((index, index), params[f"gam_{m}"])
coupling_types = self._get_coupling_types(
number_particles=spin_system.number_spins(), is_spin=True
)
for m in range(self.number_boson_modes):
for ii in range(spin_system.number_spins()):
for coup in coupling_types.get((ii, m), []):
index = mixed_systems.HermitianMixedProduct(
# Identity spin operator
[f"{ii}{coup}"],
# Create a Boson coupling operator (always a + a^dagger)
[bosons.BosonProduct([], [m])],
[],
)
value = params.get(f"g_{m}_{ii}{coup}_re", 0.0) + 1j * params.get(
f"g_{m}_{ii}{coup}_im", 0.0
)
if value != complex(0.0):
model.system_set(index, value)
for key in spin_system.keys():
index = mixed_systems.HermitianMixedProduct(
# Identity spin operator
[key],
# Create a Boson occupation operator
[bosons.BosonProduct([], [])],
[],
)
model.system_set(index, spin_system.get(key))
return model
def _create_fermionic_model(
self,
params: Parameters,
fermionic_system: fermions.FermionHamiltonianSystem,
) -> mixed_systems.MixedLindbladOpenSystem:
"""Construct model and extract spectrum.
Args:
params (Parameters): The parameters of the model.
fermionic_system (FermionHamiltonianSystem): The coherent part of the system.
Returns:
MixedLindbladOpenSystem: created based on the model parameters.
Raises:
ValueError: Coupling type not supported
"""
model = mixed_systems.MixedLindbladOpenSystem(
[], [self.number_boson_modes], [fermionic_system.number_modes()]
)
for m in range(self.number_boson_modes):
index = mixed_systems.HermitianMixedProduct(
[],
# Create a Boson occupation operator
[bosons.BosonProduct([m], [m])],
# Identity spin operator
[fermions.FermionProduct([], [])],
)
model.system_set(index, params[f"x0_{m}"])
index = mixed_systems.MixedDecoherenceProduct(
# Identity spin operator
[],
# Create a Boson occupation operator
[bosons.BosonProduct([], [m])],
# Identity spin operator
[fermions.FermionProduct([], [])],
)
model.noise_set((index, index), params[f"gam_{m}"])
for m in range(self.number_boson_modes):
for ii in range(fermionic_system.number_modes()):
for jj in range(ii, fermionic_system.number_modes()):
if self.coupling_indices is None or (ii, jj) in self.coupling_indices:
index = mixed_systems.HermitianMixedProduct(
[],
# Create a Boson coupling operator (always a + a^dagger)
[bosons.BosonProduct([], [m])],
[f"c{ii}a{jj}"],
)
value = params.get(f"g_{m}_{ii}_{jj}_re", 0.0) + 1j * params.get(
f"g_{m}_{ii}_{jj}_im", 0.0
)
if value != complex(0.0):
model.system_set(index, value)
for key in fermionic_system.keys():
index = mixed_systems.HermitianMixedProduct(
[],
# Create a Boson occupation operator
[bosons.BosonProduct([], [])],
[key],
)
model.system_set(index, fermionic_system.get(key))
model = model.truncate(1e-14)
return model