diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..767a925 --- /dev/null +++ b/.gitignore @@ -0,0 +1,132 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +pip-wheel-metadata/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +.python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + + +.DS_STORE \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..4d6bf4c --- /dev/null +++ b/README.md @@ -0,0 +1,4 @@ +# IC-POVM +Informationally complete POVM + +This is a first try diff --git a/base_povm.py b/base_povm.py new file mode 100644 index 0000000..36a9feb --- /dev/null +++ b/base_povm.py @@ -0,0 +1,137 @@ +import numpy as np + +from qiskit.quantum_info import Operator, DensityMatrix + + +class Povm: + """Abstract base class that collects all information that any POVM should specifiy.""" + + def __init__(self, povm_ops: np.ndarray): + """Initialize from explicit POVM operators. + + Args: + povm_operators: np.ndarray that contains the explicit list of POVM operators""" + + if not (len(povm_ops.shape) == 3 and povm_ops.shape[1] == povm_ops.shape[1]): + raise ValueError( + f"POVM operators need to be square instead of {povm_ops.shape[1:]}" + ) + + self.n_outcomes = povm_ops.shape[0] + self.dimension = povm_ops.shape[1] + self.povm_operators = [Operator(op) for op in povm_ops] + self.array_ops = None + + self.dual_operators = None + self.frame_superop = None + self.informationlly_complete = None + + def check_validity(self) -> bool: + """Checks if POVM axioms are fulfilled.""" + + summed_op = np.zeros((self.dimension, self.dimension), dtype=complex) + + for k in range(len(self)): + + if not np.allclose(self.povm_operators[k].data.conj().T - self.povm_operators[k].data, 0.0, atol=1e-5): + raise ValueError(f"POVM operator {k} is not hermitian.") + + for eigval in np.linalg.eigvalsh(self.povm_operators[k].data): + if eigval.real < -1e-6 or np.abs(eigval.imag) > 1e-5: + raise ValueError( + f"Negative eigenvalue {eigval} in POVM operator {k}." + ) + + summed_op += self.povm_operators[k].data + + if not np.allclose(summed_op - np.identity(self.dimension, dtype=complex), 0.0, atol=1e-5): + raise ValueError(f"POVM operators not summing up to the identity : \n{summed_op}") + + return True + + def __getitem__(self, index: slice) -> np.ndarray: + """Return a povm operator or a list of povm operators.""" + return self.povm_operators[index] + + def __len__(self): + return len(self.povm_operators) + + def get_prob(self, rho: DensityMatrix) -> np.ndarray: + return np.array([ + np.real(np.trace(rho.data @ povm_op.data)) + for povm_op in self.povm_operators]) + + @classmethod + def from_vectors(cls, povm_vectors: np.ndarray): + """Initialize a POVM from the bloch vectors |psi> (not normalized!) such that Pi = |psi> np.ndarray: + # """Return a numpy array of shape (n_outcomes, d, d) that includes all povm operators.""" + # if isinstance(index, int) : + # return self.povm_operators[index].data + # elif isinstance(index, slice) : + # return np.array([op.data for op in self.povm_operators[index]]) + # else: + # raise TypeError("Invalid Argument Type") + + def get_ops(self, idx:slice=...) -> np.ndarray: + """Return a numpy array of shape (n_outcomes, d, d) that includes all povm operators.""" + + if self.array_ops is None: + self.array_ops = np.zeros((self.n_outcomes, self.dimension, self.dimension), dtype=complex) + for k, op in enumerate(self.povm_operators): + self.array_ops[k] = op.data + + return self.array_ops[idx] + +''' diff --git a/povm_implementation.py b/povm_implementation.py new file mode 100644 index 0000000..8b54514 --- /dev/null +++ b/povm_implementation.py @@ -0,0 +1,96 @@ +import numpy as np +from collections import Counter +from qiskit.circuit import QuantumCircuit, ParameterVector +from base_povm import Povm +from single_qubit_povm import SingleQubitPOVM +from product_povm import ProductPOVM + + +class POVMImplementation: + + def __init__( + self, + n_qubit: int, + ) -> None: + + self.n_qubit = n_qubit + self.parametrized_qc = self._build_qc() + + def _build_qc(self) -> QuantumCircuit: + raise NotImplementedError("The subclass of POVMImplementation must implement `_build_from_param` method.") + + def get_parameter_and_shot(self, shot: int) -> QuantumCircuit: + raise NotImplementedError("The subclass of POVMImplementation must implement `get_parameter_and_shot` method.") + + def to_povm(self) -> Povm: + raise NotImplementedError("The subclass of POVMImplementation must implement `to_povm` method.") + + +class ProductPVMSimPOVMImplementation(POVMImplementation): + + def __init__( + self, + n_qubit: int, + parameters: np.ndarray | None = None, + ) -> None: + + super().__init__(n_qubit) + self._set_parameters(parameters) + + def _set_parameters(self, parameters: np.ndarray) -> None: + # n_param = n_qubit*(3*self.n_PVM-1) + if len(parameters) % self.n_qubit != 0: + raise ValueError('The length of the parameter array is expected to be multiple of the number of qubits') + elif (len(parameters) / self.n_qubit + 1) % 3 != 0: + raise ValueError('The number of parameters per qubit is expected to be of the form 3*n_PVM-1') + else: + self.n_PVM = int((len(parameters) // self.n_qubit + 1) // 3) + parameters = parameters.reshape((self.n_qubit, self.n_PVM * 3 - 1)) + self.angles = parameters[:, :2 * self.n_PVM].reshape((self.n_qubit, self.n_PVM, 2)) + self.PVM_distributions = np.concatenate((parameters[:, 2 * self.n_PVM:], np.ones((self.n_qubit, 1))), axis=1) + if np.any(self.PVM_distributions < 0.): + raise ValueError('There should not be any negative values in the probability distribution parameters.') + else: + self.PVM_distributions = self.PVM_distributions / self.PVM_distributions.sum(axis=1)[:, np.newaxis] + + def _build_qc(self) -> QuantumCircuit: + + theta = ParameterVector('theta', length=self.n_qubit) + phi = ParameterVector('phi', length=self.n_qubit) + + qc = QuantumCircuit(self.n_qubit) + for i in range(self.n_qubit): + qc.u(theta=theta[i], phi=phi[i], lam=0, qubit=i) + + return qc + + def get_parameter_and_shot(self, shot: int) -> QuantumCircuit: + """ + Returns a list with concrete parameter values and associated number of shots. + """ + + PVM_idx = np.zeros((shot, self.n_qubit), dtype=int) + + for i in range(self.n_qubit): + PVM_idx[:, i] = np.random.choice(self.n_PVM, size=int(shot), replace=True, p=self.PVM_distributions[i]) + counts = Counter(tuple(x) for x in PVM_idx) + + return [tuple(([self.angles[i, combination[i]] for i in range(self.n_qubit)], counts[combination])) for combination in counts] + + def to_povm(self) -> Povm: + + stabilizers = np.zeros((self.n_qubit, self.n_PVM, 2, 2), dtype=complex) + + stabilizers[:, :, 0, 0] = np.cos(self.angles[:, :, 0] / 2.) + stabilizers[:, :, 0, 1] = (np.cos(self.angles[:, :, 1]) + 1.j * np.sin(self.angles[:, :, 1])) * np.sin(self.angles[:, :, 0] / 2.) + stabilizers[:, :, 1, 0] = stabilizers[:, :, 0, 1].conjugate() + stabilizers[:, :, 1, 1] = -stabilizers[:, :, 0, 0] + + stabilizers = np.multiply(stabilizers.T, np.sqrt(self.PVM_distributions).T).T + stabilizers = stabilizers.reshape((self.n_qubit, 2 * self.n_PVM, 2)) + + sq_povms = [] + for vecs in stabilizers: + sq_povms.append(SingleQubitPOVM.from_vectors(vecs)) + + return ProductPOVM(sq_povms) diff --git a/product_povm.py b/product_povm.py new file mode 100644 index 0000000..83c6e00 --- /dev/null +++ b/product_povm.py @@ -0,0 +1,30 @@ +from typing import List + +# from qiskit.quantum_info import DensityMatrix, SparsePauliOp + +from base_povm import Povm + + +class ProductPOVM(Povm): + """Class to represent a set of product POVM operators.""" + + def __init__(self, povm_list: List[Povm]): + """Initialize from explicit list of POVMs.""" + self.dimension = 1 + self.n_outcomes = 1 + self.n_operators = 0 + for povm in povm_list: + self.dimension *= povm.dimension + self.n_outcomes *= povm.n_outcomes + self.n_operators += povm.n_outcomes + + self.povm_list = povm_list + +# def get_prob(self, rho: DensityMatrix): +# return get_p_from_paulis(SparsePauliOp.from_operator(rho), self.list_povm).ravel() + +# def plot_bloch_sphere(self, dual=False, colors=None): +# list_fig = [] +# for sqpovm in self.list_povm: +# list_fig.append(sqpovm.plot_bloch_sphere(dual, colors)) +# return list_fig diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..dbe64a4 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +python +numpy +qiskit>=1.0 \ No newline at end of file diff --git a/single_qubit_povm.py b/single_qubit_povm.py new file mode 100644 index 0000000..bddd679 --- /dev/null +++ b/single_qubit_povm.py @@ -0,0 +1,47 @@ +import numpy as np + +from qiskit.quantum_info import SparsePauliOp + +from base_povm import Povm + + +class SingleQubitPOVM(Povm): + """Class to represent a set of IC single-qubit POVM operators.""" + + def __init__(self, povm_ops: np.ndarray): + """Initialize from explicit POVM operators.""" + + super().__init__(povm_ops) + + # self.check_validity() + + self.povm_pauli_ops = [SparsePauliOp.from_operator(op) for op in self.povm_operators] + self.povm_pauli_decomp = [ + SingleQubitPOVM.pauli_op_to_array(pauli_op) for pauli_op in self.povm_pauli_ops + ] + + def check_validity(self): + if not self.dimension == 2: + raise ValueError(f"Dimension of Single Qubit POVM operator space should be 2, not {self.dimension}.") + return Povm.check_validity(self) + + @staticmethod + def pauli_op_to_array(pauli_op: SparsePauliOp): + """Convert a single-qubit SparsePauliOp into an array [c0, c1, c2, c3], + where the indices represent paulis as {"I": 0, "X": 1, "Y": 2, "Z": 3}.""" + + labels = np.zeros(4) + 0j + + for pauli_idx, coeff in pauli_op.label_iter(): + if pauli_idx == "I": + labels[0] = coeff + elif pauli_idx == "X": + labels[1] = coeff + elif pauli_idx == "Y": + labels[2] = coeff + elif pauli_idx == "Z": + labels[3] = coeff + else: + raise ValueError(f"Unexpected pauli string {pauli_idx}") + + return np.real_if_close(labels) diff --git a/test.ipynb b/test.ipynb new file mode 100644 index 0000000..dca02e7 --- /dev/null +++ b/test.ipynb @@ -0,0 +1,1401 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.append('test/')\n", + "\n", + "import numpy as np\n", + "from scipy.stats import unitary_group\n", + "from base_povm import Povm\n", + "from single_qubit_povm import SingleQubitPOVM\n", + "from povm_implementation import ProductPVMSimPOVMImplementation\n", + "from test_base_povm import TestBasePovm\n", + "\n", + "from qiskit.quantum_info import random_density_matrix, Operator" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2\n" + ] + }, + { + "data": { + "text/plain": [ + "[array([0.25, 0. , 0. , 0.25]),\n", + " array([ 0.25, 0. , 0. , -0.25]),\n", + " array([0.15, 0.15, 0. , 0. ]),\n", + " array([ 0.15, -0.15, 0. , 0. ]),\n", + " array([0.1, 0. , 0.1, 0. ]),\n", + " array([ 0.1, 0. , -0.1, 0. ])]" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ops = np.arange(24).reshape((6,2,2))\n", + "povm1 = Povm(povm_ops=ops)\n", + "print(povm1.dimension)\n", + "\n", + "pauli_X = np.array([[0,1],[1,0]])\n", + "\n", + "basis_0 = np.array([1.,0], dtype=complex)\n", + "basis_1 = np.array([0,1.], dtype=complex)\n", + "basis_plus = 1.0/np.sqrt(2) * (basis_0 + basis_1)\n", + "basis_minus = 1.0/np.sqrt(2) * (basis_0 - basis_1)\n", + "basis_plus_i = 1.0/np.sqrt(2) * (basis_0 + 1.j * basis_1)\n", + "basis_minus_i = 1.0/np.sqrt(2) * (basis_0 - 1.j * basis_1)\n", + "\n", + "q = [0.5, 0.3, 0.2]\n", + "\n", + "Z0 = np.outer(basis_0, basis_0.conj())\n", + "Z1 = np.outer(basis_1, basis_1.conj())\n", + "X0 = np.outer(basis_plus, basis_plus.conj())\n", + "X1 = np.outer(basis_minus, basis_minus.conj())\n", + "Y0 = np.outer(basis_plus_i, basis_plus_i.conj())\n", + "Y1 = np.outer(basis_minus_i, basis_minus_i.conj())\n", + "\n", + "sqpovm = SingleQubitPOVM(np.array([q[0]* Z0,q[0]* Z1,q[1]* X0,q[1]* X1,q[2]* Y0,q[2]* Y1]))\n", + "\n", + "sqpovm.povm_pauli_decomp\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Operator([[0.5+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n" + ] + } + ], + "source": [ + "sqpovm[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([5.01299530e-06+1.11641635e-24j, 3.02104770e-05+8.59363292e-25j,\n", + " 4.19594048e-05+4.14111853e-26j, 1.24877387e-02-2.07842297e-21j,\n", + " 9.89677783e-01-9.92563858e-21j, 9.97757296e-01-1.98538787e-18j])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dim = 2 \n", + "n_out = 6\n", + "n_param = n_out*(2*dim-1) - dim**2\n", + "\n", + "np.random.seed(4)\n", + "param = np.random.uniform(low=0, high=1, size=n_param)\n", + "povm = Povm.from_param(param, dim)\n", + "#povm = sqpovm\n", + "\n", + "np.unique(np.array([np.trace(op.data) for op in povm.povm_operators]))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DensityMatrix([[ 0.46571833-8.48656244e-18j, -0.02617055-4.60556496e-01j],\n", + " [-0.02617055+4.60556496e-01j, 0.53428167+8.48656244e-18j]],\n", + " dims=(2,))\n" + ] + } + ], + "source": [ + "rho = random_density_matrix(dims=2, seed = 4)\n", + "rho\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([4.68781168e-01, 5.24874638e-01, 6.31079159e-03, 1.17603086e-05,\n", + " 1.96265836e-06, 1.96791650e-05])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p = povm.get_prob(rho)\n", + "p" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([4.68781168e-01, 5.24874638e-01, 6.31079159e-03, 1.17603086e-05,\n", + " 1.96265836e-06, 1.96791650e-05])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.real(np.trace(np.dot(povm.get_ops(),rho.data),axis1=1,axis2=2))" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "84 µs ± 356 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "povm = Povm.from_param(param, dim)\n", + "p = np.real(np.trace(np.dot(povm.get_ops(),rho.data).T))" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "84 µs ± 184 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "povm = Povm.from_param(param, dim)\n", + "p = np.real(np.trace(np.dot(povm.get_ops(),rho.data),axis1=1,axis2=2))" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "94.1 µs ± 158 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "povm = Povm.from_param(param, dim)\n", + "p = povm.get_prob(rho)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "337 ns ± 1.07 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "\n", + "for k in range(povm.n_outcomes) :\n", + " y=povm.povm_operators[k].data" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "515 ns ± 1.79 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "\n", + "for op in povm :\n", + " y=op.data" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "627 ns ± 0.665 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "\n", + "for k in range(povm.n_outcomes) :\n", + " y=povm.povm_operators[k].data" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "833 ns ± 1.09 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "\n", + "for op in povm :\n", + " y=op" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "989 ns ± 4.8 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "\n", + "for op in povm :\n", + " y=op" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "620 ns ± 2.12 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "\n", + "for k in range(povm.n_outcomes) :\n", + " y=povm.povm_operators[k].data" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "ename": "TypeError", + "evalue": "unsupported operand type(s) for *: 'float' and 'SingleQubitPOVM'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[6], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m (\u001b[38;5;241;43m0.3\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43msqpovm\u001b[49m \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m0.7\u001b[39m\u001b[38;5;241m*\u001b[39msqpovm)\u001b[38;5;241m.\u001b[39mcheck_validity()\n", + "\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for *: 'float' and 'SingleQubitPOVM'" + ] + } + ], + "source": [ + "(0.3*povm + 0.7*sqpovm).check_validity()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True\n" + ] + }, + { + "data": { + "text/plain": [ + "array([ 1.00000000e+00+0.j, -1.11022302e-16+0.j, 0.00000000e+00+0.j,\n", + " 0.00000000e+00+0.j])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dim = 2\n", + "eival = np.random.uniform(low=-5, high=5, size=dim)\n", + "x = unitary_group.rvs(dim)#, random_state=seed_obs[i])\n", + "obs = x @ np.diag(eival) @ x.T.conj()\n", + "\n", + "_, V = np.linalg.eigh(obs)\n", + "\n", + "sqpovm = SingleQubitPOVM.from_vectors()\n", + "print(sqpovm.check_validity())\n", + "\n", + "sum = np.zeros(4, dtype=complex)\n", + "for coef in sqpovm.povm_pauli_decomp:\n", + " sum += coef\n", + "\n", + "sum" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3\n", + "2\n" + ] + } + ], + "source": [ + "ops = np.arange(54).reshape((6,3,3))\n", + "povm2 = Povm(povm_ops=ops)\n", + "print(povm2.dimension)\n", + "print(povm1.dimension)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3\n" + ] + }, + { + "data": { + "text/plain": [ + "3" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "povm2.print(n=3)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Point(lat=1, long=2)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from dataclasses import dataclass\n", + "\n", + "@dataclass\n", + "\n", + "class Point:\n", + "\n", + " lat: float\n", + "\n", + " long: float\n", + "\n", + "\n", + "p = Point(lat=1,long=2)\n", + "p\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ops = np.arange(24).reshape((6,2,2))\n", + "povm1 = SingleQubitPOVM(povm_ops=ops)\n", + "print(povm1.dimension)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "Dimension of Single Qubit POVM operator space should be 2, not 3.", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[6], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m ops \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m54\u001b[39m)\u001b[38;5;241m.\u001b[39mreshape((\u001b[38;5;241m6\u001b[39m,\u001b[38;5;241m3\u001b[39m,\u001b[38;5;241m3\u001b[39m))\n\u001b[0;32m----> 2\u001b[0m povm2 \u001b[38;5;241m=\u001b[39m \u001b[43mSingleQubitPOVM\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpovm_ops\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mops\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 3\u001b[0m \u001b[38;5;28mprint\u001b[39m(povm2\u001b[38;5;241m.\u001b[39mdimension)\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28mprint\u001b[39m(povm1\u001b[38;5;241m.\u001b[39mdimension)\n", + "File \u001b[0;32m~/Documents/Internship/IC-POVM/single_qubit_povm.py:18\u001b[0m, in \u001b[0;36mSingleQubitPOVM.__init__\u001b[0;34m(self, povm_ops)\u001b[0m\n\u001b[1;32m 14\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Initialize from explicit POVM operators.\"\"\"\u001b[39;00m\n\u001b[1;32m 16\u001b[0m \u001b[38;5;28msuper\u001b[39m()\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__init__\u001b[39m(povm_ops)\n\u001b[0;32m---> 18\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcheck_validity\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 20\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpovm_pauli_ops \u001b[38;5;241m=\u001b[39m [SparsePauliOp\u001b[38;5;241m.\u001b[39mfrom_operator(op) \u001b[38;5;28;01mfor\u001b[39;00m op \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpovm_operators]\n\u001b[1;32m 21\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpovm_pauli_decomp \u001b[38;5;241m=\u001b[39m [\n\u001b[1;32m 22\u001b[0m SingleQubitPOVM\u001b[38;5;241m.\u001b[39mpauli_op_to_array(pauli_op) \u001b[38;5;28;01mfor\u001b[39;00m pauli_op \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpovm_pauli_ops\n\u001b[1;32m 23\u001b[0m ]\n", + "File \u001b[0;32m~/Documents/Internship/IC-POVM/single_qubit_povm.py:27\u001b[0m, in \u001b[0;36mSingleQubitPOVM.check_validity\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mcheck_validity\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m 26\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdimension\u001b[38;5;241m==\u001b[39m\u001b[38;5;241m2\u001b[39m:\n\u001b[0;32m---> 27\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mDimension of Single Qubit POVM operator space should be 2, not \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdimension\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 28\u001b[0m Povm\u001b[38;5;241m.\u001b[39mcheck_validity(\u001b[38;5;28mself\u001b[39m)\n", + "\u001b[0;31mValueError\u001b[0m: Dimension of Single Qubit POVM operator space should be 2, not 3." + ] + } + ], + "source": [ + "ops = np.arange(54).reshape((6,3,3))\n", + "povm2 = SingleQubitPOVM(povm_ops=ops)\n", + "print(povm2.dimension)\n", + "print(povm1.dimension)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[[ 0. 1.]\n", + " [ 2. 3.]]\n", + "\n", + " [[ 4. 5.]\n", + " [ 6. 7.]]\n", + "\n", + " [[ 8. 9.]\n", + " [10. 11.]]]\n", + "\n", + "\n", + " [[[12. 13.]\n", + " [14. 15.]]\n", + "\n", + " [[16. 17.]\n", + " [18. 19.]]\n", + "\n", + " [[20. 21.]\n", + " [22. 23.]]]\n", + "\n", + "\n", + " [[[24. 25.]\n", + " [26. 27.]]\n", + "\n", + " [[28. 29.]\n", + " [30. 31.]]\n", + "\n", + " [[32. 33.]\n", + " [34. 35.]]]]\n" + ] + } + ], + "source": [ + "a = np.arange(36.0).reshape((3,3,2,2))\n", + "print(a)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[0. 1. 2.]\n", + " [3. 4. 5.]\n", + " [6. 7. 8.]]\n" + ] + } + ], + "source": [ + "\n", + "b = np.arange(9.0).reshape((3,3))\n", + "print(b)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[[ 0., 0.],\n", + " [ 0., 0.],\n", + " [ 4., 5.],\n", + " [ 6., 7.],\n", + " [ 16., 18.],\n", + " [ 20., 22.]],\n", + "\n", + " [[ 36., 39.],\n", + " [ 42., 45.],\n", + " [ 64., 68.],\n", + " [ 72., 76.],\n", + " [100., 105.],\n", + " [110., 115.]],\n", + "\n", + " [[144., 150.],\n", + " [156., 162.],\n", + " [196., 203.],\n", + " [210., 217.],\n", + " [256., 264.],\n", + " [272., 280.]]])" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(np.multiply(a.T,b.T).T).reshape((3,6,2))" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
     ┌──────────────────────┐\n",
