hardware experiment plots

.gitignore

@@ -14,4 +14,5 @@ __pypackages__/
 .mvsp/
 *.venv/
 pyenv.cfg
-docs/source/api_reference/*.rst
+docs/source/api_reference/*.rst
+data/hardware_experiment/density_data/

mvsp/preprocessing/gauss.py

@@ -0,0 +1,39 @@
+import numpy as np
+from scipy.stats import multivariate_normal
+
+
+def gaussian_function(mean=None, cov=None):
+    rv = multivariate_normal(mean=mean, cov=cov)
+    return rv.pdf
+
+
+def gaussian_function_old(
+    vec, mean, cov_inv
+):  # note, takes inverse covariance in input
+    factor = np.sqrt(np.linalg.det(cov_inv)) / (2 * np.pi)
+    value = factor * np.exp(
+        -0.5 * (vec - mean).transpose().conjugate() @ cov_inv @ (vec - mean)
+    )
+    return value[0, 0]
+
+
+def coefficient_function_old(vec, mean, cov):  # note, takes covariance in input
+    factor = 2 ** (-2)
+    value = factor * np.exp(
+        -np.pi * 1j * mean.transpose().conjugate() @ vec
+        - 0.5 * np.pi**2 * vec.transpose().conjugate() @ cov @ vec
+    )
+    return value[0, 0]
+
+
+def gauss_coefficient_function(mean, cov):  # note, takes covariance in input
+    factor = 2 ** (-2)
+
+    def fun(vec):
+        return factor * np.exp(
+            -np.pi * 1j * mean.conjugate() @ vec
+            - 0.5 * np.pi**2 * vec.conjugate() @ cov @ vec
+        )
+
+    fun = np.vectorize(fun, signature="(a)->()")
+    return fun

mvsp/utils/discrete_density.py

@@ -0,0 +1,181 @@
+"""Implementation of the discrete density module."""
+
+import numpy as np
+from numpy.typing import NDArray
+from scipy.stats import multivariate_normal
+
+
+def gaussian_function(mean=None, cov=None):
+    rv = multivariate_normal(mean=mean, cov=cov)
+    return rv.pdf
+
+
+class DiscreteDensity:
+    def __init__(
+        self,
+        discrete_positive_function: NDArray,
+        coordinates: list[NDArray] | None = None,
+        dx: list[float] | None = None,
+        domain: list[tuple[float, float]] | None = None,
+        normalize: bool = True,
+        positivize: str | bool = False,
+    ):
+        if dx is not None:
+            assert len(dx) == len(discrete_positive_function.shape)
+            if domain is None:
+                raise Exception("If dx is defined you also need to define domain.")
+        elif coordinates is not None:
+            assert len(coordinates) == len(discrete_positive_function.shape)
+        else:
+            raise Exception("Either coordinates or dx must be defined.")
+
+        self.coordinates = coordinates
+        self.dx = dx
+        self.domain = domain
+        self.normalize = normalize
+        self.positivize = positivize
+        if positivize is True:
+            self.positivize = "absolute"
+        if self.positivize == "absolute":
+            self.unnormalized_density = np.abs(discrete_positive_function)
+        elif self.positivize == "absolute square":
+            self.unnormalized_density = np.abs(discrete_positive_function) ** 2
+        else:
+            self.unnormalized_density = discrete_positive_function
+
+        self.normalization = self._get_normalization_constant()
+        if normalize:
+            self.data = self.unnormalized_density / self.normalization
+        else:
+            self.data = self.unnormalized_density
+        self.dim = len(self.unnormalized_density.shape)
+
+    def _get_normalization_constant(self):
+        integrator = self.unnormalized_density
+        for i in range(len(self.unnormalized_density.shape) - 1, -1, -1):
+            if self.dx is not None:
+                integrator = np.trapz(integrator, dx=self.dx[i], axis=i)
+            else:
+                integrator = np.trapz(integrator, x=self.coordinates[i], axis=i)
+        return integrator
+
+    def get_marginal(self, axes: tuple, normalize: bool | None = None):
+        marginal_discrete_density = self.data.copy()
+
+        if normalize is None:
+            normalize = self.normalize
+
+        for i in range(len(axes) - 1, -1, -1):
+            axis = axes[i]
+            if self.dx is not None:
+                marginal_discrete_density = np.trapz(
+                    marginal_discrete_density, dx=self.dx[axis], axis=axis
+                )
+            else:
+                marginal_discrete_density = np.trapz(
+                    marginal_discrete_density, x=self.coordinates[axis], axis=axis
+                )
+
+        new_axes = list(set(range(len(self.data.shape))) - set(axes))
+        if self.dx is None:
+            new_dx = None
+            new_domain = None
+        else:
+            new_dx = [self.dx[axis] for axis in new_axes]
+            new_domain = [self.domain[axis] for axis in new_axes]
+        if self.coordinates is None:
+            new_coordinates = None
+        else:
+            new_coordinates = [self.coordinates[axis] for axis in new_axes]
+        return DiscreteDensity(
+            marginal_discrete_density,
+            coordinates=new_coordinates,
+            dx=new_dx,
+            domain=new_domain,
+            normalize=normalize,
+        )
+
+    def integrate_discrete_function(self, discrete_function):
+        assert discrete_function.shape == self.data.shape
+        integrator = discrete_function * self.data
+        for i in range(len(discrete_function.shape) - 1, -1, -1):
+            if self.dx is not None:
+                integrator = np.trapz(integrator, dx=self.dx[i], axis=i)
+            else:
+                integrator = np.trapz(integrator, x=self.coordinates[i], axis=i)
+
+        return integrator
+
+    def _single_axis_meshgrid(self, axis):
+        assert (
+            len(self.data.shape) < 53
+        ), "Dimension of density must be smaller than 53."
+
+        from string import ascii_lowercase, ascii_uppercase
+
+        alphabet = list(ascii_lowercase) + list(ascii_uppercase)
+
+        other_axes = list(set(range(len(self.data.shape))) - set([axis]))
+
+        einsum_input_string = (
+            "".join([alphabet[i] for i in other_axes]) + f",{alphabet[axis]}"
+        )
+        einsum_output_string = "".join(
+            [alphabet[i] for i in range(len(self.data.shape))]
+        )
+        einsum_string = einsum_input_string + "->" + einsum_output_string
+
+        if self.dx is not None:
+            array_list = [np.ones(np.array(self.data.shape)[other_axes])] + [
+                np.linspace(
+                    self.domain[axis][0],
+                    self.domain[axis][1],
+                    round((self.domain[axis][1] - self.domain[axis][0]) / self.dx[axis])
+                    + 1,
+                )
+            ]
+        else:
+            array_list = [np.ones(np.array(self.data.shape)[other_axes])] + [
+                self.coordinates[axis]
+            ]
+
+        return np.einsum(einsum_string, *array_list)
+
+    def mean_value_of_axis(self, axis):
+        integrator = self._single_axis_meshgrid(axis=axis)
+
+        return self.integrate_discrete_function(integrator)
+
+    def mean_value(self):
+        return [self.mean_value_of_axis(axis) for axis in range(len(self.data.shape))]
+
+    def covariance_matrix(self):
+        mean = self.mean_value()
+
+        X = [
+            self._single_axis_meshgrid(axis=axis)
+            for axis in range(len(self.data.shape))
+        ]
+
+        integrator = [
+            [(X[i] - mean[i]) * (X[j] - mean[j]) for j in range(len(mean))]
+            for i in range(len(mean))
+        ]
+
+        covariance_matrix = np.array(
+            [
+                [
+                    self.integrate_discrete_function(integrator[i][j])
+                    for j in range(len(mean))
+                ]
+                for i in range(len(mean))
+            ]
+        )
+
+        return covariance_matrix
+
+    def correlation_matrix(self):
+        covariance_matrix = self.covariance_matrix()
+        standard_deviations = np.sqrt(covariance_matrix.diagonal())
+        std_outer = np.outer(standard_deviations, standard_deviations)
+        return covariance_matrix / std_outer

mvsp/utils/shot_distribution.py

@@ -0,0 +1,102 @@
+from itertools import product
+
+import numpy as np
+from numpy.typing import NDArray
+from scipy.ndimage import gaussian_filter
+from scipy.stats import moment as scipy_moment
+
+
+class ShotDistribution:
+    def __init__(
+        self,
+        shots: list[tuple] | NDArray,
+        outcome_dict: dict,
+    ):
+        self.shots = shots
+        self.outcome_dict = outcome_dict
+        self._outcomes = None
+        self._shot_density = None
+        self._outcome_density = None
+
+    @property
+    def outcomes(self) -> list:
+        if self._outcomes is None:
+            self._outcomes = [self.outcome_dict[tuple(shot)] for shot in self.shots]
+        return self._outcomes
+
+    @property
+    def shot_density(self):
+        if self._shot_density is None:
+            bit_tuples = list(product([0, 1], repeat=len(self.shots[0])))
+            shots_tuples = list(map(tuple, self.shots))
+            shot_density = {}
+            for tup in bit_tuples:
+                shot_density[tup] = shots_tuples.count(tup) / len(shots_tuples)
+            self._shot_density = shot_density
+        return self._shot_density
+
+    @property
+    def outcome_density(self):
+        if self._outcome_density is None:
+            shot_density = self.shot_density
+            outcomes_ordered = list(self.outcome_dict.values())
+            if not isinstance(outcomes_ordered[0], (int, float)):
+                print("Transform outcomes to tuples.")
+                outcomes_ordered = list(map(tuple, outcomes_ordered))
+            self._outcome_density = dict(zip(outcomes_ordered, shot_density.values()))
+        return self._outcome_density
+
+    def filtered_outcome_density(self, sigma=1.0, dimension=1, shape=None):
+        keys = list(self.outcome_density.keys())
+        vals = list(self.outcome_density.values())
+        if dimension == 1:
+            filtered_vals = gaussian_filter(vals, sigma=sigma)
+            res = {}
+            res["outcomes"] = keys
+            res["values"] = filtered_vals
+        else:
+            keys_reshape = np.array(keys).reshape(tuple(list(shape) + [dimension]))
+            vals_reshape = np.array(vals).reshape(shape)
+            filtered_vals_reshape = gaussian_filter(vals_reshape, sigma=sigma)
+            res = {}
+            res["outcomes"] = keys_reshape
+            res["values"] = filtered_vals_reshape
+        return res
+
+    def statistical_moment(self, moment: int, center: float | None = None):
+        return scipy_moment(
+            self.outcomes,
+            moment=moment,
+            center=center,
+        )
+
+
+def generate_outcome_dict(n_qubits, shots_postselect, dim=1):
+    if dim == 1:
+        bit_tuples = np.array(list(product([0, 1], repeat=n_qubits)))
+        n_all_outcome = 2**n_qubits
+
+        x_min = 0
+        x_max = 1
+        outcomes = np.linspace(x_min, x_max, n_all_outcome)
+
+        assert len(outcomes) == len(bit_tuples)
+        outcome_dict = dict(zip(map(tuple, bit_tuples), outcomes))
+    elif dim == 2:
+        assert not n_qubits % 2, "n_qubits must be even for dim = 2."
+
+        bit_tuples = np.array(list(product([0, 1], repeat=len(shots_postselect[0]))))
+        n_all_outcome = 2 ** (n_qubits // 2)
+        print(n_all_outcome)
+
+        x_min = 0
+        x_max = 1
+        x_space = np.linspace(x_min, x_max, n_all_outcome)
+        xx, yy = np.meshgrid(x_space, x_space)
+        pos = np.stack((xx, yy), axis=-1).transpose((1, 0, 2))
+        outcomes = pos.reshape((np.prod(pos.shape[:2]), 2))
+
+        assert len(outcomes) == len(bit_tuples)
+        outcome_dict = dict(zip(map(tuple, bit_tuples), outcomes))
+
+    return outcome_dict

requirements.txt

@@ -22,4 +22,5 @@ sympy==1.12
 stim
 nbconvert
 nbformat
-matplotlib
+matplotlib
+scikit-learn