errors.py

import numpy as np
from matplotlib import pyplot as plt
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error, pauli_error, thermal_relaxation_error
from sympy.physics.quantum.qexpr import QuantumError

from grover import GroverSolver
from plot_results import plot_results
from subset_sum_oracle import SubsetSumOracle
from qiskit_ibm_runtime import QiskitRuntimeService

from verify import solve_subset_sum_with_bitstrings, check_counts_statistically_accurate

DEFAULT_BASIS_GATES_1Q: list[str] = ['h', 's', 't', 'x', 'y', 'z', 'id', 'rz', 'sx']
DEFAULT_BASIS_GATES_2Q: list[str] = ['cx']
DEFAULT_BASIS_GATES_3Q: list[str] = []

"""
Check out https://quantum.cloud.ibm.com/docs/en/guides/build-noise-models#combining-quantum-errors for the reference

Do NOT use 0 error rate if you can, pass None instead. 
I have included checks against that, but in case of some omission it yields really weird results
"""


def _apply_errors(
        noise_model: NoiseModel, error: QuantumError,
        noisy_gates: str | list[str],
        qubits: list[int] | None = None
):
    """
    Apply the errors to all qubits or just the given list of qubits, depending on whether qubits is None or not.
    :param noise_model: Noise model to update
    :param error: Error to apply
    :param noisy_gates: Gates to apply the noise to
    :param qubits: List of qubits to apply the error to. If none, error is applied to all qubits
    :return: No return value. Noise model is modified in-place.
    """
    if qubits is None:
        noise_model.add_all_qubit_quantum_error(error, noisy_gates)
    else:
        noise_model.add_quantum_error(error, noisy_gates, qubits)


def _check_basis_gates(
        basis_gates_1q: list[str] | None,
        basis_gates_2q: list[str] | None,
        basis_gates_3q: list[str] | None
) -> tuple[list[str], list[str], list[str]]:
    """
    Sanitize the basis gates input - if none, apply the default sets
    :param basis_gates_1q: Basis 1-qubit gates
    :param basis_gates_2q: Basis 2-qubit gates
    :param basis_gates_3q: Basis 3-qubit gates
    :return: Sanitized basis gates in the order of: 1-qubit, 2-qubit, 3-qubit
    """
    if basis_gates_1q is None:
        basis_gates_1q = DEFAULT_BASIS_GATES_1Q
    if basis_gates_2q is None:
        basis_gates_2q = DEFAULT_BASIS_GATES_2Q
    if basis_gates_3q is None:
        basis_gates_3q = DEFAULT_BASIS_GATES_3Q

    return basis_gates_1q, basis_gates_2q, basis_gates_3q


def model_depolarizing_error(
        err: float | None = 0.05,
        basis_gates_1q: list[str] | None = None,
        basis_gates_2q: list[str] | None = None,
        basis_gates_3q: list[str] | None = None,
        qubits: list[int] | None = None,
        noise_model: NoiseModel = None
) -> NoiseModel:
    """
    Applies depolarizing error to the noise model.

    Depolarizing error causes information loss in the qubit, creating a mixed state.
    Selecting the depolarizing channel error model introduces a random error inbetween each operation on a qubit.
    These errors are injected in a form of bit-flips (x error), phase-flips (z error) or both at the same time (y error)
    (https://www.quantum-inspire.com/kbase/cqasm-error-models)

    Simplest usage: model_depolarizing_error(err=0.05)

    :param err: Error rate, the probability of occurrence of a depolarizing error. If None, no noise is applied.
    :param basis_gates_1q: 1-qubit gates to apply the noise to. If provided with an empty list, these gates are omitted.
        The default list is defined by {DEFAULT_BASIS_GATES_1Q}
    :param basis_gates_2q: 2-qubit gates to apply the noise to. If provided with an empty list, these gates are omitted.
        The default list is defined by {DEFAULT_BASIS_GATES_2Q}
    :param basis_gates_3q: 3-qubit gates to apply the noise to. If provided with an empty list, these gates are omitted.
        The default list is defined by {DEFAULT_BASIS_GATES_3Q}
    :param qubits: List of qubits to apply the error to. If None, error is applied to all qubits
    :param noise_model: Noise model to update. If None, a new noise model is created.
    :return: an updated noise model with depolarizing error applied.
    """
    basis_gates_1q, basis_gates_2q, basis_gates_3q = _check_basis_gates(basis_gates_1q, basis_gates_2q, basis_gates_3q)
    if noise_model is None:
        noise_model = NoiseModel()

    if err is None or err == 0.0:
        return noise_model

    error = depolarizing_error(err, 1)

    if basis_gates_1q:
        _apply_errors(noise_model, error, basis_gates_1q, qubits)
    if basis_gates_2q:
        error_2q = error.tensor(error)
        _apply_errors(noise_model, error_2q, basis_gates_2q, qubits)
    if basis_gates_3q:
        error_3q = error.tensor(error).tensor(error)
        _apply_errors(noise_model, error_3q, basis_gates_3q, qubits)

    return noise_model


