Merge pull request #56 from cirKITers/sampling

Sampling

CONTRIBUTING.md

@@ -1,12 +1,11 @@
 # Contributing to QML-Essentials
 
-:tada: Welcome! :tada:
+Contributions are highly welcome! :hugging_face:
 
-Contributions are highly welcome!
 Start of by..
 1. Creating an issue using one of the templates (Bug Report, Feature Request)
    - let's discuss what's going wrong or what should be added
-   - can you contribute with code? Great! Go ahead! :rocket:
+   - can you contribute with code? Great! Go ahead! :rocket: 
 2. Forking the repository and working on your stuff. See the sections below for details on how to set things up.
 3. Creating a pull request to the main repository
 
@@ -46,17 +45,6 @@ Publishing includes
 - creating a release with the current git tag and automatically generated release notes
 - publishing the package to PyPI using the stored credentials
 
-## Re-Installing Package
-
-If you want to overwrite the latest build you can simply push to the package index with the recent changes.
-Updating from this index however is then a bit tricky, because poetry keeps a cache of package metadata.
-To overwrite an already installed package with the same version (but different content) follow these steps:
-1. `poetry remove qml_essentials`
-2. `poetry cache clear quantum --all`
-3. `poetry add qml_essentials@latest` (or a specific version)
-
-Note that this will also re-evaluate parts of other dependencies, and thus may change the `poetry.lock` file significantly.
-
 ## Documentation
 
 We use MkDocs for our documentation. To run a server locally, run:

README.md

@@ -1,6 +1,6 @@
 # QML Essentials
 
