cirq/linalg/predicates.py

# Copyright 2018 The Cirq Developers
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Utility methods for checking properties of matrices."""
from typing import cast, List, Optional, Sequence, Union, Tuple

import numpy as np

from cirq.linalg import tolerance, transformations
from cirq import value


def is_diagonal(matrix: np.ndarray, *, atol: float = 1e-8) -> np.bool_:
    """Determines if a matrix is a approximately diagonal.

    A matrix is diagonal if i!=j implies m[i,j]==0.

    Args:
        matrix: The matrix to check.
        atol: The per-matrix-entry absolute tolerance on equality.

    Returns:
        Whether the matrix is diagonal within the given tolerance.
    """
    matrix = np.copy(matrix)
    for i in range(min(matrix.shape)):
        matrix[i, i] = 0
    return tolerance.all_near_zero(matrix, atol=atol)


def is_hermitian(matrix: np.ndarray, *, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
    """Determines if a matrix is approximately Hermitian.

    A matrix is Hermitian if it's square and equal to its adjoint.

    Args:
        matrix: The matrix to check.
        rtol: The per-matrix-entry relative tolerance on equality.
        atol: The per-matrix-entry absolute tolerance on equality.

    Returns:
        Whether the matrix is Hermitian within the given tolerance.
    """
    return matrix.shape[0] == matrix.shape[1] and np.allclose(
        matrix, np.conj(matrix.T), rtol=rtol, atol=atol
    )


def is_orthogonal(matrix: np.ndarray, *, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
    """Determines if a matrix is approximately orthogonal.

    A matrix is orthogonal if it's square and real and its transpose is its
    inverse.

    Args:
        matrix: The matrix to check.
        rtol: The per-matrix-entry relative tolerance on equality.
        atol: The per-matrix-entry absolute tolerance on equality.

    Returns:
        Whether the matrix is orthogonal within the given tolerance.
    """
    return (
        matrix.shape[0] == matrix.shape[1]
        and np.all(np.imag(matrix) == 0).item()
        and np.allclose(matrix.dot(matrix.T), np.eye(matrix.shape[0]), rtol=rtol, atol=atol)
    )


def is_special_orthogonal(matrix: np.ndarray, *, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
    """Determines if a matrix is approximately special orthogonal.

    A matrix is special orthogonal if it is square and real and its transpose
    is its inverse and its determinant is one.

    Args:
        matrix: The matrix to check.
        rtol: The per-matrix-entry relative tolerance on equality.
        atol: The per-matrix-entry absolute tolerance on equality.

    Returns:
        Whether the matrix is special orthogonal within the given tolerance.
    """
    return is_orthogonal(matrix, rtol=rtol, atol=atol) and (
        matrix.shape[0] == 0 or np.allclose(np.linalg.det(matrix), 1, rtol=rtol, atol=atol)
    )


def is_unitary(matrix: np.ndarray, *, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
    """Determines if a matrix is approximately unitary.

    A matrix is unitary if it's square and its adjoint is its inverse.

    Args:
        matrix: The matrix to check.
        rtol: The per-matrix-entry relative tolerance on equality.
        atol: The per-matrix-entry absolute tolerance on equality.

    Returns:
        Whether the matrix is unitary within the given tolerance.
    """
    return matrix.shape[0] == matrix.shape[1] and np.allclose(
        matrix.dot(np.conj(matrix.T)), np.eye(matrix.shape[0]), rtol=rtol, atol=atol
    )


def is_special_unitary(matrix: np.ndarray, *, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
    """Determines if a matrix is approximately unitary with unit determinant.

    A matrix is special-unitary if it is square and its adjoint is its inverse
    and its determinant is one.

    Args:
        matrix: The matrix to check.
        rtol: The per-matrix-entry relative tolerance on equality.
        atol: The per-matrix-entry absolute tolerance on equality.
    Returns:
        Whether the matrix is unitary with unit determinant within the given
        tolerance.
    """
    return is_unitary(matrix, rtol=rtol, atol=atol) and (
        matrix.shape[0] == 0 or np.allclose(np.linalg.det(matrix), 1, rtol=rtol, atol=atol)
    )


def is_normal(matrix: np.ndarray, *, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
    """Determines if a matrix is approximately normal.

    A matrix is normal if it's square and commutes with its adjoint.

    Args:
        matrix: The matrix to check.
        rtol: The per-matrix-entry relative tolerance on equality.
        atol: The per-matrix-entry absolute tolerance on equality.

    Returns:
        Whether the matrix is normal within the given tolerance.
    """
    return matrix_commutes(matrix, matrix.T.conj(), rtol=rtol, atol=atol)


def is_cptp(*, kraus_ops: Sequence[np.ndarray], rtol: float = 1e-5, atol: float = 1e-8):
    """Determines if a channel is completely positive trace preserving (CPTP).

    A channel composed of Kraus operators K[0:n] is a CPTP map if the sum of
    the products `adjoint(K[i]) * K[i])` is equal to 1.

    Args:
        kraus_ops: The Kraus operators of the channel to check.
        rtol: The relative tolerance on equality.
        atol: The absolute tolerance on equality.
    """
    sum_ndarray = cast(np.ndarray, sum(matrix.T.conj() @ matrix for matrix in kraus_ops))
    return np.allclose(sum_ndarray, np.eye(*sum_ndarray.shape), rtol=rtol, atol=atol)


def matrix_commutes(
    m1: np.ndarray, m2: np.ndarray, *, rtol: float = 1e-5, atol: float = 1e-8
) -> bool:
    """Determines if two matrices approximately commute.

    Two matrices A and B commute if they are square and have the same size and
    AB = BA.

    Args:
        m1: One of the matrices.
        m2: The other matrix.
        rtol: The per-matrix-entry relative tolerance on equality.
        atol: The per-matrix-entry absolute tolerance on equality.

    Returns:
        Whether the two matrices have compatible sizes and a commutator equal
        to zero within tolerance.
    """
    return (
        m1.shape[0] == m1.shape[1]
        and m1.shape == m2.shape
        and np.allclose(m1.dot(m2), m2.dot(m1), rtol=rtol, atol=atol)
    )


def allclose_up_to_global_phase(
    a: np.ndarray,
    b: np.ndarray,
    *,
    rtol: float = 1.0e-5,
    atol: float = 1.0e-8,
    equal_nan: bool = False,
) -> bool:
    """Determines if a ~= b * exp(i t) for some t.

    Args:
        a: A numpy array.
        b: Another numpy array.
        rtol: Relative error tolerance.
        atol: Absolute error tolerance.
        equal_nan: Whether or not NaN entries should be considered equal to
            other NaN entries.
    """

    if a.shape != b.shape:
        return False
    a, b = transformations.match_global_phase(a, b)

    # Should now be equivalent.
    return np.allclose(a=a, b=b, rtol=rtol, atol=atol, equal_nan=equal_nan)


# TODO(#3388) Add documentation for Raises.
# pylint: disable=missing-raises-doc
def slice_for_qubits_equal_to(
    target_qubit_axes: Sequence[int],
    little_endian_qureg_value: int = 0,
    *,  # Forces keyword args.
    big_endian_qureg_value: int = 0,
    num_qubits: Optional[int] = None,
    qid_shape: Optional[Tuple[int, ...]] = None,
) -> Tuple[Union[slice, int, 'ellipsis'], ...]:
    """Returns an index corresponding to a desired subset of an np.ndarray.

    It is assumed that the np.ndarray's shape is of the form (2, 2, 2, ..., 2).

    Example:
        ```python
        # A '4 qubit' tensor with values from 0 to 15.
        r = np.array(range(16)).reshape((2,) * 4)

        # We want to index into the subset where qubit #1 and qubit #3 are ON.
        s = cirq.slice_for_qubits_equal_to([1, 3], 0b11)
        print(s)
        # (slice(None, None, None), 1, slice(None, None, None), 1, Ellipsis)

        # Get that subset. It corresponds to numbers of the form 0b*1*1.
        # where here '*' indicates any possible value.
        print(r[s])
        # [[ 5  7]
        #  [13 15]]
        ```

    Args:
        target_qubit_axes: The qubits that are specified by the index bits. All
            other axes of the slice are unconstrained.
        little_endian_qureg_value: An integer whose bits specify what value is
            desired for of the target qubits. The integer is little endian
            w.r.t. the target qubit axes, meaning the low bit of the integer
            determines the desired value of the first targeted qubit, and so
            forth with the k'th targeted qubit's value set to
            bool(qureg_value & (1 << k)).
        big_endian_qureg_value: Same as `little_endian_qureg_value` but big
            endian w.r.t. to target qubit axes, meaning the low bit of the
            integer dertemines the desired value of the last target qubit, and
            so forth.  Specify exactly one of the `*_qureg_value` arguments.
        num_qubits: If specified the slices will extend all the way up to
            this number of qubits, otherwise if it is None, the final element
            return will be Ellipsis. Optional and defaults to using Ellipsis.
        qid_shape: The qid shape of the state vector being sliced.  Specify this
            instead of `num_qubits` when using qids with dimension != 2.  The
            qureg value is interpreted to store digits with corresponding bases
            packed into an int.

    Returns:
        An index object that will slice out a mutable view of the desired subset
        of a tensor.
    """
    qid_shape_specified = qid_shape is not None
    if qid_shape is not None or num_qubits is not None:
        if num_qubits is None:
            num_qubits = len(cast(Tuple[int, ...], qid_shape))
        elif qid_shape is None:
            qid_shape = (2,) * num_qubits
        if num_qubits != len(cast(Tuple[int, ...], qid_shape)):
            raise ValueError('len(qid_shape) != num_qubits')
    if little_endian_qureg_value and big_endian_qureg_value:
        raise ValueError(
            'Specify exactly one of the arguments little_endian_qureg_value '
            'or big_endian_qureg_value.'
        )
    out_size_specified = num_qubits is not None
    out_size = (
        cast(int, num_qubits) if out_size_specified else max(target_qubit_axes, default=-1) + 1
    )
    result = cast(List[Union[slice, int, 'ellipsis']], [slice(None)] * out_size)
    if not out_size_specified:
        result.append(Ellipsis)
    if qid_shape is None:
        qid_shape = (2,) * out_size
    target_shape = tuple(qid_shape[i] for i in target_qubit_axes)
    if big_endian_qureg_value:
        digits = value.big_endian_int_to_digits(big_endian_qureg_value, base=target_shape)
    else:
        if little_endian_qureg_value < 0 and not qid_shape_specified:
            # Allow negative binary numbers
            little_endian_qureg_value &= (1 << len(target_shape)) - 1
        digits = value.big_endian_int_to_digits(little_endian_qureg_value, base=target_shape[::-1])[
            ::-1
        ]
    for axis, digit in zip(target_qubit_axes, digits):
        result[axis] = digit
    return tuple(result)


# pylint: enable=missing-raises-doc