def _infer_qubits(
        qubits: list[int] | None,
        t1s: list[np.double] | np.double | None,
        t2s: list[np.double] | np.double | None
) -> tuple[int, list[int]]:
    """
    Infer the list of relevant qubits and their number based on the given information.
    Relevant for thermal noise methods.

    :param qubits: The list of relevant qubits. If an actual list, the list is left unchanged
        and n_qubits is computed trivially.
        Otherwise, if qubits is an int, it is treated as a number of qubits to apply the noise to.
        If it is none, the qubit number is inferred from other arguments.
    :param t1s: Thermal relaxation time T1. If an int, it means the time is the same
        for all qubits and as such is not a relevant information source for this function.
    :param t2s: Dephasing time T2. If an int, it means the time is the same
        for all qubits and as such is not a relevant information source for this function.
    :return: Number of qubits to apply the noise to, list of qubits to apply the noise to.
    :exception ValueError: If qubits is None and neither t1s nor t2s is a list (cannot infer the relevant qubits)
    """
    if qubits is None:
        n_qubits = 0
        if isinstance(t1s, list):
            n_qubits = max(n_qubits, len(t1s))
        if isinstance(t2s, list):
            n_qubits = max(n_qubits, len(t2s))
        if n_qubits == 0:
            raise ValueError('Cannot infer the list of qubits')
        qubits = list(range(n_qubits))
    elif isinstance(qubits, int):
        n_qubits = qubits
        qubits = list(range(n_qubits))
    else:
        n_qubits = len(qubits)

    return n_qubits, qubits


def _infer_or_sample_relaxation_and_dephasing_times(
        t1s: list[np.double] | np.double | None,
        t2s: list[np.double] | np.double | None,
        mean_t1: np.double,
        std_t1: np.double,
        mean_t2: np.double,
        std_t2: np.double,
        n_qubits: int
) -> tuple[list[np.double], list[np.double]]:
    """
    Infer if the relaxation and dephasing times are given or should they be sampled.

    t1s and t2s is truncated and sanitized in such a manner to ensure that T2 <= 2*T1 and 0 < T1, 0 < T2.

    :param t1s: Relaxation times T1 for each qubit. If an int, it means the same T1 is applied to all qubits.
        If none, it is sampled for each qubit from Normal(mean_t1, std_t1^2)
    :param t2s: Dephasing times T2 for each qubit. If an int, it means the same T2 is applied to all qubits.
        If none, it is sampled for each qubit from Normal(mean_t2, std_t2^2)
    :param mean_t1: Mean for the normal distribution to sample T1 from if t1s is None.
    :param std_t1: Standard deviation for the normal distribution to sample T1 from if t1s is None.
    :param mean_t2: Mean for the normal distribution to sample T2 from if t2s is None.
    :param std_t2: Standard deviation for the normal distribution to sample T2 from if t2s is None.
    :param n_qubits: The number of qubits to apply thermal noise to.
    :return: NumPy arrays of relaxation times T1 and dephasing times T2 for each qubit we want to apply noise to.
    """
    if isinstance(t1s, int):
        t1s = [t1s] * n_qubits
    if isinstance(t2s, int):
        t2s = [t2s] * n_qubits
    if t1s is None:
        t1s = np.random.normal(mean_t1, std_t1, size=n_qubits)
    if t2s is None:
        t2s = np.random.normal(mean_t2, std_t2, size=n_qubits)

    # truncate t2 to satisfy the condition T2 <= 2*T1
    t2s = np.array([min(t2s[i], 2 * t1s[i]) for i in range(n_qubits)])

    # ensure t1 > 0 and t2 > 0
    t2s = np.array([max(1, t2s[i]) for i in range(n_qubits)])
    t1s = np.array([max(1, t1s[i]) for i in range(n_qubits)])

    return t1s, t2s


