🎨 re-implemented AE without qiskit application modules

src/mqt/bench/benchmarks/ae.py

@@ -4,66 +4,76 @@
 
 import numpy as np
 from qiskit import QuantumCircuit
-from qiskit_algorithms import AmplitudeEstimation, EstimationProblem
 
 
-def create_circuit(num_qubits: int) -> QuantumCircuit:
+def create_circuit(num_qubits: int, probability: float = 0.2) -> QuantumCircuit:
     """Returns a quantum circuit implementing Quantum Amplitude Estimation.
 
     Arguments:
-        num_qubits: number of qubits of the returned quantum circuit
-    """
-    ae = AmplitudeEstimation(
-        num_eval_qubits=num_qubits - 1,  # -1 because of the to be estimated qubit
-    )
-    problem = get_estimation_problem()
+        num_qubits: Total number of qubits, including evaluation and target qubits.
+        probability: Probability of the "good" state.
 
-    qc = ae.construct_circuit(problem)
-    qc.name = "ae"
-    qc.measure_all()
+    Returns:
+        QuantumCircuit: The constructed amplitude estimation circuit.
+    """
+    if num_qubits < 2:
+        msg = "Number of qubits must be at least 2 (1 evaluation + 1 target)."
+        raise ValueError(msg)
 
-    return qc
+    num_eval_qubits = num_qubits - 1  # Number of evaluation qubits
+    qc = QuantumCircuit(num_qubits, num_eval_qubits)
 
+    # Define the Bernoulli A operator
+    theta_p = 2 * np.arcsin(np.sqrt(probability))
+    a = QuantumCircuit(1)
+    a.ry(theta_p, 0)
 
-class BernoulliQ(QuantumCircuit):  # type: ignore[misc]
-    """A circuit representing the Bernoulli Q operator."""
+    # Define the Grover operator Q = -A S_0 A† S_f
+    def grover_operator() -> QuantumCircuit:
+        """Construct the Grover operator."""
+        q = QuantumCircuit(1)
+        # Apply A
+        q.ry(theta_p, 0)
+        # Apply S0 (reflection around |0>)
+        q.z(0)
+        # Apply A†
+        q.ry(-theta_p, 0)
+        return q
 
-    def __init__(self, probability: float) -> None:
-        """Initialize the Bernoulli Q operator."""
-        super().__init__(1)  # circuit on 1 qubit
+    q = grover_operator()
 
-        self._theta_p = 2 * np.arcsin(np.sqrt(probability))
-        self.ry(2 * self._theta_p, 0)
+    # Apply Hadamard gates to the evaluation qubits
+    qc.h(range(num_eval_qubits))
 
-    def __eq__(self, other: object) -> bool:
-        """Return if the operators are equal."""
-        return isinstance(other, BernoulliQ) and self._theta_p == other._theta_p
+    # Controlled applications of A and powers of Q
+    for i in range(num_eval_qubits):
+        qc.append(a.control(1), [i, num_eval_qubits])  # Controlled A
+        qc.append(q.control(1), [i, num_eval_qubits])  # Controlled powers of Q
 
-    def __hash__(self) -> int:
-        """Return a hash of the operator."""
-        return hash(self._theta_p)
+    # Apply the inverse QFT to the evaluation qubits
+    qc.append(inverse_qft(num_eval_qubits), range(num_eval_qubits))
 
-    def power(self, power: float, _matrix_power: bool = True) -> QuantumCircuit:
-        """Return a circuit implementing the power of the operator."""
-        q_k = QuantumCircuit(1)
-        q_k.ry(2 * power * self._theta_p, 0)
-        return q_k
+    # Measure the evaluation qubits
+    qc.measure(range(num_eval_qubits), range(num_eval_qubits))
 
+    qc.name = "ae"
+    return qc
 
-def get_estimation_problem() -> EstimationProblem:
-    """Returns a estimation problem instance for a fixed p value."""
-    p = 0.2
 
-    """A circuit representing the Bernoulli A operator."""
-    a = QuantumCircuit(1)
-    theta_p = 2 * np.arcsin(np.sqrt(p))
-    a.ry(theta_p, 0)
+def inverse_qft(num_qubits: int) -> QuantumCircuit:
+    """Constructs the inverse Quantum Fourier Transform circuit.
 
-    """A circuit representing the Bernoulli Q operator."""
-    q = BernoulliQ(p)
+    Arguments:
+        num_qubits: Number of qubits.
 
-    return EstimationProblem(
-        state_preparation=a,  # A operator
-        grover_operator=q,  # Q operator
-        objective_qubits=[0],  # the "good" state Psi1 is identified as measuring |1> in qubit 0
-    )
+    Returns:
+        QuantumCircuit: The inverse QFT circuit.
+    """
+    qc = QuantumCircuit(num_qubits)
+    for i in range(num_qubits // 2):
+        qc.swap(i, num_qubits - i - 1)
+    for i in range(num_qubits):
+        qc.h(i)
+        for j in range(i + 1, num_qubits):
+            qc.cp(-np.pi / (2 ** (j - i)), j, i)
+    return qc