Add quantum phase estimation in the list of algorithms

src/qutip_qip/algorithms/__init__.py

@@ -1 +1,2 @@
 from .qft import *
+from .qpe import *

src/qutip_qip/algorithms/qpe.py

@@ -0,0 +1,124 @@
+import numpy as np
+from qutip import Qobj
+from qutip_qip.operations import Gate, ControlledGate
+from qutip_qip.circuit import QubitCircuit
+from qutip_qip.algorithms import qft_gate_sequence
+
+__all__ = ["qpe"]
+
+
+class CustomGate(Gate):
+    """
+    Custom gate that wraps an arbitrary quantum operator.
+    """
+
+    def __init__(self, targets, U, **kwargs):
+        super().__init__(targets=targets, **kwargs)
+        self.targets = targets if isinstance(targets, list) else [targets]
+        self.U = U
+        self.kwargs = kwargs
+        self.latex_str = r"U"
+
+    def get_compact_qobj(self):
+        return self.U
+
+
+def create_controlled_unitary(controls, targets, U, control_value=1):
+    """
+    Create a controlled unitary gate.
+
+    Parameters
+    ----------
+    controls : list
+        Control qubits
+    targets : list
+        Target qubits
+    U : Qobj
+        Unitary operator to apply on target qubits
+    control_value : int, optional
+        Control value (default: 1)
+
+    Returns
+    -------
+    ControlledGate
+        The controlled unitary gate
+    """
+    return ControlledGate(
+        controls=controls,
+        targets=targets,
+        control_value=control_value,
+        target_gate=CustomGate,
+        U=U,
+    )
+
+
+def qpe(U, num_counting_qubits, target_qubits=None, to_cnot=False):
+    """
+    Quantum Phase Estimation circuit implementation for QuTiP.
+
+    Parameters
+    ----------
+    U : Qobj
+        Unitary operator whose eigenvalue we want to estimate.
+        Should be a unitary quantum operator.
+    num_counting_qubits : int
+        Number of counting qubits to use for the phase estimation.
+        More qubits provide higher precision.
+    target_qubits : int or list, optional
+        Index or indices of the target qubit(s) where the eigenstate is prepared.
+        If None, target_qubits is set automatically based on U's dimension.
+    to_cnot : bool, optional
+        Flag to decompose controlled phase gates to CNOT gates (default: False)
+
+    Returns
+    -------
+    qc : instance of QubitCircuit
+        Gate sequence implementing Quantum Phase Estimation.
+    """
+    if num_counting_qubits < 1:
+        raise ValueError("Minimum value of counting qubits must be 1")
+
+    # Handle target qubits specification
+    if target_qubits is None:
+        dim = U.dims[0][0]
+        num_target_qubits = int(np.log2(dim))
+        if 2**num_target_qubits != dim:
+            raise ValueError(
+                f"Unitary operator dimension {dim} is not a power of 2"
+            )
+        target_qubits = list(
+            range(num_counting_qubits, num_counting_qubits + num_target_qubits)
+        )
+    elif isinstance(target_qubits, int):
+        target_qubits = [target_qubits]
+        num_target_qubits = 1
+    else:
+        num_target_qubits = len(target_qubits)
+
+    total_qubits = num_counting_qubits + num_target_qubits
+    qc = QubitCircuit(total_qubits)
+
+    # Apply Hadamard gates to all counting qubits
+    for i in range(num_counting_qubits):
+        qc.add_gate("SNOT", targets=[i])
+
+    # Apply controlled-U gates with increasing powers
+    for i in range(num_counting_qubits):
+        power = 2 ** (num_counting_qubits - i - 1)
+        # Create U^power
+        U_power = U if power == 1 else U**power
+
+        # Add controlled-U^power gate
+        controlled_u = create_controlled_unitary(
+            controls=[i], targets=target_qubits, U=U_power
+        )
+        qc.add_gate(controlled_u)
+
+    # Add inverse QFT on counting qubits
+    inverse_qft_circuit = qft_gate_sequence(
+        num_counting_qubits, swapping=True, to_cnot=to_cnot
+    ).reverse_circuit()
+    for gate in inverse_qft_circuit.gates:
+        qc.add_gate(gate)
+
+    return qc

src/qutip_qip/operations/gateclass.py

@@ -73,6 +73,7 @@
     "R",
     "QASMU",
     "SWAP",
+    "ControlledGate",
     "ISWAP",
     "CNOT",
     "SQRTSWAP",

tests/test_qpe.py

