First commit

.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

README.md

@@ -0,0 +1,4 @@
+# IC-POVM
+Informationally complete POVM 
+
+This is a first try

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><psi|."""
+        povm_operators = np.zeros((povm_vectors.shape[0], povm_vectors.shape[1], povm_vectors.shape[1]), dtype=complex)
+        for i, vec in enumerate(povm_vectors):
+            povm_operators[i] = np.outer(vec, vec.conj())
+        return cls(povm_operators)
+
+
+'''
+    @classmethod
+    def from_dilation_unitary(cls, U, dim):
+        """Initialize a POVM from dilation unitary"""
+        return cls.from_vectors(U[:,0:dim].conj())
+
+    @classmethod
+    def from_param(cls, param_raw: np.ndarray, dim: int):
+        """Initialize a POVM from the list of parameters"""
+
+        assert (
+            (len(param_raw)+dim**2)%(2*dim-1) == 0
+        ), f"size of the parameters ({len(param_raw)}) does not match expectation."
+
+        n_out = (len(param_raw)+dim**2)//(2*dim-1)
+
+        param = []
+        param.append(param_raw[0:(n_out-1)])
+        count = n_out-1
+        for i in range(1,dim):
+            l = 2*(n_out-i)-1
+            param.append(param_raw[count:count+l])
+            count += l
+
+        u = np.zeros((n_out,n_out), dtype=complex)
+
+        k=0
+        u[:,k] = n_sphere(param[k])
+        u_gs = gs(u) #Gram-Schmidt
+
+        for k in range(1,dim):
+            x=n_sphere(param[k])
+            # construct k'th vector of u
+            for i in range(len(x)//2):
+                u[:,k] += (x[2*i] + x[2*i+1]*1j) * u_gs[:,k+i]
+            u_gs = gs(u)
+
+        for i in range(len(param)):
+            u_gs[:,i] *= np.sign(u[0,i])*np.sign(u_gs[0,i])
+
+        return cls.from_dilation_unitary(u_gs, dim)
+'''
+
+
+'''
+    #def __getitem__(self, index:slice) -> 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]
+
+'''

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)