+       "q_0: ┤ U(theta[0],phi[0],0) ├\n",
+       "     ├──────────────────────┤\n",
+       "q_1: ┤ U(theta[1],phi[1],0) ├\n",
+       "     └──────────────────────┘
" + ], + "text/plain": [ + " ┌──────────────────────┐\n", + "q_0: ┤ U(theta[0],phi[0],0) ├\n", + " ├──────────────────────┤\n", + "q_1: ┤ U(theta[1],phi[1],0) ├\n", + " └──────────────────────┘" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import numpy as np\n", + "from povm_implementation import ProductPVMSimPOVMImplementation\n", + "povm_implementation = ProductPVMSimPOVMImplementation(2, np.array([0,0, 1,0.5, 1,1, 1,1, 0,0, 1,0.5, 0.3,0.2, 2,1]))\n", + "povm_implementation._build_qc().draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[Operator([[0.5+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 0. +0.j, -0. +0.j],\n", + " [-0. +0.j, 0.5+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.19253779+0.00000000e+00j, 0.09230753-5.04278350e-02j],\n", + " [0.09230753+5.04278350e-02j, 0.05746221-9.13716138e-19j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 0.05746221+9.13716138e-19j, -0.09230753+5.04278350e-02j],\n", + " [-0.09230753-5.04278350e-02j, 0.19253779+0.00000000e+00j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.24441706+0.00000000e+00j, 0.03620368-7.33885021e-03j],\n", + " [0.03620368+7.33885021e-03j, 0.00558294-9.66574794e-20j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 0.00558294+9.66574794e-20j, -0.03620368+7.33885021e-03j],\n", + " [-0.03620368-7.33885021e-03j, 0.24441706+0.00000000e+00j]],\n", + " input_dims=(2,), output_dims=(2,))]" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "povm = povm_implementation.to_povm()\n", + "povm.povm_list[1][:]" + ] + }, + { + "cell_type": "code", + "execution_count": 87, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[([array([1.57079633, 1.57079633])], 18),\n", + " ([array([1.57079633, 0. ])], 35),\n", + " ([array([0., 0.])], 47)]" + ] + }, + "execution_count": 87, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "n_qubit = 1\n", + "parameters = np.array(n_qubit * [0., 0., 0.5 * np.pi, 0., 0.5 * np.pi, 0.5 * np.pi, 2.5, 1.5])\n", + "cs_implementation = ProductPVMSimPOVMImplementation(n_qubit=n_qubit, parameters=parameters)\n", + "cs_implementation.get_parameter_and_shot(shot=100)" + ] + }, + { + "cell_type": "code", + "execution_count": 88, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[Operator([[0.5+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 0. +0.j, -0. +0.j],\n", + " [-0. +0.j, 0.5+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.15+0.j, 0.15+0.j],\n", + " [0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 0.15+0.j, -0.15+0.j],\n", + " [-0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[1.000000e-01+0.00000000e+00j, 6.123234e-18-1.00000000e-01j],\n", + " [6.123234e-18+1.00000000e-01j, 1.000000e-01-1.23299735e-34j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 1.000000e-01+1.23299735e-34j, -6.123234e-18+1.00000000e-01j],\n", + " [-6.123234e-18-1.00000000e-01j, 1.000000e-01+0.00000000e+00j]],\n", + " input_dims=(2,), output_dims=(2,))]" + ] + }, + "execution_count": 88, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cs_povm = cs_implementation.to_povm()\n", + "cs_povm.povm_list[0][:]" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[Operator([[0.5+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0.5+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.15+0.j, 0.15+0.j],\n", + " [0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 0.15+0.j, -0.15-0.j],\n", + " [-0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.1+0.j , 0. -0.1j],\n", + " [0. +0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.1+0.j , 0. +0.1j],\n", + " [0. -0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,))]" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sqpovm[:]" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True\n", + "True\n", + "True\n", + "True\n", + "True\n", + "True\n" + ] + } + ], + "source": [ + "n_qubit = 1\n", + "q = np.random.uniform(0,5,size=3 * n_qubit).reshape((n_qubit, 3))\n", + "q /= q.sum(axis=1)[:, np.newaxis]\n", + "\n", + "parameters = np.array(n_qubit * [0., 0., 0.5 * np.pi, 0., 0.5 * np.pi, 0.5 * np.pi, 1, 1])\n", + "for i in range(n_qubit):\n", + " parameters[i * 8 + 6] = q[i,0]/q[i,2]\n", + " parameters[i * 8 + 7] = q[i,1]/q[i,2]\n", + "\n", + "cs_implementation = ProductPVMSimPOVMImplementation(n_qubit=n_qubit, parameters=parameters)\n", + "cs_povm = cs_implementation.to_povm()\n", + "for i in range(n_qubit):\n", + " sqpovm = SingleQubitPOVM(np.array([q[i,0]* Z0,q[i,0]* Z1,q[i,1]* X0,q[i,1]* X1,q[i,2]* Y0,q[i,2]* Y1]))\n", + " for k in range(sqpovm.n_outcomes):\n", + " print(np.allclose(cs_povm.povm_list[i][k].data, sqpovm[k].data))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 3\n", + "False\n", + "0 5\n", + "False\n", + "0 4\n", + "False\n", + "0 1\n", + "False\n", + "0 0\n", + "True\n", + "[3 5 4 1 2]\n", + "1 3\n", + "False\n", + "1 5\n", + "False\n", + "1 4\n", + "False\n", + "1 1\n", + "True\n", + "[3 5 4 2]\n", + "2 3\n", + "False\n", + "2 5\n", + "False\n", + "2 4\n", + "False\n", + "2 2\n", + "True\n", + "[3 5 4]\n", + "3 3\n", + "True\n", + "[5 4]\n", + "4 5\n", + "False\n", + "4 4\n", + "True\n", + "[5]\n", + "5 5\n", + "True\n", + "[]\n" + ] + } + ], + "source": [ + "idx = np.arange(cs_povm.povm_list[0].n_outcomes)\n", + "np.random.shuffle(idx)\n", + "\n", + "for k1 in range(cs_povm.povm_list[0].n_outcomes):\n", + " for i, k2 in enumerate(idx) : \n", + " print(k1,k2)\n", + " print(np.allclose(cs_povm.povm_list[0][k1], sqpovm[k2]))\n", + " if np.allclose(cs_povm.povm_list[0][k1], sqpovm[k2]):\n", + " idx = np.delete(idx, i)\n", + " print(idx)\n", + " break\n" + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 0 9\n", + "0 3\n", + "1 0 8\n", + "0 7\n", + "2 0 7\n", + "3 1 7\n", + "4 2 7\n", + "5 3 7\n", + "6 4 7\n", + "7 5 7\n", + "5 8\n", + "8 5 6\n", + "Operator([[0.15+0.j, 0.15+0.j],\n", + " [0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0.5+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0.5+0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[ 0.15+0.j, -0.15-0.j],\n", + " [-0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0.1+0.j , 0. -0.1j],\n", + " [0. +0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0.1+0.j , 0. +0.1j],\n", + " [0. -0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "[3, 7, 8]\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 67, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "q = [0.5, 0.1, 0.2, 0.2]\n", + "\n", + "sqpovm = SingleQubitPOVM(np.array([0.15* X0, # 0\n", + " q[0]* Z0, # 1\n", + " q[0]* Z1, # 2\n", + " 0.1 * X0, # 3\n", + " 0.3 * X1, # 4\n", + " q[2]* Y0, # 5\n", + " 0.1 * Y1, # 6\n", + " 0.05* X0, # 7\n", + " 0.1 * Y1])) # 8\n", + "\n", + "n_del = 0\n", + "to_delete=[]\n", + "\n", + "for k in range(sqpovm.n_outcomes):\n", + " k1 = k - n_del\n", + " print(k,k1,sqpovm.n_outcomes)\n", + " for k2 in range(k1+1,sqpovm.n_outcomes):\n", + " if np.allclose(sqpovm[k1] / np.trace(sqpovm[k1]), sqpovm[k2] / np.trace(sqpovm[k2])):\n", + " sqpovm.povm_operators[k1] = Operator(sqpovm[k1]+sqpovm[k2])\n", + " sqpovm.povm_operators.pop(k2)\n", + " sqpovm.n_outcomes-=1\n", + "\n", + " to_delete.append(k2+n_del)\n", + " print(k1,k2+n_del)\n", + "\n", + " n_del+=1\n", + "\n", + " break\n", + "for op in sqpovm:\n", + " print(op)\n", + "\n", + "print(to_delete)\n", + "\n", + "sqpovm.check_validity()" + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " k1: 0 and #outcomes: 12\n", + "0 3\n", + "[0 0 0 1 1 1 1 1 1 1 1]\n", + "0 4\n", + "[0 0 0 2 2 2 2 2 2 2]\n", + "0 9\n", + "[0 0 0 2 2 2 2 3 3]\n", + " k1: 1 and #outcomes: 9\n", + " k1: 2 and #outcomes: 9\n", + " k1: 3 and #outcomes: 9\n", + "5 9\n", + "[0 0 0 2 2 2 4 4]\n", + "5 10\n", + "[0 0 0 2 2 2 5]\n", + "5 11\n", + "[0 0 0 2 2 2]\n", + " k1: 4 and #outcomes: 6\n", + " k1: 5 and #outcomes: 6\n", + "Operator([[0.145+0.j, 0.145+0.j],\n", + " [0.145+0.j, 0.145+0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0.51+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0.51+0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0.1+0.j , 0. +0.1j],\n", + " [0. -0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[ 0.145+0.j, -0.145-0.j],\n", + " [-0.145+0.j, 0.145+0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "Operator([[0.1+0.j , 0. -0.1j],\n", + " [0. +0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,))\n", + "[3, 4, 9, 9, 10, 11]\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 75, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "q = [0.51, 0.1, 0.2, 0.2]\n", + "\n", + "sqpovm = SingleQubitPOVM(np.array([0.15* X0, # 0\n", + " q[0]* Z0, # 1\n", + " q[0]* Z1, # 2\n", + " 0.01* X0, # 3\n", + " 0.08* X0, # 4\n", + " 0.01* Y1, # 5\n", + " 0.29 * X1, # 6\n", + " q[2]* Y0, # 7\n", + " 0.09* Y1, # 8\n", + " 0.05* X0, # 9\n", + " 0.05* Y1, # 10\n", + " 0.05* Y1])) # 11\n", + "k1=0\n", + "n_del = 0\n", + "to_delete=[]\n", + "\n", + "weird_list = np.zeros(sqpovm.n_outcomes, dtype = int)\n", + "\n", + "while k1 < sqpovm.n_outcomes:\n", + " print(f' k1: {k1} and #outcomes: {sqpovm.n_outcomes}')\n", + " k2 = k1+1\n", + " while k2 < sqpovm.n_outcomes:\n", + " if np.allclose(sqpovm[k1] / np.trace(sqpovm[k1]), sqpovm[k2] / np.trace(sqpovm[k2])):\n", + " sqpovm.povm_operators[k1] = Operator(sqpovm[k1]+sqpovm[k2])\n", + " sqpovm.povm_operators.pop(k2)\n", + " sqpovm.n_outcomes-=1\n", + "\n", + " to_delete.append(k2+n_del)\n", + " print(k1+weird_list[k1],k2+n_del)\n", + " weird_list[k2:]+=1\n", + " weird_list = np.delete(weird_list,k2)\n", + " print(weird_list)\n", + " n_del+=1\n", + " k2-=1\n", + " k2+=1\n", + " k1+=1\n", + "\n", + "for op in sqpovm:\n", + " print(op)\n", + "\n", + "print(to_delete)\n", + "\n", + "sqpovm.check_validity()" + ] + }, + { + "cell_type": "code", + "execution_count": 90, + "metadata": {}, + "outputs": [], + "source": [ + "def merge_povm_effects(povm: Povm) -> Povm:\n", + "\n", + " k1=0\n", + " n_del = 0\n", + "\n", + " while k1 < povm.n_outcomes:\n", + " print(f' k1: {k1} and #outcomes: {povm.n_outcomes}')\n", + " k2 = k1+1\n", + " while k2 < povm.n_outcomes:\n", + " if np.allclose(povm[k1] / np.trace(povm[k1]), povm[k2] / np.trace(povm[k2])):\n", + " povm.povm_operators[k1] = Operator(povm[k1]+povm[k2])\n", + " povm.povm_operators.pop(k2)\n", + " povm.n_outcomes-=1\n", + " n_del+=1\n", + " k2-=1\n", + " k2+=1\n", + " k1+=1\n", + "\n", + " assert povm.check_validity(), 'Error in cleaning the POVM'\n", + "\n", + " return povm\n", + "\n", + "def sort_povm_effects(povm: Povm) -> Povm:\n", + " sorting_values = np.array([(np.trace(op.data), op.data[0,0]) for op in povm], dtype=[('tr', 'float'), ('m00', 'float')])\n", + " idx_sort = np.argsort(sorting_values, order=('tr','m00'))[::-1]\n", + " povm.povm_operators = [povm[idx] for idx in idx_sort]\n", + "\n", + " return povm\n", + "\n", + "def clean_povm(povm: Povm) -> Povm:\n", + " povm = merge_povm_effects(povm)\n", + " povm = sort_povm_effects(povm)\n", + " return povm\n", + "\n", + "def compare_povms(povm1: Povm, povm2: Povm, tol:float=1e-5) -> bool:\n", + " povm1 = clean_povm(povm1)\n", + " povm2 = clean_povm(povm2)\n", + "\n", + " idx = np.arange(povm2.n_outcomes)\n", + "\n", + " for k1 in range(povm1.n_outcomes):\n", + " continue_check = False\n", + " for i, k2 in enumerate(idx) : \n", + " print(k1,k2)\n", + " print(np.allclose(povm1[k1], povm2[k2]))\n", + " if np.allclose(povm1[k1], povm2[k2]):\n", + " idx = np.delete(idx, i)\n", + " print(idx)\n", + " continue_check = True\n", + " break\n", + " if not continue_check :\n", + " return False\n", + "\n", + " if len(idx)==0:\n", + " return True\n", + " else:\n", + " return False" + ] + }, + { + "cell_type": "code", + "execution_count": 91, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " k1: 0 and #outcomes: 6\n", + " k1: 1 and #outcomes: 6\n", + " k1: 2 and #outcomes: 6\n", + " k1: 3 and #outcomes: 6\n", + " k1: 4 and #outcomes: 6\n", + " k1: 5 and #outcomes: 6\n", + " k1: 0 and #outcomes: 6\n", + " k1: 1 and #outcomes: 6\n", + " k1: 2 and #outcomes: 6\n", + " k1: 3 and #outcomes: 6\n", + " k1: 4 and #outcomes: 6\n", + " k1: 5 and #outcomes: 6\n", + "0 0\n", + "False\n", + "0 1\n", + "False\n", + "0 2\n", + "False\n", + "0 3\n", + "False\n", + "0 4\n", + "False\n", + "0 5\n", + "False\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/d7/fx92klcd7fncls676l53b83h0000gn/T/ipykernel_63558/1908604622.py:24: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " sorting_values = np.array([(np.trace(op.data), op.data[0,0]) for op in povm], dtype=[('tr', 'float'), ('m00', 'float')])\n" + ] + }, + { + "data": { + "text/plain": [ + "False" + ] + }, + "execution_count": 91, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "compare_povms(sqpovm, cs_povm.povm_list[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Operator([[0.51+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n" + ] + } + ], + "source": [ + "sqpovm[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Operator([[0.5+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,))\n" + ] + } + ], + "source": [ + "cs_povm.povm_list[0][0]" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/d7/fx92klcd7fncls676l53b83h0000gn/T/ipykernel_63558/2815594634.py:1: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " idx_sort = np.argsort(np.array([(np.trace(op.data), op.data[0,0]) for op in sqpovm], dtype=[('tr', 'float'), ('m00', 'float')]), order=('tr','m00'))[::-1]\n" + ] + } + ], + "source": [ + "idx_sort = np.argsort(np.array([(np.trace(op.data), op.data[0,0]) for op in sqpovm], dtype=[('tr', 'float'), ('m00', 'float')]), order=('tr','m00'))[::-1]" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[Operator([[0.5+0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0.5+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[ 0.15+0.j, -0.15-0.j],\n", + " [-0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.15+0.j, 0.15+0.j],\n", + " [0.15+0.j, 0.15+0.j]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.1+0.j , 0. -0.1j],\n", + " [0. +0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,)),\n", + " Operator([[0.1+0.j , 0. +0.1j],\n", + " [0. -0.1j, 0.1+0.j ]],\n", + " input_dims=(2,), output_dims=(2,))]" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "[sqpovm[idx] for idx in idx_sort]" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/d7/fx92klcd7fncls676l53b83h0000gn/T/ipykernel_63558/2810929492.py:1: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " np.array([(np.trace(op.data), op.data[0,0]) for op in sqpovm], dtype=[('tr', 'float'), ('m00', 'float')])\n" + ] + }, + { + "data": { + "text/plain": [ + "array([(0.3, 0.15), (0.5, 0.5 ), (0.5, 0. ), (0.2, 0.1 ), (0.3, 0.15),\n", + " (0.2, 0.1 )], dtype=[('tr', ' None: + super().__init__(methodName) + + basis_0 = np.array([1.,0], dtype=complex) + basis_1 = np.array([0,1.], dtype=complex) + basis_plus = 1.0/np.sqrt(2) * (basis_0 + basis_1) + basis_minus = 1.0/np.sqrt(2) * (basis_0 - basis_1) + basis_plus_i = 1.0/np.sqrt(2) * (basis_0 + 1.j * basis_1) + basis_minus_i = 1.0/np.sqrt(2) * (basis_0 - 1.j * basis_1) + + self.Z0 = np.outer(basis_0, basis_0.conj()) + self.Z1 = np.outer(basis_1, basis_1.conj()) + self.X0 = np.outer(basis_plus, basis_plus.conj()) + self.X1 = np.outer(basis_minus, basis_minus.conj()) + self.Y0 = np.outer(basis_plus_i, basis_plus_i.conj()) + self.Y1 = np.outer(basis_minus_i, basis_minus_i.conj()) + + + def test_CS_build(self): + """Test if we can build a standard Classical Shadow POVM from the generic class""" + + q = [1./3., 1./3., 1./3.] + sqpovm = SingleQubitPOVM(np.array([q[0]* self.Z0,q[0]* self.Z1,q[1]* self.X0,q[1]* self.X1,q[2]* self.Y0,q[2]* self.Y1])) + + for n_qubit in range(1,11): + parameters = np.array(n_qubit * [0., 0., 0.5 * np.pi, 0., 0.5 * np.pi, 0.5 * np.pi, 1, 1]) + cs_implementation = ProductPVMSimPOVMImplementation(n_qubit=n_qubit, parameters=parameters) + self.assertEqual(n_qubit, cs_implementation.n_qubit) + cs_povm = cs_implementation.to_povm() + for i in range(n_qubit): + self.assertEqual(cs_povm.povm_list[i].n_outcomes, sqpovm.n_outcomes) + for k in range(sqpovm.n_outcomes): + self.assertAlmostEqual(cs_povm.povm_list[i][k], sqpovm[k]) + + def test_LBCS_build(self): + """Test if we can build a LB Classical Shadow POVM from the generic class""" + + for n_qubit in range(1,11): + q = np.random.uniform(0,5,size=3 * n_qubit).reshape((n_qubit, 3)) + q /= q.sum(axis=1)[:, np.newaxis] + + parameters = np.array(n_qubit * [0., 0., 0.5 * np.pi, 0., 0.5 * np.pi, 0.5 * np.pi, 1, 1]) + for i in range(n_qubit): + parameters[i * 8 + 6] = q[i,0]/q[i,2] + parameters[i * 8 + 7] = q[i,1]/q[i,2] + + cs_implementation = ProductPVMSimPOVMImplementation(n_qubit=n_qubit, parameters=parameters) + self.assertEqual(n_qubit, cs_implementation.n_qubit) + cs_povm = cs_implementation.to_povm() + for i in range(n_qubit): + sqpovm = SingleQubitPOVM(np.array([q[i,0]* self.Z0,q[i,0]* self.Z1,q[i,1]* self.X0,q[i,1]* self.X1,q[i,2]* self.Y0,q[i,2]* self.Y1])) + self.assertEqual(cs_povm.povm_list[i].n_outcomes, sqpovm.n_outcomes) + for k in range(sqpovm.n_outcomes): + self.assertTrue(np.allclose(cs_povm.povm_list[i][k], sqpovm[k])) + diff --git a/test/test_sq_povm.py b/test/test_sq_povm.py new file mode 100644 index 0000000..4aecef3 --- /dev/null +++ b/test/test_sq_povm.py @@ -0,0 +1,45 @@ +"""Tests for Base POVM utils""" + +from unittest import TestCase +import numpy as np +from scipy.stats import unitary_group + + +from single_qubit_povm import SingleQubitPOVM + + +class TestSingleQubitPovm(TestCase): + """Test that we can create valid single qubit POVM and get warnings if invalid.""" + + def test_random_operators(self): + """Test """ + + ops = np.random.uniform(-1, 1, (6, 2, 2)) + 1.j * np.random.uniform(-1, 1, (6, 2, 2)) + + while np.abs(ops[0, 0, 0].imag) < 1e-6: + ops = np.random.uniform(-1, 1, (6, 2, 2)) + 1.j * np.random.uniform(-1, 1, (6, 2, 2)) + + with self.assertRaises(ValueError): + povm1 = SingleQubitPOVM(povm_ops=ops) + povm1.check_validity() + + def test_pauli_decomposition(self): + """Test """ + + # TODO : select a random POVM ... + dim = 2 + eival = np.random.uniform(low=-5, high=5, size=dim) + x = unitary_group.rvs(dim) # , random_state=seed_obs[i]) + obs = x @ np.diag(eival) @ x.T.conj() + + _, V = np.linalg.eigh(obs) + + sqpovm = SingleQubitPOVM.from_vectors(V) + + sum = np.zeros(4, dtype=complex) + for coef in sqpovm.povm_pauli_decomp: + sum += coef + + self.assertTrue(np.allclose(sum, np.array([1., 0., 0., 0.], dtype=complex))) + + # also check that the decomposition is correct TODO diff --git a/utilities.py b/utilities.py new file mode 100644 index 0000000..adc3949 --- /dev/null +++ b/utilities.py @@ -0,0 +1,27 @@ +import numpy as np + + +# Gram-Schmidt +def gs(X: np.ndarray) -> np.ndarray: + """Returns the orthonormal basis resulting from Gram-Schmidt process of X""" + + Q, _ = np.linalg.qr(X) + return Q + + +def n_sphere(param: np.ndarray) -> np.ndarray: + """Returns a unit vector on the n-sphere""" + + n = len(param) + x = np.ones(n + 1) + for i in range(n - 1): + x[i] *= np.cos(np.pi * param[i]) + x[i + 1:] *= np.sin(np.pi * param[i]) + x[-2] *= np.cos(2 * np.pi * param[-1]) + x[-1] *= np.sin(2 * np.pi * param[-1]) + + return x + + +def povms_union(): + return None