@@ -0,0 +1,145 @@
+import numpy as np
+from numpy.testing import assert_, assert_equal, assert_allclose
+import unittest
+from qutip import Qobj, sigmax, sigmay, sigmaz, identity, basis, tensor
+from qutip_qip.circuit import QubitCircuit
+from qutip_qip.operations import gate_sequence_product, ControlledGate
+
+from qutip_qip.algorithms.qpe import qpe, CustomGate, create_controlled_unitary
+
+
+class TestQPE(unittest.TestCase):
+    """
+    A test class for the Quantum Phase Estimation implementation
+    """
+
+    def test_custom_gate(self):
+        """
+        Test if CustomGate correctly stores and returns the quantum object
+        """
+        U = Qobj([[0, 1], [1, 0]])
+
+        custom = CustomGate(targets=[0], U=U)
+
+        qobj = custom.get_compact_qobj()
+        assert_((qobj - U).norm() < 1e-12)
+
+        assert_equal(custom.targets, [0])
+
+    def test_controlled_unitary(self):
+        """
+        Test if create_controlled_unitary creates a correct controlled gate
+        """
+        U = Qobj([[0, 1], [1, 0]])
+
+        controlled_u = create_controlled_unitary(
+            controls=[0], targets=[1], U=U
+        )
+
+        assert_equal(controlled_u.controls, [0])
+        assert_equal(controlled_u.targets, [1])
+        assert_equal(controlled_u.control_value, 1)
+
+        assert_(controlled_u.target_gate == CustomGate)
+
+        assert_("U" in controlled_u.kwargs)
+        assert_((controlled_u.kwargs["U"] - U).norm() < 1e-12)
+
+    def test_qpe_validation(self):
+        """
+        Test input validation in qpe function
+        """
+        U = sigmaz()
+
+        with self.assertRaises(ValueError):
+            qpe(U, num_counting_qubits=0)
+
+        invalid_U = Qobj([[1, 0, 0], [0, 1, 0], [0, 0, 1]])  # 3x3 matrix
+        with self.assertRaises(ValueError):
+            qpe(invalid_U, num_counting_qubits=3)
+
+    def test_qpe_circuit_structure(self):
+        """
+        Test if qpe creates circuit with correct structure
+        """
+        U = sigmaz()
+
+        num_counting = 3
+        circuit = qpe(
+            U, num_counting_qubits=num_counting, target_qubits=num_counting
+        )
+
+        assert_equal(circuit.N, num_counting + 1)
+
+        for i in range(num_counting):
+            assert_equal(circuit.gates[i].targets, [i])
+
+        for i in range(num_counting):
+            gate = circuit.gates[num_counting + i]
+            assert_(isinstance(gate, ControlledGate))
+            assert_equal(gate.controls, [i])
+            assert_equal(gate.targets, [num_counting])
+
+    def test_qpe_different_target_specifications(self):
+        """
+        Test QPE with different ways of specifying target qubits
+        """
+        U = sigmaz()
+        num_counting = 2
+
+        circuit1 = qpe(
+            U, num_counting_qubits=num_counting, target_qubits=num_counting
+        )
+        assert_equal(circuit1.N, num_counting + 1)
+
+        circuit2 = qpe(
+            U, num_counting_qubits=num_counting, target_qubits=[num_counting]
+        )
+        assert_equal(circuit2.N, num_counting + 1)
+
+        circuit3 = qpe(U, num_counting_qubits=num_counting, target_qubits=None)
+        assert_equal(circuit3.N, num_counting + 1)
+
+        U2 = tensor(sigmaz(), sigmaz())
+        circuit4 = qpe(
+            U2,
+            num_counting_qubits=num_counting,
+            target_qubits=[num_counting, num_counting + 1],
+        )
+        assert_equal(circuit4.N, num_counting + 2)
+
+    def test_qpe_controlled_gate_powers(self):
+        """
+        Test if QPE correctly applies powers of U
+        """
+        phase = 0.25
+        U = Qobj([[1, 0], [0, np.exp(1j * np.pi * phase)]])
+
+        num_counting = 3
+        circuit = qpe(
+            U, num_counting_qubits=num_counting, target_qubits=num_counting
+        )
+
+        for i in range(num_counting):
+            gate = circuit.gates[num_counting + i]
+            power = 2 ** (num_counting - i - 1)
+
+            u_power = gate.kwargs["U"]
+            expected_u_power = U if power == 1 else U**power
+
+            assert_((u_power - expected_u_power).norm() < 1e-12)
+
+    def test_qpe_to_cnot_flag(self):
+        """
+        Test if to_cnot flag works correctly in QPE
+        """
+        U = sigmaz()
+        num_counting = 2
+
+        circuit1 = qpe(U, num_counting_qubits=num_counting, to_cnot=False)
+
+        circuit2 = qpe(U, num_counting_qubits=num_counting, to_cnot=True)
+
+        has_cnot = any(gate.name == "CNOT" for gate in circuit2.gates)
+        assert_(has_cnot)
+        assert_(len(circuit2.gates) > len(circuit1.gates))