def model_thermal_error_1q(
        t_reset: np.double = 1000,
        t_measure: np.double = 1000,
        t_basis_gates: np.double | list[np.double] = 100,
        t1s: list[np.double] | np.double | None = None,
        t2s: list[np.double] | np.double | None = None,
        mean_t1: np.double = 50e3,
        std_t1: np.double = 10e3,
        mean_t2: np.double = 70e3,
        std_t2: np.double = 10e3,
        basis_gates_1q: list[str] | None = None,
        qubits: list[int] | int | None = None,
        noise_model: NoiseModel = None
) -> NoiseModel:
    """
    Applies thermal error for 1-qubit gates to the noise model.

    Thermal noise depends on the environment, which may cause qubit decoherence.
    Relaxation time T1 refers to qubits contracting to the |0> state regardless of their initial state.
    (https://quantumcomputing.stackexchange.com/questions/8352/what-is-the-difference-between-relaxation-dephasing-and-decoherence)
    Dephasing time T2 refers to loss of information on the relative phase.
    Dephasing can happen when a system becomes correlated to an external environment, instead of the system of interest
    (https://en.wikipedia.org/wiki/Dephasing)

    Simplest usage: model_thermal_error_1q(qubits=7, mean_t1=2e3, std_t1=1e2, mean_t2=3e3, std_t2=1e2,basis_gates_1q=['h','s','t','x','y','z','id'])

    :param t_reset: Time it takes for a qubit to be to |0>
    :param t_measure: Time it takes for a qubit to be measured.
    :param t_basis_gates: Time it takes to perform the 1-qubit gates we want to apply the noise to.
        If an int, the same T is applied to all gates
    :param t1s: Relaxation times T1 for each qubit. If an int, it means the same T1 is applied to all qubits.
        If none, it is sampled for each qubit from Normal(mean_t1, std_t1^2)
    :param t2s: Dephasing times T2 for each qubit. If an int, it means the same T2 is applied to all qubits.
        If none, it is sampled for each qubit from Normal(mean_t2, std_t2^2)
    :param mean_t1: Mean for the normal distribution to sample T1 from if t1s is None.
    :param std_t1: Standard deviation for the normal distribution to sample T1 from if t1s is None.
    :param mean_t2: Mean for the normal distribution to sample T2 from if t2s is None.
    :param std_t2: Standard deviation for the normal distribution to sample T2 from if t2s is None.
    :param basis_gates_1q: 1-qubit gates we want to apply the noise to.
    :param qubits: The list of qubits to apply thermal noise to. Note that the order matters.
        If an int, it is treated as a number of qubits to apply the noise to (counting from the 0th qubit).
        If none, the number of qubits is attempted to be inferred from t1s and t2s
    :param noise_model: Noise model to update. If None, a new noise model is created.
    :return: an updated noise model with thermal error applied.
    :exception ValueError: If qubits is None and neither t1s nor t2s is a list (cannot infer the relevant qubits)
    """
    if basis_gates_1q is None:
        basis_gates_1q = DEFAULT_BASIS_GATES_1Q
    if noise_model is None:
        noise_model = NoiseModel()

    n_qubits, qubits = _infer_qubits(qubits, t1s, t2s)

    t1s, t2s = _infer_or_sample_relaxation_and_dephasing_times(t1s, t2s, mean_t1, std_t1, mean_t2, std_t2, n_qubits)

    errors_reset = [
        thermal_relaxation_error(t1, t2, t_reset) for t1, t2 in zip(t1s, t2s)
    ]
    errors_meas = [
        thermal_relaxation_error(t1, t2, t_measure) for t1, t2 in zip(t1s, t2s)
    ]

    if isinstance(t_basis_gates, int):
        errors_gates_1q = [
            thermal_relaxation_error(t1, t2, t_basis_gates) for t1, t2 in zip(t1s, t2s)
        ]
    else:
        errors_gates_1q = [
            [
                thermal_relaxation_error(t1, t2, t_gate) for t1, t2 in zip(t1s, t2s)
            ] for t_gate in t_basis_gates
        ]

    if isinstance(t_basis_gates, int):
        for i, q in enumerate(qubits):
            noise_model.add_quantum_error(errors_reset[i], "reset", [q])
            noise_model.add_quantum_error(errors_meas[i], "measure", [q])
            noise_model.add_quantum_error(errors_gates_1q[i], basis_gates_1q, [q])
    else:
        for i, q in enumerate(qubits):
            noise_model.add_quantum_error(errors_reset[i], "reset", [q])
            noise_model.add_quantum_error(errors_meas[i], "measure", [q])
            for j, gate in enumerate(basis_gates_1q):
                noise_model.add_quantum_error(errors_gates_1q[j], gate, [q])

    return noise_model


