Allowing arbitrary initialization of the input nodes (TeamGraphix#135)

Co-authored-by: Thierry Martinez <thierry.martinez@inria.fr>

.gitignore

@@ -14,4 +14,7 @@ graphix/sim/graphix.code-workspace
 graphix/graphix.code-workspace
 *~
 .vscode/settings.json
+graphix/sim/graphix.code-workspace
+.vscode/settings.json
+.pre-commit-config.yaml
 graphix/_version.py

CHANGELOG.md

@@ -24,6 +24,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Pauli-flow finding algorithm (#117)
 - workflow for isort, codecov (#148, #147)
 
+- Allow arbitrary states for initializing input nodes in state vector
+  and density matrix backends.
+
 ### Fixed
 - Fix output node order sorting bug in Pauli preprocessing `measure_pauli` (#145)
 
@@ -40,6 +43,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   empty output set.
 - Completely migrated to pytest, no `unittest` usage remains (#134)
 
+- Basic states are now defined in `states.BasicStates` and no longer
+  in `ops.States`.
+
 ## [0.2.11] - 2024-03-16
 
 ### Added

examples/MBQCvqe.py

@@ -22,9 +22,11 @@
 import networkx as nx
 import numpy as np
 from scipy.optimize import minimize
+
 from graphix import Circuit
 from graphix.simulator import PatternSimulator
 
+
 # %%
 # Define the Hamiltonian for the VQE problem (Example: H = Z0Z1 + X0 + X1)
 def create_hamiltonian():
@@ -33,6 +35,7 @@ def create_hamiltonian():
     H = np.kron(Z, Z) + np.kron(X, np.eye(2)) + np.kron(np.eye(2), X)
     return H
 
+
 # %%
 # Function to build the VQE circuit
 def build_vqe_circuit(n_qubits, params):
@@ -45,6 +48,7 @@ def build_vqe_circuit(n_qubits, params):
         circuit.cnot(i, i + 1)
     return circuit
 
+
 # %%
 class MBQCVQE:
     def __init__(self, n_qubits, hamiltonian):
@@ -85,6 +89,7 @@ def compute_energy(self, params):
         energy = tn.expectation_value(self.hamiltonian, qubit_indices=range(self.n_qubits))
         return energy
 
+
 # %%
 # Set parameters for VQE
 n_qubits = 2
@@ -94,18 +99,20 @@ def compute_energy(self, params):
 # Instantiate the MBQCVQE class
 mbqc_vqe = MBQCVQE(n_qubits, hamiltonian)
 
+
 # %%
 # Define the cost function
 def cost_function(params):
     return mbqc_vqe.compute_energy(params)
 
+
 # %%
 # Random initial parameters
 initial_params = np.random.rand(n_qubits * 3)
 
 # %%
 # Perform the optimization using COBYLA
-result = minimize(cost_function, initial_params, method='COBYLA', options={'maxiter': 100})
+result = minimize(cost_function, initial_params, method="COBYLA", options={"maxiter": 100})
 
 print(f"Optimized parameters: {result.x}")
 print(f"Optimized energy: {result.fun}")

graphix/linalg_validations.py

@@ -18,24 +18,29 @@ def check_square(matrix: np.ndarray) -> bool:
     return True
 
 
+def truncate(s: str, max_length: int = 80, ellipsis: str = "...") -> str:
+    "Auxilliary function to truncate a long string for formatting error messages."
+    if len(s) <= max_length:
+        return s
+    return s[: max_length - len(ellipsis)] + ellipsis
+
+
 def check_psd(matrix: np.ndarray, tol: float = 1e-15) -> bool:
     """
     check if a density matrix is positive semidefinite by diagonalizing.
-    After check_square and check_hermitian (osef) so that it already is square with power of 2 dimension.
-
 
     Parameters
     ----------
     matrix : np.ndarray
-        matrix to check. Normally already square and 2**n x 2**n
+        matrix to check
     tol : float
         tolerance on the small negatives. Default 1e-15.
     """
 
     evals = np.linalg.eigvalsh(matrix)
 
     if not all(evals >= -tol):
-        raise ValueError("The matrix is not positive semi-definite.")
+        raise ValueError("The matrix {truncate(str(matrix))} is not positive semi-definite.")
 
     return True
 
@@ -70,7 +75,6 @@ def check_data_normalization(data: Union[list, tuple, np.ndarray]) -> bool:
 
 
 def check_data_dims(data: Union[list, tuple, np.ndarray]) -> bool:
-
     # convert to set to remove duplicates
     dims = set([i["operator"].shape for i in data])
 
@@ -85,7 +89,6 @@ def check_data_dims(data: Union[list, tuple, np.ndarray]) -> bool:
 
 
 def check_data_values_type(data: Union[list, tuple, np.ndarray]) -> bool:
-
     if not all(
         isinstance(i, dict) for i in data
     ):  # ni liste ni ensemble mais iterable (lazy) pas stocké, executé au besoin

graphix/ops.py

@@ -8,16 +8,6 @@
 import numpy as np
 
 
-class States:
-    plus = np.array([1.0 / np.sqrt(2), 1.0 / np.sqrt(2)])  # plus
-    minus = np.array([1.0 / np.sqrt(2), -1.0 / np.sqrt(2)])  # minus
-    zero = np.array([1.0, 0.0])  # zero
-    one = np.array([0.0, 1.0])  # one
-    iplus = np.array([1.0 / np.sqrt(2), 1.0j / np.sqrt(2)])  # +1 eigenstate of Pauli Y
-    iminus = np.array([1.0 / np.sqrt(2), -1.0j / np.sqrt(2)])  # -1 eigenstate of Pauli Y
-    vec = [plus, minus, zero, one, iplus, iminus]
-
-
 class Ops:
     """Basic single- and two-qubits operators"""
 

graphix/pattern.py

@@ -100,7 +100,7 @@ def add(self, cmd):
         cmd : list
             MBQC command.
         """
-        assert type(cmd) == list
+        assert isinstance(cmd, list)
         assert cmd[0] in ["N", "E", "M", "X", "Z", "S", "C"]
         if cmd[0] == "N":
             if cmd[1] in self.__output_nodes:
@@ -775,7 +775,7 @@ def get_layers(self):
         not_measured = set(self.__input_nodes)
         for cmd in self.__seq:
             if cmd[0] == "N":
-                if not cmd[1] in self.output_nodes:
+                if cmd[1] not in self.output_nodes:
                     not_measured = not_measured | {cmd[1]}
         depth = 0
         l_k = dict()
@@ -814,6 +814,7 @@ def connected_edges(self, node, edges):
         connected: set of tuple
                 set of connected edges
         """
+
         connected = set()
         for edge in edges:
             if edge[0] == node:
@@ -1125,10 +1126,10 @@ def connected_nodes(self, node, prepared=None):
         if not ind == "end":  # end -> 'node' is isolated
             while self.__seq[ind][0] == "E":
                 if self.__seq[ind][1][0] == node:
-                    if not self.__seq[ind][1][1] in prepared:
+                    if self.__seq[ind][1][1] not in prepared:
                         node_list.append(self.__seq[ind][1][1])
                 elif self.__seq[ind][1][1] == node:
-                    if not self.__seq[ind][1][0] in prepared:
+                    if self.__seq[ind][1][0] not in prepared:
                         node_list.append(self.__seq[ind][1][0])
                 ind += 1
         return node_list
@@ -1217,7 +1218,7 @@ def _reorder_pattern(self, meas_commands):
         # add isolated nodes
         for cmd in self.__seq:
             if cmd[0] == "N":
-                if not cmd[1] in prepared:
+                if cmd[1] not in prepared:
                     new.append(["N", cmd[1]])
         for cmd in self.__seq:
             if cmd[0] == "E":

graphix/random_objects.py

@@ -1,5 +1,6 @@
+from __future__ import annotations
+
 import numpy as np
-import numpy.typing as npt
 import scipy.linalg
 from scipy.stats import unitary_group
 
@@ -8,48 +9,48 @@
 from graphix.sim.density_matrix import DensityMatrix
 
 
-def rand_herm(l: int):
+def rand_herm(sz: int):
     """
-    generate random hermitian matrix of size l*l
+    generate random hermitian matrix of size sz*sz
     """
-    tmp = np.random.rand(l, l) + 1j * np.random.rand(l, l)
+    tmp = np.random.rand(sz, sz) + 1j * np.random.rand(sz, sz)
     return tmp + tmp.conj().T
 
 
-def rand_unit(l: int):
+def rand_unit(sz: int):
     """
-    generate haar random unitary matrix of size l*l
+    generate haar random unitary matrix of size sz*sz
     """
-    if l == 1:
+    if sz == 1:
         return np.array([np.exp(1j * np.random.rand(1) * 2 * np.pi)])
     else:
-        return unitary_group.rvs(l)
+        return unitary_group.rvs(sz)
 
 
 UNITS = np.array([1, 1j])
 
 
-def rand_dm(dim: int, rank: int = None, dm_dtype=True) -> DensityMatrix:
-    """Returns a "density matrix" as a DensityMatrix object ie a positive-semidefinite (hence Hermitian) matrix with unit trace
-    Note, not a proper DM since its dim can be something else than a power of 2.
-    The rank is random between 1 (pure) and dim if not specified
-    Thanks to Ulysse Chabaud.
+def rand_dm(dim: int, rank: int | None = None, dm_dtype=True) -> DensityMatrix | np.ndarray:
+    """Utility to generate random density matrices (positive semi-definite matrices with unit trace).
+    Returns either a :class:`graphix.sim.density_matrix.DensityMatrix` or a :class:`np.ndarray` depending on the parameter `dm_dtype`.
 
     :param dim: Linear dimension of the (square) matrix
     :type dim: int
-    :param rank: If rank not specified then random between 1 and matrix dimension
-        If rank is one : then pure state else mixed state. Defaults to None
+    :param rank: Rank of the density matrix (1 = pure state). If not specified then sent to dim (maximal rank).
+        Defaults to None
     :type rank: int, optional
-    :param dm_dtype: If True returns a :class:`graphix.sim.density_matrix.DensityMatrix` or a numpy.ndarray if False. Defaults to True.
+    :param dm_dtype: If `True` returns a :class:`graphix.sim.density_matrix.DensityMatrix` object. If `False`returns a :class:`np.ndarray`
     :type dm_dtype: bool, optional
-    :return: Random density matrix as a :class:`graphix.sim.density_matrix.DensityMatrix` object or a numpy.ndarray.
-    :rtype: :class:`graphix.sim.density_matrix.DensityMatrix` or numpy.ndarray.
-
+    :return: the density matrix in the specified format.
+    :rtype: DensityMatrix | np.ndarray
+    .. note::
+        Thanks to Ulysse Chabaud.
+    .. warning::
+        Note that setting `dm_dtype=False` allows to generate "density matrices" inconsistent with qubits i.e. with dimensions not being powers of 2.
     """
 
-    # if not provided, use a random value.
     if rank is None:
-        rank = np.random.randint(1, dim + 1)
+        rank = dim
 
     evals = np.random.rand(rank)
 
@@ -68,7 +69,7 @@ def rand_dm(dim: int, rank: int = None, dm_dtype=True) -> DensityMatrix:
         return dm
 
 
-def rand_gauss_cpx_mat(dim: int, sig: float = 1 / np.sqrt(2)) -> npt.NDArray:
+def rand_gauss_cpx_mat(dim: int, sig: float = 1 / np.sqrt(2)) -> np.ndarray:
     """
     Returns a square array of standard normal complex random variates.
     Code from QuTiP: https://qutip.org/docs/4.0.2/modules/qutip/random_objects.html