Repository: TeamGraphix/graphix

Issue:
Density matrix backend
**Is your feature request related to a problem? Please describe.**
We cannot simulate with noise/error at the moment.

**Describe the solution you'd like**
Add density matrix backend so that we can simulate various quantum channels (ones described by Kraus ops).

**Describe alternatives you've considered**
na

**Additional context**
Consider how this could be done in TN backend.


Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index 6400c6f..f8203c1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,3 +13,4 @@ docs/source/benchmarks/**
 graphix/sim/graphix.code-workspace
 graphix/graphix.code-workspace
 *~
+.vscode/settings.json
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0dde497..a155c44 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Added flow and gflow verifiers ([#99](https://github.com/TeamGraphix/graphix/issues/99)).
 - Added `gflow.flow_from_pattern` method.
+- Added noisy MBQC simulation backend.
+    - `sim.density_matrix` module added for density matrix simulation backend, which is incorporated into the `simulator.PatternSimulator` interface.
+    - `noise_models` module, containing abstractclass `NoiseModel` and a simplified model (no noise) `NoiseLessNoiseModel`, to define operaion-specfic channels (e.g. 'N' and 'E' commands have separate noise models expressed by Kraus channels).
+    - `channels` module, defining `KrausChannel` class.
+    - `random_objects` and `linalg_validations` module for math support: random state, random unitary, random maps, matrix validations for channel definition.
+
 
 ### Fixed
 
@@ -42,6 +48,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   Both search for gflow or flow, and if found, plot them. If not found,
   in the case of from the graph, only the graph is drawn, and in the case
   of from the pattern, both the graph and all correction sets are drawn.
+- Removed `paddle` from benchmarks following github dependabot alert.
+- `PatternSimulator` takes optional argument noise_model during init, to specify noise model for `densitymatrix` simualtion.
 
 ## [0.2.10] - 2024-01-03
 
diff --git a/benchmarks/requirements.txt b/benchmarks/requirements.txt
index b01c053..0285b5d 100644
--- a/benchmarks/requirements.txt
+++ b/benchmarks/requirements.txt
@@ -1,4 +1,2 @@
-paddlepaddle==2.4.0rc0
-paddle-quantum==2.2.0
 networkx==3.2.1
 rustworkx==0.13.2
diff --git a/docs/source/channels.rst b/docs/source/channels.rst
new file mode 100644
index 0000000..35a3726
--- /dev/null
+++ b/docs/source/channels.rst
@@ -0,0 +1,36 @@
+
+Quantum channels and noise models
++++++++++++++++++++++++++++++++++
+
+Kraus channel
+-------------
+
+.. currentmodule:: graphix.channels
+
+.. autoclass:: KrausChannel
+    :members:
+
+.. autofunction:: dephasing_channel
+
+.. autofunction:: depolarising_channel
+
+.. autofunction:: pauli_channel
+
+.. autofunction:: two_qubit_depolarising_channel
+
+.. autofunction:: two_qubit_depolarising_tensor_channel
+
+
+Noise model classes
+-------------------
+
+
+.. currentmodule:: graphix.noise_models.noise_model
+
+.. autoclass:: NoiseModel
+    :members:
+
+.. currentmodule:: graphix.noise_models.noiseless_noise_model
+
+.. autoclass:: NoiselessNoiseModel
+    :members:
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 917bde4..75c2949 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -57,9 +57,9 @@ html_theme_options = {
 
 sphinx_gallery_conf = {
     # path to your example scripts
-    "examples_dirs": ["../../examples", "../../benchmarks"],
+    "examples_dirs": ["../../examples"], 
     # path to where to save gallery generated output
-    "gallery_dirs": ["gallery", "benchmarks"],
+    "gallery_dirs": ["gallery"], 
     "filename_pattern": "/",
     "thumbnail_size": (800, 550),
 }
diff --git a/docs/source/references.rst b/docs/source/references.rst
index 6c5f388..59dd127 100644
--- a/docs/source/references.rst
+++ b/docs/source/references.rst
@@ -13,4 +13,5 @@ Module reference
     clifford
     device_interface
     visualization
+    channels
 
diff --git a/docs/source/simulator.rst b/docs/source/simulator.rst
index c1c94f0..c83180f 100644
--- a/docs/source/simulator.rst
+++ b/docs/source/simulator.rst
@@ -42,3 +42,14 @@ Statevector
 
 .. autoclass:: Statevec
     :members:
+
+Density Matrix
+--------------
+
+.. currentmodule:: graphix.sim.density_matrix
+
+.. autoclass:: DensityMatrixBackend
+    :members:
+
+.. autoclass:: DensityMatrix
+    :members:
diff --git a/examples/noisy_mbqc.py b/examples/noisy_mbqc.py
new file mode 100644
index 0000000..d97a0fc
--- /dev/null
+++ b/examples/noisy_mbqc.py
@@ -0,0 +1,109 @@
+"""
+Simulating noisy MBQC
+=====================
+
+:class:`~graphix.density_matrix.DensityMatrix` class (through :class:`~graphix.simulator.PatternSimulator`) 
+allows the simulation of MBQC with customizable noise model.
+
+In this example, we simulate a simple MBQC pattern with various noise models to see their effects.
+First, let us import relevant modules and define a pattern
+"""
+# %%
+import numpy as np
+from graphix import Circuit
+
+circuit = Circuit(2)
+theta = np.random.rand(2)
+circuit.rz(0, theta[0])
+circuit.rz(1, theta[1])
+circuit.cnot(0,1)
+
+# %%
+# Now we transpile into measurement pattern using :meth:`~graphix.transpiler.Circuit.transpile` method.
+# This returns :class:`~graphix.pattern.Pattern` object containing measurement pattern:
+
+pattern = circuit.transpile()
+pattern.print_pattern(lim=30)
+# pattern.draw_graph()
+
+# %%
+# simulate with statevector backend:
+out_state = pattern.simulate_pattern(backend="statevector")
+print(out_state.flatten())
+
+# %%
+# Now let us define a noise model. We specify Kraus channels for each of the command executions.
+# Here, we apply dephasing noise to the qubit preparation.
+from graphix.noise_models.noiseless_noise_model import NoiselessNoiseModel
+from graphix.noise_models.noise_model import NoiseModel
+from graphix.channels import (
+    KrausChannel,
+    dephasing_channel,
+)
+
+class NoisyGraphState(NoiseModel):
+
+    def __init__(self, p_z=0.1):
+        self.p_z = p_z
+
+    def prepare_qubit(self):
+        """return the channel to apply after clean single-qubit preparation. Here just we prepare dephased qubits."""
+        return dephasing_channel(self.p_z)
+
+    def entangle(self):
+        """return noise model to qubits that happens after the CZ gate. just identity no noise for this noise model."""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(4)}])
+
+    def measure(self):
+        """apply noise to qubit to be measured."""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def confuse_result(self, cmd):
+        """imperfect measurement effect. here we do nothing (no error).
+        cmd = "M"
+        """
+        pass
+
+    def byproduct_x(self):
+        """apply noise to qubits after X gate correction. here no error (identity)."""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def byproduct_z(self):
+        """apply noise to qubits after Z gate correction. here no error (identity)."""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def clifford(self):
+        """apply noise to qubits that happens in the Clifford gate process. here no error (identity)."""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def tick_clock(self):
+        """notion of time in real devices - this is where we apply effect of T1 and T2.
+        we assume commands that lie between 'T' commands run simultaneously on the device.
+        
+        here we assume no idle error.
+        """
+        pass
+
+#%%
+# simulate with the noise model
+from graphix.simulator import PatternSimulator
+
+simulator = PatternSimulator(pattern, backend="densitymatrix", noise_model=NoisyGraphState(p_z=0.01))
+dm_result = simulator.run()
+print(dm_result.fidelity(out_state.psi.flatten()))
+
+# %%
+import matplotlib.pyplot as plt
+
+err_arr = np.logspace(-4, -1, 10)
+fidelity = np.zeros(10)
+for i in range(10):
+    simulator = PatternSimulator(pattern, backend="densitymatrix", noise_model=NoisyGraphState(p_z=err_arr[i]))
+    dm_result = simulator.run()
+    fidelity[i] = dm_result.fidelity(out_state.psi.flatten())
+
+plt.semilogx(err_arr, fidelity, "o:")
+plt.xlabel("dephasing error of qubit preparation")
+plt.ylabel("Final fidelity")
+plt.show()
+# %%
diff --git a/graphix/channels.py b/graphix/channels.py
new file mode 100644
index 0000000..0653d4c
--- /dev/null
+++ b/graphix/channels.py
@@ -0,0 +1,212 @@
+import numpy as np
+
+from graphix.linalg_validations import check_data_dims, check_data_normalization, check_data_values_type, check_rank
+from graphix.ops import Ops
+
+
+class KrausChannel:
+    """quantum channel class in the Kraus representation.
+    Defined by Kraus operators :math:`K_i` with scalar prefactors :code:`coef`) :math:`c_i`,
+    where the channel act on density matrix as :math:`\\rho'  = \sum_i K_i^\dagger \\rho K_i`.
+    The data should satisfy :math:`\sum K_i^\dagger K_i = I`
+
+    Attributes
+    ----------
+    nqubit : int
+        number of qubits acted on by the Kraus operators
+    size : int
+        number of Kraus operators (== Choi rank)
+    kraus_ops : array_like(dict())
+        the data in format
+        array_like(dict): [{coef: scalar, operator: array_like}, {coef: scalar, operator: array_like}, ...]
+
+    Returns
+    -------
+    Channel object
+        containing the corresponding Kraus operators
+
+    """
+
+    def __init__(self, kraus_data):
+        """
+        Parameters
+        ----------
+        kraus_data : array_like
+            array of Kraus operator data.
+            array_like(dict): [{coef: scalar, operator: array_like}, {parameter: scalar, operator: array_like}, ...]
+            only works for square Kraus operators
+
+        Raises
+        ------
+        ValueError
+            If empty array_like is provided.
+        """
+
+        # check there is data
+        if not kraus_data:
+            raise ValueError("Cannot instantiate the channel with empty data.")
+
+        if not isinstance(kraus_data, (list, np.ndarray, tuple)):
+            raise TypeError(f"The data must be a list, a numpy.ndarray or a tuple not a {type(kraus_data)}.")
+
+        # check that data is correctly formatted before assigning it to the object.
+        assert check_data_values_type(kraus_data)
+        assert check_data_dims(kraus_data)
+
+        # check that the channel is properly normalized i.e
+        # \sum_K_i^\dagger K_i = Identity
+        assert check_data_normalization(kraus_data)
+        assert check_rank(kraus_data)
+
+        self.nqubit = int(np.log2(kraus_data[0]["operator"].shape[0]))
+        self.kraus_ops = kraus_data
+
+        # np.asarray(data, dtype=np.complex128)
+        # number of Kraus operators in the Channel
+
+        self.size = len(kraus_data)
+
+    def __repr__(self):
+        return f"KrausChannel object with {self.size} Kraus operators of dimension {self.nqubit}."
+
+    def is_normalized(self):
+        return check_data_normalization(self.kraus_ops)
+
+
+def dephasing_channel(prob: float) -> KrausChannel:
+    """single-qubit dephasing channel, :math:`(1-p) \\rho + p Z  \\rho Z`
+
+    Parameters
+    ----------
+    prob : float
+        The probability associated to the channel
+
+    Returns
+    -------
+    :class:`graphix.channel.KrausChannel` object
+        containing the corresponding Kraus operators
+    """
+    return KrausChannel(
+        [{"coef": np.sqrt(1 - prob), "operator": np.eye(2)}, {"coef": np.sqrt(prob), "operator": Ops.z}]
+    )
+
+
+def depolarising_channel(prob: float) -> KrausChannel:
+    """single-qubit depolarizing channel
+
+    .. math::
+        (1-p) \\rho + \\frac{p}{3} (X \\rho X + Y \\rho Y + Z \\rho Z) = (1 - 4 \\frac{p}{3}) \\rho + 4 \\frac{p}{3} id
+
+    Parameters
+    ----------
+    prob : float
+        The probability associated to the channel
+
+    """
+    return KrausChannel(
+        [
+            {"coef": np.sqrt(1 - prob), "operator": np.eye(2)},
+            {"coef": np.sqrt(prob / 3.0), "operator": Ops.x},
+            {"coef": np.sqrt(prob / 3.0), "operator": Ops.y},
+            {"coef": np.sqrt(prob / 3.0), "operator": Ops.z},
+        ]
+    )
+
+
+def pauli_channel(px: float, py: float, pz: float) -> KrausChannel:
+    """single-qubit pauli channel,
+
+    .. math::
+        (1-p_X-p_Y-p_Z) \\rho + p_X X \\rho X + p_Y Y \\rho Y + p_Z Z \\rho Z)
+
+    """
+    if px + py + pz > 1:
+        raise ValueError("The sum of probabilities must not exceed 1.")
+    pI = 1 - px - py - pz
+    return KrausChannel(
+        [
+            {"coef": np.sqrt(1 - pI), "operator": np.eye(2)},
+            {"coef": np.sqrt(px / 3.0), "operator": Ops.x},
+            {"coef": np.sqrt(py / 3.0), "operator": Ops.y},
+            {"coef": np.sqrt(pz / 3.0), "operator": Ops.z},
+        ]
+    )
+
+
+def two_qubit_depolarising_channel(prob: float) -> KrausChannel:
+    """two-qubit depolarising channel.
+
+    .. math::
+        \mathcal{E} (\\rho) = (1-p) \\rho + \\frac{p}{15}  \sum_{P_i \in \{id, X, Y ,Z\}^{\otimes 2}/(id \otimes id)}P_i \\rho P_i
+
+    Parameters
+    ----------
+    prob : float
+        The probability associated to the channel
+
+    Returns
+    -------
+    :class:`graphix.channel.KrausChannel` object
+        containing the corresponding Kraus operators
+    """
+
+    return KrausChannel(
+        [
+            {"coef": np.sqrt(1 - prob), "operator": np.kron(np.eye(2), np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, Ops.x)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(np.eye(2), Ops.x)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(np.eye(2), Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(np.eye(2), Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, Ops.x)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, Ops.x)},
+        ]
+    )
+
+
+def two_qubit_depolarising_tensor_channel(prob: float) -> KrausChannel:
+    """two-qubit tensor channel of single-qubit depolarising channels with same probability.
+    Kraus operators:
+
+    .. math::
+        \Big\{ \sqrt{(1-p)} id, \sqrt{(p/3)} X, \sqrt{(p/3)} Y , \sqrt{(p/3)} Z \Big\} \otimes \Big\{ \sqrt{(1-p)} id, \sqrt{(p/3)} X, \sqrt{(p/3)} Y , \sqrt{(p/3)} Z \Big\}
+
+    Parameters
+    ----------
+    prob : float
+        The probability associated to the channel
+
+    Returns
+    -------
+    :class:`graphix.channel.KrausChannel` object
+        containing the corresponding Kraus operators
+    """
+
+    return KrausChannel(
+        [
+            {"coef": 1 - prob, "operator": np.kron(np.eye(2), np.eye(2))},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.x, Ops.x)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.y, Ops.y)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.z, Ops.z)},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(Ops.x, np.eye(2))},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(Ops.y, np.eye(2))},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(Ops.z, np.eye(2))},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(np.eye(2), Ops.x)},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(np.eye(2), Ops.y)},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(np.eye(2), Ops.z)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.x, Ops.y)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.x, Ops.z)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.y, Ops.x)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.y, Ops.z)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.z, Ops.x)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.z, Ops.y)},
+        ]
+    )
diff --git a/graphix/graphsim/basegraphstate.py b/graphix/graphsim/basegraphstate.py
index b9340ac..527a3ea 100644
--- a/graphix/graphsim/basegraphstate.py
+++ b/graphix/graphsim/basegraphstate.py
@@ -4,8 +4,7 @@ from abc import ABC, abstractmethod
 from typing import Union
 
 import networkx as nx
-from networkx import Graph
-from networkx.classes.reportviews import EdgeView, NodeView
+import networkx.classes.reportviews as nx_reportviews
 
 from graphix.clifford import CLIFFORD_HSZ_DECOMPOSITION, CLIFFORD_MUL
 from graphix.ops import Ops
@@ -23,9 +22,9 @@ except ModuleNotFoundError:
     rx = None
     PyGraph = None
 
-NodesObject = Union[NodeView, NodeList]
-EdgesObject = Union[EdgeView, EdgeList]
-GraphObject = Union[Graph, PyGraph]
+NodesObject = Union[nx_reportviews.NodeView, NodeList]
+EdgesObject = Union[nx_reportviews.EdgeView, EdgeList]
+GraphObject = Union[nx.Graph, PyGraph]
 
 
 class BaseGraphState(ABC):
diff --git a/graphix/linalg_validations.py b/graphix/linalg_validations.py
new file mode 100644
index 0000000..277b505
--- /dev/null
+++ b/graphix/linalg_validations.py
@@ -0,0 +1,116 @@
+import numpy as np
+from typing import Union
+
+
+def check_square(matrix: np.ndarray) -> bool:
+    """
+    check if matrix is a square matrix with a power of 2 dimension.
+    """
+
+    if len(matrix.shape) != 2:
+        raise ValueError(f"The object has {len(matrix.shape)} axes but must have 2 to be a matrix.")
+    if matrix.shape[0] != matrix.shape[1]:
+        raise ValueError(f"Matrix must be square but has different dimensions {matrix.shape}.")
+    size = matrix.shape[0]
+    if size & (size - 1) != 0:
+        raise ValueError(f"Matrix size must be a power of two but is {size}.")
+    return True
+
+
+def check_psd(matrix: np.ndarray, tol: float = 1e-15) -> bool:
+    """
+    check if a density matrix is positive semidefinite by diagonalizing.
+    After check_square and check_hermitian (osef) so that it already is square with power of 2 dimension.
+
+
+    Parameters
+    ----------
+    matrix : np.ndarray
+        matrix to check. Normally already square and 2**n x 2**n
+    tol : float
+        tolerance on the small negatives. Default 1e-15.
+    """
+
+    evals = np.linalg.eigvalsh(matrix)
+
+    if not all(evals >= -tol):
+        raise ValueError("The matrix is not positive semi-definite.")
+
+    return True
+
+
+def check_hermitian(matrix: np.ndarray) -> bool:
+    """
+    check if matrix is hermitian. After check_square.
+    """
+
+    if not np.allclose(matrix, matrix.transpose().conjugate()):
+        raise ValueError("The matrix is not Hermitian.")
+    return True
+
+
+def check_unit_trace(matrix: np.ndarray) -> bool:
+    """
+    check if matrix has trace 1. After check_square.
+    """
+
+    if not np.allclose(matrix.trace(), 1.0):
+        raise ValueError("The matrix does not have unit trace.")
+    return True
+
+
+def check_data_normalization(data: Union[list, tuple, np.ndarray]) -> bool:
+    # NOTE use np.conjugate() instead of object.conj() to certify behaviour when using non-numpy float/complex types
+    opsu = np.array([i["coef"] * np.conj(i["coef"]) * i["operator"].conj().T @ i["operator"] for i in data])
+
+    if not np.allclose(np.sum(opsu, axis=0), np.eye(2 ** int(np.log2(len(data[0]["operator"]))))):
+        raise ValueError(f"The specified channel is not normalized {np.sum(opsu, axis=0)}.")
+    return True
+
+
+def check_data_dims(data: Union[list, tuple, np.ndarray]) -> bool:
+
+    # convert to set to remove duplicates
+    dims = set([i["operator"].shape for i in data])
+
+    # check all the same dimensions and that they are square matrices
+    # TODO replace by using array.ndim
+    if len(dims) != 1:
+        raise ValueError(f"All provided Kraus operators do not have the same dimension {dims}!")
+
+    assert check_square(data[0]["operator"])
+
+    return True
+
+
+def check_data_values_type(data: Union[list, tuple, np.ndarray]) -> bool:
+
+    if not all(
+        isinstance(i, dict) for i in data
+    ):  # ni liste ni ensemble mais iterable (lazy) pas stocké, executé au besoin
+        raise TypeError("All values are not dictionaries.")
+
+    if not all(set(i.keys()) == {"coef", "operator"} for i in data):
+        raise KeyError("The keys of the indivudal Kraus operators must be coef and operator.")
+
+    if not all(isinstance(i["operator"], np.ndarray) for i in data):
+        raise TypeError("All operators don't have the same type and must be np.ndarray.")
+
+    for i in data:
+        if not i["operator"].dtype in (int, float, complex, np.float64, np.complex128):
+            raise TypeError(f"All operators dtype must be scalar and not {i['operator'].dtype}.")
+
+    if not all(isinstance(i["coef"], (int, float, complex, np.float64, np.complex128)) for i in data):
+        raise TypeError("All coefs dtype must be scalar.")
+
+    return True
+
+
+def check_rank(data: Union[list, tuple, np.ndarray]) -> bool:
+    # already checked that the data is list of square matrices
+    if len(data) > data[0]["operator"].shape[0] ** 2:
+        raise ValueError(
+            f"Incorrect number of Kraus operators in the expansion. This number must be an integer between 1 and the dimension squared."
+        )
+
+    return True
diff --git a/graphix/noise_models/__init__.py b/graphix/noise_models/__init__.py
new file mode 100644
index 0000000..587c959
--- /dev/null
+++ b/graphix/noise_models/__init__.py
@@ -0,0 +1,2 @@
+from .noiseless_noise_model import NoiselessNoiseModel
+from .noise_model import NoiseModel
diff --git a/graphix/noise_models/noise_model.py b/graphix/noise_models/noise_model.py
new file mode 100644
index 0000000..22f54ec
--- /dev/null
+++ b/graphix/noise_models/noise_model.py
@@ -0,0 +1,57 @@
+import abc
+
+
+class NoiseModel(abc.ABC):
+    """Abstract base class for all noise models."""
+
+    # define the noise model parameters in there
+    @abc.abstractmethod
+    def __init__(self):
+        pass
+
+    # shared by all objects of the child class.
+    def assign_simulator(self, simulator):
+        self.simulator = simulator
+
+    @abc.abstractmethod
+    def prepare_qubit(self):
+        """return qubit to be added with preparation errors."""
+        pass
+
+    @abc.abstractmethod
+    def entangle(self):
+        """apply noise to qubits that happens in the CZ gate process"""
+        pass
+
+    @abc.abstractmethod
+    def measure(self):
+        """apply noise to qubits that happens in the measurement process"""
+        pass
+
+    @abc.abstractmethod
+    def confuse_result(self):
+        """assign wrong measurement result"""
+        pass
+
+    @abc.abstractmethod
+    def byproduct_x(self):
+        """apply noise to qubits that happens in the X gate process"""
+        pass
+
+    @abc.abstractmethod
+    def byproduct_z(self):
+        """apply noise to qubits that happens in the Z gate process"""
+        pass
+
+    @abc.abstractmethod
+    def clifford(self):
+        """apply noise to qubits that happens in the Clifford gate process"""
+        # NOTE might be different depending on the gate.
+        pass
+
+    @abc.abstractmethod
+    def tick_clock(self):
+        """notion of time in real devices - this is where we apply effect of T1 and T2.
+        we assume commands that lie between 'T' commands run simultaneously on the device.
+        """
+        pass
diff --git a/graphix/noise_models/noiseless_noise_model.py b/graphix/noise_models/noiseless_noise_model.py
new file mode 100644
index 0000000..77e00b5
--- /dev/null
+++ b/graphix/noise_models/noiseless_noise_model.py
@@ -0,0 +1,51 @@
+import numpy as np
+from graphix.noise_models.noise_model import NoiseModel
+from graphix.channels import KrausChannel
+
+
+class NoiselessNoiseModel(NoiseModel):
+    """Noiseless noise model for testing.
+    Only return the identity channel.
+
+    :param NoiseModel: Parent abstract class class:`graphix.noise_model.NoiseModel`
+    :type NoiseModel: class
+    """
+
+    def __init__(self):
+        pass
+
+    def prepare_qubit(self):
+        """return the channel to apply after clean single-qubit preparation. Here just identity."""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def entangle(self):
+        """return noise model to qubits that happens after the CZ gats"""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(4)}])
+
+    def measure(self):
+        """apply noise to qubit to be measured."""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def confuse_result(self, cmd):
+        """assign wrong measurement result"""
+        p = 0.0
+        if np.random.rand() < p:
+            self.simulator.results[cmd[1]] = 1 - self.simulator.results[cmd[1]]
+
+    def byproduct_x(self):
+        """apply noise to qubits after X gate correction"""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def byproduct_z(self):
+        """apply noise to qubits after Z gate correction"""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def clifford(self):
+        """apply noise to qubits that happens in the Clifford gate process"""
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def tick_clock(self):
+        """notion of time in real devices - this is where we apply effect of T1 and T2.
+        we assume commands that lie between 'T' commands run simultaneously on the device.
+        """
+        pass
diff --git a/graphix/ops.py b/graphix/ops.py
index 305a0a5..b21d7a4 100644
--- a/graphix/ops.py
+++ b/graphix/ops.py
@@ -3,6 +3,8 @@ quantum states and operators
 """
 
 import numpy as np
+from itertools import product
+from functools import reduce
 
 
 class States:
@@ -38,6 +40,7 @@ class Ops:
             [0, 0, 0, 0, 0, 0, 1, 0],
         ]
     )
+    Pauli_ops = [np.eye(2), x, y, z]
 
     @staticmethod
     def Rx(theta):
@@ -52,6 +55,7 @@ class Ops:
         ----------
         operator : 2*2 np.array
         """
+
         return np.array([[np.cos(theta / 2), -1j * np.sin(theta / 2)], [-1j * np.sin(theta / 2), np.cos(theta / 2)]])
 
     @staticmethod
@@ -102,3 +106,23 @@ class Ops:
         operator : 4*4 np.array
         """
         return Ops.cnot @ np.kron(np.eye(2), Ops.Rz(theta)) @ Ops.cnot
+
+    @staticmethod
+    def build_tensor_Pauli_ops(n_qubits: int):
+        """Method to build all the 4^n tensor Pauli operators {I, X, Y, Z}^{\otimes n}
+
+        :param n_qubits: number of copies (qubits) to consider
+        :type n_qubits: int
+        :return: the array of the 4^n operators of shape (2^n, 2^n)
+        :rtype: np.ndarray
+        """
+
+        if isinstance(n_qubits, int):
+            if not 1 <= n_qubits:
+                raise ValueError(f"The number of qubits must be an integer <= 1 and not {n_qubits}.")
+        else:
+            raise TypeError(f"The number of qubits must be an integer and not {n_qubits}.")
+
+        tensor_Pauli_ops = [reduce(lambda x, y: np.kron(x, y), i) for i in product(Ops.Pauli_ops, repeat=n_qubits)]
+
+        return np.array(tensor_Pauli_ops)
diff --git a/graphix/pattern.py b/graphix/pattern.py
index 5a9f6c0..1f683ed 100644
--- a/graphix/pattern.py
+++ b/graphix/pattern.py
@@ -158,6 +158,9 @@ class Pattern:
         """iterate over commands"""
         return iter(self.__seq)
 
+    def __getitem__(self, index):
+        return self.__seq[index]
+
     @property
     def Nnode(self):
         """count of nodes that are either `input_nodes` or prepared with `N` commands"""
diff --git a/graphix/random_objects.py b/graphix/random_objects.py
new file mode 100644
index 0000000..0f917f0
--- /dev/null
+++ b/graphix/random_objects.py
@@ -0,0 +1,174 @@
+import numpy as np
+import numpy.typing as npt
+import scipy.linalg
+from scipy.stats import unitary_group
+
+from graphix.channels import KrausChannel
+from graphix.sim.density_matrix import DensityMatrix
+from graphix.ops import Ops
+
+
+def rand_herm(l: int):
+    """
+    generate random hermitian matrix of size l*l
+    """
+    tmp = np.random.rand(l, l) + 1j * np.random.rand(l, l)
+    return tmp + tmp.conj().T
+
+
+def rand_unit(l: int):
+    """
+    generate haar random unitary matrix of size l*l
+    """
+    if l == 1:
+        return np.array([np.exp(1j * np.random.rand(1) * 2 * np.pi)])
+    else:
+        return unitary_group.rvs(l)
+
+
+UNITS = np.array([1, 1j])
+
+
+def rand_dm(dim: int, rank: int = None, dm_dtype=True) -> DensityMatrix:
+    """Returns a "density matrix" as a DensityMatrix object ie a positive-semidefinite (hence Hermitian) matrix with unit trace
+    Note, not a proper DM since its dim can be something else than a power of 2.
+    The rank is random between 1 (pure) and dim if not specified
+    Thanks to Ulysse Chabaud.
+
+    :param dim: Linear dimension of the (square) matrix
+    :type dim: int
+    :param rank: If rank not specified then random between 1 and matrix dimension
+        If rank is one : then pure state else mixed state. Defaults to None
+    :type rank: int, optional
+    :param dm_dtype: If True returns a :class:`graphix.sim.density_matrix.DensityMatrix` or a numpy.ndarray if False. Defaults to True.
+    :type dm_dtype: bool, optional
+    :return: Random density matrix as a :class:`graphix.sim.density_matrix.DensityMatrix` object or a numpy.ndarray.
+    :rtype: :class:`graphix.sim.density_matrix.DensityMatrix` or numpy.ndarray.
+
+    """
+
+    # if not provided, use a random value.
+    if rank is None:
+        rank = np.random.randint(1, dim + 1)
+
+    evals = np.random.rand(rank)
+
+    padded_evals = np.zeros(dim)
+    padded_evals[: len(evals)] = evals
+
+    dm = np.diag(padded_evals / np.sum(padded_evals))
+
+    randU = rand_unit(dim)
+    dm = randU @ dm @ randU.transpose().conj()
+
+    if dm_dtype:
+        # will raise an error if incorrect dimension
+        return DensityMatrix(data=dm)
+    else:
+        return dm
+
+
+def rand_gauss_cpx_mat(dim: int, sig: float = 1 / np.sqrt(2)) -> npt.NDArray:
+    """
+    Returns a square array of standard normal complex random variates.
+    Code from QuTiP: https://qutip.org/docs/4.0.2/modules/qutip/random_objects.html
+
+    Parameters
+    ----------
+    dim : int
+        Linear dimension of the (square) matrix
+    sig : float
+        standard deviation of random variates.
+        ``sig = 'ginibre`` draws from the Ginibre ensemble ie  sig = 1 / sqrt(2 * dim).
+
+    """
+
+    if sig == "ginibre":
+        sig = 1.0 / np.sqrt(2 * dim)
+
+    return np.sum(np.random.normal(loc=0.0, scale=sig, size=((dim,) * 2 + (2,))) * UNITS, axis=-1)
+
+
+def rand_channel_kraus(dim: int, rank: int = None, sig: float = 1 / np.sqrt(2)) -> KrausChannel:
+    """
+    Returns a random :class:`graphix.sim.channels.KrausChannel`object of given dimension and rank following the method of
+    [KNPPZ21] Kukulski, Nechita, Pawela, Puchała, Życzkowsk https://arxiv.org/pdf/2011.02994.pdf
+
+    Parameters
+    ----------
+    dim : int
+        Linear dimension of the (square) matrix of each Kraus operator.
+        Only square operators so far.
+
+    rank : int (default to full rank dimension**2)
+        Choi rank ie the number of Kraus operators. Must be between one and dim**2.
+
+    sig : see rand_cpx
+
+    """
+
+    if rank is None:
+        rank = dim**2
+
+    if sig == "ginibre":
+        sig = 1.0 / np.sqrt(2 * dim)
+
+    if not isinstance(rank, int):
+        raise TypeError("The rank of a Kraus expansion must be an integer.")
+
+    if not 1 <= rank:
+        raise ValueError("The rank of a Kraus expansion must be greater or equal than 1.")
+
+    pre_kraus_list = [rand_gauss_cpx_mat(dim=dim, sig=sig) for _ in range(rank)]
+    Hmat = np.sum([m.transpose().conjugate() @ m for m in pre_kraus_list], axis=0)
+    kraus_list = np.array(pre_kraus_list) @ scipy.linalg.inv(scipy.linalg.sqrtm(Hmat))
+
+    return KrausChannel([{"coef": 1.0 + 0.0 * 1j, "operator": kraus_list[i]} for i in range(rank)])
+
+
+# or merge with previous with a "pauli" kwarg?
+### continue here
+def rand_Pauli_channel_kraus(dim: int, rank: int = None) -> KrausChannel:
+    if not isinstance(dim, int):
+        raise ValueError(f"The dimension must be an integer and not {dim}.")
+
+    if not dim & (dim - 1) == 0:
+        raise ValueError(f"The dimension must be a power of 2 and not {dim}.")
+
+    nqb = int(np.log2(dim))
+
+    # max number of ops (Choi rank) is d**2
+    # default is full rank.
+    if rank is None:
+        rank = dim**2
+    else:
+        if not isinstance(rank, int):
+            raise TypeError("The rank of a Kraus expansion must be an integer.")
+        if not 1 <= rank:
+            raise ValueError("The rank of a Kraus expansion must be an integer greater or equal than 1.")
+
+    # full probability has to have dim**2 operators.
+    prob_list = np.zeros(dim**2)
+    # generate rank random numbers and normalize
+    tmp_list = np.random.uniform(size=rank)
+    tmp_list /= tmp_list.sum()
+
+    # fill the list and shuffle
+    prob_list[:rank] = tmp_list
+    np.random.shuffle(prob_list)
+
+    tensor_Pauli_ops = Ops.build_tensor_Pauli_ops(nqb)
+    target_indices = np.nonzero(prob_list)
+
+    params = prob_list[target_indices]
+    ops = tensor_Pauli_ops[target_indices]
+
+    # TODO see how to use zip and dict to convert from tuple to dict
+    # https://www.tutorialspoint.com/How-I-can-convert-a-Python-Tuple-into-Dictionary
+
+    data = [{"coef": np.sqrt(params[i]), "operator": ops[i]} for i in range(0, rank)]
+
+    # NOTE retain a strong probability on the identity or not?
+    # think we don't really care
+
+    return KrausChannel(data)
diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py
new file mode 100644
index 0000000..737e617
--- /dev/null
+++ b/graphix/sim/base_backend.py
@@ -0,0 +1,53 @@
+import numpy as np
+import graphix.clifford
+import graphix.pauli
+
+
+class Backend:
+    def __init__(self, pr_calc: bool = True):
+        """
+        Parameters
+        ----------
+            pr_calc : bool
+                whether or not to compute the probability distribution before choosing the measurement result.
+                if False, measurements yield results 0/1 with 50% probabilities each.
+        """
+        # whether to compute the probability
+        self.pr_calc = pr_calc
+
+    def _perform_measure(self, cmd):
+        s_signal = np.sum([self.results[j] for j in cmd[4]])
+        t_signal = np.sum([self.results[j] for j in cmd[5]])
+        angle = cmd[3] * np.pi
+        if len(cmd) == 7:
+            vop = cmd[6]
+        else:
+            vop = 0
+        measure_update = graphix.pauli.MeasureUpdate.compute(
+            graphix.pauli.Plane[cmd[2]], s_signal % 2 == 1, t_signal % 2 == 1, graphix.clifford.TABLE[vop]
+        )
+        angle = angle * measure_update.coeff + measure_update.add_term
+        vec = measure_update.new_plane.polar(angle)
+        loc = self.node_index.index(cmd[1])
+
+        def op_mat_from_result(result: bool) -> np.ndarray:
+            op_mat = np.eye(2, dtype=np.complex128) / 2
+            sign = (-1) ** result
+            for i in range(3):
+                op_mat += sign * vec[i] * graphix.clifford.CLIFFORD[i + 1] / 2
+            return op_mat
+
+        if self.pr_calc:
+            op_mat = op_mat_from_result(False)
+            prob_0 = self.state.expectation_single(op_mat, loc)
+            result = np.random.rand() > prob_0
+            if result:
+                op_mat = op_mat_from_result(True)
+        else:
+            # choose the measurement result randomly
+            result = np.random.choice([0, 1])
+            op_mat = op_mat_from_result(result)
+        self.results[cmd[1]] = result
+        self.state.evolve_single(op_mat, loc)
+        self.node_index.remove(cmd[1])
+        return loc
diff --git a/graphix/sim/density_matrix.py b/graphix/sim/density_matrix.py
new file mode 100644
index 0000000..6317fbe
--- /dev/null
+++ b/graphix/sim/density_matrix.py
@@ -0,0 +1,406 @@
+"""Density matrix simulator.
+
+Simulate MBQC with density matrix representation.
+"""
+
+from copy import deepcopy
+
+import numpy as np
+
+from graphix.linalg_validations import check_square, check_hermitian, check_unit_trace
+from graphix.channels import KrausChannel
+from graphix.ops import Ops
+from graphix.clifford import CLIFFORD
+from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, meas_op
+import graphix.sim.base_backend
+
+
+class DensityMatrix:
+    """DensityMatrix object."""
+
+    def __init__(self, data=None, plus_state=True, nqubit=1):
+        """
+        Parameters
+        ----------
+            data : DensityMatrix, list, tuple, np.ndarray or None
+                Density matrix of shape (2**nqubits, 2**nqubits).
+            nqubit : int
+                Number of qubits. Default is 1. If both `data` and `nqubit` are specified, `nqubit` is ignored.
+        """
+        if data is None:
+            assert nqubit >= 0
+            self.Nqubit = nqubit
+            if plus_state:
+                self.rho = np.ones((2**nqubit, 2**nqubit)) / 2**nqubit
+            else:
+                self.rho = np.zeros((2**nqubit, 2**nqubit))
+                self.rho[0, 0] = 1.0
+        else:
+            if isinstance(data, DensityMatrix):
+                data = data.rho
+            elif isinstance(data, (list, tuple)):
+                data = np.asarray(data, dtype=complex)
+            elif isinstance(data, np.ndarray):
+                pass
+            else:
+                raise TypeError("data must be DensityMatrix, list, tuple, or np.ndarray.")
+
+            assert check_square(data)
+            self.Nqubit = len(data).bit_length() - 1
+
+            self.rho = data
+        assert check_hermitian(self.rho)
+        assert check_unit_trace(self.rho)
+
+    def __repr__(self):
+        return f"DensityMatrix, data={self.rho}, shape={self.dims()}"
+
+    def evolve_single(self, op, i):
+        """Single-qubit operation.
+
+        Parameters
+        ----------
+            op : np.ndarray
+                2*2 matrix.
+            i : int
+                Index of qubit to apply operator.
+        """
+        assert i >= 0 and i < self.Nqubit
+        if op.shape != (2, 2):
+            raise ValueError("op must be 2*2 matrix.")
+
+        rho_tensor = self.rho.reshape((2,) * self.Nqubit * 2)
+        rho_tensor = np.tensordot(np.tensordot(op, rho_tensor, axes=[1, i]), op.conj().T, axes=[i + self.Nqubit, 0])
+        rho_tensor = np.moveaxis(rho_tensor, (0, -1), (i, i + self.Nqubit))
+        self.rho = rho_tensor.reshape((2**self.Nqubit, 2**self.Nqubit))
+
+    def evolve(self, op, qargs):
+        """Multi-qubit operation
+
+        Args:
+            op (np.array): 2^n*2^n matrix
+            qargs (list of ints): target qubits' indexes
+        """
+
+        d = op.shape
+        # check it is a matrix.
+        if len(d) == 2:
+            # check it is square
+            if d[0] == d[1]:
+                pass
+            else:
+                raise ValueError(f"The provided operator has shape {op.shape} and is not a square matrix.")
+        else:
+            raise ValueError(f"The provided data has incorrect shape {op.shape}.")
+
+        nqb_op = np.log2(len(op))
+        if not np.isclose(nqb_op, int(nqb_op)):
+            raise ValueError("Incorrect operator dimension: not consistent with qubits.")
+        nqb_op = int(nqb_op)
+
+        if nqb_op != len(qargs):
+            raise ValueError("The dimension of the operator doesn't match the number of targets.")
+
+        if not all(0 <= i < self.Nqubit for i in qargs):
+            raise ValueError("Incorrect target indices.")
+        if len(set(qargs)) != nqb_op:
+            raise ValueError("A repeated target qubit index is not possible.")
+
+        op_tensor = op.reshape((2,) * 2 * nqb_op)
+
+        rho_tensor = self.rho.reshape((2,) * self.Nqubit * 2)
+
+        rho_tensor = np.tensordot(
+            np.tensordot(op_tensor, rho_tensor, axes=[tuple(nqb_op + i for i in range(len(qargs))), tuple(qargs)]),
+            op.conj().T.reshape((2,) * 2 * nqb_op),
+            axes=[tuple(i + self.Nqubit for i in qargs), tuple(i for i in range(len(qargs)))],
+        )
+        rho_tensor = np.moveaxis(
+            rho_tensor,
+            [i for i in range(len(qargs))] + [-i for i in range(1, len(qargs) + 1)],
+            [i for i in qargs] + [i + self.Nqubit for i in reversed(list(qargs))],
+        )
+        self.rho = rho_tensor.reshape((2**self.Nqubit, 2**self.Nqubit))
+
+    def expectation_single(self, op, i):
+        """Expectation value of single-qubit operator.
+
+        Args:
+            op (np.array): 2*2 Hermite operator
+            loc (int): Index of qubit on which to apply operator.
+        Returns:
+            complex: expectation value (real for hermitian ops!).
+        """
+
+        if not (0 <= i < self.Nqubit):
+            raise ValueError(f"Wrong target qubit {i}. Must between 0 and {self.Nqubit-1}.")
+
+        if op.shape != (2, 2):
+            raise ValueError("op must be 2x2 matrix.")
+
+        st1 = deepcopy(self)
+        st1.normalize()
+
+        rho_tensor = st1.rho.reshape((2,) * st1.Nqubit * 2)
+        rho_tensor = np.tensordot(op, rho_tensor, axes=[1, i])
+        rho_tensor = np.moveaxis(rho_tensor, 0, i)
+        st1.rho = rho_tensor.reshape((2**self.Nqubit, 2**self.Nqubit))
+
+        return np.trace(st1.rho)
+
+    def dims(self):
+        return self.rho.shape
+
+    def tensor(self, other):
+        r"""Tensor product state with other density matrix.
+        Results in self :math:`\otimes` other.
+
+        Parameters
+        ----------
+            other : :class: `DensityMatrix` object
+                DensityMatrix object to be tensored with self.
+        """
+        if not isinstance(other, DensityMatrix):
+            other = DensityMatrix(other)
+        self.rho = np.kron(self.rho, other.rho)
+        self.Nqubit += other.Nqubit
+
+    def cnot(self, edge):
+        """Apply CNOT gate to density matrix.
+
+        Parameters
+        ----------
+            edge : (int, int) or [int, int]
+                Edge to apply CNOT gate.
+        """
+
+        self.evolve(CNOT_TENSOR.reshape(4, 4), edge)
+
+    def swap(self, edge):
+        """swap qubits
+
+        Parameters
+        ----------
+            edge : (int, int) or [int, int]
+                (control, target) qubits indices.
+        """
+
+        self.evolve(SWAP_TENSOR.reshape(4, 4), edge)
+
+    def entangle(self, edge):
+        """connect graph nodes
+
+        Parameters
+        ----------
+            edge : (int, int) or [int, int]
+                (control, target) qubit indices.
+        """
+
+        self.evolve(CZ_TENSOR.reshape(4, 4), edge)
+
+    def normalize(self):
+        """normalize density matrix"""
+        self.rho /= np.trace(self.rho)
+
+    def ptrace(self, qargs):
+        """partial trace
+
+        Parameters
+        ----------
+            qargs : list of ints or int
+                Indices of qubit to trace out.
+        """
+        n = int(np.log2(self.rho.shape[0]))
+        if isinstance(qargs, int):
+            qargs = [qargs]
+        assert isinstance(qargs, (list, tuple))
+        qargs_num = len(qargs)
+        nqubit_after = n - qargs_num
+        assert n > 0
+        assert all([qarg >= 0 and qarg < n for qarg in qargs])
+
+        rho_res = self.rho.reshape((2,) * n * 2)
+        # ket, bra indices to trace out
+        trace_axes = list(qargs) + [n + qarg for qarg in qargs]
+        rho_res = np.tensordot(
+            np.eye(2**qargs_num).reshape((2,) * qargs_num * 2), rho_res, axes=(list(range(2 * qargs_num)), trace_axes)
+        )
+
+        self.rho = rho_res.reshape((2**nqubit_after, 2**nqubit_after))
+        self.Nqubit = nqubit_after
+
+    def fidelity(self, statevec):
+        """calculate the fidelity against reference statevector.
+
+        Parameters
+        ----------
+            statevec : numpy array
+                statevector (flattened numpy array) to compare with
+        """
+        return np.abs(statevec.transpose().conj() @ self.rho @ statevec)
+
+    def apply_channel(self, channel: KrausChannel, qargs):
+        """Applies a channel to a density matrix.
+
+        Parameters
+        ----------
+        :rho: density matrix.
+        channel: :class:`graphix.channel.KrausChannel` object
+            KrausChannel to be applied to the density matrix
+        qargs: target qubit indices
+
+        Returns
+        -------
+        nothing
+
+        Raises
+        ------
+        ValueError
+            If the final density matrix is not normalized after application of the channel.
+            This shouldn't happen since :class:`graphix.channel.KrausChannel` objects are normalized by construction.
+        ....
+        """
+
+        result_array = np.zeros((2**self.Nqubit, 2**self.Nqubit), dtype=np.complex128)
+        tmp_dm = deepcopy(self)
+
+        if not isinstance(channel, KrausChannel):
+            raise TypeError("Can't apply a channel that is not a Channel object.")
+
+        for k_op in channel.kraus_ops:
+            tmp_dm.evolve(k_op["operator"], qargs)
+            result_array += k_op["coef"] * np.conj(k_op["coef"]) * tmp_dm.rho
+            # reinitialize to input density matrix
+            tmp_dm = deepcopy(self)
+
+        # Performance?
+        self.rho = deepcopy(result_array)
+
+        if not np.allclose(self.rho.trace(), 1.0):
+            raise ValueError("The output density matrix is not normalized, check the channel definition.")
+
+
+class DensityMatrixBackend(graphix.sim.base_backend.Backend):
+    """MBQC simulator with density matrix method."""
+
+    def __init__(self, pattern, max_qubit_num=12, pr_calc=True):
+        """
+        Parameters
+        ----------
+            pattern : :class:`graphix.pattern.Pattern` object
+                Pattern to be simulated.
+            max_qubit_num : int
+                optional argument specifying the maximum number of qubits
+                to be stored in the statevector at a time.
+            pr_calc : bool
+                whether or not to compute the probability distribution before choosing the measurement result.
+                if False, measurements yield results 0/1 with 50% probabilities each.
+        """
+        # check that pattern has output nodes configured
+        assert len(pattern.output_nodes) > 0
+        self.pattern = pattern
+        self.results = deepcopy(pattern.results)
+        self.state = None
+        self.node_index = []
+        self.Nqubit = 0
+        self.max_qubit_num = max_qubit_num
+        if pattern.max_space() > max_qubit_num:
+            raise ValueError("Pattern.max_space is larger than max_qubit_num. Increase max_qubit_num and try again.")
+        super().__init__(pr_calc)
+
+    def add_nodes(self, nodes, qubit_to_add=None):
+        """add new qubit to the internal density matrix
+        and asign the corresponding node number to list self.node_index.
+
+        Parameters
+        ----------
+        nodes : list
+            list of node indices
+        qubit_to_add : DensityMatrix object
+            qubit to be added to the graph states
+        """
+        if not self.state:
+            self.state = DensityMatrix(nqubit=0)
+        n = len(nodes)
+        if qubit_to_add is None:
+            dm_to_add = DensityMatrix(nqubit=n)
+        else:
+            assert isinstance(qubit_to_add, DensityMatrix)
+            assert qubit_to_add.nqubit == 1
+            dm_to_add = qubit_to_add
+        self.state.tensor(dm_to_add)
+        self.node_index.extend(nodes)
+        self.Nqubit += n
+
+    def entangle_nodes(self, edge):
+        """Apply CZ gate to the two connected nodes.
+
+        Parameters
+        ----------
+            edge : tuple (int, int)
+                a pair of node indices
+        """
+        target = self.node_index.index(edge[0])
+        control = self.node_index.index(edge[1])
+        self.state.entangle((target, control))
+
+    def measure(self, cmd):
+        """Perform measurement on the specified node and trase out the qubit.
+
+        Parameters
+        ----------
+            cmd : list
+                measurement command : ['M', node, plane, angle, s_domain, t_domain]
+        """
+        loc = self._perform_measure(cmd)
+        self.state.normalize()
+        # perform ptrace right after the measurement (destructive measurement).
+        self.state.ptrace(loc)
+
+    def correct_byproduct(self, cmd):
+        """Byproduct correction
+        correct for the X or Z byproduct operators,
+        by applying the X or Z gate.
+        """
+        if np.mod(np.sum([self.results[j] for j in cmd[2]]), 2) == 1:
+            loc = self.node_index.index(cmd[1])
+            if cmd[0] == "X":
+                op = Ops.x
+            elif cmd[0] == "Z":
+                op = Ops.z
+            self.state.evolve_single(op, loc)
+
+    def sort_qubits(self):
+        """sort the qubit order in internal density matrix"""
+        for i, ind in enumerate(self.pattern.output_nodes):
+            if not self.node_index[i] == ind:
+                move_from = self.node_index.index(ind)
+                self.state.swap((i, move_from))
+                self.node_index[i], self.node_index[move_from] = (
+                    self.node_index[move_from],
+                    self.node_index[i],
+                )
+
+    def apply_clifford(self, cmd):
+        """backend version of apply_channel
+        Parameters
+        ----------
+            qargs : list of ints. Target qubits
+        """
+        loc = self.node_index.index(cmd[1])
+        self.state.evolve_single(CLIFFORD[cmd[2]], loc)
+
+    def apply_channel(self, channel: KrausChannel, qargs):
+        """backend version of apply_channel
+        Parameters
+        ----------
+            qargs : list of ints. Target qubits
+        """
+
+        indices = [self.node_index.index(i) for i in qargs]
+        self.state.apply_channel(channel, indices)
+
+    def finalize(self):
+        """To be run at the end of pattern simulation."""
+        self.sort_qubits()
+        self.state.normalize()
diff --git a/graphix/sim/statevec.py b/graphix/sim/statevec.py
index d50842d..1aa929c 100644
--- a/graphix/sim/statevec.py
+++ b/graphix/sim/statevec.py
@@ -4,14 +4,13 @@ import numpy as np
 
 from graphix.clifford import CLIFFORD, CLIFFORD_CONJ, CLIFFORD_MUL
 from graphix.ops import Ops
-import graphix.clifford
-import graphix.pauli
+import graphix.sim.base_backend
 
 
-class StatevectorBackend:
+class StatevectorBackend(graphix.sim.base_backend.Backend):
     """MBQC simulator with statevector method."""
 
-    def __init__(self, pattern, max_qubit_num=20):
+    def __init__(self, pattern, max_qubit_num=20, pr_calc=True):
         """
         Parameters
         -----------
@@ -22,6 +21,9 @@ class StatevectorBackend:
         max_qubit_num : int
             optional argument specifying the maximum number of qubits
             to be stored in the statevector at a time.
+        pr_calc : bool
+            whether or not to compute the probability distribution before choosing the measurement result.
+            if False, measurements yield results 0/1 with 50% probabilities each.
         """
         # check that pattern has output nodes configured
         assert len(pattern.output_nodes) > 0
@@ -35,6 +37,7 @@ class StatevectorBackend:
         self.max_qubit_num = max_qubit_num
         if pattern.max_space() > max_qubit_num:
             raise ValueError("Pattern.max_space is larger than max_qubit_num. Increase max_qubit_num and try again")
+        super().__init__(pr_calc)
 
     def qubit_dim(self):
         """Returns the qubit number in the internal statevector
@@ -82,32 +85,8 @@ class StatevectorBackend:
         cmd : list
             measurement command : ['M', node, plane angle, s_domain, t_domain]
         """
-        # choose the measurement result randomly
-        result = np.random.choice([0, 1])
-        self.results[cmd[1]] = result
-
-        # extract signals for adaptive angle
-        s_signal = np.sum([self.results[j] for j in cmd[4]])
-        t_signal = np.sum([self.results[j] for j in cmd[5]])
-        angle = cmd[3] * np.pi
-        if len(cmd) == 7:
-            vop = cmd[6]
-        else:
-            vop = 0
-        measure_update = graphix.pauli.MeasureUpdate.compute(
-            graphix.pauli.Plane[cmd[2]], s_signal % 2 == 1, t_signal % 2 == 1, graphix.clifford.TABLE[vop]
-        )
-        angle = angle * measure_update.coeff + measure_update.add_term
-        vec = measure_update.new_plane.polar(angle)
-        op_mat = np.eye(2, dtype=np.complex128) / 2
-        for i in range(3):
-            op_mat += (-1) ** (result) * vec[i] * CLIFFORD[i + 1] / 2
-
-        loc = self.node_index.index(cmd[1])
-        self.state.evolve_single(op_mat, loc)
-
+        loc = self._perform_measure(cmd)
         self.state.remove_qubit(loc)
-        self.node_index.remove(cmd[1])
         self.Nqubit -= 1
 
     def correct_byproduct(self, cmd):
@@ -246,6 +225,7 @@ class Statevec:
             target qubits' indices
         """
         op_dim = int(np.log2(len(op)))
+        # TODO shape = (2,)* 2 * op_dim
         shape = [2 for _ in range(2 * op_dim)]
         op_tensor = op.reshape(shape)
         self.psi = np.tensordot(
@@ -404,7 +384,7 @@ class Statevec:
         """
         st1 = deepcopy(self)
         st1.normalize()
-        st2 = st1.deepcopy(st1)
+        st2 = deepcopy(st1)
         st1.evolve_single(op, loc)
         return np.dot(st2.psi.flatten().conjugate(), st1.psi.flatten())
 
diff --git a/graphix/simulator.py b/graphix/simulator.py
index 6cf7f52..b1f03ea 100644
--- a/graphix/simulator.py
+++ b/graphix/simulator.py
@@ -4,8 +4,13 @@ Simulates MBQC by executing the pattern.
 
 """
 
-from graphix.sim.tensornet import TensorNetworkBackend
+import numpy as np
+
+from graphix.sim.density_matrix import DensityMatrixBackend
 from graphix.sim.statevec import StatevectorBackend
+from graphix.sim.tensornet import TensorNetworkBackend
+from graphix.noise_models import NoiseModel
+import warnings
 
 
 class PatternSimulator:
@@ -14,32 +19,55 @@ class PatternSimulator:
     Executes the measurement pattern.
     """
 
-    def __init__(self, pattern, backend="statevector", **kwargs):
+    def __init__(self, pattern, backend="statevector", noise_model=None, **kwargs):
         """
         Parameters
         -----------
         pattern: :class:`graphix.pattern.Pattern` object
             MBQC pattern to be simulated.
-        backend: str, 'statevector' or 'tensornetwork'
+        backend: str, 'statevector', 'densitymatrix or 'tensornetwork'
             simulation backend (optional), default is 'statevector'.
+        noise_model:
         kwargs: keyword args for specified backend.
 
         .. seealso:: :class:`graphix.sim.statevec.StatevectorBackend`\
-            :class:`graphix.sim.tensornet.TensorNetworkBackend`
+            :class:`graphix.sim.tensornet.TensorNetworkBackend`\
+            :class:`graphix.sim.density_matrix.DensityMatrixBackend`\
         """
         # check that pattern has output nodes configured
         assert len(pattern.output_nodes) > 0
-        if backend == "statevector":
+
+        if backend == "statevector" and noise_model is None:
+            self.noise_model = None
             self.backend = StatevectorBackend(pattern, **kwargs)
-        elif backend in {"tensornetwork", "mps"}:
+        elif backend == "densitymatrix":
+            if noise_model is None:
+                self.noise_model = None
+                # no noise: no need to compute probabilities
+                self.backend = DensityMatrixBackend(pattern, **kwargs)
+                warnings.warn(
+                    "Simulating using densitymatrix backend with no noise. To add noise to the simulation, give an object of `graphix.noise_models.Noisemodel` to `noise_model` keyword argument."
+                )
+            if noise_model is not None:
+                self.set_noise_model(noise_model)
+                # if noise: have to compute the probabilities
+                self.backend = DensityMatrixBackend(pattern, pr_calc=True, **kwargs)
+        elif backend in {"tensornetwork", "mps"} and noise_model is None:
+            self.noise_model = None
             self.backend = TensorNetworkBackend(pattern, **kwargs)
+        # TODO or just do the noiseless sim with a warning?
+        elif backend in {"statevector", "tensornetwork", "mps"} and noise_model is not None:
+            raise ValueError(f"The backend {backend} doesn't support noise but noisemodel was provided.")
         else:
-            raise ValueError("unknown backend")
+            raise ValueError("Unknown backend.")
         self.pattern = pattern
         self.results = self.backend.results
         self.state = self.backend.state
         self.node_index = []
 
+    def set_noise_model(self, model):
+        self.noise_model = model
+
     def run(self):
         """Perform the simulation.
 
@@ -49,22 +77,57 @@ class PatternSimulator:
             the output quantum state,
             in the representation depending on the backend used.
         """
+
         self.backend.add_nodes(self.pattern.input_nodes)
-        for cmd in self.pattern:
-            if cmd[0] == "N":
-                self.backend.add_nodes([cmd[1]])
-            elif cmd[0] == "E":
-                self.backend.entangle_nodes(cmd[1])
-            elif cmd[0] == "M":
-                self.backend.measure(cmd)
-            elif cmd[0] == "X":
-                self.backend.correct_byproduct(cmd)
-            elif cmd[0] == "Z":
-                self.backend.correct_byproduct(cmd)
-            elif cmd[0] == "C":
-                self.backend.apply_clifford(cmd)
-            else:
-                raise ValueError("invalid commands")
-        self.backend.finalize()
+        if self.noise_model is None:
+            for cmd in self.pattern:
+                if cmd[0] == "N":
+                    self.backend.add_nodes([cmd[1]])
+                elif cmd[0] == "E":
+                    self.backend.entangle_nodes(cmd[1])
+                elif cmd[0] == "M":
+                    self.backend.measure(cmd)
+                elif cmd[0] == "X":
+                    self.backend.correct_byproduct(cmd)
+                elif cmd[0] == "Z":
+                    self.backend.correct_byproduct(cmd)
+                elif cmd[0] == "C":
+                    self.backend.apply_clifford(cmd)
+                else:
+                    raise ValueError("invalid commands")
+            self.backend.finalize()
+        else:
+            self.noise_model.assign_simulator(self)
+            for node in self.pattern.input_nodes:
+                self.backend.apply_channel(self.noise_model.prepare_qubit(), [node])
+            for cmd in self.pattern:
+                if cmd[0] == "N":  # prepare clean qubit and apply channel
+                    self.backend.add_nodes([cmd[1]])
+                    self.backend.apply_channel(self.noise_model.prepare_qubit(), [cmd[1]])
+                elif cmd[0] == "E":  # for "E" cmd[1] is already a tuyple
+                    self.backend.entangle_nodes(cmd[1])  # for some reaon entangle doesn't get the whole command
+                    self.backend.apply_channel(self.noise_model.entangle(), cmd[1])
+                elif cmd[0] == "M":  # apply channel before measuring, then measur and confuse_result
+                    self.backend.apply_channel(self.noise_model.measure(), [cmd[1]])
+                    self.backend.measure(cmd)
+                    self.noise_model.confuse_result(cmd)
+                elif cmd[0] == "X":
+                    self.backend.correct_byproduct(cmd)
+                    if np.mod(np.sum([self.results[j] for j in cmd[2]]), 2) == 1:
+                        self.backend.apply_channel(self.noise_model.byproduct_x(), [cmd[1]])
+                elif cmd[0] == "Z":
+                    self.backend.correct_byproduct(cmd)
+                    if np.mod(np.sum([self.results[j] for j in cmd[2]]), 2) == 1:
+                        self.backend.apply_channel(self.noise_model.byproduct_z(), [cmd[1]])
+                elif cmd[0] == "C":
+                    self.backend.apply_clifford(cmd)
+                    self.backend.apply_channel(self.noise_model.clifford(), [cmd[1]])
+                elif cmd[0] == "T":
+                    # T command is a flag for one clock cycle in simulated experiment,
+                    # to be added via hardware-agnostic pattern modifier
+                    self.noise_model.tick_clock()
+                else:
+                    raise ValueError("Invalid commands.")
+            self.backend.finalize()
 
         return self.backend.state
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index ec63049..abda951 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -3,6 +3,7 @@
 accepts desired gate operations and transpile into MBQC measurement patterns.
 
 """
+
 from __future__ import annotations
 
 from copy import deepcopy
diff --git a/setup.py b/setup.py
index d05870b..a8b7dc9 100644
--- a/setup.py
+++ b/setup.py
@@ -37,7 +37,7 @@ info = {
     ],
     "python_requires": ">=3.8,<3.12",
     "install_requires": requirements,
-    "extras_require": {"extra": ["graphix-ibmq", "graphix-perceval"]},
+    "extras_require": {"extra": ["graphix-ibmq", "graphix-perceval"], "test": ["graphix-ibmq", "scipy"]},
 }
 
 setup(**(info))


Test patch:
diff --git a/tests/test_clifford.py b/tests/test_clifford.py
index 4426df5..4ce0a36 100644
--- a/tests/test_clifford.py
+++ b/tests/test_clifford.py
@@ -1,6 +1,8 @@
 import unittest
+
 import numpy as np
-from graphix.clifford import CLIFFORD, CLIFFORD_CONJ, CLIFFORD_MEASURE, CLIFFORD_MUL, CLIFFORD_HSZ_DECOMPOSITION
+
+from graphix.clifford import CLIFFORD, CLIFFORD_CONJ, CLIFFORD_HSZ_DECOMPOSITION, CLIFFORD_MEASURE, CLIFFORD_MUL
 
 
 class TestClifford(unittest.TestCase):
diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py
new file mode 100644
index 0000000..689f94d
--- /dev/null
+++ b/tests/test_density_matrix.py
@@ -0,0 +1,881 @@
+import random
+import unittest
+from copy import deepcopy
+
+import numpy as np
+
+import graphix.random_objects as randobj
+from graphix import Circuit
+from graphix.channels import KrausChannel, dephasing_channel, depolarising_channel
+from graphix.ops import Ops
+from graphix.sim.density_matrix import DensityMatrix, DensityMatrixBackend
+from graphix.sim.statevec import CNOT_TENSOR, CZ_TENSOR, SWAP_TENSOR, Statevec, StatevectorBackend
+
+
+class TestDensityMatrix(unittest.TestCase):
+    """Test for DensityMatrix class."""
+
+    def test_init_without_data_fail(self):
+        with self.assertRaises(AssertionError):
+            DensityMatrix(nqubit=-2)
+        with self.assertRaises(TypeError):
+            DensityMatrix(nqubit="hello")
+        with self.assertRaises(TypeError):
+            DensityMatrix(nqubit=[])
+
+    def test_init_with_invalid_data_fail(self):
+        with self.assertRaises(TypeError):
+            DensityMatrix("hello")
+        with self.assertRaises(TypeError):
+            DensityMatrix(1)
+        # deprecated data shape (these test might be unnecessary)
+        with self.assertRaises(ValueError):
+            DensityMatrix([1, 2, [3]])
+
+        # check with hermitian dm but not unit trace
+        with self.assertRaises(ValueError):
+            DensityMatrix(data=randobj.rand_herm(2 ** np.random.randint(2, 5)))
+
+        # check with non hermitian dm but unit trace
+        with self.assertRaises(ValueError):
+            l = 2 ** np.random.randint(2, 5)
+            tmp = np.random.rand(l, l) + 1j * np.random.rand(l, l)
+            DensityMatrix(data=tmp / np.trace(tmp))
+        # check with non hermitian dm and not unit trace
+        with self.assertRaises(ValueError):
+            l = 2 ** np.random.randint(2, 5)  # np.random.randint(2, 20)
+            DensityMatrix(data=np.random.rand(l, l) + 1j * np.random.rand(l, l))
+
+        # check not square matrix
+        with self.assertRaises(ValueError):
+            # l = 2 ** np.random.randint(2, 5) # np.random.randint(2, 20)
+            DensityMatrix(data=np.random.rand(3, 2))
+
+        # check higher dimensional matrix
+        with self.assertRaises(ValueError):
+            DensityMatrix(data=np.random.rand(2, 2, 3))
+
+        # check square and hermitian but with incorrect dimension (non-qubit type)
+        with self.assertRaises(ValueError):
+            # not really a dm since not PSD but ok.
+            data = randobj.rand_herm(5)
+            data /= np.trace(data)
+            DensityMatrix(data=data)
+
+    def test_init_without_data_success(self):
+        for n in range(3):
+            dm = DensityMatrix(nqubit=n)
+            expected_density_matrix = np.outer(np.ones((2,) * n), np.ones((2,) * n)) / 2**n
+            assert dm.Nqubit == n
+            assert dm.rho.shape == (2**n, 2**n)
+            assert np.allclose(dm.rho, expected_density_matrix)
+
+            dm = DensityMatrix(plus_state=False, nqubit=n)
+            expected_density_matrix = np.zeros((2**n, 2**n))
+            expected_density_matrix[0, 0] = 1
+            assert dm.Nqubit == n
+            assert dm.rho.shape == (2**n, 2**n)
+            assert np.allclose(dm.rho, expected_density_matrix)
+
+    def test_init_with_data_success(self):
+        # don't use rand_dm here since want to check
+        for n in range(3):
+            data = randobj.rand_herm(2**n)
+
+            data /= np.trace(data)
+            dm = DensityMatrix(data=data)
+            assert dm.Nqubit == n
+            assert dm.rho.shape == (2**n, 2**n)
+            assert np.allclose(dm.rho, data)
+
+    def test_evolve_single_fail(self):
+        dm = DensityMatrix(nqubit=2)
+        # generate random 4 x 4 unitary matrix
+        op = randobj.rand_unit(4)
+
+        with self.assertRaises(AssertionError):
+            dm.evolve_single(op, 2)
+        with self.assertRaises(ValueError):
+            dm.evolve_single(op, 1)
+
+    def test_evolve_single_success(self):
+        # generate random 2 x 2 unitary matrix
+        op = randobj.rand_unit(2)
+
+        n = 10
+        for i in range(n):
+            sv = Statevec(nqubit=n)
+            sv.evolve_single(op, i)
+            expected_density_matrix = np.outer(sv.psi, sv.psi.conj())
+            dm = DensityMatrix(nqubit=n)
+            dm.evolve_single(op, i)
+            assert np.allclose(dm.rho, expected_density_matrix)
+
+    def test_expectation_single_fail(self):
+        nqb = 3
+        dm = DensityMatrix(nqubit=nqb)
+
+        # wrong dimensions
+        # generate random 4 x 4 unitary matrix
+        op = randobj.rand_unit(4)
+
+        with self.assertRaises(ValueError):
+            dm.expectation_single(op, 2)
+        with self.assertRaises(ValueError):
+            dm.expectation_single(op, 1)
+
+        # wrong qubit indices
+        op = randobj.rand_unit(2)
+        with self.assertRaises(ValueError):
+            dm.expectation_single(op, -3)
+        with self.assertRaises(ValueError):
+            dm.expectation_single(op, nqb + 3)
+
+    def test_expectation_single_success(self):
+        """compare to pure state case
+        hence only pure states
+        but by linearity ok"""
+
+        nqb = np.random.randint(1, 4)
+        # NOTE a statevector object so can't use its methods
+        target_qubit = np.random.randint(0, nqb)
+
+        psi = np.random.rand(2**nqb) + 1j * np.random.rand(2**nqb)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+
+        op = randobj.rand_unit(2)
+
+        dm.expectation_single(op, target_qubit)
+
+        # by hand: copy paste from SV backend
+
+        psi1 = np.tensordot(op, psi.reshape((2,) * nqb), (1, target_qubit))
+        psi1 = np.moveaxis(psi1, 0, target_qubit)
+        psi1 = psi1.reshape(2**nqb)
+
+        # watch out ordering. Expval unitary is cpx so psi1 on the right to match DM.
+        np.testing.assert_allclose(np.dot(psi.conjugate(), psi1), dm.expectation_single(op, target_qubit))
+
+    def test_tensor_fail(self):
+        dm = DensityMatrix(nqubit=1)
+        with self.assertRaises(TypeError):
+            dm.tensor("hello")
+        with self.assertRaises(TypeError):
+            dm.tensor(1)
+
+    def test_tensor_without_data_success(self):
+        for n in range(3):
+            dm_a = DensityMatrix(nqubit=n)
+            dm_b = DensityMatrix(nqubit=n + 1)
+            dm_a.tensor(dm_b)
+            assert dm_a.Nqubit == 2 * n + 1
+            assert dm_a.rho.shape == (2 ** (2 * n + 1), 2 ** (2 * n + 1))
+
+    def test_tensor_with_data_success(self):
+        for n in range(3):
+            data_a = randobj.rand_dm(2**n, dm_dtype=False)
+            dm_a = DensityMatrix(data=data_a)
+
+            data_b = randobj.rand_dm(2 ** (n + 1), dm_dtype=False)
+            dm_b = DensityMatrix(data=data_b)
+            dm_a.tensor(dm_b)
+            assert dm_a.Nqubit == 2 * n + 1
+            assert dm_a.rho.shape == (2 ** (2 * n + 1), 2 ** (2 * n + 1))
+            assert np.allclose(dm_a.rho, np.kron(data_a, data_b))
+
+    def test_cnot_fail(self):
+        dm = DensityMatrix(nqubit=2)
+        with self.assertRaises(ValueError):
+            dm.cnot((1, 1))
+        with self.assertRaises(ValueError):
+            dm.cnot((-1, 1))
+        with self.assertRaises(ValueError):
+            dm.cnot((1, -1))
+        with self.assertRaises(ValueError):
+            dm.cnot((1, 2))
+        with self.assertRaises(ValueError):
+            dm.cnot((2, 1))
+
+    def test_cnot_success(self):
+        dm = DensityMatrix(nqubit=2)
+        original_matrix = dm.rho.copy()
+        dm.cnot((0, 1))
+        expected_matrix = np.array([[1, 1, 1, 1] * 4]).reshape((4, 4)) / 4
+        assert np.allclose(dm.rho, expected_matrix)
+        dm.cnot((0, 1))
+        assert np.allclose(dm.rho, original_matrix)
+
+        # test on 2 qubits only
+        psi = np.random.rand(4) + 1j * np.random.rand(4)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+        edge = (0, 1)
+        dm.cnot(edge)
+        rho = dm.rho.copy()
+        psi = psi.reshape((2, 2))
+        psi = np.tensordot(CNOT_TENSOR, psi, ((2, 3), edge))
+        psi = np.moveaxis(psi, (0, 1), edge)
+        expected_matrix = np.outer(psi, psi.conj())
+        assert np.allclose(rho, expected_matrix)
+
+        # test on arbitrary number of qubits and random pair
+        n = np.random.randint(2, 4)
+        psi = np.random.rand(2**n) + 1j * np.random.rand(2**n)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+
+        # for test only. Sample distinct pairs (== without replacement).
+        # https://docs.python.org/2/library/random.html#random.sample
+
+        edge = tuple(random.sample(range(n), 2))
+        dm.cnot(edge)
+        rho = dm.rho.copy()
+        psi = psi.reshape((2,) * n)
+        psi = np.tensordot(CNOT_TENSOR, psi, ((2, 3), edge))
+        psi = np.moveaxis(psi, (0, 1), edge)
+        expected_matrix = np.outer(psi, psi.conj())
+        np.testing.assert_allclose(rho, expected_matrix)
+
+    def test_swap_fail(self):
+        dm = DensityMatrix(nqubit=2)
+        with self.assertRaises(ValueError):
+            dm.swap((1, 1))
+        with self.assertRaises(ValueError):
+            dm.swap((-1, 1))
+        with self.assertRaises(ValueError):
+            dm.swap((1, -1))
+        with self.assertRaises(ValueError):
+            dm.swap((1, 2))
+        with self.assertRaises(ValueError):
+            dm.swap((2, 1))
+
+    def test_swap_success(self):
+        dm = DensityMatrix(nqubit=2)
+        original_matrix = dm.rho.copy()
+        dm.swap((0, 1))
+        expected_matrix = np.array([[1, 1, 1, 1] * 4]).reshape((4, 4)) / 4
+        assert np.allclose(dm.rho, expected_matrix)
+        dm.swap((0, 1))
+        assert np.allclose(dm.rho, original_matrix)
+
+        psi = np.random.rand(4) + 1j * np.random.rand(4)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+        edge = (0, 1)
+        dm.swap(edge)
+        rho = dm.rho
+        psi = psi.reshape((2, 2))
+        psi = np.tensordot(SWAP_TENSOR, psi, ((2, 3), edge))
+        psi = np.moveaxis(psi, (0, 1), edge)
+        expected_matrix = np.outer(psi, psi.conj())
+        assert np.allclose(rho, expected_matrix)
+
+    def test_entangle_fail(self):
+        dm = DensityMatrix(nqubit=3)
+        with self.assertRaises(ValueError):
+            dm.entangle((1, 1))
+        with self.assertRaises(ValueError):
+            dm.entangle(((1, 3)))
+        with self.assertRaises(ValueError):
+            dm.entangle((0, 1, 2))
+
+    def test_entangle_success(self):
+        dm = DensityMatrix(nqubit=2)
+        original_matrix = dm.rho.copy()
+        dm.entangle((0, 1))
+        expected_matrix = np.array([[1, 1, 1, -1], [1, 1, 1, -1], [1, 1, 1, -1], [-1, -1, -1, 1]]) / 4
+        assert np.allclose(dm.rho, expected_matrix)
+        dm.entangle((0, 1))
+        assert np.allclose(dm.rho, original_matrix)
+
+        dm = DensityMatrix(nqubit=3)
+        dm.entangle((0, 1))
+        dm.entangle((2, 1))
+        dm1 = dm.rho.copy()
+
+        dm = DensityMatrix(nqubit=3)
+        dm.entangle((2, 1))
+        dm.entangle((0, 1))
+        dm2 = dm.rho
+        assert np.allclose(dm1, dm2)
+
+        psi = np.random.rand(4) + 1j * np.random.rand(4)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+        edge = (0, 1)
+        dm.entangle(edge)
+        rho = dm.rho
+        psi = psi.reshape((2, 2))
+        psi = np.tensordot(CZ_TENSOR, psi, ((2, 3), edge))
+        psi = np.moveaxis(psi, (0, 1), edge)
+        expected_matrix = np.outer(psi, psi.conj())
+        np.testing.assert_allclose(rho, expected_matrix)
+
+    def test_evolve_success(self):
+        # single-qubit gate
+        # check against evolve_single
+
+        N_qubits = np.random.randint(2, 4)
+        N_qubits_op = 1
+
+        # random statevector
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+
+        # density matrix calculation
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+        dm_single = deepcopy(dm)
+
+        op = randobj.rand_unit(2**N_qubits_op)
+        i = np.random.randint(0, N_qubits)
+
+        # need a list format for a single target
+        dm.evolve(op, [i])
+        dm_single.evolve_single(op, i)
+
+        np.testing.assert_allclose(dm.rho, dm_single.rho)
+
+        # 2-qubit gate
+
+        N_qubits = np.random.randint(2, 4)
+        N_qubits_op = 2
+
+        # random unitary
+        op = randobj.rand_unit(2**N_qubits_op)
+
+        # random pair of indices
+        edge = tuple(random.sample(range(N_qubits), 2))
+
+        # random statevector to compare to
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+
+        # density matrix calculation
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+        dm.evolve(op, edge)
+        rho = dm.rho
+
+        psi = psi.reshape((2,) * N_qubits)
+        psi = np.tensordot(op.reshape((2,) * 2 * N_qubits_op), psi, ((2, 3), edge))
+        psi = np.moveaxis(psi, (0, 1), edge)
+        expected_matrix = np.outer(psi, psi.conj())
+        np.testing.assert_allclose(rho, expected_matrix)
+
+        # 3-qubit gate
+        N_qubits = np.random.randint(3, 5)
+        N_qubits_op = 3
+
+        # random unitary
+        op = randobj.rand_unit(2**N_qubits_op)
+
+        # 3 random indices
+        targets = tuple(random.sample(range(N_qubits), 3))
+
+        # random statevector to compare to
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+
+        # density matrix calculation
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+        dm.evolve(op, targets)
+        rho = dm.rho
+
+        psi = psi.reshape((2,) * N_qubits)
+        psi = np.tensordot(op.reshape((2,) * 2 * N_qubits_op), psi, ((3, 4, 5), targets))
+        psi = np.moveaxis(psi, (0, 1, 2), targets)
+        expected_matrix = np.outer(psi, psi.conj())
+        np.testing.assert_allclose(rho, expected_matrix)
+
+    def test_evolve_fail(self):
+        # test on 3-qubit gate just in case.
+        N_qubits = np.random.randint(3, 5)
+        N_qubits_op = 3
+
+        # random unitary
+        op = randobj.rand_unit(2**N_qubits_op)
+        # 3 random indices
+        dm = DensityMatrix(nqubit=N_qubits)
+
+        # dimension mismatch
+        with self.assertRaises(ValueError):
+            dm.evolve(op, (1, 1))
+        with self.assertRaises(ValueError):
+            dm.evolve(op, (0, 1, 2, 3))
+        # incorrect range
+        with self.assertRaises(ValueError):
+            dm.evolve(op, (-1, 0, 1))
+        # repeated index
+        with self.assertRaises(ValueError):
+            dm.evolve(op, (0, 1, 1))
+
+        # check not square matrix
+        with self.assertRaises(ValueError):
+            dm.evolve(np.random.rand(2, 3), (0, 1))
+
+        # check higher dimensional matrix
+        with self.assertRaises(ValueError):
+            dm.evolve(np.random.rand(2, 2, 3), (0, 1))
+
+        # check square but with incorrect dimension (non-qubit type)
+        with self.assertRaises(ValueError):
+            dm.evolve(np.random.rand(5, 5), (0, 1))
+
+    # TODO the test for normalization is done at initialization with data. Now check that all operations conserve the norm.
+    def test_normalize(self):
+        #  tmp = np.random.rand(4, 4) + 1j * np.random.rand(4, 4)
+
+        data = randobj.rand_herm(2 ** np.random.randint(2, 4))
+
+        dm = DensityMatrix(data / data.trace())
+        dm.normalize()
+        assert np.allclose(np.trace(dm.rho), 1)
+
+    def test_ptrace_fail(self):
+        dm = DensityMatrix(nqubit=0)
+        with self.assertRaises(AssertionError):
+            dm.ptrace((0,))
+        dm = DensityMatrix(nqubit=2)
+        with self.assertRaises(AssertionError):
+            dm.ptrace((2,))
+
+    def test_ptrace(self):
+        psi = np.kron(np.array([1, 0]), np.array([1, 1]) / np.sqrt(2))
+        data = np.outer(psi, psi)
+        dm = DensityMatrix(data=data)
+        dm.ptrace((0,))
+        expected_matrix = np.array([[0.5, 0.5], [0.5, 0.5]])
+        assert np.allclose(dm.rho, expected_matrix)
+
+        dm = DensityMatrix(data=data)
+        dm.ptrace((1,))
+        expected_matrix = np.array([[1, 0], [0, 0]])
+        assert np.allclose(dm.rho, expected_matrix)
+
+        dm = DensityMatrix(data=data)
+        dm.ptrace((0, 1))
+        expected_matrix = np.array([1])
+        assert np.allclose(dm.rho, expected_matrix)
+
+        dm = DensityMatrix(nqubit=4)
+        dm.ptrace((0, 1, 2))
+        expected_matrix = np.array([[1, 1], [1, 1]]) / 2
+        assert np.allclose(dm.rho, expected_matrix)
+
+        dm = DensityMatrix(nqubit=4)
+        dm.ptrace((0, 1, 2, 3))
+        expected_matrix = np.array([1])
+        assert np.allclose(dm.rho, expected_matrix)
+
+        psi = np.kron(np.kron(np.array([1, np.sqrt(2)]) / np.sqrt(3), np.array([1, 0])), np.array([0, 1]))
+        data = np.outer(psi, psi)
+        dm = DensityMatrix(data=data)
+        dm.ptrace((2,))
+        expected_matrix = np.array(
+            [[1 / 3, 0, np.sqrt(2) / 3, 0], [0, 0, 0, 0], [np.sqrt(2) / 3, 0, 2 / 3, 0], [0, 0, 0, 0]]
+        )
+        assert np.allclose(dm.rho, expected_matrix)
+
+    def test_apply_dephasing_channel(self):
+        # check on single qubit first
+        # # create random density matrix
+        # data = randobj.rand_herm(2 ** np.random.randint(2, 4))
+        data = randobj.rand_herm(2)
+        data /= np.trace(data)
+        dm = DensityMatrix(data=data)
+
+        # copy of initial dm
+        rho_test = deepcopy(dm.rho)
+
+        # create dephasing channel
+        prob = np.random.rand()
+        dephase_channel = dephasing_channel(prob)
+
+        # useless since checked in apply_channel method.
+        assert isinstance(dephase_channel, KrausChannel)
+        # useless since checked in the constructor.
+        assert dephase_channel.is_normalized()
+
+        # apply channel. list with single element needed.
+        # if Channel.nqubit == 1 use list with single element.
+        dm.apply_channel(dephase_channel, [0])
+        id = np.array([[1.0, 0.0], [0.0, 1.0]])
+
+        # compare
+        expected_dm = (
+            np.sqrt(1 - prob) ** 2 * id @ rho_test @ id.conj().T
+            + np.sqrt(prob) ** 2 * Ops.z @ rho_test @ Ops.z.conj().T
+        )
+
+        np.testing.assert_allclose(expected_dm.trace(), 1.0)
+        np.testing.assert_allclose(dm.rho, expected_dm)
+
+        N_qubits = np.random.randint(2, 5)
+
+        i = np.random.randint(0, N_qubits)
+
+        # create random density matrix from statevector
+
+        # random statevector to compare to
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+
+        # build DensityMatrix
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+
+        # create dephasing channel
+        prob = np.random.rand()
+        dephase_channel = dephasing_channel(prob)
+
+        # useless since checked in apply_channel method.
+        assert isinstance(dephase_channel, KrausChannel)
+        # useless since checked in the constructor.
+        assert dephase_channel.is_normalized()
+
+        # apply channel. list with single element needed.
+        # if Channel.nqubit == 1 use list with single element.
+        dm.apply_channel(dephase_channel, [i])
+
+        # compute on the statevector
+        # psi.reshape((2,) * N_qubits)
+        # tmp = np.zeros(psi.shape)
+
+        id = np.array([[1.0, 0.0], [0.0, 1.0]])
+
+        # by hand: operator list and gate application
+        psi_evolved = np.tensordot(id, psi.reshape((2,) * N_qubits), (1, i))
+        psi_evolved = np.moveaxis(psi_evolved, 0, i)
+
+        psi_evolvedb = np.tensordot(Ops.z, psi.reshape((2,) * N_qubits), (1, i))
+        psi_evolvedb = np.moveaxis(psi_evolvedb, 0, i)
+
+        # compute final density matrix
+        psi_evolved = np.reshape(psi_evolved, (2**N_qubits))
+        psi_evolvedb = np.reshape(psi_evolvedb, (2**N_qubits))
+        expected_dm = np.sqrt(1 - prob) ** 2 * np.outer(psi_evolved, psi_evolved.conj()) + np.sqrt(
+            prob
+        ) ** 2 * np.outer(psi_evolvedb, psi_evolvedb.conj())
+
+        # compare
+        np.testing.assert_allclose(expected_dm.trace(), 1.0)
+        np.testing.assert_allclose(dm.rho, expected_dm)
+
+    def test_apply_depolarising_channel(self):
+        # check on single qubit first
+        # # create random density matrix
+        # data = randobj.rand_herm(2 ** np.random.randint(2, 4))
+        data = randobj.rand_herm(2)
+        data /= np.trace(data)
+        dm = DensityMatrix(data=data)
+
+        # copy of initial dm
+        rho_test = deepcopy(dm.rho)
+
+        # create dephasing channel
+        prob = np.random.rand()
+        depol_channel = depolarising_channel(prob)
+
+        # useless since checked in apply_channel method.
+        assert isinstance(depol_channel, KrausChannel)
+        # useless since checked in the constructor.
+        assert depol_channel.is_normalized()
+
+        # apply channel. list with single element needed.
+        # if Channel.nqubit == 1 use list with single element.
+        dm.apply_channel(depol_channel, [0])
+        id = np.array([[1.0, 0.0], [0.0, 1.0]])
+
+        # compare
+        expected_dm = (
+            np.sqrt(1 - prob) ** 2 * id @ rho_test @ id.conj().T
+            + np.sqrt(prob / 3.0) ** 2 * Ops.x @ rho_test @ Ops.x.conj().T
+            + np.sqrt(prob / 3.0) ** 2 * Ops.y @ rho_test @ Ops.y.conj().T
+            + np.sqrt(prob / 3.0) ** 2 * Ops.z @ rho_test @ Ops.z.conj().T
+        )
+
+        np.testing.assert_allclose(expected_dm.trace(), 1.0)
+        np.testing.assert_allclose(dm.rho, expected_dm)
+
+        # chek against statevector backend by hand for now.
+        # create random density matrix
+
+        N_qubits = np.random.randint(2, 5)
+
+        # target qubit index
+        i = np.random.randint(0, N_qubits)
+
+        # create random density matrix from statevector
+
+        # random statevector to compare to
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+
+        # build DensityMatrix
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+
+        # create dephasing channel
+        prob = np.random.rand()
+        depol_channel = depolarising_channel(prob)
+
+        # useless since checked in apply_channel method.
+        assert isinstance(depol_channel, KrausChannel)
+        # useless since checked in the constructor.
+        assert depol_channel.is_normalized()
+
+        # apply channel. list with single element needed.
+        # if Channel.nqubit == 1 use list with single element.
+        dm.apply_channel(depol_channel, [i])
+
+        # compute on the statevector
+        # psi.reshape((2,) * N_qubits)
+        # tmp = np.zeros(psi.shape)
+
+        id = np.array([[1.0, 0.0], [0.0, 1.0]])
+
+        # by hand: operator list and gate application
+        psi_evolved = np.tensordot(id, psi.reshape((2,) * N_qubits), (1, i))
+        psi_evolved = np.moveaxis(psi_evolved, 0, i)
+
+        psi_evolvedb = np.tensordot(Ops.x, psi.reshape((2,) * N_qubits), (1, i))
+        psi_evolvedb = np.moveaxis(psi_evolvedb, 0, i)
+
+        psi_evolvedc = np.tensordot(Ops.y, psi.reshape((2,) * N_qubits), (1, i))
+        psi_evolvedc = np.moveaxis(psi_evolvedc, 0, i)
+
+        psi_evolvedd = np.tensordot(Ops.z, psi.reshape((2,) * N_qubits), (1, i))
+        psi_evolvedd = np.moveaxis(psi_evolvedd, 0, i)
+
+        # compute final density matrix
+        psi_evolved = np.reshape(psi_evolved, (2**N_qubits))
+        psi_evolvedb = np.reshape(psi_evolvedb, (2**N_qubits))
+        psi_evolvedc = np.reshape(psi_evolvedc, (2**N_qubits))
+        psi_evolvedd = np.reshape(psi_evolvedd, (2**N_qubits))
+
+        expected_dm = (
+            np.sqrt(1 - prob) ** 2 * np.outer(psi_evolved, psi_evolved.conj())
+            + np.sqrt(prob / 3.0) ** 2 * np.outer(psi_evolvedb, psi_evolvedb.conj())
+            + np.sqrt(prob / 3.0) ** 2 * np.outer(psi_evolvedc, psi_evolvedc.conj())
+            + np.sqrt(prob / 3.0) ** 2 * np.outer(psi_evolvedd, psi_evolvedd.conj())
+        )
+
+        # compare
+        np.testing.assert_allclose(expected_dm.trace(), 1.0)
+        np.testing.assert_allclose(dm.rho, expected_dm)
+
+    def test_apply_random_channel_one_qubit(self):
+        """
+        test random 1-qubit channel.
+        Especially checks for complex parameters.
+        """
+
+        # check against statevector backend by hand for now.
+        # create random density matrix
+
+        N_qubits = np.random.randint(2, 5)
+        # id = np.array([[1.0, 0.0], [0.0, 1.0]])
+
+        # target qubit index
+        i = np.random.randint(0, N_qubits)
+
+        # create random density matrix from statevector
+
+        # random statevector to compare to
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+
+        # build DensityMatrix
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+
+        # create random channel
+        # random_channel utility already checked for type and such
+        # here dim = 2 (single qubit) and rank is between 1 and 4
+        nqb = 1
+        dim = 2**nqb
+        rk = np.random.randint(1, dim**2 + 1)
+        channel = randobj.rand_channel_kraus(dim=dim, rank=rk)
+
+        # apply channel. list with single element needed.
+        # if Channel.nqubit == 1 use list with single element.
+        dm.apply_channel(channel, [i])
+
+        # compute on the statevector
+        # psi.reshape((2,) * N_qubits)
+        # tmp = np.zeros(psi.shape)
+
+        # initialize. NOT a DM object, just a matrix.
+        expected_dm = np.zeros((2**N_qubits, 2**N_qubits), dtype=np.complex128)
+
+        for elem in channel.kraus_ops:  # kraus_ops is a list of dicts
+            psi_evolved = np.tensordot(elem["operator"], psi.reshape((2,) * N_qubits), (1, i))
+            psi_evolved = np.moveaxis(psi_evolved, 0, i)
+            expected_dm += elem["coef"] * np.conj(elem["coef"]) * np.outer(psi_evolved, np.conj(psi_evolved))
+
+        # compare
+        np.testing.assert_allclose(expected_dm.trace(), 1.0)
+        np.testing.assert_allclose(dm.rho, expected_dm)
+
+    def test_apply_random_channel_two_qubits(self):
+        """
+        test random 2-qubit channel on a rank 1 dm (pure state). Generalizes by linearity.
+        Especially checks for complex parameters.
+        """
+
+        N_qubits = np.random.randint(2, 5)
+
+        # target qubits indices
+        qubits = tuple(random.sample(range(N_qubits), 2))
+
+        # create random density matrix from statevector
+        # random statevector to compare to
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+        # build DensityMatrix
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+
+        # create random channel
+        # for 2 qubits, rank between 1 and 16
+        # number of qubits it acts on
+        nqb = 2
+        dim = 2**nqb
+        rk = np.random.randint(1, dim**2 + 1)
+        channel = randobj.rand_channel_kraus(dim=dim, rank=rk)
+
+        dm.apply_channel(channel, qubits)
+
+        # initialize. NOT a DM object, just a matrix.
+        expected_dm = np.zeros((2**N_qubits, 2**N_qubits), dtype=np.complex128)
+        # reshape statevec since not in tensor format
+        for elem in channel.kraus_ops:  # kraus_ops is a list of dicts
+            psi_evolved = np.tensordot(
+                elem["operator"].reshape((2,) * 2 * nqb), psi.reshape((2,) * N_qubits), ((2, 3), qubits)
+            )
+            psi_evolved = np.moveaxis(psi_evolved, (0, 1), qubits)
+            expected_dm += elem["coef"] * np.conj(elem["coef"]) * np.outer(psi_evolved, np.conj(psi_evolved))
+
+        np.testing.assert_allclose(expected_dm.trace(), 1.0)
+        np.testing.assert_allclose(dm.rho, expected_dm)
+
+    def test_apply_channel_fail(self):
+        """
+        test apply a channel that is not a Channel object
+        """
+        N_qubits = np.random.randint(2, 5)
+        i = np.random.randint(0, N_qubits)
+
+        psi = np.random.rand(2**N_qubits) + 1j * np.random.rand(2**N_qubits)
+        psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+
+        # build DensityMatrix
+        dm = DensityMatrix(data=np.outer(psi, psi.conj()))
+
+        with self.assertRaises(TypeError):
+            dm.apply_channel("a", [i])
+
+
+class DensityMatrixBackendTest(unittest.TestCase):
+    """Test for DensityMatrixBackend class."""
+
+    def test_init_fail(self):
+        with self.assertRaises(TypeError):
+            DensityMatrixBackend()
+
+    def test_init_success(self):
+        circ = Circuit(1)
+        circ.rx(0, np.pi / 2)
+        pattern = circ.transpile()
+        backend = DensityMatrixBackend(pattern)
+        assert backend.pattern == pattern
+        assert backend.results == pattern.results
+        assert backend.state is None
+        assert backend.node_index == []
+        assert backend.Nqubit == 0
+        assert backend.max_qubit_num == 12
+
+    def test_add_nodes(self):
+        circ = Circuit(1)
+        pattern = circ.transpile()
+        backend = DensityMatrixBackend(pattern)
+        backend.add_nodes([0, 1])
+        expected_matrix = np.array([0.25] * 16).reshape(4, 4)
+        np.testing.assert_allclose(backend.state.rho, expected_matrix)
+
+    def test_entangle_nodes(self):
+        circ = Circuit(1)
+        pattern = circ.transpile()
+        backend = DensityMatrixBackend(pattern)
+        backend.add_nodes([0, 1])
+        backend.entangle_nodes((0, 1))
+        expected_matrix = np.array([[1, 1, 1, -1], [1, 1, 1, -1], [1, 1, 1, -1], [-1, -1, -1, 1]]) / 4
+        np.testing.assert_allclose(backend.state.rho, expected_matrix)
+
+        backend.entangle_nodes((0, 1))
+        np.testing.assert_allclose(backend.state.rho, np.array([0.25] * 16).reshape(4, 4))
+
+    def test_measure(self):
+        circ = Circuit(1)
+        circ.rx(0, np.pi / 2)
+        pattern = circ.transpile()
+
+        backend = DensityMatrixBackend(pattern)
+        backend.add_nodes([0, 1, 2])
+        backend.entangle_nodes((0, 1))
+        backend.entangle_nodes((1, 2))
+        backend.measure(backend.pattern[-4])
+
+        expected_matrix_1 = np.kron(np.array([[1, 0], [0, 0]]), np.ones((2, 2)) / 2)
+        expected_matrix_2 = np.kron(np.array([[0, 0], [0, 1]]), np.array([[0.5, -0.5], [-0.5, 0.5]]))
+        assert np.allclose(backend.state.rho, expected_matrix_1) or np.allclose(backend.state.rho, expected_matrix_2)
+
+    def test_measure_pr_calc(self):
+        # circuit there just to provide a measurement command to try out. Weird.
+        circ = Circuit(1)
+        circ.rx(0, np.pi / 2)
+        pattern = circ.transpile()
+
+        backend = DensityMatrixBackend(pattern, pr_calc=True)
+        backend.add_nodes([0, 1, 2])
+        backend.entangle_nodes((0, 1))
+        backend.entangle_nodes((1, 2))
+        backend.measure(backend.pattern[-4])
+
+        # 3-qubit linear graph state: |+0+> + |-1->
+        expected_matrix_1 = np.kron(np.array([[1, 0], [0, 0]]), np.ones((2, 2)) / 2)
+        expected_matrix_2 = np.kron(np.array([[0, 0], [0, 1]]), np.array([[0.5, -0.5], [-0.5, 0.5]]))
+
+        assert np.allclose(backend.state.rho, expected_matrix_1) or np.allclose(backend.state.rho, expected_matrix_2)
+
+    def test_correct_byproduct(self):
+        np.random.seed(0)
+
+        circ = Circuit(1)
+        circ.rx(0, np.pi / 2)
+        pattern = circ.transpile()
+
+        backend = DensityMatrixBackend(pattern)
+        backend.add_nodes([0, 1, 2])
+        backend.entangle_nodes((0, 1))
+        backend.entangle_nodes((1, 2))
+        backend.measure(backend.pattern[-4])
+        backend.measure(backend.pattern[-3])
+        backend.correct_byproduct(backend.pattern[-2])
+        backend.correct_byproduct(backend.pattern[-1])
+        backend.finalize()
+        rho = backend.state.rho
+
+        backend = StatevectorBackend(pattern)
+        backend.add_nodes([0, 1, 2])
+        backend.entangle_nodes((0, 1))
+        backend.entangle_nodes((1, 2))
+        backend.measure(backend.pattern[-4])
+        backend.measure(backend.pattern[-3])
+        backend.correct_byproduct(backend.pattern[-2])
+        backend.correct_byproduct(backend.pattern[-1])
+        backend.finalize()
+        psi = backend.state.psi
+
+        np.testing.assert_allclose(rho, np.outer(psi, psi.conj()))
+
+
+if __name__ == "__main__":
+    np.random.seed(32)
+    unittest.main()
diff --git a/tests/test_kraus.py b/tests/test_kraus.py
new file mode 100644
index 0000000..6518150
--- /dev/null
+++ b/tests/test_kraus.py
@@ -0,0 +1,232 @@
+import unittest
+
+import numpy as np
+
+import graphix.random_objects as randobj
+from graphix.channels import (
+    KrausChannel,
+    two_qubit_depolarising_tensor_channel,
+    two_qubit_depolarising_channel,
+    dephasing_channel,
+    depolarising_channel,
+)
+from graphix.ops import Ops
+
+
+class TestChannel(unittest.TestCase):
+    """Tests for Channel class"""
+
+    def test_init_with_data_success(self):
+        "test for successful intialization"
+
+        prob = np.random.rand()
+        mychannel = KrausChannel(
+            [
+                {"coef": np.sqrt(1 - prob), "operator": np.array([[1.0, 0.0], [0.0, 1.0]])},
+                {"coef": np.sqrt(prob), "operator": np.array([[1.0, 0.0], [0.0, -1.0]])},
+            ]
+        )
+        assert isinstance(mychannel.nqubit, int)
+        assert mychannel.nqubit == 1
+        assert mychannel.size == 2
+        assert isinstance(mychannel.kraus_ops, (list, np.ndarray, tuple))
+        assert mychannel.is_normalized
+
+    def test_init_with_data_fail(self):
+        "test for unsuccessful intialization"
+
+        prob = np.random.rand()
+
+        # empty data
+        with self.assertRaises(ValueError):
+            mychannel = KrausChannel([])
+
+        # incorrect parameter type
+        with self.assertRaises(TypeError):
+            mychannel = KrausChannel("a")
+
+        # incorrect "parameter" key
+        with self.assertRaises(KeyError):
+            mychannel = KrausChannel(
+                [
+                    {"coefficients": np.sqrt(1 - prob), "operator": np.array([[1.0, 0.0], [0.0, 1.0]])},
+                    {"coef": np.sqrt(prob), "operator": np.array([[1.0, 0.0], [0.0, -1.0]])},
+                ]
+            )
+
+        # incorrect "operator" key
+        with self.assertRaises(KeyError):
+            mychannel = KrausChannel(
+                [
+                    {"coef": np.sqrt(1 - prob), "oertor": np.array([[1.0, 0.0], [0.0, 1.0]])},
+                    {"coef": np.sqrt(prob), "operator": np.array([[1.0, 0.0], [0.0, -1.0]])},
+                ]
+            )
+
+        # incorrect parameter type
+        with self.assertRaises(TypeError):
+            mychannel = KrausChannel(
+                [
+                    {"coef": "a", "operator": np.array([[1.0, 0.0], [0.0, 1.0]])},
+                    {"coef": np.sqrt(prob), "operator": np.array([[1.0, 0.0], [0.0, -1.0]])},
+                ]
+            )
+
+        # incorrect operator type
+        with self.assertRaises(TypeError):
+            mychannel = KrausChannel(
+                [
+                    {"coef": np.sqrt(1 - prob), "operator": "a"},
+                    {"coef": np.sqrt(prob), "operator": np.array([[1.0, 0.0], [0.0, -1.0]])},
+                ]
+            )
+
+        # incorrect operator dimension
+        with self.assertRaises(ValueError):
+            mychannel = KrausChannel(
+                [
+                    {"coef": np.sqrt(1 - prob), "operator": np.array([1.0, 0.0])},
+                    {"coef": np.sqrt(prob), "operator": np.array([[1.0, 0.0], [0.0, -1.0]])},
+                ]
+            )
+
+        # incorrect operator dimension: square but not qubits
+        with self.assertRaises(ValueError):
+            mychannel = KrausChannel(
+                [
+                    {"coef": np.sqrt(1 - prob), "operator": np.random.rand(3, 3)},
+                    {"coef": np.sqrt(prob), "operator": np.random.rand(3, 3)},
+                ]
+            )
+
+        # doesn't square to 1. Not normalized. Parameter.
+        with self.assertRaises(ValueError):
+            mychannel = KrausChannel(
+                [
+                    {"coef": 2 * np.sqrt(1 - prob), "operator": np.array([[1.0, 0.0], [0.0, 1.0]])},
+                    {"coef": np.sqrt(prob), "operator": np.array([[1.0, 0.0], [0.0, -1.0]])},
+                ]
+            )
+
+        # doesn't square to 1. Not normalized. Operator.
+        with self.assertRaises(ValueError):
+            mychannel = KrausChannel(
+                [
+                    {"coef": np.sqrt(1 - prob), "operator": np.array([[1.0, 0.0], [0.0, 1.0]])},
+                    {"coef": np.sqrt(prob), "operator": np.array([[1.0, 3.0], [0.0, -1.0]])},
+                ]
+            )
+
+        # incorrect rank (number of kraus_operators)
+        # use a random channel to do that.
+        with self.assertRaises(ValueError):
+            randobj.rand_channel_kraus(dim=2**2, rank=20)
+
+    def test_dephasing_channel(self):
+
+        prob = np.random.rand()
+        data = [
+            {"coef": np.sqrt(1 - prob), "operator": np.array([[1.0, 0.0], [0.0, 1.0]])},
+            {"coef": np.sqrt(prob), "operator": Ops.z},
+        ]
+        dephase_channel = dephasing_channel(prob)
+        assert isinstance(dephase_channel, KrausChannel)
+        assert dephase_channel.nqubit == 1
+        assert dephase_channel.size == 2
+        assert dephase_channel.is_normalized
+
+        for i in range(len(dephase_channel.kraus_ops)):
+            np.testing.assert_allclose(dephase_channel.kraus_ops[i]["coef"], data[i]["coef"])
+            np.testing.assert_allclose(dephase_channel.kraus_ops[i]["operator"], data[i]["operator"])
+
+    def test_depolarising_channel(self):
+
+        prob = np.random.rand()
+        data = [
+            {"coef": np.sqrt(1 - prob), "operator": np.eye(2)},
+            {"coef": np.sqrt(prob / 3.0), "operator": Ops.x},
+            {"coef": np.sqrt(prob / 3.0), "operator": Ops.y},
+            {"coef": np.sqrt(prob / 3.0), "operator": Ops.z},
+        ]
+
+        depol_channel = depolarising_channel(prob)
+
+        assert isinstance(depol_channel, KrausChannel)
+        assert depol_channel.nqubit == 1
+        assert depol_channel.size == 4
+        assert depol_channel.is_normalized
+
+        for i in range(len(depol_channel.kraus_ops)):
+            np.testing.assert_allclose(depol_channel.kraus_ops[i]["coef"], data[i]["coef"])
+            np.testing.assert_allclose(depol_channel.kraus_ops[i]["operator"], data[i]["operator"])
+
+    def test_2_qubit_depolarising_channel(self):
+
+        prob = np.random.rand()
+        data = [
+            {"coef": np.sqrt(1 - prob), "operator": np.kron(np.eye(2), np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, Ops.x)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, np.eye(2))},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(np.eye(2), Ops.x)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(np.eye(2), Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(np.eye(2), Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.x, Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, Ops.x)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.y, Ops.z)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, Ops.y)},
+            {"coef": np.sqrt(prob / 15.0), "operator": np.kron(Ops.z, Ops.x)},
+        ]
+
+        depol_channel_2_qubit = two_qubit_depolarising_channel(prob)
+
+        assert isinstance(depol_channel_2_qubit, KrausChannel)
+        assert depol_channel_2_qubit.nqubit == 2
+        assert depol_channel_2_qubit.size == 16
+        assert depol_channel_2_qubit.is_normalized
+
+        for i in range(len(depol_channel_2_qubit.kraus_ops)):
+            np.testing.assert_allclose(depol_channel_2_qubit.kraus_ops[i]["coef"], data[i]["coef"])
+            np.testing.assert_allclose(depol_channel_2_qubit.kraus_ops[i]["operator"], data[i]["operator"])
+
+    def test_2_qubit_depolarising_tensor_channel(self):
+
+        prob = np.random.rand()
+        data = [
+            {"coef": 1 - prob, "operator": np.kron(np.eye(2), np.eye(2))},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.x, Ops.x)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.y, Ops.y)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.z, Ops.z)},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(Ops.x, np.eye(2))},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(Ops.y, np.eye(2))},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(Ops.z, np.eye(2))},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(np.eye(2), Ops.x)},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(np.eye(2), Ops.y)},
+            {"coef": np.sqrt(1 - prob) * np.sqrt(prob / 3.0), "operator": np.kron(np.eye(2), Ops.z)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.x, Ops.y)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.x, Ops.z)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.y, Ops.x)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.y, Ops.z)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.z, Ops.x)},
+            {"coef": prob / 3.0, "operator": np.kron(Ops.z, Ops.y)},
+        ]
+
+        depol_tensor_channel_2_qubit = two_qubit_depolarising_tensor_channel(prob)
+
+        assert isinstance(depol_tensor_channel_2_qubit, KrausChannel)
+        assert depol_tensor_channel_2_qubit.nqubit == 2
+        assert depol_tensor_channel_2_qubit.size == 16
+        assert depol_tensor_channel_2_qubit.is_normalized
+
+        for i in range(len(depol_tensor_channel_2_qubit.kraus_ops)):
+            np.testing.assert_allclose(depol_tensor_channel_2_qubit.kraus_ops[i]["coef"], data[i]["coef"])
+            np.testing.assert_allclose(depol_tensor_channel_2_qubit.kraus_ops[i]["operator"], data[i]["operator"])
+
+
+if __name__ == "__main__":
+    np.random.seed(2)
+    unittest.main()
diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
new file mode 100644
index 0000000..41fde7d
--- /dev/null
+++ b/tests/test_noisy_density_matrix.py
@@ -0,0 +1,457 @@
+import unittest
+
+import numpy as np
+
+from graphix import Circuit
+from graphix.noise_models.noise_model import NoiseModel
+from graphix.channels import (
+    KrausChannel,
+    depolarising_channel,
+    two_qubit_depolarising_tensor_channel,
+    two_qubit_depolarising_channel,
+)
+from graphix.noise_models.noiseless_noise_model import NoiselessNoiseModel
+from graphix.ops import Ops
+
+
+class TestNoiseModel(NoiseModel):
+    """Test noise model for testing.
+    Only return the identity channel.
+
+    :param NoiseModel: Parent abstract class class:`graphix.noise_model.NoiseModel`
+    :type NoiseModel: class
+    """
+
+    def __init__(
+        self,
+        prepare_error_prob=0.0,
+        x_error_prob=0.0,
+        z_error_prob=0.0,
+        entanglement_error_prob=0.0,
+        measure_channel_prob=0.0,
+        measure_error_prob=0.0,
+    ):
+        self.prepare_error_prob = prepare_error_prob
+        self.x_error_prob = x_error_prob
+        self.z_error_prob = z_error_prob
+        self.entanglement_error_prob = entanglement_error_prob
+        self.measure_error_prob = measure_error_prob
+        self.measure_channel_prob = measure_channel_prob
+
+    def prepare_qubit(self):
+        """return the channel to apply after clean single-qubit preparation. Here just identity."""
+        return depolarising_channel(self.prepare_error_prob)
+
+    def entangle(self):
+        """return noise model to qubits that happens after the CZ gate"""
+        # return two_qubit_depolarising_tensor_channel(self.entanglement_error_prob)
+        return two_qubit_depolarising_channel(self.entanglement_error_prob)
+
+    def measure(self):
+        """apply noise to qubit to be measured."""
+        return depolarising_channel(self.measure_channel_prob)
+
+    def confuse_result(self, cmd):
+        """assign wrong measurement result
+        cmd = "M"
+        """
+
+        if np.random.rand() < self.measure_error_prob:
+            self.simulator.results[cmd[1]] = 1 - self.simulator.results[cmd[1]]
+
+    def byproduct_x(self):
+        """apply noise to qubits after X gate correction"""
+        return depolarising_channel(self.x_error_prob)
+
+    def byproduct_z(self):
+        """apply noise to qubits after Z gate correction"""
+        return depolarising_channel(self.z_error_prob)
+
+    def clifford(self):
+        """apply noise to qubits that happens in the Clifford gate process"""
+        # TODO list separate different Cliffords to allow customization
+        return KrausChannel([{"coef": 1.0, "operator": np.eye(2)}])
+
+    def tick_clock(self):
+        """notion of time in real devices - this is where we apply effect of T1 and T2.
+        we assume commands that lie between 'T' commands run simultaneously on the device.
+        """
+        pass
+
+
+class NoisyDensityMatrixBackendTest(unittest.TestCase):
+    """Test for Noisy DensityMatrixBackend simultation."""
+
+    def setUp(self):
+        # set up the random numbers
+        self.rng = np.random.default_rng()  # seed=422
+        np.random.seed()  # 42 vs 422
+
+        # Hadamard pattern
+        ncirc = Circuit(1)
+        ncirc.h(0)
+        self.hadamardpattern = ncirc.transpile()
+
+        # Rz(alpha) pattern
+        self.alpha = 2 * np.pi * self.rng.random()
+        print(f"alpha is {self.alpha}")
+        circ = Circuit(1)
+        circ.rz(0, self.alpha)
+        self.rzpattern = circ.transpile()
+        self.rz_exact_res = 0.5 * np.array([[1.0, np.exp(-1j * self.alpha)], [np.exp(1j * self.alpha), 1.0]])
+
+    # test noiseless noisy vs noiseless
+    def test_noiseless_noisy_hadamard(self):
+
+        noiselessres = self.hadamardpattern.simulate_pattern(backend="densitymatrix")
+        # noiseless noise model
+        noisynoiselessres = self.hadamardpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=NoiselessNoiseModel()
+        )
+        np.testing.assert_allclose(noiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]]))
+        # result should be |0>
+        np.testing.assert_allclose(noisynoiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]]))
+
+    # test measurement confuse outcome
+    def test_noisy_measure_confuse_hadamard(self):
+
+        res = self.hadamardpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(measure_error_prob=1.0)
+        )
+        # result should be |1>
+        np.testing.assert_allclose(res.rho, np.array([[0.0, 0.0], [0.0, 1.0]]))
+
+        # arbitrary probability
+        measure_error_pr = self.rng.random()
+        print(f"measure_error_pr = {measure_error_pr}")
+        res = self.hadamardpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(measure_error_prob=measure_error_pr)
+        )
+        # result should be |1>
+        assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) or np.allclose(
+            res.rho, np.array([[0.0, 0.0], [0.0, 1.0]])
+        )
+
+    def test_noisy_measure_channel_hadamard(self):
+
+        measure_channel_pr = self.rng.random()
+        print(f"measure_channel_pr = {measure_channel_pr}")
+        # measurement error only
+        res = self.hadamardpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(measure_channel_prob=measure_channel_pr)
+        )
+        # just TP the depolarizing channel
+        np.testing.assert_allclose(
+            res.rho, np.array([[1 - 2 * measure_channel_pr / 3.0, 0.0], [0.0, 2 * measure_channel_pr / 3.0]])
+        )
+
+    # test Pauli X error
+    def test_noisy_X_hadamard(self):
+
+        # x error only
+        x_error_pr = self.rng.random()
+        print(f"x_error_pr = {x_error_pr}")
+        res = self.hadamardpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(x_error_prob=x_error_pr)
+        )
+        # analytical result since deterministic pattern output is |0>.
+        # if no X applied, no noise. If X applied X noise on |0><0|
+
+        assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) or np.allclose(
+            res.rho, np.array([[1 - 2 * x_error_pr / 3.0, 0.0], [0.0, 2 * x_error_pr / 3.0]])
+        )
+
+    # test entanglement error
+    def test_noisy_entanglement_hadamard(self):
+
+        entanglement_error_pr = np.random.rand()
+        res = self.hadamardpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(entanglement_error_prob=entanglement_error_pr)
+        )
+        # analytical result for tensor depolarizing channel
+        # np.testing.assert_allclose(
+        #     res.rho,
+        #     np.array(
+        #         [
+        #             [1 - 4 * entanglement_error_pr / 3.0 + 8 * entanglement_error_pr**2 / 9, 0.0],
+        #             [0.0, 4 * entanglement_error_pr / 3.0 - 8 * entanglement_error_pr**2 / 9],
+        #         ]
+        #     ),
+        # )
+
+        # analytical result for true 2-qubit depolarizing channel
+        np.testing.assert_allclose(
+            res.rho,
+            np.array(
+                [
+                    [1 - 8 * entanglement_error_pr / 15.0, 0.0],
+                    [0.0, 8 * entanglement_error_pr / 15.0],
+                ]
+            ),
+        )
+
+    # test preparation error
+    def test_noisy_preparation_hadamard(self):
+
+        prepare_error_pr = self.rng.random()
+        print(f"prepare_error_pr = {prepare_error_pr}")
+        res = self.hadamardpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(prepare_error_prob=prepare_error_pr)
+        )
+        # analytical result
+        np.testing.assert_allclose(
+            res.rho, np.array([[1 - 2 * prepare_error_pr / 3.0, 0.0], [0.0, 2 * prepare_error_pr / 3.0]])
+        )
+
+    ###  Test rz gate
+
+    # test noiseless noisy vs noiseless
+    def test_noiseless_noisy_rz(self):
+
+        noiselessres = self.rzpattern.simulate_pattern(backend="densitymatrix")
+        # noiseless noise model or TestNoiseModel() since all probas are 0
+        noisynoiselessres = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel()
+        )  # NoiselessNoiseModel()
+        np.testing.assert_allclose(
+            noiselessres.rho, 0.5 * np.array([[1.0, np.exp(-1j * self.alpha)], [np.exp(1j * self.alpha), 1.0]])
+        )
+        # result should be |0>
+        np.testing.assert_allclose(
+            noisynoiselessres.rho, 0.5 * np.array([[1.0, np.exp(-1j * self.alpha)], [np.exp(1j * self.alpha), 1.0]])
+        )
+
+    # test preparation error
+    def test_noisy_preparation_rz(self):
+
+        prepare_error_pr = self.rng.random()
+        print(f"prepare_error_pr = {prepare_error_pr}")
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(prepare_error_prob=prepare_error_pr)
+        )
+        # analytical result
+        np.testing.assert_allclose(
+            res.rho,
+            0.5
+            * np.array(
+                [
+                    [
+                        1.0,
+                        (3 - 4 * prepare_error_pr) ** 2
+                        * (3 * np.cos(self.alpha) + 1j * (-3 + 4 * prepare_error_pr) * np.sin(self.alpha))
+                        / 27,
+                    ],
+                    [
+                        (3 - 4 * prepare_error_pr) ** 2
+                        * (3 * np.cos(self.alpha) - 1j * (-3 + 4 * prepare_error_pr) * np.sin(self.alpha))
+                        / 27,
+                        1.0,
+                    ],
+                ]
+            ),
+        )
+
+    # test entanglement error
+    def test_noisy_entanglement_rz(self):
+
+        entanglement_error_pr = np.random.rand()
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(entanglement_error_prob=entanglement_error_pr)
+        )
+        # analytical result for tensor depolarizing channel
+        # np.testing.assert_allclose(
+        #     res.rho,
+        #     0.5
+        #     * np.array(
+        #         [
+        #             [
+        #                 1.0,
+        #                 (-3 + 4 * entanglement_error_pr) ** 3
+        #                 * (-3 * np.cos(self.alpha) + 1j * (3 - 4 * entanglement_error_pr) * np.sin(self.alpha))
+        #                 / 81,
+        #             ],
+        #             [
+        #                 (-3 + 4 * entanglement_error_pr) ** 3
+        #                 * (-3 * np.cos(self.alpha) - 1j * (3 - 4 * entanglement_error_pr) * np.sin(self.alpha))
+        #                 / 81,
+        #                 1.0,
+        #             ],
+        #         ]
+        #     ),
+        # )
+
+        # analytical result for true 2-qubit depolarizing channel
+        np.testing.assert_allclose(
+            res.rho,
+            0.5
+            * np.array(
+                [
+                    [
+                        1.0,
+                        np.exp(-1j * self.alpha) * (15 - 16 * entanglement_error_pr) ** 2 / 225,
+                    ],
+                    [
+                        np.exp(1j * self.alpha) * (15 - 16 * entanglement_error_pr) ** 2 / 225,
+                        1.0,
+                    ],
+                ]
+            ),
+        )
+
+    def test_noisy_measure_channel_rz(self):
+
+        measure_channel_pr = self.rng.random()
+        print(f"measure_channel_pr = {measure_channel_pr}")
+        # measurement error only
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(measure_channel_prob=measure_channel_pr)
+        )
+
+        np.testing.assert_allclose(
+            res.rho,
+            0.5
+            * np.array(
+                [
+                    [
+                        1.0,
+                        (-3 + 4 * measure_channel_pr)
+                        * (-3 * np.cos(self.alpha) + 1j * (3 - 4 * measure_channel_pr) * np.sin(self.alpha))
+                        / 9,
+                    ],
+                    [
+                        (-3 + 4 * measure_channel_pr)
+                        * (-3 * np.cos(self.alpha) - 1j * (3 - 4 * measure_channel_pr) * np.sin(self.alpha))
+                        / 9,
+                        1.0,
+                    ],
+                ]
+            ),
+        )
+
+    def test_noisy_X_rz(self):
+
+        # x error only
+        x_error_pr = self.rng.random()
+        print(f"x_error_pr = {x_error_pr}")
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(x_error_prob=x_error_pr)
+        )
+
+        # only two cases: if no X correction, Z or no Z correction but exact result.
+        # If X correction the noise result is the same with or without the PERFECT Z correction.
+        assert np.allclose(
+            res.rho, 0.5 * np.array([[1.0, np.exp(-1j * self.alpha)], [np.exp(1j * self.alpha), 1.0]])
+        ) or np.allclose(
+            res.rho,
+            0.5
+            * np.array(
+                [
+                    [1.0, np.exp(-1j * self.alpha) * (3 - 4 * x_error_pr) / 3],
+                    [np.exp(1j * self.alpha) * (3 - 4 * x_error_pr) / 3, 1.0],
+                ]
+            ),
+        )
+
+    def test_noisy_Z_rz(self):
+
+        # z error only
+        z_error_pr = self.rng.random()
+        print(f"z_error_pr = {z_error_pr}")
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(z_error_prob=z_error_pr)
+        )
+
+        # only two cases: if no Z correction, X or no X correction but exact result.
+        # If Z correction the noise result is the same with or without the PERFECT X correction.
+        assert np.allclose(
+            res.rho, 0.5 * np.array([[1.0, np.exp(-1j * self.alpha)], [np.exp(1j * self.alpha), 1.0]])
+        ) or np.allclose(
+            res.rho,
+            0.5
+            * np.array(
+                [
+                    [1.0, np.exp(-1j * self.alpha) * (3 - 4 * z_error_pr) / 3],
+                    [np.exp(1j * self.alpha) * (3 - 4 * z_error_pr) / 3, 1.0],
+                ]
+            ),
+        )
+
+    def test_noisy_XZ_rz(self):
+
+        # x and z errors
+        x_error_pr = self.rng.random()
+        print(f"x_error_pr = {x_error_pr}")
+        z_error_pr = self.rng.random()
+        print(f"z_error_pr = {z_error_pr}")
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(x_error_prob=x_error_pr, z_error_prob=z_error_pr)
+        )
+
+        # 4 cases : no corr, noisy X, noisy Z, noisy XZ.
+        assert (
+            np.allclose(res.rho, 0.5 * np.array([[1.0, np.exp(-1j * self.alpha)], [np.exp(1j * self.alpha), 1.0]]))
+            or np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * self.alpha) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9],
+                        [np.exp(1j * self.alpha) * (3 - 4 * x_error_pr) * (3 - 4 * z_error_pr) / 9, 1.0],
+                    ]
+                ),
+            )
+            or np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * self.alpha) * (3 - 4 * z_error_pr) / 3],
+                        [np.exp(1j * self.alpha) * (3 - 4 * z_error_pr) / 3, 1.0],
+                    ]
+                ),
+            )
+            or np.allclose(
+                res.rho,
+                0.5
+                * np.array(
+                    [
+                        [1.0, np.exp(-1j * self.alpha) * (3 - 4 * x_error_pr) / 3],
+                        [np.exp(1j * self.alpha) * (3 - 4 * x_error_pr) / 3, 1.0],
+                    ]
+                ),
+            )
+        )
+
+    # test measurement confuse outcome
+    def test_noisy_measure_confuse_rz(self):
+
+        # probability 1 to shift both outcome
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(measure_error_prob=1.0)
+        )
+        # result X, XZ or Z
+
+        assert (
+            np.allclose(res.rho, Ops.x @ self.rz_exact_res @ Ops.x)
+            or np.allclose(res.rho, Ops.z @ self.rz_exact_res @ Ops.z)
+            or np.allclose(res.rho, Ops.z @ Ops.x @ self.rz_exact_res @ Ops.x @ Ops.z)
+        )
+
+        # arbitrary probability
+        measure_error_pr = self.rng.random()
+        print(f"measure_error_pr = {measure_error_pr}")
+        res = self.rzpattern.simulate_pattern(
+            backend="densitymatrix", noise_model=TestNoiseModel(measure_error_prob=measure_error_pr)
+        )
+        # just add the case without readout errors
+        assert (
+            np.allclose(res.rho, self.rz_exact_res)
+            or np.allclose(res.rho, Ops.x @ self.rz_exact_res @ Ops.x)
+            or np.allclose(res.rho, Ops.z @ self.rz_exact_res @ Ops.z)
+            or np.allclose(res.rho, Ops.z @ Ops.x @ self.rz_exact_res @ Ops.x @ Ops.z)
+        )
+
+
+if __name__ == "__main__":
+    np.random.seed(32)
+    unittest.main()
diff --git a/tests/test_random_utilities.py b/tests/test_random_utilities.py
new file mode 100644
index 0000000..9f68be2
--- /dev/null
+++ b/tests/test_random_utilities.py
@@ -0,0 +1,195 @@
+import unittest
+
+import numpy as np
+
+import graphix.random_objects as randobj
+from graphix.channels import KrausChannel
+from graphix.linalg_validations import check_data_dims, check_hermitian, check_psd, check_square, check_unit_trace
+from graphix.sim.density_matrix import DensityMatrix
+from graphix.ops import Ops
+
+
+class TestUtilities(unittest.TestCase):
+    def test_rand_herm(self):
+        tmp = randobj.rand_herm(np.random.randint(2, 20))
+        np.testing.assert_allclose(tmp, tmp.conj().T)
+
+    # TODO : work on that. Verify an a random vector and not at the operator level...
+
+    def test_rand_unit(self):
+        d = np.random.randint(2, 20)
+        tmp = randobj.rand_unit(d)
+        print(type(tmp), tmp.dtype)
+
+        # different default values for testing.assert_allclose and all_close!
+        np.testing.assert_allclose(tmp @ tmp.conj().T, np.eye(d), atol=1e-15)
+        np.testing.assert_allclose(tmp.conj().T @ tmp, np.eye(d), atol=1e-15)
+
+    def test_random_channel_success(self):
+
+        nqb = np.random.randint(1, 5)
+        dim = 2**nqb  # np.random.randint(2, 8)
+
+        # no rank feature
+        channel = randobj.rand_channel_kraus(dim=dim)
+
+        assert isinstance(channel, KrausChannel)
+        assert check_data_dims(channel.kraus_ops)
+        assert channel.kraus_ops[0]["operator"].shape == (dim, dim)
+        assert channel.nqubit == nqb
+        assert channel.size == dim**2
+        assert channel.is_normalized
+
+        rk = np.random.randint(1, dim**2 + 1)
+        channel = randobj.rand_channel_kraus(dim=dim, rank=rk)
+
+        assert isinstance(channel, KrausChannel)
+        assert check_data_dims(channel.kraus_ops)
+        assert channel.kraus_ops[0]["operator"].shape == (dim, dim)
+        assert channel.nqubit == nqb
+        assert channel.size == rk
+        assert channel.is_normalized
+
+    def test_random_channel_fail(self):
+
+        # incorrect rank type
+        with self.assertRaises(TypeError):
+            mychannel = randobj.rand_channel_kraus(dim=2**2, rank=3.0)
+
+        # null rank
+        with self.assertRaises(ValueError):
+            mychannel = randobj.rand_channel_kraus(dim=2**2, rank=0)
+
+    def test_rand_gauss_cpx(self):
+
+        nsample = int(1e4)
+
+        dim = np.random.randint(2, 20)
+        tmp = [randobj.rand_gauss_cpx_mat(dim=dim) for _ in range(nsample)]
+
+        dimset = {i.shape for i in tmp}
+        assert len(dimset) == 1
+        assert list(dimset)[0] == (dim, dim)
+
+    def test_check_psd_success(self):
+
+        # Generate a random mixed state from state vectors with same probability
+        # We know this is PSD
+
+        nqb = np.random.randint(2, 7)
+
+        dim = 2**nqb
+        m = np.random.randint(1, dim)
+
+        dm = np.zeros((dim,) * 2, dtype=np.complex128)
+
+        # TODO optimize that
+        for _ in range(m):
+            psi = np.random.rand(dim) + 1j * np.random.rand(dim)
+            psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
+            dm += np.outer(psi, psi.conj()) / m
+
+        assert check_psd(dm)
+
+    def test_check_psd_fail(self):
+
+        # not hermitian
+        # don't use dim = 2, too easy to have a PSD matrix.
+        # NOTE useless test since eigvalsh treats the matrix as hermitian and takes only the L or U part
+
+        l = np.random.randint(5, 20)
+
+        mat = np.random.rand(l, l) + 1j * np.random.rand(l, l)
+
+        # eigvalsh doesn't raise a LinAlgError since just use upper or lower part of the matrix.
+        # instead Value error
+        with self.assertRaises(ValueError):
+            check_psd(mat)
+
+        # hermitian but not positive eigenvalues
+        mat = randobj.rand_herm(l)
+
+        with self.assertRaises(ValueError):
+            check_psd(mat)
+
+    def test_rand_dm(self):
+        # needs to be power of 2 dimension since builds a DM object
+        dm = randobj.rand_dm(2 ** np.random.randint(2, 5))
+
+        assert isinstance(dm, DensityMatrix)
+        assert check_square(dm.rho)
+        assert check_hermitian(dm.rho)
+        assert check_psd(dm.rho)
+        assert check_unit_trace(dm.rho)
+
+    # try with incorrect dimension
+    def test_rand_dm_fail(self):
+        with self.assertRaises(ValueError):
+            dm = randobj.rand_dm(2 ** np.random.randint(2, 5) + 1)
+
+    def test_rand_dm_rank(self):
+
+        rk = 3
+        dm = randobj.rand_dm(2 ** np.random.randint(2, 5), rank=rk)
+
+        assert isinstance(dm, DensityMatrix)
+        assert check_square(dm.rho)
+        assert check_hermitian(dm.rho)
+        assert check_psd(dm.rho)
+        assert check_unit_trace(dm.rho)
+
+        evals = np.linalg.eigvalsh(dm.rho)
+
+        evals[np.abs(evals) < 1e-15] = 0
+
+        assert rk == np.count_nonzero(evals)
+
+    # TODO move that somewhere else?
+    def test_pauli_tensor_ops(self):
+        nqb = np.random.randint(2, 6)
+        Pauli_tensor_ops = Ops.build_tensor_Pauli_ops(nqb)
+
+        assert len(Pauli_tensor_ops) == 4**nqb
+
+        dims = np.array([i.shape for i in Pauli_tensor_ops])
+        # or np.apply_along_axis ?
+        assert np.all(dims == (2**nqb, 2**nqb))
+
+    def test_pauli_tensor_ops_fail(self):
+
+        with self.assertRaises(TypeError):
+            Pauli_tensor_ops = Ops.build_tensor_Pauli_ops(np.random.randint(2, 6) + 0.5)
+
+        with self.assertRaises(ValueError):
+            Pauli_tensor_ops = Ops.build_tensor_Pauli_ops(0)
+
+    def test_random_pauli_channel_success(self):
+
+        nqb = np.random.randint(2, 6)
+        rk = np.random.randint(1, 2**nqb + 1)
+        Pauli_channel = randobj.rand_Pauli_channel_kraus(dim=2**nqb, rank=rk)  # default is full rank
+
+        assert isinstance(Pauli_channel, KrausChannel)
+        assert Pauli_channel.nqubit == nqb
+        assert Pauli_channel.size == rk
+        assert Pauli_channel.is_normalized
+
+    def test_random_pauli_channel_fail(self):
+        nqb = 3
+        rk = 2
+        with self.assertRaises(TypeError):
+            dm = randobj.rand_Pauli_channel_kraus(dim=2**nqb, rank=rk + 0.5)
+
+        with self.assertRaises(ValueError):
+            dm = randobj.rand_Pauli_channel_kraus(dim=2**nqb + 0.5, rank=rk)
+
+        with self.assertRaises(ValueError):
+            dm = randobj.rand_Pauli_channel_kraus(dim=2**nqb, rank=-3)
+
+        with self.assertRaises(ValueError):
+            dm = randobj.rand_Pauli_channel_kraus(dim=2**nqb + 1, rank=rk)
+
+
+if __name__ == "__main__":
+    np.random.seed(2)
+    unittest.main()
diff --git a/tests/test_runner.py b/tests/test_runner.py
index 881b77e..bb3d194 100644
--- a/tests/test_runner.py
+++ b/tests/test_runner.py
@@ -17,6 +17,9 @@ from graphix.device_interface import PatternRunner
 gx_ibmq_mock = MagicMock()
 sys.modules["graphix_ibmq.runner"] = gx_ibmq_mock
 
+import tests.random_circuit as rc
+from graphix.device_interface import PatternRunner
+
 
 def modify_statevector(statevector, output_qubit):
     statevector = np.asarray(statevector)