Simulating noisy MBQC

DensityMatrix class (through 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 transpile() method. This returns Pattern object containing measurement pattern:

pattern = circuit.transpile().pattern
pattern.print_pattern(lim=30)
# pattern.draw_graph()

N, node = 2
N, node = 3
E, nodes = (0, 2)
E, nodes = (2, 3)
M, node = 0, plane = XY, angle(pi) = -0.17831963662165098, s-domain = [], t_domain = []
M, node = 2, plane = XY, angle(pi) = 0, s-domain = [], t_domain = []
X byproduct, node = 3, domain = [2]
Z byproduct, node = 3, domain = [0]
N, node = 4
N, node = 5
E, nodes = (1, 4)
E, nodes = (4, 5)
M, node = 1, plane = XY, angle(pi) = -0.09111772772805722, s-domain = [], t_domain = []
M, node = 4, plane = XY, angle(pi) = 0, s-domain = [], t_domain = []
X byproduct, node = 5, domain = [4]
Z byproduct, node = 5, domain = [1]
N, node = 6
N, node = 7
E, nodes = (5, 6)
E, nodes = (3, 6)
E, nodes = (6, 7)
M, node = 5, plane = XY, angle(pi) = 0, s-domain = [], t_domain = []
M, node = 6, plane = XY, angle(pi) = 0, s-domain = [], t_domain = []
X byproduct, node = 7, domain = [6]
Z byproduct, node = 7, domain = [5]
Z byproduct, node = 3, domain = [5]

simulate with statevector backend:

out_state = pattern.simulate_pattern(backend="statevector")
print(out_state.flatten())

from graphix.channels import KrausChannel, dephasing_channel
from graphix.noise_models.noise_model import NoiseModel

[0.5       +0.j         0.47965405+0.1411807j  0.33131835+0.37447049j
 0.42357239+0.26568106j]

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


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()))

0.9470903893460423

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()