diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst index e857eacaf8..64d385a6ab 100644 --- a/doc/source/api-reference/qibo.rst +++ b/doc/source/api-reference/qibo.rst @@ -1952,11 +1952,22 @@ For more details, see G. Chiribella *et al.*, *Theoretical framework for quantum `Physical Review A 80.2 (2009): 022339 `_. + .. autoclass:: qibo.quantum_info.quantum_networks.QuantumNetwork :members: :member-order: bysource +.. autoclass:: qibo.quantum_info.quantum_networks.QuantumComb + :members: + :member-order: bysource + + +.. autoclass:: qibo.quantum_info.quantum_networks.QuantumChannel + :members: + :member-order: bysource + + Random Ensembles ^^^^^^^^^^^^^^^^ diff --git a/doc/source/code-examples/tutorials/quantum_networks/README.md b/doc/source/code-examples/tutorials/quantum_networks/README.md deleted file mode 100644 index e0b074640c..0000000000 --- a/doc/source/code-examples/tutorials/quantum_networks/README.md +++ /dev/null @@ -1,122 +0,0 @@ -# Quantum Networks - -## The Quantum Network Model - -The quantum network model is a mathematical framework that allows us to uniquely describe quantum information processing that involves multiple points in time and space. -Each distinguished point in time and space is treated as a linear system $\mathcal{H}_ i$. -A quantum network involving $n$ points in time and space is a Hermitian operator $\mathcal{N}$ that acts on the tensor product of the linear systems $\mathcal{H}_ 0 \otimes \mathcal{H}_ 1 \otimes \cdots \otimes \mathcal{H}_ {n-1}$. -Each system $\mathcal{H}_ {i}$ is either an input or an output of the network. - -A physically implementable quantum network is described by a semi-positive definite operator $\mathcal{N}$ that satisfies the causal constraints. - -A simple example is a quantum channel $\Gamma: \mathcal{H}_ 0 \to \mathcal{H}_ 1$, where $\mathcal{H}_ 0$ is the input system and $\mathcal{H}_ 1$ is the output system. -The quantum channel is a linear map, such that it maps any input quantum state to an output quantum state, which is a sufficient and necessary condition for the map to be physical. -A Hermitian operator $J^\Gamma$ acting on $\mathcal{H}_ 0\otimes \mathcal{H}_ 1$ is associated with a quantum channel $\Gamma$, if $J^\Gamma$ satisfies the following conditions: -$$J^\Gamma \geq 0, \quad \text{and} \quad \text{Tr}_ {\mathcal{H}_ 1} J^\Gamma = \mathbb{I}_ {\mathcal{H} _0} .$$ - -The first condition is called *complete positivity*, and the second condition is called *trace-preserving*. -In particular, the second condition ensures that the information of the input system is only accessible through the output system. - -In particular, a quantum state $\rho$ may be also considered as a quantum network, where the input system is the trivial system $\mathbb{C}$, and the output system is the quantum system $\mathcal{H}$. -The constraints on the quantum channels are then equivalent to the constraints on the quantum states: -$$\rho \geq 0, \quad \text{and} \quad \text{Tr} \rho = \mathbb{I}_ \mathbb{C} = 1\ .$$ - -> For more details, see G. Chiribella *et al.*, *Theoretical framework for quantum networks*, -> [Physical Review A 80.2 (2009): 022339](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.80.022339). - -## Quantum Network in `qibo` - -The manipulation of quantum networks in `qibo` is done through the `QuantumNetwork` class. - -```python -from qibo.quantum_info.quantum_networks import QuantumNetwork -``` - -A quantum state is a quantum network with a single input system and a single output system, where the input system is the trivial 1-dimensional system. -We need to specify the dimensions of each system in the `partition` argument. - -```python -from qibo.quantum_info import random_density_matrix, random_unitary - -state = random_density_matrix(2) -state_choi = QuantumNetwork(state, (1,2)) -print(f'A quantum state is a quantum netowrk of the form {state_choi}') -``` - -``` ->>> A quantum state is a quantum netowrk of the form J[1 -> 2] -``` - -A general quantum channel can be created in a similar way. - -```python -from qibo.gates import DepolarizingChannel - -test_ch = DepolarizingChannel(0,0.5) -N = len(test_ch.target_qubits) -partition = (2**N, 2**N) - -depolar_choi = QuantumNetwork(test_ch.to_choi(), partition) -print(f'A quantum channel is a quantum netowrk of the form {depolar_choi}') -``` - -``` ->>> A quantum channel is a quantum netowrk of the form J[2 -> 2] -``` - -One may apply a quantum channel to a quantum state, or compose two quantum channels, using the `@` operator. - -```python -new_state = depolar_choi @ depolar_choi @ state_choi -``` - -## Example - -For 3-dimensional systems, an unital channel may not be a mixed unitary channel. - -> Example 4.3 in (Watrous, John. The theory of quantum information. Cambridge university press, 2018.) - -```python -A1 = np.array([ - [0,0,0], - [0,0,1/np.sqrt(2)], - [0,-1/np.sqrt(2),0], -]) -A2 = np.array([ - [0,0,1/np.sqrt(2)], - [0,0,0], - [-1/np.sqrt(2),0,0], -]) -A3 = np.array([ - [0,1/np.sqrt(2),0], - [-1/np.sqrt(2),0,0], - [0,0,0], -]) - -Choi1 = QuantumNetwork(A1, (3,3), pure=True) * 3 -Choi2 = QuantumNetwork(A2, (3,3), pure=True)*3 -Choi3 = QuantumNetwork(A3, (3,3), pure=True)*3 -``` - -The three channels are pure but not unital. Which means they are not unitary. - -```python -print(f"Choi1 is unital: {Choi1.unital()}") -print(f"Choi2 is unital: {Choi2.unital()}") -print(f"Choi3 is unital: {Choi3.unital()}") -``` - -``` ->>> Choi1 is unital: False -Choi2 is unital: False -Choi3 is unital: False -``` - -However, the mixture of the three operators is unital. -As the matrices are orthogonal, they are the extreme points of the convex set of the unital channels. -Therefore, this mixed channel is not a mixed unitary channel. - -```python -Choi = Choi1/3 + Choi2/3 + Choi3/3 -print(f"The mixed channel is unital: {Choi.unital()}") -``` diff --git a/src/qibo/quantum_info/quantum_networks.py b/src/qibo/quantum_info/quantum_networks.py index 48b34a0ff2..691f40a515 100644 --- a/src/qibo/quantum_info/quantum_networks.py +++ b/src/qibo/quantum_info/quantum_networks.py @@ -1,7 +1,7 @@ """Module defining the `QuantumNetwork` class and adjacent functions.""" -import re from functools import reduce +from logging import warning from operator import mul from typing import List, Optional, Tuple, Union @@ -12,30 +12,30 @@ class QuantumNetwork: - """This class stores the Choi operator of the quantum network as a tensor, - which is an unique representation of the quantum network. + """This class stores the representation of the quantum network as a tensor. + This is a unique representation of the quantum network. A minimum quantum network is a quantum channel, which is a quantum network of the form :math:`J[n \\to m]`, where :math:`n` is the dimension of the input system , and :math:`m` is the dimension of the output system. - A quantum state is a quantum network of the form :math:`J[1 \\to n]`, + A quantum state is a quantum network of the form :math:`J: 1 \\to n`, such that the input system is trivial. - An observable is a quantum network of the form :math:`J[n \\to 1]`, + An observable is a quantum network of the form :math:`J: n \\to 1`, such that the output system is trivial. A quantum network may contain multiple input and output systems. - For example, a "quantum comb" is a quantum network of the form :math:`J[n', n \\to m, m']`, - which convert a quantum channel of the form :math:`J[n \\to m]` - to a quantum channel of the form :math:`J[n' \\to m']`. + For example, a "quantum comb" is a quantum network of the form :math:`J: n', n \\to m, m'`, + which convert a quantum channel of the form :math:`J: n \\to m` + to a quantum channel of the form :math:`J: n' \\to m'`. Args: - matrix (ndarray): input Choi operator. - partition (List[int] or Tuple[int]): partition of ``matrix``. - system_output (List[bool] or Tuple[bool], optional): mask on the output system of the + tensor (ndarray): input Choi operator. + partition (List[int] or Tuple[int]): partition of ``tensor``. + system_input (List[bool] or Tuple[bool], optional): mask on the output system of the Choi operator. If ``None``, defaults to - ``(False,True,False,True,...)``, where ``len(system_output)=len(partition)``. + ``(True,False,True,False,...)``, where ``len(system_input)=len(partition)``. Defaults to ``None``. - pure (bool, optional): ``True`` when ``matrix`` is a "pure" representation (e.g. a pure + pure (bool, optional): ``True`` when ``tensor`` is a "pure" representation (e.g. a pure state, a unitary operator, etc.), ``False`` otherwise. Defaults to ``False``. backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. @@ -44,162 +44,237 @@ class QuantumNetwork: def __init__( self, - matrix, - partition: Union[List[int], Tuple[int]], - system_output: Optional[Union[List[bool], Tuple[bool]]] = None, + tensor, + partition: Optional[Union[List[int], Tuple[int]]] = None, + system_input: Optional[Union[List[bool], Tuple[bool]]] = None, pure: bool = False, backend=None, ): - self._run_checks(partition, system_output, pure) - - self._matrix = matrix - self.partition = partition - self.system_output = system_output + self._tensor = tensor + self.partition = tuple(partition) + self.system_input = system_input self._pure = pure self._backend = backend - self.dims = reduce(mul, self.partition) - self._set_tensor_and_parameters() + self._run_checks(self.partition, self.system_input, self._pure) - def matrix(self, backend=None): - """Returns the Choi operator of the quantum network in matrix form. + self._set_parameters() - Args: - backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used - to return the Choi operator. If ``None``, defaults to the backend defined - when initializing the :class:`qibo.quantum_info.quantum_networks.QuantumNetwork` - object. Defaults to ``None``. + self.dims = reduce(mul, self.partition) if len(self.partition) > 0 else 1 - Returns: - ndarray: Choi matrix of the quantum network. - """ - if backend is None: # pragma: no cover - backend = self._backend + @staticmethod + def _order_tensor_to_operator(dims: int): + """Returns the order to reshape a tensor into an operator. - return backend.cast(self._matrix, dtype=self._matrix.dtype) + Given a tenosr of ``2 * dims`` leads, the order is + :math:`[0, 2, 4, ..., 1, 3, 5, ...]`. - def is_pure(self): - """Returns bool indicading if the Choi operator of the network is pure.""" - return self._pure + Args: + dims (int): dimension. - def is_hermitian( - self, order: Optional[Union[int, str]] = None, precision_tol: float = 1e-8 - ): - """Returns bool indicating if the Choi operator :math:`\\mathcal{E}` of the network is Hermitian. + Returns: + list: order to reshape tensor into an operator. + """ + return list(range(0, 2 * dims, 2)) + list(range(1, 2 * dims, 2)) - Hermicity is calculated as distance between :math:`\\mathcal{E}` and - :math:`\\mathcal{E}^{\\dagger}` with respect to a given norm. - Default is the ``Hilbert-Schmidt`` norm (also known as ``Frobenius`` norm). + @staticmethod + def _order_operator_to_tensor(nsystems: int): + """Returns the order to reshape an operator to a tensor. - For specifications on the other possible values of the - parameter ``order`` for the ``tensorflow`` backend, please refer to - `tensorflow.norm `_. - For all other backends, please refer to - `numpy.linalg.norm `_. + Given a operator of :math:`2n` systems, the order is + :math:`[0, n, 1, n+1, 2, n+2, ...]`. Args: - order (str or int, optional): order of the norm. Defaults to ``None``. - precision_tol (float, optional): threshold :math:`\\epsilon` that defines if - Choi operator of the network is :math:`\\epsilon`-close to Hermicity in - the norm given by ``order``. Defaults to :math:`10^{-8}`. + nsystems (int): number of systems. Returns: - bool: Hermiticity condition. + list: order to reshape operator into tensor. """ - if precision_tol < 0.0: + return list( + sum(zip(list(range(0, nsystems)), list(range(nsystems, nsystems * 2))), ()) + ) + + @classmethod + def _operator_to_tensor(cls, operator, partition: List[int]): + + n = len(partition) + order = cls._order_operator_to_tensor(n) + + # Check if the `partition` matches the shape of the input matrix + if np.prod(tuple(operator.shape)) != np.prod( + tuple(dim**2 for dim in partition) + ): raise_error( ValueError, - f"``precision_tol`` must be non-negative float, but it is {precision_tol}", + "``partition`` does not match the shape of the input matrix. " + + f"Cannot reshape matrix of size {operator.shape} to partition {partition}", ) - if order is None and self._backend.__class__.__name__ == "TensorflowBackend": - order = "euclidean" + # Check if `operator` is a pytourch tensor + tensor = operator.reshape(list(partition) * 2) + if operator.__class__.__name__ == "Tensor": + tensor = tensor.permute(order) + else: + tensor = tensor.transpose(order) + return tensor.reshape([dim**2 for dim in partition]) + + @classmethod + def from_operator( + cls, + operator, + partition: Optional[Union[List[int], Tuple[int]]] = None, + system_input: Optional[Union[List[bool], Tuple[bool]]] = None, + pure: bool = False, + backend=None, + ): + """Construct a :class:`qibo.quantum_info.QuantumNetwork` object from a ndarray. - self._matrix = self._full() - self._pure = False + This method converts a Choi operator to the internal representation of + :class:`qibo.quantum_info.quantum_networks.QuantumNetwork`. + The input array can be a pure state, a Choi operator, a unitary operator, etc. - reshaped = self._backend.cast( - self._backend.np.reshape(self._matrix, (self.dims, self.dims)), - dtype=self._matrix.dtype, - ) - reshaped = self._backend.cast( - self._backend.np.conj(reshaped).T - reshaped, dtype=reshaped.dtype - ) - norm = self._backend.calculate_norm_density_matrix(reshaped, order=order) + Args: + arr (ndarray): input numpy array. + partition (List[int] or Tuple[int], optional): partition of ``arr``. If ``None``, + defaults to the shape of ``arr``. Defaults to ``None``. + system_input (List[bool] or Tuple[bool], optional): mask on the input system of the + Choi operator. If ``None``, defaults to + ``(True,False,True,False...)``, where ``len(system_input)=len(partition)``. + Defaults to ``None``. + pure (bool, optional): ``True`` when ``arr`` is a "pure" representation (e.g. a pure + state, a unitary operator, etc.), ``False`` otherwise. Defaults to ``False``. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in + calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. - return float(norm) <= precision_tol + Returns: + :class:`qibo.quantum_info.quantum_networks.QuantumNetwork`: + quantum network constructed from the input Choi operator. + """ - def is_unital( - self, order: Optional[Union[int, str]] = None, precision_tol: float = 1e-8 - ): - """Returns bool indicating if the Choi operator :math:`\\mathcal{E}` of the network is unital. + if pure: + if partition is None: + partition = tuple(operator.shape) + tensor = operator + else: + if tuple(partition) not in [ + tuple(operator.shape), + tuple(int(np.sqrt(dim)) for dim in operator.shape) * 2, + ]: + raise_error( + ValueError, + "``partition`` does not match the shape of the input matrix. " + + f"Cannot reshape matrix of size {operator.shape} " + + f"to partition {partition}", + ) + + tensor = operator.reshape(partition) + else: + # check if arr is a valid choi operator + len_sys = len(operator.shape) + if (len_sys % 2 != 0) or ( + operator.shape[: len_sys // 2] != operator.shape[len_sys // 2 :] + ): + raise_error( + ValueError, + "The opertor must be a square operator where the first half of the shape " + + "is the same as the second half of the shape. " + + f"However, the shape of the input is {operator.shape}. " + + "If the input is pure, set `pure=True`.", + ) - Unitality is calculated as distance between the partial trace of :math:`\\mathcal{E}` - and the Identity operator :math:`I`, with respect to a given norm. - Default is the ``Hilbert-Schmidt`` norm (also known as ``Frobenius`` norm). + if partition is None: + partition = operator.shape[: len_sys // 2] - For specifications on the other possible values of the - parameter ``order`` for the ``tensorflow`` backend, please refer to - `tensorflow.norm `_. - For all other backends, please refer to - `numpy.linalg.norm `_. + tensor = cls._operator_to_tensor(operator, partition) + + return cls( + tensor, + partition=partition, + system_input=system_input, + pure=pure, + backend=backend, + ) + + def operator(self, full: bool = False, backend=None): + """Returns the Choi operator of the quantum network. + + The shape of the returned operator is :math:`(*self.partition, *self.partition)`. Args: - order (str or int, optional): order of the norm. Defaults to ``None``. - precision_tol (float, optional): threshold :math:`\\epsilon` that defines - if Choi operator of the network is :math:`\\epsilon`-close to unitality - in the norm given by ``order``. Defaults to :math:`10^{-8}`. + full (bool, optional): If this is ``False``, and the network is pure, the method + will only return the eigenvector (unique when the network is pure). + If ``True``, returns the full tensor of the quantum network. Defaults to ``False``. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used + to return the Choi operator. If ``None``, defaults to the backend defined + when initializing the :class:`qibo.quantum_info.quantum_networks.QuantumNetwork` + object. Defaults to ``None``. Returns: - bool: Unitality condition. + ndarray: Choi operator of the quantum network. """ - if precision_tol < 0.0: - raise_error( - ValueError, - f"``precision_tol`` must be non-negative float, but it is {precision_tol}", - ) + if backend is None: # pragma: no cover + backend = self._backend - if order is None and self._backend.__class__.__name__ == "TensorflowBackend": - order = "euclidean" + if self.is_pure() and not full: + return backend.cast(self._tensor, dtype=self._tensor.dtype) - self._matrix = self._full() - self._pure = False + tensor = self.full(backend) if self.is_pure() else self._tensor - partial_trace = self._einsum("jkjl -> kl", self._matrix) - identity = self._backend.cast( - np.eye(partial_trace.shape[0]), dtype=partial_trace.dtype - ) + n = len(self.partition) + order = self._order_tensor_to_operator(n) - norm = self._backend.calculate_norm_density_matrix( - partial_trace - identity, - order=order, + operator = self._backend.np.transpose( + tensor.reshape(tuple(np.repeat(self.partition, 2))), order ) - return float(norm) <= precision_tol + return backend.cast(operator, dtype=self._tensor.dtype) - def is_causal( + def matrix(self, backend=None): + """Returns the Choi operator of the quantum network in the matrix form. + The shape of the returned operator is :math:`(self.dims, self.dims)`. + + Args: + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used + to return the Choi operator. If ``None``, defaults to the backend defined + when initializing the :class:`qibo.quantum_info.quantum_networks.QuantumNetwork` + object. Defaults to ``None``. + + Returns: + ndarray: Choi operator of the quantum network. + """ + return self.operator(full=True, backend=backend).reshape((self.dims, self.dims)) + + def is_pure(self): + """Returns bool indicading if the Choi operator of the network is pure.""" + return self._pure + + def is_hermitian( self, order: Optional[Union[int, str]] = None, precision_tol: float = 1e-8 ): - """Returns bool indicating if the Choi operator :math:`\\mathcal{E}` of the network satisfies the causal order condition. + """Returns bool indicating if the Choi operator :math:`\\mathcal{J}` is Hermitian. - Causality is calculated as distance between partial trace of :math:`\\mathcal{E}` - and the Identity operator :math:`I`, with respect to a given norm. + Hermicity is calculated as distance between :math:`\\mathcal{J}` and + :math:`\\mathcal{J}^{\\dagger}` with respect to a given norm. Default is the ``Hilbert-Schmidt`` norm (also known as ``Frobenius`` norm). For specifications on the other possible values of the parameter ``order`` for the ``tensorflow`` backend, please refer to `tensorflow.norm `_. For all other backends, please refer to - `numpy.linalg.norm `_. + `numpy.linalg.norm + `_. Args: order (str or int, optional): order of the norm. Defaults to ``None``. - precision_tol (float, optional): threshold :math:`\\epsilon` that defines - if Choi operator of the network is :math:`\\epsilon`-close to causality - in the norm given by ``order``. Defaults to :math:`10^{-8}`. + precision_tol (float, optional): threshold :math:`\\epsilon` that defines if + Choi operator of the network is :math:`\\epsilon`-close to Hermicity in + the norm given by ``order``. Defaults to :math:`10^{-8}`. Returns: - bool: Causal order condition. + bool: Hermiticity condition. If the adjoint of the Choi operator is equal to the + Choi operator, the method returns ``True``. + If the input is pure, the its always Hermitian. """ if precision_tol < 0.0: raise_error( @@ -210,101 +285,61 @@ def is_causal( if order is None and self._backend.__class__.__name__ == "TensorflowBackend": order = "euclidean" - self._matrix = self._full() - self._pure = False + if self.is_pure(): # if the input is pure, it is always hermitian + return True - partial_trace = self._einsum("jklk -> jl", self._matrix) - identity = self._backend.cast( - np.eye(partial_trace.shape[0]), dtype=partial_trace.dtype + reshaped = self._backend.cast( + self.matrix(), + dtype=self._tensor.dtype, ) + if self._backend.__class__.__name__ == "PyTorchBackend": + adjoint = self._backend.np.transpose(reshaped, (1, 0)) + else: + adjoint = self._backend.np.transpose(reshaped) - norm = self._backend.calculate_norm_density_matrix( - partial_trace - identity, - order=order, - ) + mat_diff = self._backend.np.conj(adjoint) - reshaped + norm = self._backend.calculate_norm_density_matrix(mat_diff, order=order) return float(norm) <= precision_tol def is_positive_semidefinite(self, precision_tol: float = 1e-8): - """Returns bool indicating if Choi operator :math:`\\mathcal{E}` of the network is positive-semidefinite. + """Returns bool indicating if Choi operator :math:`\\mathcal{J}` is positive-semidefinite. Args: precision_tol (float, optional): threshold value used to check if eigenvalues of - the Choi operator :math:`\\mathcal{E}` are such that - :math:`\\textup{eigenvalues}(\\mathcal{E}) >= - \\textup{precision_tol}`. + the Choi operator :math:`\\mathcal{J}` are such that + :math:`\\textup{eigenvalues}(\\mathcal{J}) >= - \\textup{precision_tol}`. Note that this parameter can be set to negative values. Defaults to :math:`0.0`. Returns: bool: Positive-semidefinite condition. """ - self._matrix = self._full() - self._pure = False + if precision_tol < 0.0: + raise_error( + ValueError, + f"``precision_tol`` must be non-negative float, but it is {precision_tol}", + ) + + if self.is_pure(): # if the input is pure, it is always positive semidefinite + return True + + reshaped = self._backend.cast( + self.matrix(), + dtype=self._tensor.dtype, + ) - reshaped = self._backend.np.reshape(self._matrix, (self.dims, self.dims)) if self.is_hermitian(): eigenvalues = self._backend.calculate_eigenvalues(reshaped) else: - if self._backend.__class__.__name__ in [ - "CupyBackend", - "CuQuantumBackend", - ]: # pragma: no cover - reshaped = np.array(reshaped.tolist(), dtype=reshaped.dtype) - eigenvalues = self._backend.calculate_eigenvalues(reshaped, hermitian=False) + return False return all( self._backend.np.real(eigenvalue) >= -precision_tol for eigenvalue in eigenvalues ) - def is_channel( - self, - order: Optional[Union[int, str]] = None, - precision_tol_causal: float = 1e-8, - precision_tol_psd: float = 1e-8, - ): - """Returns bool indicating if Choi operator :math:`\\mathcal{E}` is a channel. - - Args: - order (int or str, optional): order of the norm used to calculate causality. - Defaults to ``None``. - precision_tol_causal (float, optional): threshold :math:`\\epsilon` that defines if - Choi operator of the network is :math:`\\epsilon`-close to causality in the norm - given by ``order``. Defaults to :math:`10^{-8}`. - precision_tol_psd (float, optional): threshold value used to check if eigenvalues of - the Choi operator :math:`\\mathcal{E}` are such that - :math:`\\textup{eigenvalues}(\\mathcal{E}) >= \\textup{precision_tol_psd}`. - Note that this parameter can be set to negative values. - Defaults to :math:`0.0`. - - Returns: - bool: Channel condition. - """ - return self.is_causal( - order, precision_tol_causal - ) and self.is_positive_semidefinite(precision_tol_psd) - - def apply(self, state): - """Apply the Choi operator :math:`\\mathcal{E}` to ``state`` :math:`\\varrho`. - - It is assumed that ``state`` :math:`\\varrho` is a density matrix. - - Args: - state (ndarray): density matrix of a ``state``. - - Returns: - ndarray: Resulting state :math:`\\mathcal{E}(\\varrho)`. - """ - matrix = self._backend.cast(self._matrix, copy=True) - - if self.is_pure(): - return self._einsum( - "kj,ml,jl -> km", matrix, self._backend.np.conj(matrix), state - ) - - return self._einsum("jklm,km -> jl", matrix, state) - - def link_product(self, second_network, subscripts: str = "ij,jk -> ik"): + def link_product(self, subscripts: str, second_network): """Link product between two quantum networks. The link product is not commutative. Here, we assume that @@ -313,97 +348,40 @@ def link_product(self, second_network, subscripts: str = "ij,jk -> ik"): in order to simplify notation. Args: - second_network (:class:`qibo.quantum_info.quantum_networks.QuantumNetwork`): Quantum - network to be applied to the original network. subscripts (str, optional): Specifies the subscript for summation using the Einstein summation convention. For more details, please refer to - `numpy.einsum `_. + `numpy.einsum + `_. + second_network (:class:`qibo.quantum_info.quantum_networks.QuantumNetwork`): Quantum + network to be applied to the original network. Returns: :class:`qibo.quantum_info.quantum_networks.QuantumNetwork`: Quantum network resulting from the link product between two quantum networks. """ - if not isinstance(second_network, QuantumNetwork): - raise_error( - TypeError, - "It is not possible to implement link product of a " - + "``QuantumNetwork`` with a non-``QuantumNetwork``.", - ) - if not isinstance(subscripts, str): - raise_error( - TypeError, - f"subscripts must be type str, but it is type {type(subscripts)}.", - ) - - subscripts = subscripts.replace(" ", "") - pattern_two, pattern_four = self._check_subscript_pattern(subscripts) - channel_subscripts = pattern_two and subscripts[1] == subscripts[3] - inv_subscripts = pattern_two and subscripts[0] == subscripts[4] - super_subscripts = ( - pattern_four - and subscripts[1] == subscripts[5] - and subscripts[2] == subscripts[6] - ) - - if not channel_subscripts and not inv_subscripts and not super_subscripts: - raise_error( - NotImplementedError, - "Subscripts do not match any implemented pattern.", - ) + return link_product(subscripts, self, second_network, backend=self._backend) - first_matrix = self._full() - second_matrix = second_network._full() # pylint: disable=W0212 - - if super_subscripts: - cexpr = "jklmnopq,klop->jmnq" - return QuantumNetwork( - self._einsum(cexpr, first_matrix, second_matrix), - [self.partition[0] + self.partition[-1]], - backend=self._backend, - ) - - cexpr = "jkab,klbc->jlac" - - if inv_subscripts: - return QuantumNetwork( - self._einsum(cexpr, second_matrix, first_matrix), - [second_network.partition[0], self.partition[1]], - backend=self._backend, - ) - - return QuantumNetwork( - self._einsum(cexpr, first_matrix, second_matrix), - [self.partition[0], second_network.partition[1]], + def copy(self): + """Returns a copy of the :class:`qibo.quantum_info.QuantumNetwork` object.""" + return self.__class__( + self._backend.np.copy(self._tensor), + partition=self.partition, + system_input=self.system_input, + pure=self._pure, backend=self._backend, ) - def copy(self): - """Returns a copy of the :class:`qibo.quantum_info.quantum_networks.QuantumNetwork` object.""" + def conj(self): + """Returns the conjugate of the quantum network.""" return self.__class__( - self._backend.np.copy(self._matrix), + self._backend.np.conj(self._tensor), partition=self.partition, - system_output=self.system_output, + system_input=self.system_input, pure=self._pure, backend=self._backend, ) - def to_full(self, backend=None): - """Convert the internal representation to the full Choi operator of the network. - - Returns: - (:class:`qibo.quantum_info.quantum_networks.QuantumNetwork`): The full representation - of the Quantum network. - """ - if backend is None: # pragma: no cover - backend = self._backend - - if self.is_pure(): - self._matrix = self._full() - self._pure = False - - return self.matrix(backend) - def __add__(self, second_network): """Add two Quantum Networks by adding their Choi operators. @@ -424,23 +402,23 @@ def __add__(self, second_network): + f"and and object of type ``{type(second_network)}``.", ) - if self._full().shape != second_network._full().shape: + if self.full().shape != second_network.full().shape: raise_error( ValueError, - f"The Choi operators must have the same shape, but {self._matrix.shape} != " - + f"{second_network.matrix(second_network._backend).shape}.", + f"The Choi operators must have the same shape, but {self.full().shape} != " + + f"{second_network.full().shape}.", ) - if self.system_output != second_network.system_output: - raise_error(ValueError, "The networks must have the same output system.") + if self.system_input != second_network.system_input: + raise_error(ValueError, "The networks must have the same input systems.") - new_first_matrix = self._full() - new_second_matrix = second_network._full() + new_first_tensor = self.full() + new_second_tensor = second_network.full() return QuantumNetwork( - new_first_matrix + new_second_matrix, + new_first_tensor + new_second_tensor, self.partition, - self.system_output, + self.system_input, pure=False, backend=self._backend, ) @@ -468,19 +446,19 @@ def __mul__(self, number: Union[float, int]): if self.is_pure() and number > 0.0: return QuantumNetwork( - np.sqrt(number) * self.matrix(backend=self._backend), + np.sqrt(number) * self._tensor, partition=self.partition, - system_output=self.system_output, + system_input=self.system_input, pure=True, backend=self._backend, ) - matrix = self._full() + tensor = self.full() return QuantumNetwork( - number * matrix, + number * tensor, partition=self.partition, - system_output=self.system_output, + system_input=self.system_input, pure=False, backend=self._backend, ) @@ -513,9 +491,9 @@ def __truediv__(self, number: Union[float, int]): number = np.sqrt(number) if self.is_pure() and number > 0.0 else number return QuantumNetwork( - self.matrix(backend=self._backend) / number, + self._tensor / number, partition=self.partition, - system_output=self.system_output, + system_input=self.system_input, pure=self.is_pure(), backend=self._backend, ) @@ -579,29 +557,21 @@ def __matmul__(self, second_network): + "Use `link_product` method to specify the subscript.", ) - return self.link_product(second_network, subscripts=subscripts) + return self.link_product(subscripts, second_network) # pylint: disable=E0606 def __str__(self): """Method to define how to print relevant information of the quantum network.""" - string_in = ", ".join( - [ - str(self.partition[k]) - for k in range(len(self.partition)) - if not self.system_output[k] - ] - ) + systems = [] - string_out = ", ".join( - [ - str(self.partition[k]) - for k in range(len(self.partition)) - if self.system_output[k] - ] - ) + for i, dim in enumerate(self.partition): + if self.system_input[i]: + systems.append(f"┍{dim}┑") + else: + systems.append(f"┕{dim}┙") - return f"J[{string_in} -> {string_out}]" + return f"J[{', '.join(systems)}]" - def _run_checks(self, partition, system_output, pure): + def _run_checks(self, partition, system_input, pure): """Checks if all inputs are correct in type and value.""" if not isinstance(partition, (list, tuple)): raise_error( @@ -624,11 +594,11 @@ def _run_checks(self, partition, system_output, pure): + "but contains non-positive integers.", ) - if system_output is not None and len(system_output) != len(partition): + if system_input is not None and len(system_input) != len(partition): raise_error( ValueError, - "``len(system_output)`` must be the same as ``len(partition)``, " - + f"but {len(system_output)} != {len(partition)}.", + "``len(system_input)`` must be the same as ``len(partition)``, " + + f"but {len(system_input)} != {len(partition)}.", ) if not isinstance(pure, bool): @@ -637,55 +607,512 @@ def _run_checks(self, partition, system_output, pure): f"``pure`` must be type ``bool``, but it is type ``{type(pure)}``.", ) - def _set_tensor_and_parameters(self): - """Sets tensor based on inputs.""" + @staticmethod + def _check_system_input(system_input, partition) -> Tuple[bool]: + """ + If `system_input` not defined, assume the network follows the order of a quantum Comb. + """ + + if system_input is None: + system_input = [ + False, + ] * len(partition) + for k in range(len(partition) // 2): + system_input[k * 2] = True + return tuple(system_input) + + def _set_parameters(self): + """Standarize the parameters.""" self._backend = _check_backend(self._backend) + self.partition = tuple(self.partition) + + self.system_input = self._check_system_input(self.system_input, self.partition) + self._einsum = self._backend.np.einsum + self._tensordot = self._backend.np.tensordot + self._tensor = self._backend.cast(self._tensor, dtype=self._tensor.dtype) + + if self._pure: + if np.prod(tuple(self._tensor.shape)) != np.prod(tuple(self.partition)): + raise_error( + ValueError, + "``partition`` does not match the shape of the input matrix. " + + f"Cannot reshape matrix of size {self._tensor.shape} " + + f"to partition {self.partition}.", + ) + self._tensor = self._backend.np.reshape(self._tensor, self.partition) + else: + if np.prod(tuple(self._tensor.shape)) != np.prod( + tuple(dim**2 for dim in self.partition) + ): + raise_error( + ValueError, + "``partition`` does not match the shape of the input matrix. " + + f"Cannot reshape matrix of size {self._tensor.shape} " + + f"to partition {self.partition}.", + ) + matrix_partition = [dim**2 for dim in self.partition] + self._tensor = self._backend.np.reshape(self._tensor, matrix_partition) + + def full(self, update: bool = False, backend=None): + """Convert the internal representation to the full tensor of the network. + + Args: + update (bool, optional): If ``True``, updates the internal representation of the + network to the full tensor. Defaults to ``False``. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in + calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + + Returns: + ndarray: full reprentation of the quantum network. + """ + if backend is None: # pragma: no cover + backend = self._backend + tensor = self._backend.np.copy(self._tensor) + tensor = backend.cast(tensor, dtype=self._tensor.dtype) + conj = backend.np.conj + + if self.is_pure(): + # Reshapes input matrix based on purity. + tensor.reshape(self.dims) + if self._backend.__class__.__name__ == "PyTorchBackend": + tensor = self._tensordot(tensor, conj(tensor), dims=0) + else: + tensor = self._tensordot(tensor, conj(tensor), axes=0) + tensor = self._operator_to_tensor(tensor, self.partition) + + if update: + self._tensor = tensor + self._pure = False + + return tensor + - if isinstance(self.partition, list): - self.partition = tuple(self.partition) +class QuantumComb(QuantumNetwork): + """Stores a Quantum comb, which is a network in which the systems follows a sequential order. - try: - if self._pure: - self._matrix = self._backend.np.reshape(self._matrix, self.partition) + It is also called the *non-Markovian quantum process* in many literatures. + A quantum comb is a quantum network of the form :math:`J[┍i1┑,┕o1┙,┍i2┑,┕o2┙, ...]`, + where the process first take an input state from system :math:`i1`, + then output a state to system :math:`o1`, and so on. + This is a non-Markovian process as the output of the system :math:`o2` may depend on + what happened in systems :math:`i1`, and :math:`o1`. + + A quantum channel is a special case of quantum comb, where there are only one input + system and one output system. + + Args: + tensor (ndarray): the tensor representations of the quantum Comb. + partition (List[int] or Tuple[int]): partition of ``matrix``. + system_input (List[bool] or Tuple[bool], optional): mask on the input system of the + Choi operator. If ``None``, defaults to + ``(True,False,True,False,...)``, where ``len(system_input)=len(partition)``. + Defaults to ``None``. + pure (bool, optional): ``True`` when ``tensor`` is a "pure" representation (e.g. a pure + state, a unitary operator, etc.), ``False`` otherwise. Defaults to ``False``. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in + calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + """ + + def __init__( + self, + tensor, + partition: Optional[Union[List[int], Tuple[int]]] = None, + system_input: Optional[Union[List[bool], Tuple[bool]]] = None, + pure: bool = False, + backend=None, + ): + if partition is None: + if pure: + partition = tensor.shape else: - matrix_partition = self.partition * 2 - self._matrix = self._backend.np.reshape(self._matrix, matrix_partition) - except: + partition = tuple(int(np.sqrt(d)) for d in tensor.shape) + if len(partition) % 2 != 0: raise_error( ValueError, - "``partition`` does not match the shape of the input matrix. " - + f"Cannot reshape matrix of size {self._matrix.shape} to partition {self.partition}", + "A quantum comb should only contain equal number of input and output systems. " + + "For general quantum networks, one should use the ``QuantumNetwork`` class.", + ) + if system_input is not None: + warning("system_input is ignored for QuantumComb") + + super().__init__( + tensor, partition, [True, False] * (len(partition) // 2), pure, backend + ) + + def is_causal( + self, order: Optional[Union[int, str]] = None, precision_tol: float = 1e-8 + ): + """Returns bool indicating if the Choi operator :math:`\\mathcal{J}` satisfies causal order + + Causality is calculated based on a recursive constrains. + This method reduce a n-comb to a (n-1)-comb at each step, + and checks if the reduced comb is independent on the last output system. + + Args: + order (str or int, optional): order of the norm. Defaults to ``None``. + precision_tol (float, optional): threshold :math:`\\epsilon` that defines + if Choi operator of the network is :math:`\\epsilon`-close to causality + in the norm given by ``order``. Defaults to :math:`10^{-8}`. + + Returns: + bool: Causal order condition. + """ + if precision_tol < 0.0: + raise_error( + ValueError, + f"``precision_tol`` must be non-negative float, but it is {precision_tol}", ) - if self.system_output is None: - self.system_output = [ - True, - ] * len(self.partition) - for k in range(len(self.partition) // 2): - self.system_output[k * 2] = False - self.system_output = tuple(self.system_output) + if order is None and self._backend.__class__.__name__ == "TensorflowBackend": + order = "euclidean" + + backend = self._backend + + dim_out = self.partition[-1] + dim_in = self.partition[-2] + + trace_out = TraceOperation(dim_out, backend=backend).full() + trace_in = TraceOperation(dim_in, backend=backend).full() + + if self._backend.__class__.__name__ == "PyTorchBackend": + reduced = self._tensordot(self.full(), trace_out, dims=([-1], [0])) + sub_comb = self._tensordot(reduced, trace_in, dims=([-1], [0])) + expected = self._tensordot(sub_comb, trace_in / dim_in, dims=0) else: - self.system_output = tuple(self.system_output) + reduced = self._tensordot(self.full(), trace_out, axes=(-1, 0)) + sub_comb = self._tensordot(reduced, trace_in, axes=(-1, 0)) + expected = self._tensordot(sub_comb, trace_in / dim_in, axes=0) - def _full(self): - """Reshapes input matrix based on purity.""" - matrix = self._backend.cast(self._matrix, copy=True) + norm = self._backend.calculate_norm(reduced - expected, order=order) - if self.is_pure(): - matrix = self._einsum( - "jk,lm -> kjml", matrix, self._backend.np.conj(matrix) + if float(norm) > precision_tol: + return False + + if len(self.partition) == 2: + return True + + return QuantumComb( + sub_comb, self.partition[:-2], pure=False, backend=self._backend + ).is_causal(order, precision_tol) + + @classmethod + def from_operator( + cls, operator, partition=None, inverse=False, pure=False, backend=None + ): # pylint: disable=W0237 + comb = super().from_operator(operator, partition, None, pure, backend) + if inverse: + # Convert mathmetical convention of Choi operator to physical convention + comb.partition = comb.partition[::-1] + comb._tensor = comb._tensor.T # pylint: disable=W0212 + return comb + + +class QuantumChannel(QuantumComb): + """Stores a Quantum channel, which is a special case of quantum comb. + + A quantum channel is a quantum comb with only one input and one output. + This class includes all quantum channels, unitary operators, and quantum states. + + To construct a `QuantumChannel` object, one can use the `QuantumNetwork.from_nparray` method. + **Note**: if one try to construct a quantum network from a unitary operator or Choi operator, + the first system will be the output. + However, here we assume the first system is the input system. + It is important to specify `inverse=True` when constructing by `QuantumNetwork.from_nparray`. + + Args: + tensor (ndarray): the tensor representations of the quantum comb. + partition (List[int] or Tuple[int], optional): partition of ``matrix``. + If not provided and `system_input` is `None`, assume the input is a quantum state, + whose input is a trivial system. If `system_input` is set to `True`, + assume the input is an observable, whose output is a trivial system. + system_input (List[bool] or Tuple[bool], optional): mask on the input system of the + Choi operator. If ``None`` the default is ``(True,False)``. + Defaults to ``None``. + pure (bool, optional): ``True`` when ``tensor`` is a "pure" representation (e.g. a pure + state, a unitary operator, etc.), ``False`` otherwise. Defaults to ``False``. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in + calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + """ + + def __init__( + self, + tensor, + partition: Optional[Union[List[int], Tuple[int]]] = None, + system_input: Optional[Union[List[bool], Tuple[bool]]] = None, + pure: bool = False, + backend=None, + ): + if isinstance(partition, int): + partition = (partition,) + + if partition is not None: + if len(partition) > 2: + raise_error( + ValueError, + "A quantum channel should only contain one input system and one output system." + + "For general quantum networks, one should use the ``QuantumNetwork`` class.", + ) + if len(partition) == 1: + if system_input is None: # Assume the input is a quantum state + partition = (1, partition[0]) + else: + if isinstance(system_input, bool): + system_input = (system_input,) + + partition = ( + (partition[0], 1) if system_input[0] else (1, partition[0]) + ) + + super().__init__(tensor, partition, pure=pure, backend=backend) + + def is_unital( + self, order: Optional[Union[int, str]] = None, precision_tol: float = 1e-8 + ): + """Returns bool indicating if the Choi operator :math:`\\mathcal{J}` is unital. + + A map is unital if it preserves the identity operator. + Unitality is calculated as distance between the partial trace of :math:`\\mathcal{J}` + and the Identity operator :math:`I`, with respect to a given norm. + Default is the ``Hilbert-Schmidt`` norm (also known as ``Frobenius`` norm). + + For specifications on the other possible values of the + parameter ``order`` for the ``tensorflow`` backend, please refer to + `tensorflow.norm `_. + For all other backends, please refer to + `numpy.linalg.norm + `_. + + Args: + order (str or int, optional): order of the norm. Defaults to ``None``. + precision_tol (float, optional): threshold :math:`\\epsilon` that defines + if Choi operator of the network is :math:`\\epsilon`-close to unitality + in the norm given by ``order``. Defaults to :math:`10^{-8}`. + + Returns: + bool: Unitality condition. + """ + if precision_tol < 0.0: + raise_error( + ValueError, + f"``precision_tol`` must be non-negative float, but it is {precision_tol}", + ) + + if order is None and self._backend.__class__.__name__ == "TensorflowBackend": + order = "euclidean" + + backend = self._backend + + dim_out = self.partition[1] + dim_in = self.partition[0] + + trace_out = TraceOperation(dim_out, backend=backend).full() + trace_in = TraceOperation(dim_in, backend=backend).full() + + if self._backend.__class__.__name__ == "PyTorchBackend": + reduced = self._tensordot(self.full(), trace_in, dims=([0], [0])) + sub_comb = self._tensordot( + reduced, + trace_out, + dims=([0], [0]), ) + expected = self._tensordot(trace_out / dim_out, sub_comb, dims=0) + else: + reduced = self._tensordot(self.full(), trace_in, axes=(0, 0)) + sub_comb = self._tensordot(reduced, trace_out, axes=(0, 0)) + expected = self._tensordot(trace_out / dim_out, sub_comb, axes=0) + + norm = self._backend.calculate_norm((reduced - expected), order=order) + if float(norm) > precision_tol: + return False + + if len(self.partition) == 2: + return True - return matrix + # Unital is defined for quantum channels only. + # But we can extend it to quantum combs as follows: + return QuantumChannel( # pragma: no cover + sub_comb, self.partition[2:], pure=False, backend=self._backend + ).is_unital(order, precision_tol) - def _check_subscript_pattern(self, subscripts: str): - """Checks if input subscript match any implemented pattern.""" - braket = "[a-z]" - pattern_two = re.compile(braket * 2 + "," + braket * 2 + "->" + braket * 2) - pattern_four = re.compile(braket * 4 + "," + braket * 2 + "->" + braket * 2) + def is_channel( + self, + order: Optional[Union[int, str]] = None, + precision_tol_causal: float = 1e-8, + precision_tol_psd: float = 1e-8, + ): + """Returns bool indicating if Choi operator :math:`\\mathcal{E}` is a channel. - return bool(re.match(pattern_two, subscripts)), bool( - re.match(pattern_four, subscripts) + Args: + order (int or str, optional): order of the norm used to calculate causality. + Defaults to ``None``. + precision_tol_causal (float, optional): threshold :math:`\\epsilon` that defines if + Choi operator of the network is :math:`\\epsilon`-close to causality in the norm + given by ``order``. Defaults to :math:`10^{-8}`. + precision_tol_psd (float, optional): threshold value used to check if eigenvalues of + the Choi operator :math:`\\mathcal{E}` are such that + :math:`\\textup{eigenvalues}(\\mathcal{E}) >= \\textup{precision_tol_psd}`. + Note that this parameter can be set to negative values. + Defaults to :math:`0.0`. + + Returns: + bool: Channel condition. + """ + return self.is_causal( + order, precision_tol_causal + ) and self.is_positive_semidefinite(precision_tol_psd) + + def apply(self, state): + """Apply the Choi operator :math:`\\mathcal{E}` to ``state`` :math:`\\varrho`. + + It is assumed that ``state`` :math:`\\varrho` is a density matrix. + + Args: + state (ndarray): density matrix of a ``state``. + + Returns: + ndarray: Resulting state :math:`\\mathcal{E}(\\varrho)`. + """ + operator = self.copy().operator() + conj = self._backend.np.conj + + if self.is_pure(): + return self._einsum("ij,lk,il", operator, conj(operator), state) + + return self._einsum("ijkl, jl", operator, state) + + +def link_product( + subscripts: str, + *operands: QuantumNetwork, + backend=None, + surpress_warning=False, +): + """Link product between two quantum networks. + + The link product is not commutative. Here, we assume that + :math:`A.\\textup{link_product}(B)` means "applying :math:`B` to :math:`A`". + However, the ``link_product`` is associative, so we override the `@` operation + in order to simplify notation. + + Args: + subscripts (str, optional): Specifies the subscript for summation using + the Einstein summation convention. For more details, please refer to + `numpy.einsum `_. + operands (:class:`qibo.quantum_info.quantum_networks.QuantumNetwork`): Quantum + networks to be contracted. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in + calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + surpress_warning (bool, optional): If ``True``, surpresses the warning + regarding if the same index connects two input or two output + systems. Defaults to ``False``. + + Returns: + :class:`qibo.quantum_info.quantum_networks.QuantumNetwork`: Quantum network resulting + from the link product between two quantum networks. + """ + + if not isinstance(subscripts, str): + raise_error( + TypeError, + f"subscripts must be type str, but it is type {type(subscripts)}.", ) + + for i, operand in enumerate(operands): + if not isinstance(operand, QuantumNetwork): + raise_error(TypeError, f"The {i}-th operator is not a ``QuantumNetwork``.") + + if backend is None: # pragma: no cover + backend = operands[0]._backend # pylint: disable=W0212 + + tensors = [ + ( + backend.to_numpy(operand.full()) + if operand.is_pure() + else backend.to_numpy(operand._tensor) # pylint: disable=W0212 + ) + for operand in operands + ] + + # keep track of the `partition` and `system_input` of the network + _, contracrtion_list = np.einsum_path( + subscripts, *tensors, optimize=False, einsum_call=True + ) + + inds, idx_rm, einsum_str, _, _ = contracrtion_list[0] + input_str, results_index = einsum_str.split("->") + inputs = input_str.split(",") + + # Warning if the same index connects two input or two output systems + if not surpress_warning: + for ind in idx_rm: + found = 0 + for i, script in enumerate(inputs): + index = script.find(ind) + if index < 0: + continue + found += 1 + if found > 1 and is_input == operands[inds[i]].system_input[index]: + warning( + f"Index {ind} connects two {'input' if is_input else 'output'} systems." + ) + is_input = operands[inds[i]].system_input[index] + if found > 2: + warning( + f"Index {ind} appears multiple times in the input subscripts {input_str}." + ) + + # set correct order of the `partition` and `system_input` + partition = [] + system_input = [] + for ind in results_index: + for i, script in enumerate(inputs): + index = script.find(ind) + if index < 0: + continue + + partition.append(operands[inds[i]].partition[index]) + system_input.append(operands[inds[i]].system_input[index]) + + new_tensor = np.einsum(subscripts, *tensors) + + return QuantumNetwork(new_tensor, partition, system_input, backend=backend) + + +class IdentityChannel(QuantumChannel): + """The Identity channel with the given dimension. + + Args: + dim (int): Dimension of the Identity operator. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in + calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + """ + + def __init__(self, dim: int, backend=None): + + identity = np.eye(dim, dtype=complex) + identity = backend.cast(identity, dtype=identity.dtype) + super().__init__(identity, [dim, dim], pure=True, backend=backend) + + +class TraceOperation(QuantumNetwork): + """The Trace operator with the given dimension. + + Args: + dim (int): Dimension of the Trace operator. + backend (:class:`qibo.backends.abstract.Backend`, optional): Backend to be used in + calculations. If ``None``, defaults to :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + """ + + def __init__(self, dim: int, backend=None): + + identity = np.eye(dim, dtype=complex) + identity = backend.cast(identity, dtype=identity.dtype) + super().__init__(identity, [dim], [True], pure=False, backend=backend) diff --git a/tests/test_quantum_info_quantum_networks.py b/tests/test_quantum_info_quantum_networks.py index 5ebbbea3d3..daa8ab53f7 100644 --- a/tests/test_quantum_info_quantum_networks.py +++ b/tests/test_quantum_info_quantum_networks.py @@ -4,7 +4,14 @@ import pytest from qibo import gates -from qibo.quantum_info.quantum_networks import QuantumNetwork +from qibo.quantum_info.quantum_networks import ( + IdentityChannel, + QuantumChannel, + QuantumComb, + QuantumNetwork, + TraceOperation, + link_product, +) from qibo.quantum_info.random_ensembles import ( random_density_matrix, random_gaussian_matrix, @@ -18,7 +25,13 @@ def test_errors(backend): nqubits = len(channel.target_qubits) dims = 2**nqubits partition = (dims, dims) - network = QuantumNetwork( + network = QuantumNetwork.from_operator( + channel.to_choi(backend=backend), partition, backend=backend + ) + quantum_comb = QuantumComb.from_operator( + channel.to_choi(backend=backend), partition, backend=backend + ) + quantum_channel = QuantumChannel.from_operator( channel.to_choi(backend=backend), partition, backend=backend ) @@ -32,7 +45,7 @@ def test_errors(backend): comb_sys_out = (False, True) * 2 comb = random_density_matrix(2**4, backend=backend) comb_choi = QuantumNetwork( - comb, comb_partition, system_output=comb_sys_out, backend=backend + comb, comb_partition, system_input=comb_sys_out, backend=backend ) with pytest.raises(TypeError): @@ -46,20 +59,26 @@ def test_errors(backend): with pytest.raises(ValueError): QuantumNetwork( - channel.to_choi(backend=backend), partition=(1, 2), system_output=(1, 2, 3) + channel.to_choi(backend=backend), partition=(1, 2), system_input=(1, 2, 3) ) with pytest.raises(TypeError): QuantumNetwork(channel.to_choi(backend=backend), partition=(1, 2), pure="True") + with pytest.raises(TypeError): + QuantumNetwork(channel.to_choi(backend=backend), partition=1, pure=True) + with pytest.raises(ValueError): network.is_hermitian(precision_tol=-1e-8) with pytest.raises(ValueError): - network.is_unital(precision_tol=-1e-8) + network.is_positive_semidefinite(precision_tol=-1e-8) + + with pytest.raises(ValueError): + quantum_comb.is_causal(precision_tol=-1e-8) with pytest.raises(ValueError): - network.is_causal(precision_tol=-1e-8) + quantum_channel.is_unital(precision_tol=-1e-8) with pytest.raises(TypeError): network + 1 @@ -75,23 +94,20 @@ def test_errors(backend): network_2 = network.copy() with pytest.raises(ValueError): - network_2.system_output = (False,) + network_2.system_input = (False,) network += network_2 # Multiplying QuantumNetwork with non-QuantumNetwork with pytest.raises(TypeError): - network @ network.matrix(backend) + network @ network.operator(backend=backend) # Linking QuantumNetwork with non-QuantumNetwork with pytest.raises(TypeError): - network.link_product(network.matrix(backend)) + network.link_product(network.operator(backend=backend)) with pytest.raises(TypeError): network.link_product(network, subscripts=True) - with pytest.raises(NotImplementedError): - network.link_product(network, subscripts="jk,lm->no") - with pytest.raises(NotImplementedError): net @ net @@ -113,6 +129,44 @@ def test_errors(backend): with pytest.raises(ValueError): QuantumNetwork(matrix, (1, 2), backend=backend) + with pytest.raises(ValueError): + QuantumNetwork(matrix, (1, 1), pure=True, backend=backend) + + with pytest.raises(ValueError): + QuantumNetwork.from_operator(matrix, (1, 2), pure=True, backend=backend) + + vec = np.random.rand(4) + vec = backend.cast(vec, dtype=vec.dtype) + vec = backend.cast(vec, dtype=vec.dtype) + with pytest.raises(ValueError): + QuantumNetwork.from_operator(vec, backend=backend) + + with pytest.raises(ValueError): + QuantumComb.from_operator(vec, pure=True, backend=backend) + + with pytest.raises(ValueError): + QuantumChannel(matrix, partition=(2, 2, 2), pure=True, backend=backend) + + with pytest.raises(TypeError): + link_product(1, quantum_comb, backend=backend) + + with pytest.raises(TypeError): + link_product("ij, i", quantum_comb, matrix, backend=backend) + + # raise warning + link_product("ii", quantum_channel, backend=backend) + link_product("ij, kj", network_state, quantum_channel, backend=backend) + link_product("ij, jj", network_state, quantum_channel, backend=backend) + link_product( + "ij, jj, jj", network_state, quantum_channel, quantum_channel, backend=backend + ) + + +def test_class_methods(backend): + matrix = random_density_matrix(2**2, backend=backend) + with pytest.raises(ValueError): + QuantumNetwork._operator_to_tensor(matrix, (3,)) + def test_operational_logic(backend): lamb = float(np.random.rand()) @@ -120,7 +174,7 @@ def test_operational_logic(backend): nqubits = len(channel.target_qubits) dims = 2**nqubits partition = (dims, dims) - network = QuantumNetwork( + network = QuantumNetwork.from_operator( channel.to_choi(backend=backend), partition, backend=backend ) @@ -129,31 +183,40 @@ def test_operational_logic(backend): # Sum with itself has to match multiplying by int backend.assert_allclose( - (network + network).matrix(backend), (2 * network).matrix(backend) + (network + network).operator(backend=backend), + (2 * network).operator(backend=backend), ) backend.assert_allclose( - (network_state_pure + network_state_pure).matrix(backend), - (2 * network_state_pure).to_full(), + (network_state_pure + network_state_pure).operator(backend=backend), + (2 * network_state_pure).operator(full=True, backend=backend), ) # Sum with itself has to match multiplying by float backend.assert_allclose( - (network + network).matrix(backend), (2.0 * network).matrix(backend) + (network + network).operator(backend=backend), + (2.0 * network).operator(backend=backend), ) backend.assert_allclose( - (network_state_pure + network_state_pure).matrix(backend), - (2.0 * network_state_pure).to_full(), + (network_state_pure + network_state_pure).operator(backend=backend), + (2.0 * network_state_pure).operator(full=True, backend=backend), ) # Multiplying and dividing by same scalar has to bring back to original network backend.assert_allclose( - ((2.0 * network) / 2).matrix(backend), network.matrix(backend) + ((2.0 * network) / 2).operator(backend=backend), + network.operator(backend=backend), ) unitary = random_unitary(dims, backend=backend) network_unitary = QuantumNetwork(unitary, (dims, dims), pure=True, backend=backend) backend.assert_allclose( - (network_unitary / 2).matrix(backend), unitary / np.sqrt(2), atol=1e-5 + (network_unitary / 2).operator(backend=backend), unitary / np.sqrt(2), atol=1e-5 + ) + + # Complex conjugate of a network has to match the complex conjugate of the operator + backend.assert_allclose( + network.conj().operator(backend=backend), + backend.np.conj(network.operator(backend=backend)), ) @@ -165,20 +228,33 @@ def test_parameters(backend): dims = 2**nqubits partition = (dims, dims) - network = QuantumNetwork( - channel.to_choi(backend=backend), partition, backend=backend + choi = channel.to_choi(backend=backend) + + network = QuantumNetwork.from_operator(choi, partition, backend=backend) + quantum_comb = QuantumComb.from_operator(choi, partition, backend=backend) + quantum_channel = QuantumChannel.from_operator( + choi, partition, backend=backend, inverse=True + ) + + rand = random_density_matrix(dims**2, backend=backend) + non_channel = QuantumChannel.from_operator( + rand, partition, backend=backend, inverse=True ) - backend.assert_allclose(network.matrix(backend=backend).shape, (2, 2, 2, 2)) + backend.assert_allclose(network.operator(backend=backend).shape, (2, 2, 2, 2)) backend.assert_allclose(network.dims, 4) backend.assert_allclose(network.partition, partition) - backend.assert_allclose(network.system_output, (False, True)) + backend.assert_allclose(network.system_input, (True, False)) - assert network.is_causal() - assert network.is_unital() assert network.is_hermitian() assert network.is_positive_semidefinite() - assert network.is_channel() + assert quantum_comb.is_causal() + assert quantum_channel.is_unital() + assert quantum_channel.is_channel() + + # Test non-unital and non_causal + assert not non_channel.is_causal() + assert not non_channel.is_unital() def test_with_states(backend): @@ -186,24 +262,20 @@ def test_with_states(backend): dims = 2**nqubits state = random_density_matrix(dims, backend=backend) - network_state = QuantumNetwork(state, (1, 2), backend=backend) + network_state = QuantumChannel.from_operator(state, backend=backend) lamb = float(np.random.rand()) channel = gates.DepolarizingChannel(0, lamb) - network_channel = QuantumNetwork( - channel.to_choi(backend=backend), (dims, dims), backend=backend + network_channel = QuantumChannel.from_operator( + channel.to_choi(backend=backend), (dims, dims), backend=backend, inverse=True ) state_output = channel.apply_density_matrix(backend, state, nqubits) state_output_network = network_channel.apply(state) - state_output_link = network_state.link_product( - network_channel, subscripts="ij,jk -> ik" - ) + state_output_link = network_state.link_product("ij,jk -> ik", network_channel) backend.assert_allclose(state_output_network, state_output) - backend.assert_allclose( - state_output_link.matrix(backend=backend).reshape((dims, dims)), state_output - ) + backend.assert_allclose(state_output_link.matrix(backend=backend), state_output) assert network_state.is_hermitian() assert network_state.is_positive_semidefinite() @@ -217,51 +289,120 @@ def test_with_unitaries(backend, subscript): unitary_1 = random_unitary(dims, backend=backend) unitary_2 = random_unitary(dims, backend=backend) - network_1 = QuantumNetwork(unitary_1, (dims, dims), pure=True, backend=backend) - network_2 = QuantumNetwork(unitary_2, (dims, dims), pure=True, backend=backend) - network_3 = QuantumNetwork( - unitary_2 @ unitary_1, (dims, dims), pure=True, backend=backend + network_1 = QuantumComb.from_operator( + unitary_1, (dims, dims), pure=True, backend=backend, inverse=True ) - network_4 = QuantumNetwork( - unitary_1 @ unitary_2, (dims, dims), pure=True, backend=backend + network_2 = QuantumComb.from_operator( + unitary_2, (dims, dims), pure=True, backend=backend, inverse=True + ) + network_3 = QuantumComb.from_operator( + unitary_2 @ unitary_1, (dims, dims), pure=True, backend=backend, inverse=True + ) + network_4 = QuantumComb.from_operator( + unitary_1 @ unitary_2, (dims, dims), pure=True, backend=backend, inverse=True ) - test = network_1.link_product(network_2, subscript).to_full(backend=backend) + test = network_1.link_product(subscript, network_2).full( + backend=backend, update=True + ) if subscript[1] == subscript[3]: - backend.assert_allclose(test, network_3.to_full(backend), atol=1e-8) + backend.assert_allclose( + test, network_3.full(backend=backend, update=True), atol=1e-8 + ) backend.assert_allclose( - test, (network_1 @ network_2).to_full(backend=backend), atol=1e-8 + test, (network_1 @ network_2).full(backend=backend, update=True), atol=1e-8 ) if subscript[0] == subscript[4]: - backend.assert_allclose(test, network_4.to_full(backend)) + backend.assert_allclose(test, network_4.full(backend)) - backend.assert_allclose(test, (network_2 @ network_1).to_full(backend=backend)) + backend.assert_allclose(test, (network_2 @ network_1).full(backend=backend)) + + # Check properties for pure states + assert network_1.is_causal() + assert network_1.is_hermitian() + assert network_1.is_positive_semidefinite() def test_with_comb(backend): subscript = "jklm,kl->jm" comb_partition = (2,) * 4 channel_partition = (2,) * 2 - comb_sys_out = (False, True) * 2 - channel_sys_out = (False, True) + comb_sys_in = (False, True) * 2 + channel_sys_in = (False, True) + + rand_choi = random_density_matrix(4**2, backend=backend) + unitary_1 = random_unitary(4, backend=backend) + unitary_2 = random_unitary(4, backend=backend) + non_channel = QuantumNetwork.from_operator( + rand_choi, + (2, 2, 2, 2), + system_input=(True, True, False, False), + backend=backend, + ) + unitary_channel = QuantumNetwork.from_operator( + unitary_1, + (2, 2, 2, 2), + system_input=(True, True, False, False), + pure=True, + backend=backend, + ) + unitary_channel2 = QuantumNetwork.from_operator( + unitary_2, + (2, 2, 2, 2), + system_input=(True, True, False, False), + pure=True, + backend=backend, + ) + + non_comb = link_product( + "ij kl, km on -> jl mn", non_channel, unitary_channel, backend=backend + ) + non_comb = QuantumComb( + non_comb.full(backend=backend), + (2, 2, 2, 2), + system_input=(True, False, True, False), + backend=backend, + ) + two_comb = link_product( + "ij kl, km on, i, o", + unitary_channel, + unitary_channel2, + TraceOperation(2, backend=backend), + TraceOperation(2, backend=backend), + backend=backend, + ) + two_comb = QuantumComb( + two_comb.full(backend=backend), + (2, 2, 2, 2), + system_input=(True, False, True, False), + backend=backend, + ) comb = random_density_matrix(2**4, backend=backend) channel = random_density_matrix(2**2, backend=backend) - comb_choi = QuantumNetwork( - comb, comb_partition, system_output=comb_sys_out, backend=backend + comb_choi = QuantumNetwork.from_operator( + comb, comb_partition, system_input=comb_sys_in, backend=backend ) - channel_choi = QuantumNetwork( - channel, channel_partition, system_output=channel_sys_out, backend=backend + channel_choi = QuantumNetwork.from_operator( + channel, channel_partition, system_input=channel_sys_in, backend=backend ) - test = comb_choi.link_product(channel_choi, subscript).to_full(backend) + test = comb_choi.link_product(subscript, channel_choi).full( + update=True, backend=backend + ) channel_choi2 = comb_choi @ channel_choi - backend.assert_allclose(test, channel_choi2.to_full(backend), atol=1e-5) + backend.assert_allclose(test, channel_choi2.full(backend), atol=1e-5) + + assert non_comb.is_hermitian() + assert not non_comb.is_causal() + + assert two_comb.is_hermitian() + assert two_comb.is_causal() def test_apply(backend): @@ -270,7 +411,9 @@ def test_apply(backend): state = random_density_matrix(dims, backend=backend) unitary = random_unitary(dims, backend=backend) - network = QuantumNetwork(unitary, (dims, dims), pure=True, backend=backend) + network = QuantumChannel.from_operator( + unitary, (dims, dims), pure=True, backend=backend, inverse=True + ) applied = network.apply(state) target = unitary @ state @ backend.np.conj(unitary).T @@ -283,11 +426,95 @@ def test_non_hermitian_and_prints(backend): dims = 2**nqubits matrix = random_gaussian_matrix(dims**2, backend=backend) - network = QuantumNetwork(matrix, (dims, dims), pure=False, backend=backend) + network = QuantumNetwork.from_operator( + matrix, (dims, dims), pure=False, backend=backend + ) assert not network.is_hermitian() - assert not network.is_causal() assert not network.is_positive_semidefinite() - assert not network.is_channel() - assert network.__str__() == "J[4 -> 4]" + assert network.__str__() == "J[┍4┑, ┕4┙]" + + +def test_uility_function(): + # _order_tensor2operator should convert + # (a0,a1,b0,b1,...) to (a0,b0,..., a1,b1,...) + old_shape = (0, 10, 1, 11, 2, 12, 3, 13) + test_ls = np.ones(old_shape) + n = len(test_ls.shape) // 2 + + order2op = QuantumNetwork._order_tensor_to_operator(n) + order2tensor = QuantumNetwork._order_operator_to_tensor(n) + + new_shape = test_ls.transpose(order2op).shape + for i in range(n): + assert (new_shape[i] - new_shape[i + n]) == -10 + + assert tuple(test_ls.transpose(order2op).transpose(order2tensor).shape) == old_shape + + +def test_predefined(backend): + tr_ch = TraceOperation(2, backend=backend) + + id_ch = IdentityChannel(2, backend=backend) + id_mat = id_ch.matrix(backend=backend) + + backend.assert_allclose( + id_mat, + backend.cast( + np.array([[1, 0, 0, 1], [0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 1]]), + dtype=id_mat.dtype, + ), + atol=1e-8, + ) + + traced = link_product("ij,j", id_ch, tr_ch, backend=backend) + + backend.assert_allclose( + tr_ch.matrix(backend=backend), traced.matrix(backend=backend), atol=1e-8 + ) + + +def test_default_construction(backend): + vec = np.random.rand(4).reshape([4, 1]) + mat = np.random.rand(16).reshape([2, 2, 2, 2]) + tensor = np.random.rand(16).reshape([4, 4]) + vec = backend.cast(vec, dtype=vec.dtype) + mat = backend.cast(mat, dtype=mat.dtype) + tensor = backend.cast(tensor, dtype=tensor.dtype) + vec = backend.cast(vec, dtype=vec.dtype) + mat = backend.cast(mat, dtype=mat.dtype) + tensor = backend.cast(tensor, dtype=tensor.dtype) + network = QuantumNetwork.from_operator(vec, pure=True, backend=backend) + assert network.partition == (4, 1) + assert network.system_input == (True, False) + comb1 = QuantumComb.from_operator(vec, (4, 1), pure=True, backend=backend) + assert comb1.system_input == (True, False) + comb2 = QuantumComb.from_operator(vec, pure=True, backend=backend) + assert comb2.partition == (4, 1) + assert comb2.system_input == (True, False) + comb3 = QuantumComb.from_operator(mat, pure=False, backend=backend) + assert comb3.partition == (2, 2) + assert comb3.system_input == (True, False) + comb3 = QuantumComb(vec, system_input=(True, True), pure=True, backend=backend) + assert comb3.partition == (4, 1) + assert comb3.system_input == (True, False) + channel1 = QuantumChannel.from_operator(vec, pure=True, backend=backend) + assert channel1.partition == (4, 1) + assert channel1.system_input == (True, False) + channel2 = QuantumChannel( + vec, partition=4, system_input=True, pure=True, backend=backend + ) + assert channel2.partition == (4, 1) + assert channel2.system_input == (True, False) + channel3 = QuantumChannel(vec, partition=4, pure=True, backend=backend) + assert channel3.partition == (1, 4) + assert channel3.system_input == (True, False) + channel4 = QuantumChannel( + vec, partition=4, system_input=False, pure=True, backend=backend + ) + assert channel4.partition == (1, 4) + assert channel4.system_input == (True, False) + channel5 = QuantumChannel(tensor, pure=False, backend=backend) + assert channel5.partition == (2, 2) + assert channel5.system_input == (True, False)