def model_thermal_error_2q(
        t_basis_gates: int | list[int] = 100,
        t1s: list[int] | int | None = None,
        t2s: list[int] | int | None = None,
        mean_t1: int = 50e3,
        std_t1: int = 10e3,
        mean_t2: int = 70e3,
        std_t2: int = 10e3,
        basis_gates_2q: list[str] | None = None,
        qubits: list[int] | int | None = None,
        noise_model: NoiseModel = None
) -> NoiseModel:
    """
    Applies thermal error for 2-qubit gates to the noise model.

    Thermal noise depends on the environment, which may cause qubit decoherence.
    Relaxation time T1 refers to qubits contracting to the |0> state regardless of their initial state.
    (https://quantumcomputing.stackexchange.com/questions/8352/what-is-the-difference-between-relaxation-dephasing-and-decoherence)
    Dephasing time T2 refers to loss of information on the relative phase.
    Dephasing can happen when a system becomes correlated to an external environment, instead of the system of interest
    (https://en.wikipedia.org/wiki/Dephasing)

    Simplest usage: model_thermal_error_2q(qubits=7, mean_t1=2e3, std_t1=1e2, mean_t2=3e3, std_t2=1e2,basis_gates_2q=['cx'])

    :param t_basis_gates: Time it takes to perform the 2-qubit gates we want to apply the noise to.
        If an int, the same T is applied to all gates.
    :param t1s: Relaxation times T1 for each qubit. If an int, it means the same T1 is applied to all qubits.
        If none, it is sampled for each qubit from Normal(mean_t1, std_t1^2)
    :param t2s: Dephasing times T2 for each qubit. If an int, it means the same T2 is applied to all qubits.
        If none, it is sampled for each qubit from Normal(mean_t2, std_t2^2)
    :param mean_t1: Mean for the normal distribution to sample T1 from if t1s is None.
    :param std_t1: Standard deviation for the normal distribution to sample T1 from if t1s is None.
    :param mean_t2: Mean for the normal distribution to sample T2 from if t2s is None.
    :param std_t2: Standard deviation for the normal distribution to sample T2 from if t2s is None.
    :param basis_gates_2q: 2-qubit gates we want to apply the noise to.
    :param qubits: The list of qubits to apply thermal noise to. Note that the order matters.
        If an int, it is treated as a number of qubits to apply the noise to (counting from the 0th qubit).
        If none, the number of qubits is attempted to be inferred from t1s and t2s
    :param noise_model: Noise model to update. If None, a new noise model is created.
    :return: an updated noise model with thermal error applied.
    :exception ValueError: If qubits is None and neither t1s nor t2s is a list (cannot infer the relevant qubits)
    """
    if basis_gates_2q is None:
        basis_gates_2q = DEFAULT_BASIS_GATES_2Q
    if noise_model is None:
        noise_model = NoiseModel()

    n_qubits, qubits = _infer_qubits(qubits, t1s, t2s)

    t1s, t2s = _infer_or_sample_relaxation_and_dephasing_times(t1s, t2s, mean_t1, std_t1, mean_t2, std_t2, n_qubits)

    if isinstance(t_basis_gates, int):
        errors_gates_1q = [
            thermal_relaxation_error(t1, t2, t_basis_gates) for t1, t2 in zip(t1s, t2s)
        ]
        errors_gates_2q = [
            [
                err1.tensor(err2) for err1 in errors_gates_1q
            ] for err2 in errors_gates_1q
        ]
    else:
        errors_gates_1q = [
            [
                thermal_relaxation_error(t1, t2, t_gate) for t1, t2 in zip(t1s, t2s)
            ] for t_gate in t_basis_gates
        ]
        errors_gates_2q = [
            [
                [
                    err1.tensor(err2) for err1 in errors_gates_1q[i]
                ] for err2 in errors_gates_1q[i]
            ] for i in range(len(t_basis_gates))
        ]

    if isinstance(t_basis_gates, int):
        for i, q in enumerate(qubits):
            for j, q2 in enumerate(qubits):
                noise_model.add_quantum_error(errors_gates_2q[i][j], basis_gates_2q, [q, q2])
    else:
        for i, q in enumerate(qubits):
            for j, q2 in enumerate(qubits):
                noise_model.add_quantum_error(errors_gates_2q[i][j], basis_gates_2q, [q, q2])
                for k, gate in enumerate(basis_gates_2q):
                    noise_model.add_quantum_error(errors_gates_2q[k][i][j], gate, [q, q2])

    return noise_model


