This page was generated from the notebook examples/notebooks/4c_Quantum_Phase_Estimation.ipynb.

from math import floor

import numpy as np
import numpy.typing as npt

from mpqp import (
    QCircuit,
    H,
    CRk,
    Barrier,
    SWAP,
    Gate,
    CustomGate,
    Ry,
    X,
    BasisMeasure,
    IBMDevice,
    run,
)
from mpqp.qasm.qasm_to_mpqp import qasm2_parse
from mpqp.tools.maths import normalize

Generate circuit from state

from qiskit import QuantumCircuit
from qiskit import qasm2
import re


def generate_circuit(state: npt.NDArray[np.complex128], nb_qubits: int):
    circuit = QuantumCircuit(nb_qubits)
    circuit.initialize(state, list(range(nb_qubits)))
    decomposed_circuit = (
        circuit.decompose().decompose().decompose().decompose().decompose()
    ) # weird trick, we know
    qasm_code = qasm2.dumps(decomposed_circuit.reverse_bits())
    cleaned_code = re.sub(r"reset q\[\d+\];\n", "", qasm_code)
    circ = qasm2_parse(cleaned_code)

    return circ

print(generate_circuit(normalize(np.array([2,5,4,6])), 2))

       ┌───────────────┐
q_0: ──┤ U(1.8587,0,0) ├─────■────────────────────
     ┌─┴───────────────┴──┐┌─┴─┐┌────────────────┐
q_1: ┤ U(1.7783,-2π,-π/2) ├┤ X ├┤ U(0.60229,0,0) ├
     └────────────────────┘└───┘└────────────────┘

Quantum Fourier Transform

def QFT(nb_qubits: int, inverse: bool=False):

    qft = QCircuit(nb_qubits)
    for j in range(nb_qubits):
        qft.add(H(j))
        qft.add([CRk(i+1, i, j) for i in range(j+1, nb_qubits)])
        qft.add(Barrier())
    qft.add([SWAP(i, nb_qubits-1-i) for i in range(nb_qubits // 2)])

    if inverse:
        return qft.inverse()
    return qft

print(QFT(3))
print(QFT(3, inverse=True))

     ┌───┐                   ░                ░       ░
q_0: ┤ H ├─■────────■────────░────────────────░───────░──X─
     └───┘ │P(π/2)  │        ░ ┌───┐          ░       ░  │
q_1: ──────■────────┼────────░─┤ H ├─■────────░───────░──┼─
                    │P(π/4)  ░ └───┘ │P(π/4)  ░ ┌───┐ ░  │
q_2: ───────────────■────────░───────■────────░─┤ H ├─░──X─
                             ░                ░ └───┘ ░
         ░       ░                 ░                     ┌───┐
q_0: ─X──░───────░─────────────────░──■─────────■────────┤ H ├
      │  ░       ░           ┌───┐ ░  │         │P(-π/2) └───┘
q_1: ─┼──░───────░──■────────┤ H ├─░──┼─────────■─────────────
      │  ░ ┌───┐ ░  │P(-π/4) └───┘ ░  │P(-π/4)
q_2: ─X──░─┤ H ├─░──■──────────────░──■───────────────────────
         ░ └───┘ ░                 ░

Build successive control U^k gates

def c_Uk(unitary: Gate, n_1: int, n_2: int) -> CustomGate:
    matrix = np.zeros((2**(n_1 + n_2), 2**(n_1 + n_2)))

    n_1_power = 2**n_1
    for k in range(n_1_power):
        k_vec = np.zeros((n_1_power, n_1_power))
        k_vec[k,k] = 1
        matrix += np.kron(k_vec, unitary.power(k).to_matrix())
    return CustomGate(matrix, list(range(n_1+n_2)))

c_Uk(Ry(0.34,0),2,1)

CustomGate(UnitaryMatrix(array([[ 1. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 1. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0.98558477, -0.16918235, 0. , 0. , 0. , 0. ], [ 0. , 0. , 0.16918235, 0.98558477, 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0.94275467, -0.33348709, 0. , 0. ], [ 0. , 0. , 0. , 0. , 0.33348709, 0.94275467, 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0.87274451, -0.48817725], [ 0. , 0. , 0. , 0. , 0. , 0. , 0.48817725, 0.87274451]])), [0, 1, 2] )

Building the QPE circuit

def QPE(
    eigen_state: npt.NDArray[np.complex128], U: Gate, n_1: int, n_2: int, shots: int
):

    first_register = QCircuit(H.range(n_1))
    state_circuit = generate_circuit(eigen_state, n_2)
    whool_circuit = first_register @ state_circuit
    whool_circuit.add(c_Uk(U, n_1, n_2))
    whool_circuit.append(QFT(n_1, inverse=True))

    # Extract result
    whool_circuit.add(BasisMeasure(list(range(n_1)), shots=shots))
    result = run(whool_circuit, IBMDevice.AER_SIMULATOR)
    probabilities = list(result.probabilities)
    print(probabilities)
    most_probable = probabilities.index(max(probabilities))
    print(most_probable)

    # From basis state to phase
    max_count_state = format(most_probable, '0{}b'.format(n_1))[
        ::-1
    ]  # Reverse the bit order
    print(max_count_state)
    phase_estimation = int(max_count_state, 2) / (2**n_1)

    return phase_estimation

state = np.array([0, 1])
U = CustomGate(np.diag([1, -1]), [0])
print(np.linalg.eig(U.to_matrix()))
n_1 = 3
n_2 = 1
phase = QPE(state, U, n_1, n_2, 10000)
print("phase :", phase)
print("eigenvalue", np.exp(1j * 2 * np.pi * phase))

EigResult(eigenvalues=array([ 1., -1.]), eigenvectors=array([[1., 0.],
       [0., 1.]]))
[0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
4
001
phase : 0.125
eigenvalue (0.7071067811865476+0.7071067811865476j)

EigResult(eigenvalues=array([ 1., -1.]), eigenvectors=array([[1., 0.],
       [0., 1.]]))

X(0).to_matrix()

array([[0., 1.],
       [1., 0.]])