fxeb_noisy.py

def fidelity_xeb_noisy(n_qubits: int, trials: int, n_samples: int, prob_no_error: float):
    dim = 2**n_qubits

    # keep track of ideal output probabilities
    ideal_probs = []

    # identify depolarizing operators over n-qubit space
    depolarizing_ops = []
    for x in itertools.product(paulis, repeat=n_qubits):
        op = functools.reduce(lambda a, b: np.kron(a, b), x)
        depolarizing_ops.append(op)

    # loop over random programs
    for _ in range(trials):
        unitary = random_unitary(dim)

        # sample an operator according to specified probabilities
        all_ops = [unitary] + depolarizing_ops
        probabilities = [prob_no_error] + len(depolarizing_ops)*[(1-prob_no_error)/len(depolarizing_ops)]
        op_idx = np.random.choice(len(all_ops), p=probabilities)
        op = all_ops[op_idx]

        # draw samples from the resultant state
        sample_probs = [simulate_probability(op, bb) for bb in range(dim)]
        samples = np.random.choice(dim, size=n_samples, p=sample_probs)

        # collect ideal sampling probability for these samples
        for sample in samples:
            ideal_prob = simulate_probability(unitary, sample)
            ideal_probs.append(ideal_prob)

    # calculate and return the fidelity of the XEB
    return dim*np.mean(ideal_probs) - 1