Generalize effective noise channels

pulser-simulation/pulser_simulation/simconfig.py

@@ -34,7 +34,9 @@
             " `pip install typing-extensions`."
         )
 
-NOISE_TYPES = Literal["doppler", "amplitude", "SPAM", "dephasing"]
+NOISE_TYPES = Literal[
+    "doppler", "amplitude", "SPAM", "dephasing", "depolarizing", "eff_noise"
+]
 MASS = 1.45e-25  # kg
 KB = 1.38e-23  # J/K
 KEFF = 8.7  # µm^-1
@@ -54,6 +56,13 @@ class SimConfig:
             noise types:
 
             - "dephasing": Random phase (Z) flip
+            - "depolarizing": Quantum noise where the state(rho) is
+              turned into a mixed state I/2 with probability p,
+              and left unchanged with probability 1-p.
+            - "eff_noise": General effective noise channel defined by
+              the set of collapse operators **eff_noise_opers**
+              and the corresponding probability distribution
+              **eff_noise_probs**.
             - "doppler": Local atom detuning due to finite speed of the
               atoms and Doppler effect with respect to laser frequency
             - "amplitude": Gaussian damping due to finite laser waist
@@ -86,6 +95,9 @@ class SimConfig:
     epsilon: float = 0.01
     epsilon_prime: float = 0.05
     dephasing_prob: float = 0.05
+    depolarizing_prob: float = 0.05
+    eff_noise_probs: list[float] = field(default_factory=list, repr=False)
+    eff_noise_opers: list[qutip.Qobj] = field(default_factory=list, repr=False)
     solver_options: Optional[qutip.Options] = None
     spam_dict: dict[str, float] = field(
         init=False, default_factory=dict, repr=False
@@ -113,6 +125,7 @@ def __post_init__(self) -> None:
         self._check_noise_types()
         self._check_spam_dict()
         self._calc_sigma_doppler()
+        self._check_eff_noise()
 
     def __str__(self, solver_options: bool = False) -> str:
         lines = [
@@ -125,13 +138,22 @@ def __str__(self, solver_options: bool = False) -> str:
             lines.append("Noise types:           " + ", ".join(self.noise))
         if "SPAM" in self.noise:
             lines.append(f"SPAM dictionary:       {self.spam_dict}")
+        if "eff_noise" in self.noise:
+            lines.append(
+                f"General noise distribution:       {self.eff_noise_probs}"
+            )
+            lines.append(
+                f"General noise operators:       {self.eff_noise_opers}"
+            )
         if "doppler" in self.noise:
             lines.append(f"Temperature:           {self.temperature*1.e6}µK")
             lines.append(f"Amplitude standard dev.:  {self.amp_sigma}")
         if "amplitude" in self.noise:
             lines.append(f"Laser waist:           {self.laser_waist}μm")
         if "dephasing" in self.noise:
             lines.append(f"Dephasing probability: {self.dephasing_prob}")
+        if "depolarizing" in self.noise:
+            lines.append(f"Depolarizing probability: {self.depolarizing_prob}")
         if solver_options:
             lines.append(
                 "Solver Options: \n" + f"{str(self.solver_options)[10:-1]}"
@@ -167,6 +189,20 @@ def _check_noise_types(self) -> None:
                     + "Valid noise types: "
                     + ", ".join(get_args(NOISE_TYPES))
                 )
+        dephasing_on = "dephasing" in self.noise
+        depolarizing_on = "depolarizing" in self.noise
+        eff_noise_on = "eff_noise" in self.noise
+        eff_noise_conflict = (
+            (dephasing_on and depolarizing_on)
+            or (depolarizing_on and eff_noise_on)
+            or (dephasing_on and eff_noise_on)
+        )
+        if eff_noise_conflict:
+            raise NotImplementedError(
+                "Depolarizing, dephasing and eff_noise channels"
+                "cannot be activated at the same time in"
+                " one simulation."
+            )
 
     def _calc_sigma_doppler(self) -> None:
         # sigma = keff Deltav, keff = 8.7mum^-1, Deltav = sqrt(kB T / m)
@@ -176,3 +212,69 @@ def _calc_sigma_doppler(self) -> None:
 
     def _change_attribute(self, attr_name: str, new_value: Any) -> None:
         object.__setattr__(self, attr_name, new_value)
+
+    def _check_eff_noise(self) -> None:
+        # Check the validity of the distribution of probability
+        if "eff_noise" in self.noise:
+            if len(self.eff_noise_opers) != len(self.eff_noise_probs):
+                raise ValueError(
+                    f"The operators list length({len(self.eff_noise_opers)}) "
+                    "and probabilities list length"
+                    f"({len(self.eff_noise_probs)}) must be equal."
+                )
+            if self.eff_noise_opers == [] or self.eff_noise_probs == []:
+                raise ValueError("Fill the general noise parameters.")
+
+            for prob in self.eff_noise_probs:
+                if not isinstance(prob, float):
+                    raise TypeError(
+                        "eff_noise_probs is a list of floats"
+                        f" it must not contain a {type(prob)}."
+                    )
+
+            prob_distr = np.array(self.eff_noise_probs)
+            lower_bound = np.any(prob_distr < 0.0)
+            upper_bound = np.any(prob_distr > 1.0)
+            sum_p = not np.isclose(sum(prob_distr), 1.0)
+
+            if sum_p or lower_bound or upper_bound:
+                raise ValueError(
+                    "The distribution given is not a probability distribution."
+                )
+            # Check the validity of operators
+            for operator in self.eff_noise_opers:
+                # type checking
+
+                if type(operator) != qutip.qobj.Qobj:
+                    raise TypeError(f"{operator} is not a Qobj.")
+                if operator.type != "oper":
+                    raise TypeError(
+                        "Operators are supposed to be of type oper."
+                    )
+                if operator.shape != (2, 2):
+                    raise NotImplementedError(
+                        "Operator's shape must be (2,2) "
+                        f"not {operator.shape}."
+                    )
+            # Identity position
+            identity = qutip.qeye(2)
+            if self.eff_noise_opers[0] != identity:
+                raise NotImplementedError(
+                    "You must put the identity matrix at the "
+                    "beginning of the operator list."
+                )
+            # Completeness relation checking
+            sum_op = qutip.Qobj(shape=(2, 2))
+            length = len(self.eff_noise_probs)
+            for i in range(length):
+                sum_op += (
+                    self.eff_noise_probs[i]
+                    * self.eff_noise_opers[i]
+                    * self.eff_noise_opers[i].dag()
+                )
+
+            if sum_op != identity:
+                raise ValueError(
+                    "The completeness relation is not verified."
+                    f" Ended up with {sum_op} instead of {identity}."
+                )

pulser-simulation/pulser_simulation/simulation.py

@@ -40,7 +40,14 @@
 )
 
 SUPPORTED_NOISE = {
-    "ising": {"dephasing", "doppler", "amplitude", "SPAM"},
+    "ising": {
+        "dephasing",
+        "doppler",
+        "amplitude",
+        "SPAM",
+        "depolarizing",
+        "eff_noise",
+    },
     "XY": {"SPAM"},
 }
 
@@ -193,6 +200,8 @@ def set_config(self, cfg: SimConfig) -> None:
             self._doppler_detune = {qid: 0.0 for qid in self._qid_index}
         # Noise, samples and Hamiltonian update routine
         self._construct_hamiltonian()
+
+        Kraus_ops = []
         if "dephasing" in self.config.noise:
             if self.basis_name == "digital" or self.basis_name == "all":
                 # Go back to previous config
@@ -212,16 +221,81 @@ def set_config(self, cfg: SimConfig) -> None:
                     stacklevel=2,
                 )
             k = np.sqrt(prob * (1 - prob) ** (n - 1))
