feat(algorithms): implement Canonical Quantum Amplitude Estimation (QAE) using QPE

README.md

@@ -69,7 +69,7 @@ Python Program"]
 - **Bernstein-Vazirani**: Determining hidden linear functions
 - **Deutsch-Jozsa**: Distinguishing constant from balanced functions
 - **Quantum Random Number Generator (QRNG)**: True quantum randomness
-- **Iterative Amplitude Estimation**: Maximum-likelihood amplitude estimation for Monte Carlo-style workflows (canonical QPE-based variant planned)
+- **Amplitude Estimation**: Canonical QPE-based amplitude estimation, plus the iterative (maximum-likelihood) variant for Monte Carlo-style workflows
 
 
 ### Quantum Information & States

qpiai_quantum/algorithms/README.md

@@ -165,9 +165,9 @@ print(f"10 random values: {values}")
 ```
 
 ### 9. Amplitude Estimation
-Quantum speedup for Monte Carlo estimation.  The SDK provides the **iterative
-(maximum-likelihood) variant**; the canonical QPE-based variant is planned but
-not yet implemented.
+Quantum speedup for Monte Carlo estimation.  The SDK provides both the
+**canonical QPE-based variant** (`AmplitudeEstimation`) and the **iterative
+(maximum-likelihood) variant** (`IterativeAmplitudeEstimation`).
 
 ```python
 from qpiai_quantum.algorithms.amplitude_estimation import (
@@ -193,10 +193,45 @@ amplitude = iae.estimate(problem, shots=2000)
 print(f"Estimated amplitude: {amplitude:.4f}")
 ```
 
-> **Note:** The canonical QPE-based `AmplitudeEstimation` class is not yet
-> available and raises `NotImplementedError`.  Use
-> `IterativeAmplitudeEstimation` for currently supported amplitude estimation
-> workflows.
+The canonical QPE-based variant estimates the amplitude from the phase
+measured on an evaluation register:
+
+```python
+from qpiai_quantum.algorithms.amplitude_estimation import (
+    AmplitudeEstimation,
+    EstimationProblem,
+)
+from qpiai_quantum.circuit import Circuit
+
+A = Circuit(1)
+A.ry(0, 0.8)
+problem = EstimationProblem(state_preparation=A, objective_qubits=[0])
+
+# More evaluation qubits gives finer phase resolution; 5 or more is recommended.
+qae = AmplitudeEstimation(num_evaluation_qubits=6)
+amplitude = qae.estimate(problem, shots=1000, device_name="QpiAI-QSV-Local")
+print(f"Estimated amplitude: {amplitude:.4f}")  # true value: sin^2(0.4) = 0.1516
+```
+
+`AmplitudeEstimation` honours a custom `is_good_state` predicate: it is
+enumerated over the state register and synthesized into an explicit oracle.
+That enumeration is supported for up to `AmplitudeEstimation.MAX_ORACLE_QUBITS`
+state qubits; problems using the default `objective_qubits` marking have no
+such limit.  `IterativeAmplitudeEstimation` does not yet build its Grover
+operator from the predicate, so custom predicates are not supported there:
+
+```python
+# Bitstrings are MSB first, so bitstring[-1] is qubit 0.
+B = Circuit(2)
+B.ry(0, 0.8)
+B.ry(1, 1.1)
+
+problem = EstimationProblem(
+    state_preparation=B,
+    objective_qubits=[0],
+    is_good_state=lambda bitstring: bitstring.count("1") == 1,  # exactly one qubit set
+)
+```
 
 ## Common Operations
 
@@ -359,7 +394,7 @@ for n in [2, 3, 4, 5]:
 - **Complexity**: O(1/ε) queries vs O(1/ε²) classical Monte Carlo
 - **Use Cases**: Option pricing, risk analysis, counting problems
 - **Key Feature**: Quadratic speedup for estimating expectation values
-- **Note**: Only the iterative (maximum-likelihood) variant is currently implemented. The canonical QPE-based variant is planned.
+- **Note**: Both the canonical QPE-based variant (`AmplitudeEstimation`) and the iterative (maximum-likelihood) variant (`IterativeAmplitudeEstimation`) are implemented. Custom `is_good_state` predicates are currently supported by the canonical variant only.
 
 ## Performance Tips
 

qpiai_quantum/algorithms/amplitude_estimation.py

@@ -52,6 +52,7 @@ def __init__(
         self.state_preparation = state_preparation
         self.objective_qubits = objective_qubits
         self.is_good_state: Callable[[str], bool]
+        self._default_is_good_state: Callable[[str], bool] | None
         if is_good_state is None:
 
             def default_is_good_state(bitstring: str) -> bool:
@@ -65,43 +66,255 @@ def default_is_good_state(bitstring: str) -> bool:
                 return True
 
             self.is_good_state = default_is_good_state
+            self._default_is_good_state = default_is_good_state
         else:
             self.is_good_state = is_good_state
+            self._default_is_good_state = None
 
     @property
     def num_qubits(self) -> int:
         return self.state_preparation.num_qubits
 
+    @property
+    def has_custom_good_state(self) -> bool:
+        """
+        True when is_good_state is anything other than the default marking.
 
-class AmplitudeEstimation(QuantumAlgorithm):
-    """
-    Canonical Quantum Amplitude Estimation (QAE) — **not yet implemented**.
+        Compared by identity rather than recorded at construction time, so
+        reassigning is_good_state afterwards is still honoured by the circuit
+        rather than only by the result counting.
+        """
+        return self.is_good_state is not self._default_is_good_state
 
-    The canonical QPE-based variant is not yet available.  Canonical QAE
-    requires constructing controlled Grover operators, which in turn need
-    a general controlled-subcircuit compiler not yet present in the SDK.
 
-    Use :class:`IterativeAmplitudeEstimation` for currently supported
-    amplitude estimation workflows.
+class AmplitudeEstimation(QuantumAlgorithm):
+    """
+    Canonical Quantum Amplitude Estimation (QAE) using Quantum Phase Estimation (QPE).
+
+    QAE uses an evaluation register of m qubits initialized in superposition,
+    applies controlled powers of the Grover operator C-Q^(2^(m-1-j)) on
+    evaluation qubit j, and then applies Inverse QFT on the evaluation register
+    to estimate the amplitude a = sin^2(theta_a).
+
+    Custom ``is_good_state`` predicates are supported: the predicate is
+    enumerated over the state register and synthesized into an explicit oracle
+    when it differs from the default all-objective-qubits-are-1 marking.  That
+    enumeration is bounded by ``MAX_ORACLE_QUBITS``; problems relying on the
+    default marking are unaffected by the bound.
     """
 
+    # Upper bound on the state register for which an explicit oracle can be
+    # synthesized from a custom is_good_state predicate (2^n enumeration).
+    # Only reached when a caller supplies their own predicate.
+    MAX_ORACLE_QUBITS = 16
+
     def __init__(self, num_evaluation_qubits: int):
         self.num_evaluation_qubits = num_evaluation_qubits
         super().__init__(num_qubits=0, name="Amplitude Estimation")
         self.description = "Canonical Amplitude Estimation using QPE"
 
-    def build_circuit(self, problem: EstimationProblem) -> Circuit:
-        raise NotImplementedError(
-            "Canonical QAE requires native implementation of controlled circuits."
-        )
+    @staticmethod
+    def _apply_multi_controlled_z(circuit: Circuit, controls: list[int], target: int):
+        """Apply a phase flip on the all-ones state of controls + target."""
+        if not controls:
+            circuit.z(target)
+        elif len(controls) == 1:
+            circuit.cz(controls[0], target)
+        elif len(controls) == 2:
+            circuit.h(target)
+            circuit.ccx(controls[0], controls[1], target)
+            circuit.h(target)
+        else:
+            circuit.h(target)
+            circuit.mcx(controls, target)
+            circuit.h(target)
+
+    def _resolve_good_states(self, problem: EstimationProblem) -> list[int] | None:
+        """
+        Resolve ``problem.is_good_state`` into the basis states it marks.
+
+        Returns None when the predicate marks exactly the states covered by the
+        default "every objective qubit is 1" rule, so the compact
+        objective-qubit encoding can be used.  Otherwise returns the integer
+        encoding of every marked basis state, where bit q is qubit q.
+
+        Only a caller-supplied predicate has to be enumerated, so the state
+        register size is bounded on that path alone; the default marking is
+        emitted directly on the objective qubits at any register size.
+        """
+        if not problem.has_custom_good_state:
+            return None
+
+        n = problem.num_qubits
+        if n > self.MAX_ORACLE_QUBITS:
+            raise ValueError(
+                f"A custom is_good_state predicate can only be synthesized into "
+                f"an oracle for at most {self.MAX_ORACLE_QUBITS} state qubits, "
+                f"got {n}. Express the marked set through objective_qubits "
+                f"instead, which has no such limit."
+            )
+
+        objective_mask = 0
+        for q in problem.objective_qubits:
+            objective_mask |= 1 << q
+
+        good: list[int] = []
+        default: list[int] = []
+        for value in range(2**n):
+            # Same layout as measurement counts: MSB first, qubit q at index n-1-q.
+            if problem.is_good_state(format(value, f"0{n}b")):
+                good.append(value)
+            if value & objective_mask == objective_mask:
+                default.append(value)
+
+        return None if good == default else good
+
+    def _apply_controlled_s_chi(
+        self,
+        circuit: Circuit,
+        control_qubit: int,
+        state_offset: int,
+        problem: EstimationProblem,
+        good_states: list[int] | None = None,
+    ):
+        """
+        Apply controlled phase-flip S_chi on good states.
+
+        With ``good_states`` None the default "all objective qubits are 1"
+        marking is emitted directly on the objective qubits.  Otherwise one
+        phase flip is emitted per marked basis state, so custom
+        ``is_good_state`` predicates are honoured rather than silently ignored.
+        """
+        if good_states is None:
+            obj_qubits = [state_offset + q for q in problem.objective_qubits]
+            self._apply_multi_controlled_z(
+                circuit, [control_qubit] + obj_qubits[:-1], obj_qubits[-1]
+            )
+            return
+
+        n_state = problem.num_qubits
+        state_qubits = [state_offset + q for q in range(n_state)]
+        controls = [control_qubit] + state_qubits[:-1]
+        target = state_qubits[-1]
+
+        for value in good_states:
+            zeros = [state_qubits[q] for q in range(n_state) if not (value >> q) & 1]
+            for q in zeros:
+                circuit.x(q)
+            self._apply_multi_controlled_z(circuit, controls, target)
+            for q in zeros:
+                circuit.x(q)
+
+    def _apply_controlled_s_0(
+        self,
+        circuit: Circuit,
+        control_qubit: int,
+        state_offset: int,
+        n_state_qubits: int,
+    ):
+        """Apply controlled phase-flip S_0 on zero state |0...0>."""
+        state_qubits = [state_offset + q for q in range(n_state_qubits)]
 
-    def estimate(self, problem: EstimationProblem) -> float:
-        raise NotImplementedError(
-            "Canonical QAE requires native implementation of controlled sub-circuits, "
-            "which are not yet fully supported by the basic Circuit builder. "
-            "Please use IterativeAmplitudeEstimation instead."
+        for q in state_qubits:
+            circuit.x(q)
+
+        self._apply_multi_controlled_z(
+            circuit, [control_qubit] + state_qubits[:-1], state_qubits[-1]
         )
 
+        for q in state_qubits:
+            circuit.x(q)
+
+    def _apply_controlled_grover(
+        self,
+        circuit: Circuit,
+        control_qubit: int,
+        state_offset: int,
+        problem: EstimationProblem,
+        power: int,
+        good_states: list[int] | None = None,
+    ):
+        r"""
+        Apply controlled Grover operator Q^power controlled by control_qubit.
+
+        Controlled-Q = A · (Controlled-S_0) · A^\dagger · (Controlled-S_chi)
+        """
+        n_state = problem.num_qubits
+        state_qubits = [state_offset + q for q in range(n_state)]
+        A = problem.state_preparation
+        A_inv = A.inverse()
+
+        for _ in range(power):
+            circuit.z(control_qubit)
+            self._apply_controlled_s_chi(
+                circuit, control_qubit, state_offset, problem, good_states
+            )
+            circuit.compose(A_inv, state_qubits)
+            self._apply_controlled_s_0(circuit, control_qubit, state_offset, n_state)
+            circuit.compose(A, state_qubits)
+
+    def build_circuit(self, problem: EstimationProblem) -> Circuit:
+        m = self.num_evaluation_qubits
+        n = problem.num_qubits
+        total_qubits = m + n
+
+        self.num_qubits = total_qubits
+        self.circuit = Circuit(total_qubits, m)
+
+        # 1. Initialize evaluation register to |+>
+        for i in range(m):
+            self.circuit.h(i)
+
+        # 2. Initialize state register with A
+        state_qubits = [m + q for q in range(n)]
+        self.circuit.compose(problem.state_preparation, state_qubits)
+
+        # 3. Apply controlled Grover powers C-Q^(2^(m-1-j)).
+        #    Evaluation qubit j must accumulate the phase weight the inverse QFT
+        #    expects at that position, matching QuantumPhaseEstimation.
+        good_states = self._resolve_good_states(problem)
+        for j in range(m):
+            power = 2 ** (m - 1 - j)
+            self._apply_controlled_grover(
+                self.circuit, j, m, problem, power, good_states
+            )
+
+        # 4. Apply Inverse QFT on evaluation register
+        QFT.apply_inverse_qft_to_circuit(self.circuit, start=0, n=m)
+
+        # 5. Measure evaluation register
+        for i in range(m):
+            self.circuit.measure(i, i)
+
+        return self.circuit
+
+    def estimate(
+        self,
+        problem: EstimationProblem,
+        shots: int = 1024,
+        device_name: str = "QpiAI-QSV-Local",
+    ) -> float:
+        """
+        Estimate amplitude a = sin^2(theta_a) using Canonical QAE.
+        """
+        m = self.num_evaluation_qubits
+        circ = self.build_circuit(problem)
+        result = circ.run(shots=shots, device_name=device_name)
+        counts = result.get()["counts"]
+
+        if not counts:
+            return 0.0
+
+        total_shots = sum(counts.values())
+        weighted_amplitude = 0.0
+
+        for bitstring, count in counts.items():
+            y = int(bitstring, 2)
+            theta = math.pi * y / (2**m)
+            weighted_amplitude += (count / total_shots) * (math.sin(theta) ** 2)
+
+        return float(weighted_amplitude)
+
 
 class IterativeAmplitudeEstimation(QuantumAlgorithm):
     """