def model_thermal_error_1q_simple(
        t1: int,
        t2: int,
        t_reset: int = 1000,
        t_measure: int = 1000,
        t_basis_gates: list[int] | int = 100,
        basis_gates_1q: list[str] | None = None,
        qubits: list[int] | None = None,
        noise_model: NoiseModel = None
) -> NoiseModel:
    """
    Simplified method for thermal noise.
    Applies thermal error for 1-qubit gates to the noise model.
    It applies the noise to all or specified qubits, using the same relaxation T1 and dephasing T2 times for all qubits.

    Thermal noise depends on the environment, which may cause qubit decoherence.
    Relaxation time T1 refers to qubits contracting to the |0> state regardless of their initial state.
    (https://quantumcomputing.stackexchange.com/questions/8352/what-is-the-difference-between-relaxation-dephasing-and-decoherence)
    Dephasing time T2 refers to loss of information on the relative phase.
    Dephasing can happen when a system becomes correlated to an external environment, instead of the system of interest
    (https://en.wikipedia.org/wiki/Dephasing)

    :param t_reset: Time it takes for a qubit to be to |0>
    :param t_measure: Time it takes for a qubit to be measured.
    :param t_basis_gates: Time it takes to perform the 1-qubit gates we want to apply the noise to.
        If an int, the same T is applied to all gates.
    :param t1: Relaxation time T1 for each qubit
    :param t2: Dephasing time T2 for each qubit
    :param basis_gates_1q: 1-qubit gates we want to apply the noise to.
    :param qubits: The list of qubits to apply thermal noise to. If None, it is applied to all qubits.
    :param noise_model: Noise model to update. If None, a new noise model is created.
    :return: an updated noise model with thermal error applied.
    """
    if basis_gates_1q is None:
        basis_gates_1q = DEFAULT_BASIS_GATES_1Q
    if noise_model is None:
        noise_model = NoiseModel()

    errors_reset = thermal_relaxation_error(t1, t2, t_reset)
    errors_meas = thermal_relaxation_error(t1, t2, t_measure)

    if isinstance(t_basis_gates, int):
        errors_gates_1q = thermal_relaxation_error(t1, t2, t_basis_gates)
    else:
        errors_gates_1q = [
            thermal_relaxation_error(t1, t2, t_gate) for t_gate in t_basis_gates
        ]

    _apply_errors(noise_model, errors_reset, "reset", qubits)
    _apply_errors(noise_model, errors_meas, "measure", qubits)

    if isinstance(t_basis_gates, int):
        _apply_errors(noise_model, errors_gates_1q, basis_gates_1q, qubits)
    else:
        for j, gate in enumerate(basis_gates_1q):
            _apply_errors(noise_model, errors_gates_1q[j], gate, qubits)

    return noise_model