-[![version](https://img.shields.io/badge/version-0.1.12-green.svg)](https://ea3a0fbb-599f-4d83-86f1-0e71abe27513.ka.bw-cloud-instance.org/lc3267/quantum/) [![Quality](https://github.com/cirKITers/qml-essentials/actions/workflows/quality.yml/badge.svg)](https://github.com/cirKITers/qml-essentials/actions/workflows/quality.yml) [![Testing](https://github.com/cirKITers/qml-essentials/actions/workflows/test.yml/badge.svg)](https://github.com/cirKITers/qml-essentials/actions/workflows/test.yml) [![Documentation](https://github.com/cirKITers/qml-essentials/actions/workflows/docs.yml/badge.svg)](https://github.com/cirKITers/qml-essentials/actions/workflows/docs.yml)
+[![Quality](https://github.com/cirKITers/qml-essentials/actions/workflows/quality.yml/badge.svg)](https://github.com/cirKITers/qml-essentials/actions/workflows/quality.yml) [![Testing](https://github.com/cirKITers/qml-essentials/actions/workflows/test.yml/badge.svg)](https://github.com/cirKITers/qml-essentials/actions/workflows/test.yml) [![Documentation](https://github.com/cirKITers/qml-essentials/actions/workflows/docs.yml/badge.svg)](https://github.com/cirKITers/qml-essentials/actions/workflows/docs.yml)
 
 ## 📜 About
 
@@ -22,4 +22,4 @@ You can find details on how to use it and further documentation on the correspon
 
 ## 🚧 Contributing
 
-Contributions are very welcome! See [Contribution Guidelines](https://github.com/cirKITers/qml-essentials/blob/main/CONTRIBUTING.md).
+Contributions are highly welcome! Take a look at our [Contribution Guidelines](https://github.com/cirKITers/qml-essentials/blob/main/CONTRIBUTING.md).

docs/entanglement.md

@@ -20,7 +20,7 @@ ent_cap = Entanglement.meyer_wallach(
 )
 ```
 
-Here, `n_samples` is the number of samples for the input domain in $[0, 2\pi]$, and `seed` is the random number generator seed.
+Here, `n_samples` is the number of samples for the parameters, sampled according to the default initialization strategy of the model, and `seed` is the random number generator seed.
 Note, that every function in this class accepts keyword-arguments which are being passed to the model call, so you could e.g. enable caching by
 ```python
 ent_cap = Entanglement.meyer_wallach(

docs/expressibility.md

@@ -9,15 +9,16 @@ model = Model(
 )
 
 input_domain, bins, dist_circuit = Expressibility.state_fidelities(
-    n_bins=10,
+    seed=1000,
     n_samples=200,
+    n_bins=10,
     n_input_samples=5,
-    seed=1000,
+    input_domain=[0, 2*np.pi],
     model=model,
 )
 ```
 
-Here, `n_bins` is the number of bins that you want to use in the histogram, `n_samples` is the number of parameter sets to generate, `n_input_samples` is the number of samples for the input domain in $[-\pi, \pi]$, and `seed` is the random number generator seed.
+Here, `n_bins` is the number of bins that you want to use in the histogram, `n_samples` is the number of parameter sets to generate (using the default initialization strategy of the model), `n_input_samples` is the number of samples for the input domain in $[0, 2\pi]$, and `seed` is the random number generator seed.
 
 Note that `state_fidelities` accepts keyword arguments that are being passed to the model call.
 This allows you to utilize e.g. caching.

qml_essentials/entanglement.py

@@ -1,7 +1,8 @@
-from typing import Callable, Optional, Any
+from typing import Optional, Any
 import pennylane as qml
 import pennylane.numpy as np
 
+from qml_essentials.model import Model
 import logging
 
 log = logging.getLogger(__name__)
@@ -11,20 +12,14 @@ class Entanglement:
 
     @staticmethod
     def meyer_wallach(
-        model: Callable,  # type: ignore
-        n_samples: int,
-        seed: Optional[int],
-        **kwargs: Any
+        model: Model, n_samples: int, seed: Optional[int], **kwargs: Any  # type: ignore
     ) -> float:
         """
         Calculates the entangling capacity of a given quantum circuit
         using Meyer-Wallach measure.
 
         Args:
-            model (Callable): Function that models the quantum circuit. It must
-                have a `n_qubits` attribute representing the number of qubits.
-                It must accept a `params` argument representing the parameters
-                of the circuit.
+            model (Callable): Function that models the quantum circuit.
             n_samples (int): Number of samples per qubit.
             seed (Optional[int]): Seed for the random number generator.
             kwargs (Any): Additional keyword arguments for the model function.
@@ -33,73 +28,39 @@ def meyer_wallach(
             float: Entangling capacity of the given circuit. It is guaranteed
                 to be between 0.0 and 1.0.
         """
-
-        def _meyer_wallach(
-            evaluate: Callable[[np.ndarray], np.ndarray],
-            n_qubits: int,
-            samples: int,
-            params: np.ndarray,
-        ) -> float:
-            """
-            Calculates the Meyer-Wallach sampling of the entangling capacity
-            of a quantum circuit.
-
-            Args:
-                evaluate (Callable[[np.ndarray], np.ndarray]): Callable that
-                    evaluates the quantum circuit It must accept a `params`
-                    argument representing the parameters of the circuit and may
-                    accept additional keyword arguments.
-                n_qubits (int): Number of qubits in the circuit
-                samples (int): Number of samples to be taken
-                params (np.ndarray): Parameters of the instructor. Shape:
-                    (samples, *model.params.shape)
-
-            Returns:
-                float: Entangling capacity of the given circuit. It is
-                    guaranteed to be between 0.0 and 1.0
-            """
-            assert (
-                params.shape[0] == samples
-            ), "Number of samples does not match number of parameters"
-
-            mw_measure = np.zeros(samples, dtype=complex)
-            qb = list(range(n_qubits))
-
-            for i in range(samples):
-                # implicitly set input to none in case it's not needed
-                kwargs.setdefault("inputs", None)
-                # explicitly set execution type because everything else won't work
-                U = evaluate(params=params[i], execution_type="density", **kwargs)
-
-                entropy = 0
-
-                for j in range(n_qubits):
-                    density = qml.math.partial_trace(U, qb[:j] + qb[j + 1 :])
-                    entropy += np.trace((density @ density).real)
-
-                mw_measure[i] = 1 - entropy / n_qubits
-
-            mw = 2 * np.sum(mw_measure.real) / samples
-
-            # catch floating point errors
-            return min(max(mw, 0.0), 1.0)
-
+        rng = np.random.default_rng(seed)
         if n_samples > 0:
             assert seed is not None, "Seed must be provided when samples > 0"
             # TODO: maybe switch to JAX rng
-            rng = np.random.default_rng(seed)
-            params = rng.uniform(0, 2 * np.pi, size=(n_samples, *model.params.shape))
+            model.initialize_params(rng=rng, repeat=n_samples)
         else:
             if seed is not None:
                 log.warning("Seed is ignored when samples is 0")
             n_samples = 1
-            params = model.params.reshape(1, *model.params.shape)
+            model.initialize_params(rng=rng, repeat=1)
+
+        samples = model.params.shape[-1]
+        mw_measure = np.zeros(samples, dtype=complex)
+        qb = list(range(model.n_qubits))
+
+        # TODO: vectorize in future iterations
+        for i in range(samples):
+            # implicitly set input to none in case it's not needed
+            kwargs.setdefault("inputs", None)
+            # explicitly set execution type because everything else won't work
+            U = model(params=model.params[:, :, i], execution_type="density", **kwargs)
+
+            entropy = 0
+
+            for j in range(model.n_qubits):
+                density = qml.math.partial_trace(U, qb[:j] + qb[j + 1 :])
+                entropy += np.trace((density @ density).real)
+
+            mw_measure[i] = 1 - entropy / model.n_qubits
+
+        mw = 2 * np.sum(mw_measure.real) / samples
 
-        entangling_capability = _meyer_wallach(
-            evaluate=model,
-            n_qubits=model.n_qubits,
-            samples=n_samples,
-            params=params,
-        )
+        # catch floating point errors
+        entangling_capability = min(max(mw, 0.0), 1.0)
 
         return float(entangling_capability)

qml_essentials/expressibility.py

@@ -1,16 +1,16 @@
 import pennylane.numpy as np
-from typing import Tuple, Callable, Any
+from typing import Tuple, List, Any
 from scipy import integrate
 from scipy.special import rel_entr
 import os
 
+from qml_essentials.model import Model
+
 
 class Expressibility:
     @staticmethod
     def _sample_state_fidelities(
-        model: Callable[
-            [np.ndarray, np.ndarray], np.ndarray
-        ],  # type: ignore[name-defined]
+        model: Model,
         x_samples: np.ndarray,
         n_samples: int,
         seed: int,
@@ -20,10 +20,7 @@ def _sample_state_fidelities(
         Compute the fidelities for each pair of input samples and parameter sets.
 
         Args:
-            model (Callable[[np.ndarray, np.ndarray], np.ndarray]):
-                Function that evaluates the model. It must accept inputs
-                and params as arguments
-                and return an array of shape (n_samples, n_features).
+            model (Callable): Function that models the quantum circuit.
             x_samples (np.ndarray): Array of shape (n_input_samples, n_features)
                 containing the input samples.
             n_samples (int): Number of parameter sets to generate.
@@ -43,10 +40,9 @@ def _sample_state_fidelities(
         fidelities: np.ndarray = np.zeros((n_x_samples, n_samples))
 
         # Generate random parameter sets
-        w: np.ndarray = (
-            2 * np.pi * (1 - 2 * rng.random(size=[*model.params.shape, n_samples * 2]))
-        )
-
+        # We need two sets of parameters, as we are computing fidelities for a
+        # pair of random state vectors
+        model.initialize_params(rng=rng, repeat=n_samples * 2)
         # Batch input samples and parameter sets for efficient computation
         x_samples_batched: np.ndarray = x_samples.reshape(1, -1).repeat(
             n_samples * 2, axis=0
@@ -59,7 +55,7 @@ def _sample_state_fidelities(
             # Execution type is explicitly set to density
             sv: np.ndarray = model(
                 inputs=x_samples_batched[:, idx],
-                params=w,
+                params=model.params,
                 execution_type="density",
                 **kwargs,
             )
@@ -80,21 +76,23 @@ def _sample_state_fidelities(
 
     @staticmethod
     def state_fidelities(
-        n_bins: int,
+        seed: int,
         n_samples: int,
+        n_bins: int,
         n_input_samples: int,
-        seed: int,
-        model: Callable,  # type: ignore
+        input_domain: List[float],
+        model: Model,
         **kwargs: Any,
     ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
         """
         Sample the state fidelities and histogram them into a 2D array.
 
         Args:
-            n_bins (int): Number of histogram bins.
-            n_samples (int): Number of parameter sets to generate.
-            n_input_samples (int): Number of samples for the input domain in [-pi, pi]
             seed (int): Random number generator seed.
+            n_samples (int): Number of parameter sets to generate.
+            n_bins (int): Number of histogram bins.
+            n_input_samples (int): Number of input samples.
+            input_domain (List[float]): Input domain.
             model (Callable): Function that models the quantum circuit.
             kwargs (Any): Additional keyword arguments for the model function.
 
@@ -103,8 +101,7 @@ def state_fidelities(
                 input samples, bin edges, and histogram values.
         """
 
-        x_domain = [-1 * np.pi, 1 * np.pi]
-        x = np.linspace(x_domain[0], x_domain[1], n_input_samples, requires_grad=False)
+        x = np.linspace(*input_domain, n_input_samples, requires_grad=False)
 
         fidelities = Expressibility._sample_state_fidelities(
             x_samples=x,

qml_essentials/model.py

@@ -24,6 +24,7 @@ def __init__(
         circuit_type: Union[str, Circuit],
         data_reupload: bool = True,
         initialization: str = "random",
+        initialization_domain: List[float] = [0, 2 * np.pi],
         output_qubit: Union[List[int], int] = -1,
         shots: Optional[int] = None,
         random_seed: int = 1000,
@@ -100,9 +101,10 @@ def __init__(
         # this will also be re-used in the init method,
         # however, only if nothing is provided
         self._inialization_strategy = initialization
+        self._initialization_domain = initialization_domain
 
-        # ..here! where we only require a seed
-        self.initialize_params(random_seed)
+        # ..here! where we only require a rng
+        self.initialize_params(np.random.default_rng(random_seed))
 
         # Initialize two circuits, one with the default device and
         # one with the mixed device
@@ -205,9 +207,34 @@ def shots(self, value: Optional[int]) -> None:
             value = None
         self._shots = value
 
-    def initialize_params(self, random_seed, initialization: str = None) -> None:
+    def initialize_params(
+        self,
+        rng,
+        repeat: int = None,
+        initialization: str = None,
+        initialization_domain: List[float] = None,
+    ) -> None:
+        """
+        Initializes the parameters of the model.
+
+        Args:
+            rng: A random number generator to use for initialization.
+            repeat: The number of times to repeat the parameters.
+                If None, the number of layers is used.
+            initialization: The strategy to use for parameter initialization.
+                If None, the strategy specified in the constructor is used.
+            initialization_domain: The domain to use for parameter initialization.
+                If None, the domain specified in the constructor is used.
+
+        Returns:
+            None
+        """
+        params_shape = (
+            self._params_shape if repeat is None else [*self._params_shape, repeat]
+        )
         # use existing strategy if not specified
         initialization = initialization or self._inialization_strategy
+        initialization_domain = initialization_domain or self._initialization_domain
 
         def set_control_params(params: np.ndarray, value: float) -> np.ndarray:
             indices = self.pqc.get_control_indices(self.n_qubits)
@@ -225,25 +252,22 @@ def set_control_params(params: np.ndarray, value: float) -> np.ndarray:
                 )
             return params
 
-        rng = np.random.default_rng(random_seed)
         if initialization == "random":
             self.params: np.ndarray = rng.uniform(
-                0, 2 * np.pi, self._params_shape, requires_grad=True
+                *initialization_domain, params_shape, requires_grad=True
             )
         elif initialization == "zeros":
-            self.params: np.ndarray = np.zeros(self._params_shape, requires_grad=True)
+            self.params: np.ndarray = np.zeros(params_shape, requires_grad=True)
         elif initialization == "pi":
-            self.params: np.ndarray = (
-                np.ones(self._params_shape, requires_grad=True) * np.pi
-            )
+            self.params: np.ndarray = np.ones(params_shape, requires_grad=True) * np.pi
         elif initialization == "zero-controlled":
             self.params: np.ndarray = rng.uniform(
-                0, 2 * np.pi, self._params_shape, requires_grad=True
+                *initialization_domain, params_shape, requires_grad=True
             )
             self.params = set_control_params(self.params, 0)
         elif initialization == "pi-controlled":
             self.params: np.ndarray = rng.uniform(
-                0, 2 * np.pi, self._params_shape, requires_grad=True
+                *initialization_domain, params_shape, requires_grad=True
             )
             self.params = set_control_params(self.params, np.pi)
         else: