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", "gen_noise"
+]
 MASS = 1.45e-25  # kg
 KB = 1.38e-23  # J/K
 KEFF = 8.7  # µm^-1
@@ -54,6 +56,11 @@ class SimConfig:
             noise types:
 
             - "dephasing": Random phase (Z) flip
+            - "depolarizing": model of decohering where the qubit undergoes a
+              bit-flip error, phase-flip error or both errors.
+            - "gen_noise": general model of noise channels where the qubit
+            undergoes a modification according to the specified noise
+            operators, it encapsulates dephasing and depolarizing channels.
             - "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 +93,13 @@ class SimConfig:
     epsilon: float = 0.01
     epsilon_prime: float = 0.05
     dephasing_prob: float = 0.05
+    depolarizing_prob: float = 0.05
+    gen_noise_probs: list[float] = field(
+        init=True, default_factory=list, repr=False
+    )
+    gen_noise_opers: list[qutip.Qobj] = field(
+        init=True, default_factory=list, repr=False
+    )
     solver_options: Optional[qutip.Options] = None
     spam_dict: dict[str, float] = field(
         init=False, default_factory=dict, repr=False
@@ -101,6 +115,12 @@ def __post_init__(self) -> None:
                 f"{self.amp_sigma}) must be greater than or equal"
                 " to 0. and smaller than 1."
             )
+        if len(self.gen_noise_opers) != len(self.gen_noise_probs):
+            raise ValueError(
+                f"The operators list length ({len(self.gen_noise_opers)}) "
+                f"and probabilities list length ({len(self.gen_noise_probs)})"
+                "are supposed to be equal."
+            )
         self._process_temperature()
         self._change_attribute(
             "spam_dict",
@@ -113,6 +133,7 @@ def __post_init__(self) -> None:
         self._check_noise_types()
         self._check_spam_dict()
         self._calc_sigma_doppler()
+        self._check_gen_noise()
 
     def __str__(self, solver_options: bool = False) -> str:
         lines = [
@@ -125,13 +146,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 "gen_noise" in self.noise:
+            lines.append(
+                f"General noise distribution:       {self.gen_noise_probs}"
+            )
+            lines.append(
+                f"General noise operators:       {self.gen_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]}"
@@ -176,3 +206,62 @@ 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_gen_noise(self) -> None:
+        # Check the validity of the distribution of probability
+        if "gen_noise" in self.noise:
+            if self.gen_noise_opers == [] or self.gen_noise_probs == []:
+                raise ValueError("Fill the general noise parameters.")
+
+            for prob in self.gen_noise_probs:
+                if not isinstance(prob, float):
+                    raise TypeError(
+                        "gen_noise_probs is a list of floats"
+                        f" it must not contain a {type(prob)}."
+                    )
+
+            prob_distr = np.array(self.gen_noise_probs)
+            lower_bound = np.any(prob_distr < 0.0)
+            upper_bound = np.any(prob_distr > 1.0)
+            sum_p = 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.gen_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.gen_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.gen_noise_probs)
+            for i in range(length):
+                sum_op += (
+                    self.gen_noise_probs[i]
+                    * self.gen_noise_opers[i]
+                    * self.gen_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",
+        "gen_noise",
+    },
     "XY": {"SPAM"},
 }
 
@@ -202,7 +209,7 @@ def set_config(self, cfg: SimConfig) -> None:
                 )
             # Probability of phase (Z) flip:
             # First order in prob
-            prob = self.config.dephasing_prob / 2
+            prob = self.config.dephasing_prob
             n = self._size
             if prob > 0.1 and n > 1:
                 warnings.warn(
@@ -225,6 +232,97 @@ def set_config(self, cfg: SimConfig) -> None:
                 for qid in self._qid_index
             ]
 
+        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
+            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 = {2*prob}"
+                    " is too large for realistic results.",
+                    stacklevel=2,
+                )
+
+            k = np.sqrt((prob / 3) * (1 - prob) ** (n - 1))
+            self._collapse_ops = [
+                np.sqrt((1 - prob) ** n)
+                * qutip.tensor([self.op_matrix["I"] for _ in range(n)])
+            ]
+            self._collapse_ops += [
+                k
+                * (
+                    self.build_operator([("sigma_rr", [qid])])
+                    - self.build_operator([("sigma_gg", [qid])])
+                )
+                for qid in self._qid_index
+            ]
+            self._collapse_ops += [
+                k
+                * (
+                    self.build_operator([("sigma_gr", [qid])]).dag()
+                    + self.build_operator([("sigma_gr", [qid])])
+                )
+                for qid in self._qid_index
+            ]
+            self._collapse_ops += [
+                1j
+                * k
+                * (
+                    self.build_operator([("sigma_gr", [qid])]).dag()
+                    - self.build_operator([("sigma_gr", [qid])])
+                )
+                for qid in self._qid_index
+            ]
+
+        if "gen_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.gen_noise_opers)
+            if n > 1:
+                for i in range(1, m):
+                    prob_i = self.config.gen_noise_probs[i]
+                    if prob_i > 0.1:
+                        warnings.warn(
+                            "The general noise model is a first-order"
+                            " approximation in the noise probability."
+                            f"p={2*prob_i} is large for realistic results.",
+                            stacklevel=2,
+                        )
+                        break
+            # Building collapse operators
+            prob_id = self.config.gen_noise_probs[0]
+            self._collapse_ops += [
+                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.gen_noise_probs[i] * prob_id ** (n - 1)
+                )
+                c_op = k * self.config.gen_noise_opers[i]
+                self._collapse_ops += [
+                    self.build_operator([(c_op, [qid])])
+                    for qid in self._qid_index
+                ]
+
     def add_config(self, config: SimConfig) -> None:
         """Updates the current configuration with parameters of another one.
 
@@ -268,6 +366,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 "gen_noise" in diff_noise_set:
+            param_dict["gen_noise_opers"] = config.gen_noise_opers
+            param_dict["gen_noise_probs"] = config.gen_noise_probs
         param_dict["temperature"] *= 1.0e6
         # update runs:
         param_dict["runs"] = config.runs
@@ -445,7 +548,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", "gen_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 +980,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 "gen_noise" in self.config.noise
+            ):
                 result = qutip.mesolve(
                     self._hamiltonian,
                     self.initial_state,
@@ -902,7 +1011,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", "gen_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()