-            self._collapse_ops = [
+
+            self._collapse_ops += [
                 np.sqrt((1 - prob) ** n)
                 * qutip.tensor([self.op_matrix["I"] for _ in range(n)])
             ]
+            Kraus_ops.append(k * qutip.sigmaz())
+
+        if "depolarizing" in self.config.noise:
+            if self.basis_name == "digital" or self.basis_name == "all":
+                # Go back to previous config
+                self.set_config(prev_config)
+                raise NotImplementedError(
+                    "Cannot include depolarizing "
+                    + "noise in digital- or all-basis."
+                )
+            # Probability of error occurrence
+
+            prob = self.config.depolarizing_prob / 4
+            n = self._size
+            if prob > 0.1 and n > 1:
+                warnings.warn(
+                    "The depolarizing model is a first-order approximation"
+                    f" in the depolarizing probability. p = {4*prob}"
+                    " is too large for realistic results.",
+                    stacklevel=2,
+                )
+
+            k = np.sqrt((prob) * (1 - 3 * prob) ** (n - 1))
+            self._collapse_ops += [
+                np.sqrt((1 - 3 * prob) ** n)
+                * qutip.tensor([self.op_matrix["I"] for _ in range(n)])
+            ]
+            Kraus_ops.append(k * qutip.sigmax())
+            Kraus_ops.append(k * qutip.sigmay())
+            Kraus_ops.append(k * qutip.sigmaz())
+
+        if "eff_noise" in self.config.noise:
+            if self.basis_name == "digital" or self.basis_name == "all":
+                # Go back to previous config
+                self.set_config(prev_config)
+                raise NotImplementedError(
+                    "Cannot include general "
+                    + "noise in digital- or all-basis."
+                )
+            # Probability distribution of error occurences
+            n = self._size
+            m = len(self.config.eff_noise_opers)
+            if n > 1:
+                for i in range(1, m):
+                    prob_i = self.config.eff_noise_probs[i]
+                    if prob_i > 0.1:
+                        warnings.warn(
+                            "The effective noise model is a first-order"
+                            " approximation in the noise probability."
+                            f"p={prob_i} is large for realistic results.",
+                            stacklevel=2,
+                        )
+                        break
+            # Deriving Kraus operators
+            prob_id = self.config.eff_noise_probs[0]
             self._collapse_ops += [
-                k
-                * (
-                    self.build_operator([("sigma_rr", [qid])])
-                    - self.build_operator([("sigma_gg", [qid])])
+                np.sqrt(prob_id**n)
+                * qutip.tensor([self.op_matrix["I"] for _ in range(n)])
+            ]
+            for i in range(1, m):
+                k = np.sqrt(
+                    self.config.eff_noise_probs[i] * prob_id ** (n - 1)
                 )
+                K_op = k * self.config.eff_noise_opers[i]
+                Kraus_ops.append(K_op)
+
+        # Building collapse operators
+        for operator in Kraus_ops:
+            self._collapse_ops += [
+                self.build_operator([(operator, [qid])])
                 for qid in self._qid_index
             ]
 
@@ -268,6 +342,11 @@ def add_config(self, config: SimConfig) -> None:
             param_dict["laser_waist"] = config.laser_waist
         if "dephasing" in diff_noise_set:
             param_dict["dephasing_prob"] = config.dephasing_prob
+        if "depolarizing" in diff_noise_set:
+            param_dict["depolarizing_prob"] = config.depolarizing_prob
+        if "eff_noise" in diff_noise_set:
+            param_dict["eff_noise_opers"] = config.eff_noise_opers
+            param_dict["eff_noise_probs"] = config.eff_noise_probs
         param_dict["temperature"] *= 1.0e6
         # update runs:
         param_dict["runs"] = config.runs
@@ -445,7 +524,9 @@ def draw(
     def _extract_samples(self) -> None:
         """Populates samples dictionary with every pulse in the sequence."""
         local_noises = True
-        if set(self.config.noise).issubset({"dephasing", "SPAM"}):
+        if set(self.config.noise).issubset(
+            {"dephasing", "SPAM", "depolarizing", "eff_noise"}
+        ):
             local_noises = "SPAM" in self.config.noise and self.config.eta > 0
         samples = self.samples_obj.to_nested_dict(all_local=local_noises)
 
@@ -875,7 +956,11 @@ def _run_solver() -> CoherentResults:
             else:
                 raise ValueError("`progress_bar` must be a bool.")
 
-            if "dephasing" in self.config.noise:
+            if (
+                "dephasing" in self.config.noise
+                or "depolarizing" in self.config.noise
+                or "eff_noise" in self.config.noise
+            ):
                 result = qutip.mesolve(
                     self._hamiltonian,
                     self.initial_state,
@@ -902,7 +987,9 @@ def _run_solver() -> CoherentResults:
             )
 
         # Check if noises ask for averaging over multiple runs:
-        if set(self.config.noise).issubset({"dephasing", "SPAM"}):
+        if set(self.config.noise).issubset(
+            {"dephasing", "SPAM", "depolarizing", "eff_noise"}
+        ):
             # If there is "SPAM", the preparation errors must be zero
             if "SPAM" not in self.config.noise or self.config.eta == 0:
                 return _run_solver()