def model_thermal_error_2q_simple(
        t1: int, t2: int,
        t_basis_gates: list[int] | int = 100,
        basis_gates_2q: list[str] | None = None,
        qubits: list[list[int]] | None = None,
        noise_model: NoiseModel = None
) -> NoiseModel:
    """
    Simplified method for thermal noise.
    Applies thermal error for 2-qubit gates to the noise model.
    It applies the noise to all or specified qubits, using the same relaxation T1 and dephasing T2 times for all qubits.

    Thermal noise depends on the environment, which may cause qubit decoherence.
    Relaxation time T1 refers to qubits contracting to the |0> state regardless of their initial state.
    (https://quantumcomputing.stackexchange.com/questions/8352/what-is-the-difference-between-relaxation-dephasing-and-decoherence)
    Dephasing time T2 refers to loss of information on the relative phase.
    Dephasing can happen when a system becomes correlated to an external environment, instead of the system of interest
    (https://en.wikipedia.org/wiki/Dephasing)

    :param t_basis_gates: Time it takes to perform the 1=2-qubit gates we want to apply the noise to.
        If an int, the same T is applied to all gates.
    :param t1: Relaxation time T1 for each qubit
    :param t2: Dephasing time T2 for each qubit
    :param basis_gates_2q: 2-qubit gates we want to apply the noise to.
    :param qubits: The list of qubits to apply thermal noise to. If None, it is applied to all qubits.
    :param noise_model: Noise model to update. If None, a new noise model is created.
    :return: an updated noise model with thermal error applied.
    """
    if basis_gates_2q is None:
        basis_gates_2q = DEFAULT_BASIS_GATES_2Q
    if noise_model is None:
        noise_model = NoiseModel()

    if isinstance(t_basis_gates, int):
        errors_gates_1q = thermal_relaxation_error(t1, t2, t_basis_gates)
        errors_gates_2q = errors_gates_1q.tensor(errors_gates_1q)
    else:
        errors_gates_1q = [
            thermal_relaxation_error(t1, t2, t_gate) for t_gate in t_basis_gates
        ]
        errors_gates_2q = [
            err.tensor(err) for err in errors_gates_1q
        ]

    if isinstance(t_basis_gates, int):
        _apply_errors(noise_model, errors_gates_2q, basis_gates_2q, qubits)
    else:
        for j, gate in enumerate(basis_gates_2q):
            _apply_errors(noise_model, errors_gates_2q[j], gate, qubits)

    return noise_model


def model_flip_error_1q(
        err_reset: float | None = 0.03,
        err_meas: float | None = 0.1,
        err_basis_gates: float | list[float] | None = 0.05,
        basis_gates_1q: list[str] | None = None,
        qubits: list[int] | None = None,
        flip_type='bit',
        noise_model: NoiseModel = None
) -> NoiseModel:
    """
    Applies bit/phase flip error to 1-qubit gates.

    Bit-flip error means introducing an X-gate after an operation with a given probability.
    Phase-flip error means introducing a Z-gate after an operation with a given probability.

    Simplest usage: model_flip_error_1q(err_reset=1e-3, err_meas=1e-3, err_basis_gates=0.01, basis_gates_1q=['h', 's', 't', 'x', 'y', 'z', 'id'])

    :param err_reset: The error rate when resetting a qubit to |0>. Ignored if None
    :param err_meas: The error rate when measuring a qubit. Ignored if None
    :param err_basis_gates: The error rate for given 1-qubit gates we want to apply the noise to.
        If a float, the same rate is applied to all gates. Ignored if None.
    :param basis_gates_1q: The 1-qubit gates we want to apply the noise to.
    :param qubits: The list of qubits we want to apply the noise to. If None, noise is applied to all qubits.
    :param flip_type: The type of flipping to apply (bit/phase). Must be 'bit', 'phase', 'X' or 'Z'
    :param noise_model: Noise model to update. If None, a new noise model is created.
    :return: an updated noise model with flip error applied.
    """
    if flip_type not in ['bit', 'phase', 'X', 'Z']:
        raise ValueError('flip_type must be either "bit" or "phase"')
    if basis_gates_1q is None:
        basis_gates_1q = DEFAULT_BASIS_GATES_1Q
    if noise_model is None:
        noise_model = NoiseModel()

    flip_gate = "X" if flip_type == 'bit' else "Z" if flip_type == 'phase' else flip_type

    if err_reset is not None and err_reset > 0.0:
        error_reset = pauli_error([(flip_gate, err_reset), ("I", 1 - err_reset)])
        _apply_errors(noise_model, error_reset, "reset", qubits)
    if err_meas is not None and err_meas > 0.0:
        error_meas = pauli_error([(flip_gate, err_meas), ("I", 1 - err_meas)])
        _apply_errors(noise_model, error_meas, "measure", qubits)

    if err_basis_gates is not None:
        if isinstance(err_basis_gates, float) and err_basis_gates > 0.0:
            error_gates_1q = pauli_error([(flip_gate, err_basis_gates), ("I", 1 - err_basis_gates)])
            _apply_errors(noise_model, error_gates_1q, basis_gates_1q, qubits)
        elif isinstance(err_basis_gates, list):
            for err, gate in zip(err_basis_gates, basis_gates_1q):
                if err is None or err == 0.0:
                    continue
                error_gate = pauli_error([(flip_gate, err), ("I", 1 - err)])
                _apply_errors(noise_model, error_gate, basis_gates_1q, qubits)

    return noise_model


