Repository: PennyLaneAI/pennylane
Pull request #5842 (2024-06-12 03:56:43)

DefaultQutritMixed measurement error
**Context:**
This adds measurement error to DefaultQutritMixed device. It uses both amplitude damping and TritFlip error channels to add a realistic for qutrits measurement errror noise model. 

**Description of the Change:**
Added two keywords to qml.DefaultQutritMixed, damping_measurement_gammas and trit_flip_measurement_probs, if left as None, nothing changes. Both inputs, however, can have a list of three floats as inputs. If damping_measurement_gammas is inputted then [qml.QutritAmplitudeDamping](https://docs.pennylane.ai/en/latest/code/api/pennylane.QutritAmplitudeDamping.html) is applied with these inputs after measurement diagonalization. If trit_flip_measurement_probs is inputted then [qml.TritFlip](https://docs.pennylane.ai/en/latest/code/api/pennylane.TritFlip.html) is applied after measurement diagonalization. If both are inputted then amplitude damping is applied first.

**Benefits:**
Allows for easier and more clear measurement error applications, also allows for measurment error to be applied to non diagonal observables.

**Possible Drawbacks:**
This method is not fully generalizable, but a more general measurement error could be applied before measurement if that is necessary. The output depends on the diagonalizing gates of the observable.

**Related GitHub Issues:**
N/A


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5842?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.67%. Comparing base [(`4600c22`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/4600c22d98c6ee2a6543d28ff14c2ddb155e3a8e?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`6e715c8`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/6e715c838982c6d0335f774fe1835ed52efa4ef8?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> Report is 14 commits behind head on master.


<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##           master    #5842   +/-   ##
=======================================
  Coverage   99.67%   99.67%           
=======================================
  Files         417      420    +3     
  Lines       39978    39905   -73     
=======================================
- Hits        39847    39775   -72     
+ Misses        131      130    -1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5842?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

Gabriel-Bottrill: Hi @albi3ro, @glassnotes, @mudit2812, and @trbromley, sorry for the wait, this should be my last PR for a while, I implemented measurement/readout error. I added some comments on the implementation and a few questions and would love y'alls feedback as always!

P.S.
Is measurement_error ok or would you prefer readout_error?
Gabriel-Bottrill: Hi @trbromley, @glassnotes and I were discussing the best way for users to input error channels. There is the way that it is now where a user inputs a list or list of inputs for `qml.QutritAmplitudeDamping` and `qml.TritFlip`. The other way we were thinking was having a quantum function input that requires a wire as input that would be applied for each input. Finally, we talked about potentially having both inputs. Along with what the inputs are do you have an opinion about how a user should input them, especially in the case that it isn't just a quantum function.
trbromley: > Hi @trbromley, @glassnotes and I were discussing the best way for users to input error channels. There is the way that it is now where a user inputs a list or list of inputs for `qml.QutritAmplitudeDamping` and `qml.TritFlip`. The other way we were thinking was having a quantum function input that requires a wire as input that would be applied for each input. Finally, we talked about potentially having both inputs. Along with what the inputs are do you have an opinion about how a user should input them, especially in the case that it isn't just a quantum function? Do you think it would be best to have seperate inputs, a dictionary input, **kwargs like input, or something else?

@Gabriel-Bottrill here are you mainly thinking about readout/measurement noise? For our `NoiseModel` addition in Q2/Q3, we'll be adding support for:
```python
qml.NoiseModel({c0: n0, c1: n1, c2: n2}, meas={c3: n3}, t1=0.04)
```
where `n3` would be a quantum function that adds channel-based noise. You can check out our progress on `NoiseModel`s so far by going [here](https://docs.pennylane.ai/en/latest/code/qml_noise.html).
Gabriel-Bottrill: @trbromley, thank you, I will remove the change to **readout_probs then.
mlxd: Hi @Gabriel-Bottrill 
We are aiming to do a release right now, and the CI runners are needed to finalise the release.
I will temporarily close this PR to free the runners, so that we have the capacity to finish the release. Feel free to reopen the PR on Wednesday morning.
Gabriel-Bottrill: Hi @mlxd sorry about that, hope the release goes well!
mlxd: Hey @Gabriel-Bottrill no problem at all, and thanks again for understanding. The final PRs will be merged this evening, so you can reopen the PR as planned tomorrow morning without a problem. Feel free to reach out if there are any issues.
mudit2812: @Gabriel-Bottrill The issue with the failing tensorflow test is resolved. I've enabled auto merge and once you update your branch CI _should_ pass.
Gabriel-Bottrill: Thank you @mudit2812 

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 531aedd1c..e0c7fdff3 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -44,6 +44,14 @@
 * Observable validation for `default.qubit` is now based on execution mode (analytic vs. finite shots) and measurement type (sample measurement vs. state measurement).
   [(#5890)](https://github.com/PennyLaneAI/pennylane/pull/5890)
 
+<h4>Community contributions 🥳</h4>
+
+* `DefaultQutritMixed` readout error has been added using parameters `readout_relaxation_probs` and 
+  `readout_misclassification_probs` on the `default.qutrit.mixed` device. These parameters add a `~.QutritAmplitudeDamping`  and a `~.TritFlip` channel, respectively,
+  after measurement diagonalization. The amplitude damping error represents the potential for
+  relaxation to occur during longer measurements. The trit flip error represents misclassification during readout.
+  [(#5842)](https://github.com/PennyLaneAI/pennylane/pull/5842)
+
 <h3>Breaking changes 💔</h3>
 
 * ``qml.transforms.map_batch_transform`` has been removed, since transforms can be applied directly to a batch of tapes.
@@ -80,9 +88,10 @@
 
 This release contains contributions from (in alphabetical order):
 Guillermo Alonso,
-Ahmed Darwish,
 Astral Cai,
 Yushao Chen,
+Gabriel Bottrill,
+Ahmed Darwish,
 Lillian M. A. Frederiksen,
 Pietropaolo Frisoni,
 Christina Lee,
diff --git a/pennylane/devices/default_qutrit_mixed.py b/pennylane/devices/default_qutrit_mixed.py
index 0ade57c0c..cdad24647 100644
--- a/pennylane/devices/default_qutrit_mixed.py
+++ b/pennylane/devices/default_qutrit_mixed.py
@@ -14,7 +14,9 @@
 """The default.qutrit.mixed device is PennyLane's standard qutrit simulator for mixed-state
 computations."""
 import logging
+import warnings
 from dataclasses import replace
+from functools import partial
 from typing import Callable, Optional, Sequence, Tuple, Union
 
 import numpy as np
@@ -33,6 +35,7 @@ from .modifiers import simulator_tracking, single_tape_support
 from .preprocess import (
     decompose,
     no_sampling,
+    null_postprocessing,
     validate_device_wires,
     validate_measurements,
     validate_observables,
@@ -85,10 +88,67 @@ def accepted_sample_measurement(m: qml.measurements.MeasurementProcess) -> bool:
     return isinstance(m, qml.measurements.SampleMeasurement)
 
 
+@qml.transform
+def warn_readout_error_state(
+    tape: qml.tape.QuantumTape,
+) -> tuple[Sequence[qml.tape.QuantumTape], Callable]:
+    """If a measurement in the QNode is an analytic state or density_matrix, and a readout error
+    parameter is defined, warn that readout error will not be applied.
+
+    Args:
+        tape (QuantumTape, .QNode, Callable): a quantum circuit.
+
+    Returns:
+        qnode (pennylane.QNode) or quantum function (callable) or tuple[List[.QuantumTape], function]:
+        The unaltered input circuit.
+    """
+    if not tape.shots:
+        for m in tape.measurements:
+            if isinstance(m, qml.measurements.StateMP):
+                warnings.warn(f"Measurement {m} is not affected by readout error.")
+
+    return (tape,), null_postprocessing
+
+
+def get_readout_errors(readout_relaxation_probs, readout_misclassification_probs):
+    r"""Get the list of readout errors that should be applied to each measured wire.
+
+    Args:
+        readout_relaxation_probs (List[float]): Inputs for :class:`~.QutritAmplitudeDamping` channel
+            of the form :math:`[\gamma_{10}, \gamma_{20}, \gamma_{21}]`. This error models
+            amplitude damping associated with longer readout and varying relaxation times of
+            transmon-based qudits.
+        readout_misclassification_probs (List[float]): Inputs for :class:`~.TritFlip` channel
+            of the form :math:`[p_{01}, p_{02}, p_{12}]`. This error models misclassification events
+            in readout.
+
+    Returns:
+        readout_errors (List[Callable]): List of readout error channels that should be
+        applied to each measured wire.
+    """
+    measure_funcs = []
+    if readout_relaxation_probs is not None:
+        try:
+            with qml.queuing.QueuingManager.stop_recording():
+                qml.QutritAmplitudeDamping(*readout_relaxation_probs, wires=0)
+        except Exception as e:
+            raise qml.DeviceError("Applying damping readout error results in error:\n" + str(e))
+        measure_funcs.append(partial(qml.QutritAmplitudeDamping, *readout_relaxation_probs))
+    if readout_misclassification_probs is not None:
+        try:
+            with qml.queuing.QueuingManager.stop_recording():
+                qml.TritFlip(*readout_misclassification_probs, wires=0)
+        except Exception as e:
+            raise qml.DeviceError("Applying trit flip readout error results in error:\n" + str(e))
+        measure_funcs.append(partial(qml.TritFlip, *readout_misclassification_probs))
+
+    return None if len(measure_funcs) == 0 else measure_funcs
+
+
 @simulator_tracking
 @single_tape_support
 class DefaultQutritMixed(Device):
-    """A PennyLane Python-based device for mixed-state qutrit simulation.
+    r"""A PennyLane Python-based device for mixed-state qutrit simulation.
 
     Args:
         wires (int, Iterable[Number, str]): Number of wires present on the device, or iterable that
@@ -104,6 +164,12 @@ class DefaultQutritMixed(Device):
             If a ``jax.random.PRNGKey`` is passed as the seed, a JAX-specific sampling function using
             ``jax.random.choice`` and the ``PRNGKey`` will be used for sampling rather than
             ``numpy.random.default_rng``.
+        readout_relaxation_probs (List[float]): Input probabilities for relaxation errors implemented
+            with the :class:`~.QutritAmplitudeDamping` channel. The input defines the
+            channel's parameters :math:`[\gamma_{10}, \gamma_{20}, \gamma_{21}]`.
+        readout_misclassification_probs (List[float]):  Input probabilities for state readout
+            misclassification events implemented with the :class:`~.TritFlip` channel. The input defines the
+            channel's parameters :math:`[p_{01}, p_{02}, p_{12}]`.
 
     **Example:**
 
@@ -154,6 +220,33 @@ class DefaultQutritMixed(Device):
     >>> jax.grad(f)(jax.numpy.array(1.2))
     DeviceArray(-0.93203914, dtype=float32, weak_type=True)
 
+    .. details::
+        :title: Readout Error
+
+        ``DefaultQutritMixed`` includes readout error support. Two input arguments control
+        the parameters of error channels applied to each measured wire of the state after
+        it has been diagonalized for measurement:
+
+        * ``readout_relaxation_probs``:  Input parameters of a :class:`~.QutritAmplitudeDamping` channel.
+          This error models state relaxation error that occurs during readout of transmon-based qutrits.
+          The motivation for this readout error is described in [`1 <https://arxiv.org/abs/2003.03307>`_] (Sec II.A).
+        * ``readout_misclassification_probs``: Input parameters of a :class:`~.TritFlip` channel.
+          This error models misclassification events in readout. An example of this readout error
+          can be seen in [`2 <https://arxiv.org/abs/2309.11303>`_] (Fig 1a).
+
+        In the case that both parameters are defined, relaxation error is applied first then
+        misclassification error is applied.
+
+        .. note::
+            The readout errors will be applied to the state after it has been diagonalized for each
+            measurement. This may give different results depending on how the observable is defined.
+            This is because diagonalizing gates for the same observable may return eigenvalues in
+            different orders. For example, measuring :class:`~.THermitian` with a non-diagonal
+            GellMann matrix will result in a different measurement result then measuring the
+            equivalent :class:`~.GellMann` observable, as the THermitian eigenvalues are returned
+            in increasing order when explicitly diagonalized (i.e., ``[-1, 0, 1]``), while non-diagonal GellManns provided
+            in PennyLane have their eigenvalues hardcoded (i.e., ``[1, -1, 0]``).
+
     .. details::
         :title: Tracking
 
@@ -166,7 +259,6 @@ class DefaultQutritMixed(Device):
         * ``batches``: The number of times :meth:`~.execute` is called.
         * ``results``: The results of each call of :meth:`~.execute`
 
-
     """
 
     _device_options = ("rng", "prng_key")  # tuple of string names for all the device options.
@@ -177,11 +269,13 @@ class DefaultQutritMixed(Device):
         return "default.qutrit.mixed"
 
     @debug_logger_init
-    def __init__(
+    def __init__(  # pylint: disable=too-many-arguments
         self,
         wires=None,
         shots=None,
         seed="global",
+        readout_relaxation_probs=None,
+        readout_misclassification_probs=None,
     ) -> None:
         super().__init__(wires=wires, shots=shots)
         seed = np.random.randint(0, high=10000000) if seed == "global" else seed
@@ -193,6 +287,10 @@ class DefaultQutritMixed(Device):
             self._rng = np.random.default_rng(seed)
         self._debugger = None
 
+        self.readout_errors = get_readout_errors(
+            readout_relaxation_probs, readout_misclassification_probs
+        )
+
     @debug_logger
     def supports_derivatives(
         self,
@@ -284,6 +382,9 @@ class DefaultQutritMixed(Device):
         if config.gradient_method == "backprop":
             transform_program.add_transform(no_sampling, name="backprop + default.qutrit")
 
+        if self.readout_errors is not None:
+            transform_program.add_transform(warn_readout_error_state)
+
         return transform_program, config
 
     @debug_logger
@@ -292,7 +393,6 @@ class DefaultQutritMixed(Device):
         circuits: QuantumTape_or_Batch,
         execution_config: ExecutionConfig = DefaultExecutionConfig,
     ) -> Result_or_ResultBatch:
-
         interface = (
             execution_config.interface
             if execution_config.gradient_method in {"best", "backprop", None}
@@ -306,6 +406,7 @@ class DefaultQutritMixed(Device):
                 prng_key=self._prng_key,
                 debugger=self._debugger,
                 interface=interface,
+                readout_errors=self.readout_errors,
             )
             for c in circuits
         )
diff --git a/pennylane/devices/qutrit_mixed/measure.py b/pennylane/devices/qutrit_mixed/measure.py
index 9329bab6a..d93a0eeff 100644
--- a/pennylane/devices/qutrit_mixed/measure.py
+++ b/pennylane/devices/qutrit_mixed/measure.py
@@ -18,7 +18,7 @@ Code relevant for performing measurements on a qutrit mixed state.
 from string import ascii_letters as alphabet
 from typing import Callable
 
-from pennylane import math
+from pennylane import math, queuing
 from pennylane.measurements import (
     ExpectationMP,
     MeasurementProcess,
@@ -36,7 +36,10 @@ from .utils import QUDIT_DIM, get_num_wires, reshape_state_as_matrix
 
 
 def calculate_expval(
-    measurementprocess: ExpectationMP, state: TensorLike, is_state_batched: bool = False
+    measurementprocess: ExpectationMP,
+    state: TensorLike,
+    is_state_batched: bool = False,
+    readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Measure the expectation value of an observable.
 
@@ -44,18 +47,23 @@ def calculate_expval(
         measurementprocess (ExpectationMP): measurement process to apply to the state.
         state (TensorLike): the state to measure.
         is_state_batched (bool): whether the state is batched or not.
+        readout_errors (List[Callable]): List of chanels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         TensorLike: expectation value of observable wrt the state.
     """
-    probs = calculate_probability(measurementprocess, state, is_state_batched)
+    probs = calculate_probability(measurementprocess, state, is_state_batched, readout_errors)
     eigvals = math.asarray(measurementprocess.eigvals(), dtype="float64")
     # In case of broadcasting, `probs` has two axes and these are a matrix-vector products
     return math.dot(probs, eigvals)
 
 
 def calculate_reduced_density_matrix(
-    measurementprocess: StateMeasurement, state: TensorLike, is_state_batched: bool = False
+    measurementprocess: StateMeasurement,
+    state: TensorLike,
+    is_state_batched: bool = False,
+    _readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Get the state or reduced density matrix.
 
@@ -63,6 +71,8 @@ def calculate_reduced_density_matrix(
         measurementprocess (StateMeasurement): measurement to apply to the state.
         state (TensorLike): state to apply the measurement to.
         is_state_batched (bool): whether the state is batched or not.
+        _readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors. These are not applied on this type of measurement.
 
     Returns:
         TensorLike: state or reduced density matrix.
@@ -91,7 +101,10 @@ def calculate_reduced_density_matrix(
 
 
 def calculate_probability(
-    measurementprocess: StateMeasurement, state: TensorLike, is_state_batched: bool = False
+    measurementprocess: StateMeasurement,
+    state: TensorLike,
+    is_state_batched: bool = False,
+    readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Find the probability of measuring states.
 
@@ -99,6 +112,8 @@ def calculate_probability(
         measurementprocess (StateMeasurement): measurement to apply to the state.
         state (TensorLike): state to apply the measurement to.
         is_state_batched (bool): whether the state is batched or not.
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         TensorLike: the probability of the state being in each measurable state.
@@ -106,9 +121,15 @@ def calculate_probability(
     for op in measurementprocess.diagonalizing_gates():
         state = apply_operation(op, state, is_state_batched=is_state_batched)
 
+    wires = measurementprocess.wires
     num_state_wires = get_num_wires(state, is_state_batched)
     wire_order = Wires(range(num_state_wires))
-    wires = measurementprocess.wires
+
+    if readout_errors is not None:
+        with queuing.QueuingManager.stop_recording():
+            for wire in wires:
+                for m_error in readout_errors:
+                    state = apply_operation(m_error(wire), state, is_state_batched=is_state_batched)
 
     # probs are diagonal elements
     # stacking list since diagonal function axis selection parameter names
@@ -158,7 +179,10 @@ def calculate_probability(
 
 
 def calculate_variance(
-    measurementprocess: StateMeasurement, state: TensorLike, is_state_batched: bool = False
+    measurementprocess: StateMeasurement,
+    state: TensorLike,
+    is_state_batched: bool = False,
+    readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Find variance of observable.
 
@@ -166,11 +190,13 @@ def calculate_variance(
         measurementprocess (StateMeasurement): measurement to apply to the state.
         state (TensorLike): state to apply the measurement to.
         is_state_batched (bool): whether the state is batched or not.
+        readout_errors (List[Callable]): List of operators to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         TensorLike: the variance of the observable wrt the state.
     """
-    probs = calculate_probability(measurementprocess, state, is_state_batched)
+    probs = calculate_probability(measurementprocess, state, is_state_batched, readout_errors)
     eigvals = math.asarray(measurementprocess.eigvals(), dtype="float64")
     # In case of broadcasting, `probs` has two axes and these are a matrix-vector products
     return math.dot(probs, (eigvals**2)) - math.dot(probs, eigvals) ** 2
@@ -180,6 +206,7 @@ def calculate_expval_sum_of_terms(
     measurementprocess: ExpectationMP,
     state: TensorLike,
     is_state_batched: bool = False,
+    readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Measure the expectation value of the state when the measured observable is a ``Hamiltonian`` or ``Sum``
     and it must be backpropagation compatible.
@@ -188,6 +215,8 @@ def calculate_expval_sum_of_terms(
         measurementprocess (ExpectationMP): measurement process to apply to the state.
         state (TensorLike): the state to measure.
         is_state_batched (bool): whether the state is batched or not.
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         TensorLike: the expectation value of the sum of Hamiltonian observable wrt the state.
@@ -196,12 +225,23 @@ def calculate_expval_sum_of_terms(
         # Recursively call measure on each term, so that the best measurement method can
         # be used for each term
         return sum(
-            measure(ExpectationMP(term), state, is_state_batched=is_state_batched)
+            measure(
+                ExpectationMP(term),
+                state,
+                is_state_batched=is_state_batched,
+                readout_errors=readout_errors,
+            )
             for term in measurementprocess.obs
         )
     # else hamiltonian
     return sum(
-        c * measure(ExpectationMP(t), state, is_state_batched=is_state_batched)
+        c
+        * measure(
+            ExpectationMP(t),
+            state,
+            is_state_batched=is_state_batched,
+            readout_errors=readout_errors,
+        )
         for c, t in zip(*measurementprocess.obs.terms())
     )
 
@@ -209,7 +249,7 @@ def calculate_expval_sum_of_terms(
 # pylint: disable=too-many-return-statements
 def get_measurement_function(
     measurementprocess: MeasurementProcess,
-) -> Callable[[MeasurementProcess, TensorLike], TensorLike]:
+) -> Callable[[MeasurementProcess, TensorLike, bool, list[Callable]], TensorLike]:
     """Get the appropriate method for performing a measurement.
 
     Args:
@@ -238,7 +278,10 @@ def get_measurement_function(
 
 
 def measure(
-    measurementprocess: MeasurementProcess, state: TensorLike, is_state_batched: bool = False
+    measurementprocess: MeasurementProcess,
+    state: TensorLike,
+    is_state_batched: bool = False,
+    readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Apply a measurement process to a state.
 
@@ -246,8 +289,11 @@ def measure(
         measurementprocess (MeasurementProcess): measurement process to apply to the state.
         state (TensorLike): the state to measure.
         is_state_batched (bool): whether the state is batched or not.
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         Tensorlike: the result of the measurement process being applied to the state.
     """
-    return get_measurement_function(measurementprocess)(measurementprocess, state, is_state_batched)
+    measurement_function = get_measurement_function(measurementprocess)
+    return measurement_function(measurementprocess, state, is_state_batched, readout_errors)
diff --git a/pennylane/devices/qutrit_mixed/sampling.py b/pennylane/devices/qutrit_mixed/sampling.py
index 3c61614b0..c5395e80d 100644
--- a/pennylane/devices/qutrit_mixed/sampling.py
+++ b/pennylane/devices/qutrit_mixed/sampling.py
@@ -15,6 +15,7 @@
 Code relevant for sampling a qutrit mixed state.
 """
 import functools
+from typing import Callable
 
 import numpy as np
 
@@ -103,6 +104,7 @@ def _measure_with_samples_diagonalizing_gates(
     is_state_batched: bool = False,
     rng=None,
     prng_key=None,
+    readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Returns the samples of the measurement process performed on the given state,
     by rotating the state into the measurement basis using the diagonalizing gates
@@ -118,6 +120,8 @@ def _measure_with_samples_diagonalizing_gates(
             If no value is provided, a default RNG will be used.
         prng_key (Optional[jax.random.PRNGKey]): An optional ``jax.random.PRNGKey``. This is
             the key to the JAX pseudo random number generator. Only for simulation using JAX.
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         TensorLike[Any]: Sample measurement results
@@ -162,6 +166,7 @@ def _measure_with_samples_diagonalizing_gates(
                 wires=wires,
                 rng=rng,
                 prng_key=prng_key,
+                readout_errors=readout_errors,
             )
             processed_samples.append(_process_single_shot(samples))
 
@@ -174,6 +179,7 @@ def _measure_with_samples_diagonalizing_gates(
         wires=wires,
         rng=rng,
         prng_key=prng_key,
+        readout_errors=readout_errors,
     )
 
     return _process_single_shot(samples)
@@ -186,6 +192,7 @@ def _measure_sum_with_samples(
     is_state_batched: bool = False,
     rng=None,
     prng_key=None,
+    readout_errors: list[Callable] = None,
 ):
     """Compute expectation values of Sum or Hamiltonian Observables"""
     # mp.obs returns is the list of observables for Sum,
@@ -203,6 +210,7 @@ def _measure_sum_with_samples(
                     is_state_batched=is_state_batched,
                     rng=rng,
                     prng_key=prng_key,
+                    readout_errors=readout_errors,
                 )
             )
 
@@ -224,6 +232,7 @@ def _sample_state_jax(
     prng_key,
     is_state_batched: bool = False,
     wires=None,
+    readout_errors: list[Callable] = None,
 ) -> np.ndarray:
     """Returns a series of samples of a state for the JAX interface based on the PRNG.
 
@@ -234,6 +243,8 @@ def _sample_state_jax(
             the key to the JAX pseudo random number generator.
         is_state_batched (bool): whether the state is batched or not
         wires (Sequence[int]): The wires to sample
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         ndarray[int]: Sample values of the shape (shots, num_wires)
@@ -252,7 +263,7 @@ def _sample_state_jax(
     basis_states = np.arange(QUDIT_DIM**num_wires)
 
     with qml.queuing.QueuingManager.stop_recording():
-        probs = measure(qml.probs(wires=wires_to_sample), state, is_state_batched)
+        probs = measure(qml.probs(wires=wires_to_sample), state, is_state_batched, readout_errors)
 
     if is_state_batched:
         # Produce separate keys for each of the probabilities along the broadcasted axis
@@ -282,6 +293,7 @@ def sample_state(
     wires=None,
     rng=None,
     prng_key=None,
+    readout_errors: list[Callable] = None,
 ) -> np.ndarray:
     """Returns a series of computational basis samples of a state.
 
@@ -295,13 +307,20 @@ def sample_state(
             If no value is provided, a default RNG will be used
         prng_key (Optional[jax.random.PRNGKey]): An optional ``jax.random.PRNGKey``. This is
             the key to the JAX pseudo random number generator. Only for simulation using JAX.
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         ndarray[int]: Sample values of the shape (shots, num_wires)
     """
     if prng_key is not None:
         return _sample_state_jax(
-            state, shots, prng_key, is_state_batched=is_state_batched, wires=wires
+            state,
+            shots,
+            prng_key,
+            is_state_batched=is_state_batched,
+            wires=wires,
+            readout_errors=readout_errors,
         )
 
     rng = np.random.default_rng(rng)
@@ -314,7 +333,7 @@ def sample_state(
     basis_states = np.arange(QUDIT_DIM**num_wires)
 
     with qml.queuing.QueuingManager.stop_recording():
-        probs = measure(qml.probs(wires=wires_to_sample), state, is_state_batched)
+        probs = measure(qml.probs(wires=wires_to_sample), state, is_state_batched, readout_errors)
 
     if is_state_batched:
         # rng.choice doesn't support broadcasting
@@ -335,6 +354,7 @@ def measure_with_samples(
     is_state_batched: bool = False,
     rng=None,
     prng_key=None,
+    readout_errors: list[Callable] = None,
 ) -> TensorLike:
     """Returns the samples of the measurement process performed on the given state.
     This function assumes that the user-defined wire labels in the measurement process
@@ -350,6 +370,8 @@ def measure_with_samples(
             If no value is provided, a default RNG will be used.
         prng_key (Optional[jax.random.PRNGKey]): An optional ``jax.random.PRNGKey``. This is
             the key to the JAX pseudo random number generator. Only for simulation using JAX.
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         TensorLike[Any]: Sample measurement results
@@ -368,4 +390,5 @@ def measure_with_samples(
         is_state_batched=is_state_batched,
         rng=rng,
         prng_key=prng_key,
+        readout_errors=readout_errors,
     )
diff --git a/pennylane/devices/qutrit_mixed/simulate.py b/pennylane/devices/qutrit_mixed/simulate.py
index 54d4f35c2..69fd4c4c8 100644
--- a/pennylane/devices/qutrit_mixed/simulate.py
+++ b/pennylane/devices/qutrit_mixed/simulate.py
@@ -98,7 +98,9 @@ def get_final_state(circuit, debugger=None, interface=None, **kwargs):
     return state, is_state_batched
 
 
-def measure_final_state(circuit, state, is_state_batched, rng=None, prng_key=None) -> Result:
+def measure_final_state(  # pylint: disable=too-many-arguments
+    circuit, state, is_state_batched, rng=None, prng_key=None, readout_errors=None
+) -> Result:
     """
     Perform the measurements required by the circuit on the provided state.
 
@@ -115,6 +117,8 @@ def measure_final_state(circuit, state, is_state_batched, rng=None, prng_key=Non
             the key to the JAX pseudo random number generator. Only for simulation using JAX.
             If None, the default ``sample_state`` function and a ``numpy.random.default_rng``
             will be for sampling.
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         Tuple[TensorLike]: The measurement results
@@ -124,11 +128,12 @@ def measure_final_state(circuit, state, is_state_batched, rng=None, prng_key=Non
 
     if not circuit.shots:
         # analytic case
-
         if len(circuit.measurements) == 1:
-            return measure(circuit.measurements[0], state, is_state_batched)
+            return measure(circuit.measurements[0], state, is_state_batched, readout_errors)
 
-        return tuple(measure(mp, state, is_state_batched) for mp in circuit.measurements)
+        return tuple(
+            measure(mp, state, is_state_batched, readout_errors) for mp in circuit.measurements
+        )
 
     # finite-shot case
     rng = default_rng(rng)
@@ -140,6 +145,7 @@ def measure_final_state(circuit, state, is_state_batched, rng=None, prng_key=Non
             is_state_batched=is_state_batched,
             rng=rng,
             prng_key=prng_key,
+            readout_errors=readout_errors,
         )
         for mp in circuit.measurements
     )
@@ -151,8 +157,13 @@ def measure_final_state(circuit, state, is_state_batched, rng=None, prng_key=Non
     return results
 
 
-def simulate(
-    circuit: qml.tape.QuantumScript, rng=None, prng_key=None, debugger=None, interface=None
+def simulate(  # pylint: disable=too-many-arguments
+    circuit: qml.tape.QuantumScript,
+    rng=None,
+    prng_key=None,
+    debugger=None,
+    interface=None,
+    readout_errors=None,
 ) -> Result:
     """Simulate a single quantum script.
 
@@ -168,6 +179,8 @@ def simulate(
             generated. Only for simulation using JAX.
         debugger (_Debugger): The debugger to use
         interface (str): The machine learning interface to create the initial state with
+        readout_errors (List[Callable]): List of channels to apply to each wire being measured
+        to simulate readout errors.
 
     Returns:
         tuple(TensorLike): The results of the simulation
@@ -185,4 +198,11 @@ def simulate(
     state, is_state_batched = get_final_state(
         circuit, debugger=debugger, interface=interface, rng=rng, prng_key=prng_key
     )
-    return measure_final_state(circuit, state, is_state_batched, rng=rng, prng_key=prng_key)
+    return measure_final_state(
+        circuit,
+        state,
+        is_state_batched,
+        rng=rng,
+        prng_key=prng_key,
+        readout_errors=readout_errors,
+    )


Test patch:
diff --git a/tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py b/tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py
index 5f0c5c99b..a49265ce0 100644
--- a/tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py
+++ b/tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py
@@ -12,6 +12,8 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 """Tests for qutrit mixed device preprocessing."""
+import warnings
+
 import numpy as np
 import pytest
 
@@ -269,3 +271,47 @@ class TestPreprocessingIntegration:
         program, _ = DefaultQutritMixed().preprocess()
         with pytest.raises(qml.DeviceError, match="Operator NoMatNoDecompOp"):
             program(tapes)
+
+    @pytest.mark.parametrize(
+        "relaxations,misclassifications,req_warn",
+        [
+            [(0.1, 0.2, 0.3), None, True],
+            [None, (0.1, 0.2, 0.3), True],
+            [(0.1, 0.2, 0.3), (0.1, 0.2, 0.3), True],
+            [None, None, False],
+        ],
+    )
+    @pytest.mark.parametrize(
+        "measurements",
+        [
+            [qml.state()],
+            [qml.density_matrix(0)],
+            [qml.state(), qml.density_matrix([1, 2])],
+            [qml.state(), qml.expval(qml.GellMann(1))],
+        ],
+    )
+    def test_preprocess_warns_measurement_error_state(
+        self, relaxations, misclassifications, req_warn, measurements
+    ):
+        """Test that preprocess raises a warning if there is an analytic state measurement and
+        measurement error."""
+        tapes = [
+            qml.tape.QuantumScript(ops=[], measurements=measurements),
+            qml.tape.QuantumScript(
+                ops=[qml.THadamard(0), qml.TRZ(0.123, wires=1)], measurements=measurements
+            ),
+        ]
+        device = DefaultQutritMixed(
+            readout_relaxation_probs=relaxations, readout_misclassification_probs=misclassifications
+        )
+        program, _ = device.preprocess()
+
+        with warnings.catch_warnings(record=True) as warning:
+            program(tapes)
+            if req_warn:
+                assert len(warning) != 0
+                for warn in list(warning):
+                    print(type(warn.message))
+                    assert "is not affected by readout error" in warn.message.args[0]
+            else:
+                assert len(warning) == 0
diff --git a/tests/devices/test_default_qutrit_mixed.py b/tests/devices/test_default_qutrit_mixed.py
index 1203363b3..f26134219 100644
--- a/tests/devices/test_default_qutrit_mixed.py
+++ b/tests/devices/test_default_qutrit_mixed.py
@@ -13,13 +13,14 @@
 # limitations under the License.
 """Tests for default qutrit mixed."""
 
-from functools import reduce
+from functools import partial, reduce
 
 import numpy as np
 import pytest
 
 import pennylane as qml
 from pennylane import math
+from pennylane import numpy as qnp
 from pennylane.devices import DefaultQutritMixed, ExecutionConfig
 
 
@@ -757,23 +758,23 @@ class TestSumOfTermsDifferentiability:
     x = 0.52
 
     @staticmethod
-    def f(scale, coeffs, n_wires=5, offset=0.1):
+    def f(scale, coeffs, num_wires=5, offset=0.1):
         """Function to differentiate that implements a circuit with a SumOfTerms operator."""
-        ops = [qml.TRX(offset + scale * i, wires=i, subspace=(0, 2)) for i in range(n_wires)]
+        ops = [qml.TRX(offset + scale * i, wires=i, subspace=(0, 2)) for i in range(num_wires)]
         H = qml.Hamiltonian(
             coeffs,
             [
-                reduce(lambda x, y: x @ y, (qml.GellMann(i, 3) for i in range(n_wires))),
-                reduce(lambda x, y: x @ y, (qml.GellMann(i, 5) for i in range(n_wires))),
+                reduce(lambda x, y: x @ y, (qml.GellMann(i, 3) for i in range(num_wires))),
+                reduce(lambda x, y: x @ y, (qml.GellMann(i, 5) for i in range(num_wires))),
             ],
         )
         qs = qml.tape.QuantumScript(ops, [qml.expval(H)])
         return DefaultQutritMixed().execute(qs)
 
     @staticmethod
-    def expected(scale, coeffs, n_wires=5, offset=0.1, like="numpy"):
+    def expected(scale, coeffs, num_wires=5, offset=0.1, like="numpy"):
         """Gets the expected output of function TestSumOfTermsDifferentiability.f."""
-        phase = offset + scale * qml.math.asarray(range(n_wires), like=like)
+        phase = offset + scale * qml.math.asarray(range(num_wires), like=like)
         cosines = qml.math.cos(phase / 2) ** 2
         sines = -qml.math.sin(phase)
         return coeffs[0] * qml.math.prod(cosines) + coeffs[1] * qml.math.prod(sines)
@@ -1268,3 +1269,512 @@ class TestIntegration:
         grad_jit = jax.grad(qnode_jit)(x, y)
 
         assert qml.math.allclose(grad, grad_jit)
+
+
+@pytest.mark.parametrize("num_wires", [2, 3])
+class TestReadoutError:
+    """Tests for measurement readout error"""
+
+    setup_unitary = np.array(
+        [
+            [1 / np.sqrt(2), 1 / np.sqrt(3), 1 / np.sqrt(6)],
+            [np.sqrt(2 / 29), np.sqrt(3 / 29), -2 * np.sqrt(6 / 29)],
+            [-5 / np.sqrt(58), 7 / np.sqrt(87), 1 / np.sqrt(174)],
+        ]
+    ).T
+
+    def setup_state(self, num_wires):
+        """Sets up a basic state used for testing."""
+        qml.QutritUnitary(self.setup_unitary, wires=0)
+        qml.QutritUnitary(self.setup_unitary, wires=1)
+        if num_wires == 3:
+            qml.TAdd(wires=(0, 2))
+
+    @staticmethod
+    def get_expected_dm(num_wires):
+        """Gets the expected density matrix of the circuit for the first num_wires"""
+        state = np.array([2, 3, 6], dtype=complex) ** -(1 / 2)
+        if num_wires == 2:
+            state = np.kron(state, state)
+        if num_wires == 3:
+            state = sum(
+                [
+                    state[i] * reduce(np.kron, [v, state, v])
+                    for i, v in enumerate([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+                ]
+            )
+        return np.outer(state, state)
+
+    # Set up the sets of probabilities that are inputted as the readout errors.
+    relax_and_misclass = [
+        [(0, 0, 0), (0, 0, 0)],
+        [None, (1, 0, 0)],
+        [None, (0, 1, 0)],
+        [None, (0, 0, 1)],
+        [(1, 0, 0), None],
+        [(0, 1, 0), None],
+        [(1, 0, 1), None],
+        [(0, 1, 0), (0, 0, 1)],
+        [None, (0.1, 0.2, 0.4)],
+        [(0.2, 0.1, 0.3), None],
+        [(0.2, 0.1, 0.4), (0.1, 0.2, 0.5)],
+    ]
+
+    # Expected probabilities of measuring each state after the above readout errors are applied.
+    expected_probs = [
+        [1 / 2, 1 / 3, 1 / 6],
+        [1 / 3, 1 / 2, 1 / 6],
+        [1 / 6, 1 / 3, 1 / 2],
+        [1 / 2, 1 / 6, 1 / 3],
+        [5 / 6, 0, 1 / 6],
+        [2 / 3, 1 / 3, 0],
+        [5 / 6, 1 / 6, 0],
+        [2 / 3, 0, 1 / 3],
+        [5 / 12, 17 / 60, 0.3],
+        [7 / 12, 19 / 60, 0.1],
+        [11 / 24, 7 / 30, 37 / 120],
+    ]
+
+    @pytest.mark.parametrize(
+        "relax_and_misclass, expected", zip(relax_and_misclass, expected_probs)
+    )
+    def test_probs_with_readout_error(self, num_wires, relax_and_misclass, expected):
+        """Tests the measurement results for probs"""
+        dev = qml.device(
+            "default.qutrit.mixed",
+            wires=num_wires,
+            readout_relaxation_probs=relax_and_misclass[0],
+            readout_misclassification_probs=relax_and_misclass[1],
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            self.setup_state(num_wires)
+            return qml.probs(wires=0)
+
+        res = circuit()
+        assert np.allclose(res, expected)
+
+    # Expected expval list from circuit with diagonal observables after the readout errors
+    # defined by relax_and_misclass are applied.
+    expected_commuting_expvals = [
+        [1 / 6, 1 / (2 * np.sqrt(3)), 1 / 6],
+        [-1 / 6, 1 / (2 * np.sqrt(3)), -1 / 6],
+        [-1 / 6, -1 / (2 * np.sqrt(3)), -1 / 6],
+        [1 / 3, 0, 1 / 3],
+        [5 / 6, 1 / (2 * np.sqrt(3)), 5 / 6],
+        [1 / 3, 1 / np.sqrt(3), 1 / 3],
+        [2 / 3, 1 / np.sqrt(3), 2 / 3],
+        [4 / 6, 0, 4 / 6],
+        [2 / 15, 1 / (10 * np.sqrt(3)), 2 / 15],
+        [4 / 15, 7 / (10 * np.sqrt(3)), 4 / 15],
+        [9 / 40, 3 / (40 * np.sqrt(3)), 9 / 40],
+    ]
+
+    @pytest.mark.parametrize(
+        "relax_and_misclass, expected", zip(relax_and_misclass, expected_commuting_expvals)
+    )
+    def test_readout_expval_commuting(self, num_wires, relax_and_misclass, expected):
+        """Tests the measurement results for expval of diagonal GellMann observables (3 and 8)"""
+        dev = qml.device(
+            "default.qutrit.mixed",
+            wires=num_wires,
+            readout_relaxation_probs=relax_and_misclass[0],
+            readout_misclassification_probs=relax_and_misclass[1],
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            self.setup_state(num_wires)
+            return (
+                qml.expval(qml.GellMann(0, 3)),
+                qml.expval(qml.GellMann(0, 8)),
+                qml.expval(qml.GellMann(1, 3)),
+            )
+
+        res = circuit()
+        assert np.allclose(res, expected)
+
+    # Expected expval list from circuit with non-diagonal observables after the measurement errors
+    # defined by relax_and_misclass are applied. The readout error is applied to the circuit shown
+    # above so the pre measurement probabilities are the same.
+    expected_noncommuting_expvals = [
+        [-1 / 3, -7 / 6, 1 / 6],
+        [-1 / 6, -1, -1 / 6],
+        [1 / 3, -1 / 6, -1 / 6],
+        [-1 / 6, -5 / 6, 1 / 3],
+        [-2 / 3, -3 / 2, 5 / 6],
+        [-2 / 3, -5 / 3, 1 / 3],
+        [-5 / 6, -11 / 6, 2 / 3],
+        [-1 / 3, -1, 4 / 6],
+        [-7 / 60, -49 / 60, 2 / 15],
+        [-29 / 60, -83 / 60, 4 / 15],
+        [-3 / 20, -101 / 120, 9 / 40],
+    ]
+
+    @pytest.mark.parametrize(
+        "relax_and_misclass, expected", zip(relax_and_misclass, expected_noncommuting_expvals)
+    )
+    def test_readout_expval_non_commuting(self, num_wires, relax_and_misclass, expected):
+        """Tests the measurement results for expval of GellMann 1 observables"""
+        dev = qml.device(
+            "default.qutrit.mixed",
+            wires=num_wires,
+            readout_relaxation_probs=relax_and_misclass[0],
+            readout_misclassification_probs=relax_and_misclass[1],
+        )
+        # Create matrices for the observables with diagonalizing matrix :math:`THadamard^\dag`
+        inv_sqrt_3_i = 1j / np.sqrt(3)
+        non_commuting_obs_one = np.array(
+            [
+                [0, -1 + inv_sqrt_3_i, -1 - inv_sqrt_3_i],
+                [-1 - inv_sqrt_3_i, 0, -1 + inv_sqrt_3_i],
+                [-1 + inv_sqrt_3_i, -1 - inv_sqrt_3_i, 0],
+            ]
+        )
+        non_commuting_obs_one /= 2
+
+        non_commuting_obs_two = np.array(
+            [
+                [-2 / 3, -2 / 3 + inv_sqrt_3_i, -2 / 3 - inv_sqrt_3_i],
+                [-2 / 3 - inv_sqrt_3_i, -2 / 3, -2 / 3 + inv_sqrt_3_i],
+                [-2 / 3 + inv_sqrt_3_i, -2 / 3 - inv_sqrt_3_i, -2 / 3],
+            ]
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            self.setup_state(num_wires)
+
+            qml.THadamard(wires=0)
+            qml.THadamard(wires=1, subspace=(0, 1))
+
+            return (
+                qml.expval(qml.THermitian(non_commuting_obs_one, 0)),
+                qml.expval(qml.THermitian(non_commuting_obs_two, 0)),
+                qml.expval(qml.GellMann(1, 1)),
+            )
+
+        res = circuit()
+        assert np.allclose(res, expected)
+
+    state_relax_and_misclass = [
+        [(0, 0, 0), (0, 0, 0)],
+        [(0.1, 0.15, 0.25), (0.1, 0.15, 0.25)],
+        [(1, 0, 1), (1, 0, 0)],
+    ]
+
+    @pytest.mark.parametrize("relaxations, misclassifications", state_relax_and_misclass)
+    def test_readout_state(self, num_wires, relaxations, misclassifications):
+        """Tests the state output is not affected by readout error"""
+        dev = qml.device(
+            "default.qutrit.mixed",
+            wires=num_wires,
+            readout_relaxation_probs=relaxations,
+            readout_misclassification_probs=misclassifications,
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            self.setup_state(num_wires)
+            return qml.state()
+
+        res = circuit()
+        assert np.allclose(res, self.get_expected_dm(num_wires))
+
+    @pytest.mark.parametrize("relaxations, misclassifications", state_relax_and_misclass)
+    def test_readout_density_matrix(self, num_wires, relaxations, misclassifications):
+        """Tests the density matrix output is not affected by readout error"""
+        dev = qml.device(
+            "default.qutrit.mixed",
+            wires=num_wires,
+            readout_relaxation_probs=relaxations,
+            readout_misclassification_probs=misclassifications,
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            self.setup_state(num_wires)
+            return qml.density_matrix(wires=1)
+
+        res = circuit()
+        assert np.allclose(res, self.get_expected_dm(1))
+
+    @pytest.mark.parametrize(
+        "relaxations, misclassifications, expected",
+        [
+            ((0, 0, 0), (0, 0, 0), [(np.ones(2) * 2)] * 2),
+            (None, (0, 0, 1), [np.ones(2)] * 2),
+            ((0, 0, 1), None, [np.ones(2)] * 2),
+            (None, (0, 1, 0), [np.zeros(2)] * 2),
+            ((0, 1, 0), None, [np.zeros(2)] * 2),
+        ],
+    )
+    def test_readout_sample(self, num_wires, relaxations, misclassifications, expected):
+        """Tests the sample output with readout error"""
+        dev = qml.device(
+            "default.qutrit.mixed",
+            shots=2,
+            wires=num_wires,
+            readout_relaxation_probs=relaxations,
+            readout_misclassification_probs=misclassifications,
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            qml.QutritBasisState([2] * num_wires, wires=range(num_wires))
+            return qml.sample(wires=[0, 1])
+
+        res = circuit()
+        assert np.allclose(res, expected)
+
+    @pytest.mark.parametrize(
+        "relaxations, misclassifications, expected",
+        [
+            ((0, 0, 0), (0, 0, 0), {"22": 100}),
+            (None, (0, 0, 1), {"11": 100}),
+            ((0, 0, 1), None, {"11": 100}),
+            (None, (0, 1, 0), {"00": 100}),
+            ((0, 1, 0), None, {"00": 100}),
+        ],
+    )
+    def test_readout_counts(self, num_wires, relaxations, misclassifications, expected):
+        """Tests the counts output with readout error"""
+        dev = qml.device(
+            "default.qutrit.mixed",
+            shots=100,
+            wires=num_wires,
+            readout_relaxation_probs=relaxations,
+            readout_misclassification_probs=misclassifications,
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            qml.QutritBasisState([2] * num_wires, wires=range(num_wires))
+            return qml.counts(wires=[0, 1])
+
+        res = circuit()
+        assert res == expected
+
+    @pytest.mark.parametrize(
+        "relaxations, misclassifications, expected",
+        [
+            [None, (0.1, 0.2, 0.4), [5 / 12, 17 / 60, 0.3]],
+            [(0.2, 0.1, 0.3), None, [7 / 12, 19 / 60, 0.1]],
+            [(0.2, 0.1, 0.4), (0.1, 0.2, 0.5), [11 / 24, 7 / 30, 37 / 120]],
+        ],
+    )
+    def test_approximate_readout_counts(self, num_wires, relaxations, misclassifications, expected):
+        """Tests the counts output with readout error"""
+        num_shots = 10000
+        dev = qml.device(
+            "default.qutrit.mixed",
+            shots=num_shots,
+            wires=num_wires,
+            readout_relaxation_probs=relaxations,
+            readout_misclassification_probs=misclassifications,
+            seed=221349,
+        )
+
+        @qml.qnode(dev)
+        def circuit():
+            self.setup_state(num_wires)
+            return qml.counts(wires=[0])
+
+        res = circuit()
+        assert isinstance(res, dict)
+        assert len(res) == 3
+        cases = ["0", "1", "2"]
+        for case, expected_result in zip(cases, expected):
+            assert np.isclose(res[case] / num_shots, expected_result, atol=0.05)
+
+    @pytest.mark.parametrize(
+        "relaxations,misclassifications",
+        [
+            [(0.1, 0.2), None],
+            [None, (0.1, 0.2, 0.3, 0.1)],
+            [(0.1, 0.2, 0.3, 0.1), (0.1, 0.2, 0.3)],
+        ],
+    )
+    def test_measurement_error_validation(self, relaxations, misclassifications, num_wires):
+        """Ensure error is raised for wrong number of arguments inputted in readout errors."""
+        with pytest.raises(qml.DeviceError, match="results in error:"):
+            qml.device(
+                "default.qutrit.mixed",
+                wires=num_wires,
+                readout_relaxation_probs=relaxations,
+                readout_misclassification_probs=misclassifications,
+            )
+
+    def test_prob_type(self, num_wires):
+        """Tests that an error is raised for wrong data type in readout errors"""
+        with pytest.raises(qml.DeviceError, match="results in error:"):
+            qml.device(
+                "default.qutrit.mixed", wires=num_wires, readout_relaxation_probs=[0.1, 0.2, "0.3"]
+            )
+        with pytest.raises(qml.DeviceError, match="results in error:"):
+            qml.device(
+                "default.qutrit.mixed",
+                wires=num_wires,
+                readout_misclassification_probs=[0.1, 0.2, "0.3"],
+            )
+
+    diff_parameters = [
+        [
+            None,
+            [0.1, 0.2, 0.4],
+            [
+                [
+                    [1 / 3 - 1 / 2, 1 / 6 - 1 / 2, 0.0],
+                    [1 / 2 - 1 / 3, 0.0, 1 / 6 - 1 / 3],
+                    [0.0, 1 / 2 - 1 / 6, 1 / 3 - 1 / 6],
+                ]
+            ],
+        ],
+        [
+            [0.2, 0.1, 0.3],
+            None,
+            [
+                [
+                    [1 / 3, 1 / 6, 0.0],
+                    [-1 / 3, 0.0, 1 / 6],
+                    [0.0, -1 / 6, -1 / 6],
+                ]
+            ],
+        ],
+        [
+            [0.2, 0.1, 0.3],
+            [0.0, 0.0, 0.0],
+            [
+                [[1 / 3, 1 / 6, 0.0], [-1 / 3, 0.0, 1 / 6], [0.0, -1 / 6, -1 / 6]],
+                [
+                    [19 / 60 - 7 / 12, 0.1 - 7 / 12, 0.0],
+                    [7 / 12 - 19 / 60, 0.0, 0.1 - 19 / 60],
+                    [0.0, 7 / 12 - 0.1, 19 / 60 - 0.1],
+                ],
+            ],
+        ],
+    ]
+
+    def get_diff_function(self, interface, num_wires):
+        """Get the function to differentiate for following differentiability interface tests"""
+
+        def diff_func(relaxations, misclassifications):
+            dev = qml.device(
+                "default.qutrit.mixed",
+                wires=num_wires,
+                readout_relaxation_probs=relaxations,
+                readout_misclassification_probs=misclassifications,
+            )
+
+            @qml.qnode(dev, interface=interface)
+            def circuit():
+                self.setup_state(num_wires)
+                return qml.probs(0)
+
+            return circuit()
+
+        return diff_func
+
+    @pytest.mark.autograd
+    @pytest.mark.parametrize("relaxations, misclassifications, expected", diff_parameters)
+    def test_differentiation_autograd(self, num_wires, relaxations, misclassifications, expected):
+        """Tests the differentiation of readout errors using autograd"""
+
+        if misclassifications is None:
+            args_to_diff = (0,)
+            relaxations = qnp.array(relaxations)
+        elif relaxations is None:
+            args_to_diff = (1,)
+            misclassifications = qnp.array(misclassifications)
+        else:
+            args_to_diff = (0, 1)
+            relaxations = qnp.array(relaxations)
+            misclassifications = qnp.array(misclassifications)
+
+        diff_func = self.get_diff_function("autograd", num_wires)
+        jac = qml.jacobian(diff_func, args_to_diff)(relaxations, misclassifications)
+        assert np.allclose(jac, expected)
+
+    @pytest.mark.jax
+    @pytest.mark.parametrize("relaxations, misclassifications, expected", diff_parameters)
+    @pytest.mark.parametrize("use_jit", (True, False))
+    def test_differentiation_jax(  # pylint: disable=too-many-arguments
+        self, num_wires, relaxations, misclassifications, use_jit, expected
+    ):
+        """Tests the differentiation of readout errors using JAX"""
+        import jax
+
+        if misclassifications is None:
+            args_to_diff = (0,)
+            relaxations = jax.numpy.array(relaxations)
+        elif relaxations is None:
+            args_to_diff = (1,)
+            misclassifications = jax.numpy.array(misclassifications)
+        else:
+            args_to_diff = (0, 1)
+            relaxations = jax.numpy.array(relaxations)
+            misclassifications = jax.numpy.array(misclassifications)
+
+        diff_func = self.get_diff_function("jax", num_wires)
+        if use_jit:
+            diff_func = jax.jit(diff_func)
+        jac = jax.jacobian(diff_func, args_to_diff)(relaxations, misclassifications)
+        assert np.allclose(jac, expected)
+
+    @pytest.mark.torch
+    @pytest.mark.parametrize("relaxations, misclassifications, expected", diff_parameters)
+    def test_differentiation_torch(self, num_wires, relaxations, misclassifications, expected):
+        """Tests the differentiation of readout errors using PyTorch"""
+        import torch
+
+        if misclassifications is None:
+            relaxations = torch.tensor(relaxations, requires_grad=True, dtype=torch.float64)
+            diff_func = partial(self.get_diff_function("torch", num_wires), misclassifications=None)
+            diff_variables = relaxations
+        elif relaxations is None:
+            misclassifications = torch.tensor(
+                misclassifications, requires_grad=True, dtype=torch.float64
+            )
+
+            def diff_func(misclass):
+                return self.get_diff_function("torch", num_wires)(None, misclass)
+
+            diff_variables = misclassifications
+        else:
+            relaxations = torch.tensor(relaxations, requires_grad=True, dtype=torch.float64)
+            misclassifications = torch.tensor(
+                misclassifications, requires_grad=True, dtype=torch.float64
+            )
+            diff_variables = (relaxations, misclassifications)
+            diff_func = self.get_diff_function("torch", num_wires)
+
+        jac = torch.autograd.functional.jacobian(diff_func, diff_variables)
+        if isinstance(jac, tuple):
+            for j, expected_j in zip(jac, expected):
+                np.allclose(j.detach().numpy(), expected_j)
+        else:
+            assert np.allclose(jac.detach().numpy(), expected)
+
+    @pytest.mark.tf
+    @pytest.mark.parametrize("relaxations, misclassifications, expected", diff_parameters)
+    def test_differentiation_tensorflow(self, num_wires, relaxations, misclassifications, expected):
+        """Tests the differentiation of readout errors using TensorFlow"""
+        import tensorflow as tf
+
+        if misclassifications is None:
+            relaxations = tf.Variable(relaxations, dtype="float64")
+            diff_variables = [relaxations]
+        elif relaxations is None:
+            misclassifications = tf.Variable(misclassifications, dtype="float64")
+            diff_variables = [misclassifications]
+        else:
+            relaxations = tf.Variable(relaxations, dtype="float64")
+            misclassifications = tf.Variable(misclassifications, dtype="float64")
+            diff_variables = [relaxations, misclassifications]
+
+        diff_func = self.get_diff_function("tf", num_wires)
+        with tf.GradientTape() as grad_tape:
+            probs = diff_func(relaxations, misclassifications)
+        jac = grad_tape.jacobian(probs, diff_variables)
+        assert np.allclose(jac, expected)