From b63e4c0f61d3e903ddf654ef57e1882f1ebeaa42 Mon Sep 17 00:00:00 2001
From: davidedellagiustina <davide@dellagiustina.com>
Date: Sun, 4 Feb 2024 23:00:31 +0000
Subject: [PATCH] Working on amplitude estimation simulation

---
 src/qasp/problems/estimation.py | 69 ++++++++++++++++++++++++++++++++-
 1 file changed, 68 insertions(+), 1 deletion(-)

diff --git a/src/qasp/problems/estimation.py b/src/qasp/problems/estimation.py
index 382bb8f..4ebe584 100644
--- a/src/qasp/problems/estimation.py
+++ b/src/qasp/problems/estimation.py
@@ -5,7 +5,8 @@
 import math
 from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
 from qiskit.circuit.library import GroverOperator, QFT
-from ..oracle import QuantumOracle
+from ..oracle import Oracle, QuantumOracle
+from ..simul import __exec_circuit
 
 
 # +--------------------+
@@ -103,3 +104,69 @@ def circuit(
     circ.measure(list(range(t-m, t)), result)  # Only m bits
 
     return circ
+
+
+# +-----------------------+
+# | Algoroithm simulation |
+# +-----------------------+
+
+def __measure_to_count(measurements: str, num_search_qubits: int) -> tuple[float, int]:
+    '''Convert the result of a measurement to to actual phase factor and the resulting solutions \
+        count.
+
+    #### Arguments
+        measurements (str): Measured bits.
+        num_search_qubits (int): Number of search qubits.
+
+    #### Return
+        tuple[float, int]: Phase factor and solutions count.
+    '''
+    phi = 0.0
+    for (idx, bit) in zip(range(len(measurements)), measurements):
+        phi += int(bit) * 2**(-idx-1)
+    print(measurements, phi)
+    phase = 2 * math.pi * phi
+    count = 2**num_search_qubits * (math.sin(phase / 2))**2
+
+    return (phase, count)
+
+
+def exec_count(
+    algorithm: QuantumCircuit,
+    oracle: Oracle,
+    m: int,
+    eps: float,
+    aux_qubits: list[int] = None
+) -> tuple[QuantumCircuit, float, int]:
+    '''Simulate the amplitude estimation circuit to approximate the number of solutions of the \
+        problem.
+
+    #### Arguments
+        algorithm (QuantumCircuit): Ciruit that implements the initialization algorithm.
+        oracle (Oracle): Pair of classical and quantum oracles.
+        m (int): Desired number of binary digits to be estimated.
+        eps (float): Complement of the desired success probability.
+        aux_qubits (list[int]): List of indices of auxiliary qubits (e.g. used by the oracle) \
+            that should not be used for the search procedure. Defaults to the empty list.
+
+    #### Return
+        tuple[QuantumCircuit, float, int]: Used circuit, measured phase, and estimated number \
+            of solutions.
+    '''
+    aux_qubits = [] if aux_qubits is None else aux_qubits
+    (_, q_oracle) = oracle  # Classical oracle is unused
+
+    n = len(q_oracle.qubits) - len(aux_qubits) + 1  # Account for `aug`
+
+    # Build circuit
+    circ = circuit(algorithm, q_oracle, m, eps, aux_qubits)
+
+    # Run simulation
+    result = __exec_circuit(circ, shots=10000)
+    measurements = list(result.get_counts(circ).keys())[0]
+    print(result.get_counts())
+
+    # Compute results
+    (phase, count) = __measure_to_count(measurements, n)
+
+    return (circ, phase, count)