def model_flip_error_2q(
        err_basis_gates: float | list[float] | None = 0.05,
        basis_gates_2q: list[str] | None = None,
        qubits: list[list[int]] | None = None,
        flip_type='bit',
        noise_model: NoiseModel = None
) -> NoiseModel:
    """
        Applies bit/phase flip error to 1-qubit gates.

        Bit-flip error means introducing an X-gate after an operation with a given probability.
        Phase-flip error means introducing a Z-gate after an operation with a given probability.

        Simplest usage: model_flip_error_2q(err_basis_gates=0.01, basis_gates_2q=['cx'])

        :param err_basis_gates: The error rate for given 1-qubit gates we want to apply the noise to.
            If a float, the same rate is applied to all gates. Ignored if None.
        :param basis_gates_2q: The 2-qubit gates we want to apply the noise to.
        :param qubits: The list of qubits we want to apply the noise to. If None, noise is applied to all qubits.
        :param flip_type: The type of flipping to apply (bit/phase). Must be 'bit', 'phase', 'X' or 'Z'
        :param noise_model: Noise model to update. If None, a new noise model is created.
        :return: an updated noise model with flip error applied.
        """
    if flip_type not in ['bit', 'phase', 'X', 'Z']:
        raise ValueError('flip_type must be either "bit" or "phase"')
    if basis_gates_2q is None:
        basis_gates_2q = DEFAULT_BASIS_GATES_2Q
    if noise_model is None:
        noise_model = NoiseModel()
    if qubits is None:
        qubits = [None]

    flip_gate = "X" if flip_type == 'bit' else "Z" if flip_type == 'phase' else flip_type

    if err_basis_gates is not None:
        if isinstance(err_basis_gates, float) and err_basis_gates > 0.0:
            error_gates_1q = pauli_error([(flip_gate, err_basis_gates), ("I", 1 - err_basis_gates)])
            error_gates_2q = error_gates_1q.tensor(error_gates_1q)
            for q_pair in qubits:
                _apply_errors(noise_model, error_gates_2q, basis_gates_2q, q_pair)
        elif isinstance(err_basis_gates, list):
            for err, gate in zip(err_basis_gates, basis_gates_2q):
                if err is None or err == 0.0:
                    continue
                error_gate = pauli_error([(flip_gate, err), ("I", 1 - err)])
                error_gate_2q = error_gate.tensor(error_gate)
                for q_pair in qubits:
                    _apply_errors(noise_model, error_gate_2q, basis_gates_2q, q_pair)

    return noise_model


def build_combined_noise_model(
        depol_err: float | None,
        t1_mean: float,
        n_qubits: int
) -> NoiseModel:
    """Build noise model with combined depolarizing and thermal errors."""
    noise_model = NoiseModel()

    t2_mean = t1_mean * 1.5
    t_gate_1q = 100
    t_gate_2q = 300

    for q in range(n_qubits):
        # 1-qubit gates: compose depolarizing + thermal
        thermal_1q = thermal_relaxation_error(t1_mean, t2_mean, t_gate_1q)
        if depol_err is not None and depol_err > 0:
            depol_1q = depolarizing_error(depol_err, 1)
            combined_1q = depol_1q.compose(thermal_1q)
        else:
            combined_1q = thermal_1q

        noise_model.add_quantum_error(combined_1q, DEFAULT_BASIS_GATES_1Q, [q])

    # 2-qubit gates
    for q1 in range(n_qubits):
        for q2 in range(n_qubits):
            if q1 == q2:
                continue
            thermal_2q = thermal_relaxation_error(t1_mean, t2_mean, t_gate_2q).tensor(
                thermal_relaxation_error(t1_mean, t2_mean, t_gate_2q)
            )
            if depol_err is not None and depol_err > 0:
                depol_2q = depolarizing_error(depol_err, 1).tensor(depolarizing_error(depol_err, 1))
                combined_2q = depol_2q.compose(thermal_2q)
            else:
                combined_2q = thermal_2q

            noise_model.add_quantum_error(combined_2q, DEFAULT_BASIS_GATES_2Q, [q1, q2])

    return noise_model


def get_noise_from_ibm_backend(backend: str) -> NoiseModel:
    """
    Given a name of an IBM backend (e.g. "ibm_fez"), it outputs a noise model
    according to the latest noise data from the selected backend.
    :param backend: The name of the IBM backend.
    :return: Noise model.
    """
    service = QiskitRuntimeService()
    backend = service.backend(backend)
    noise_model = NoiseModel.from_backend(backend)

    return noise_model


def _construct_grover_solver(numbers, target, backend):
    oracle = SubsetSumOracle(numbers, target)

    def verifier(subset):
        return sum(w * s for w, s in zip(numbers, subset)) == target

    solver = GroverSolver(
        oracle=oracle,
        n_search_qubits=len(numbers),
        verifier=verifier,
        backend=backend,
    )

    return solver


def run_simple_sweep(numbers: list[int], target: int, shots: int = 1024):
    # errs = [None, 0.01, 0.02, 0.03, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25]
    errs = [None, 0.001, 0.002, 0.003, 0.005, 0.0075, 0.01, 0.015, 0.02, 0.025]

    for err in errs:
        noise_model = model_depolarizing_error(err=err)
        # noise_model = model_thermal_error_1q(t_measure=0.01, qubits=7, mean_t1=2e3, std_t1=1e2, mean_t2=3e3, std_t2=1e2,basis_gates_1q=['h','s','t','x','y','z','id'], noise_model=noise_model)
        # noise_model = model_thermal_error_2q(qubits=7, mean_t1=2e3, std_t1=1e2, mean_t2=3e3, std_t2=1e2,basis_gates_2q=['cx'], noise_model=noise_model)
        # noise_model = model_flip_error_1q(err_reset=1e-3, err_meas=1e-3, err_basis_gates=0.01, basis_gates_1q=['h', 's', 't', 'x', 'y', 'z', 'id'], noise_model=noise_model)

        print(noise_model)

        simulator = AerSimulator(noise_model=noise_model)

        solver = _construct_grover_solver(numbers, target, simulator)

        result = solver.run_unknown_solutions(shots=shots)
        # result = solver.run_known_solutions(n_solutions=2, shots=shots)
        if result is None:
            print("No verified solution found")
        else:
            print(f'result: {result}')
            correct_solutions_bitstrings = [bitstring for _, bitstring in
                                            solve_subset_sum_with_bitstrings(numbers, target)]
            print(f"Is correct: {check_counts_statistically_accurate(result, correct_solutions_bitstrings)}")
            plot_results(result, numbers)


if __name__ == '__main__':
    numbers = [1, 2, 3, 4]
    target = 5

    run_simple_sweep(numbers, target, shots=1024)