Source code for thewalrus._permanent

# Copyright 2019 Xanadu Quantum Technologies Inc.

# 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

#     http://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.
r"""
This submodule provides access to tools for finding the permanent of a matrix. The algorithms implemented
here was first derived in

* Ryser, Herbert John (1963).
  Combinatorial Mathematics, The Carus Mathematical Monographs, Vol. 14, Mathematical Association of America.
* Glynn, David G.
  (2010), "The permanent of a square matrix", European Journal of Combinatorics, 31 (7): 1887-1891.
  <doi:10.1016/j.ejc.2010.01.010>`_
"""
from itertools import chain

import numpy as np
from numba import jit, prange

from scipy.special import factorial

from ._hafnian import hafnian_repeated, find_kept_edges


[docs] def perm(A, method="bbfg"): """Returns the permanent of a matrix using various methods. Args: A (array[float or complex]): a square array. method (string): Set this to ``"ryser"`` to use the `Ryser formula <https://en.wikipedia.org/wiki/Computing_the_permanent#Ryser_formula>`_, or ``"bbfg"`` to use the `BBFG formula <https://en.wikipedia.org/wiki/Computing_the_permanent#Balasubramanian%E2%80%93Bax%E2%80%93Franklin%E2%80%93Glynn_formula>`_. Returns: float or complex: the permanent of matrix ``A`` """ if not isinstance(A, np.ndarray): raise TypeError("Input matrix must be a NumPy array.") matshape = A.shape if matshape[0] != matshape[1]: raise ValueError("Input matrix must be square.") if np.isnan(A).any(): raise ValueError("Input matrix must not contain NaNs.") if matshape[0] == 0: return A.dtype.type(1.0) if matshape[0] == 1: return A[0, 0] if matshape[0] == 2: return A[0, 0] * A[1, 1] + A[0, 1] * A[1, 0] if matshape[0] == 3: return ( A[0, 2] * A[1, 1] * A[2, 0] + A[0, 1] * A[1, 2] * A[2, 0] + A[0, 2] * A[1, 0] * A[2, 1] + A[0, 0] * A[1, 2] * A[2, 1] + A[0, 1] * A[1, 0] * A[2, 2] + A[0, 0] * A[1, 1] * A[2, 2] ) isRyser = bool(method != "bbfg") return perm_ryser(A) if isRyser else perm_bbfg(A)
@jit(nopython=True) def perm_ryser(M): # pragma: no cover """ Returns the permanent of a matrix using the Ryser formula in Gray ordering. The code is an re-implementation from a Python 2 code found in `Permanent code golf <https://codegolf.stackexchange.com/questions/97060/calculate-the-permanent-as-quickly-as-possible>`_ using Numba. Args: M (array) : a square array. Returns: float or complex: the permanent of matrix ``M`` """ n = len(M) if n == 0: return M.dtype.type(1.0) # row_comb keeps the sum of previous subsets. # Every iteration, it removes a term and/or adds a new term # to give the term to add for the next subset row_comb = np.zeros((n), dtype=M.dtype) total = 0 old_grey = 0 sign = +1 binary_power_dict = [2**i for i in range(n)] num_loops = 2**n for k in range(0, num_loops): bin_index = (k + 1) % num_loops reduced = np.prod(row_comb) total += sign * reduced new_grey = bin_index ^ (bin_index // 2) grey_diff = old_grey ^ new_grey grey_diff_index = binary_power_dict.index(grey_diff) new_vector = M[grey_diff_index] direction = (old_grey > new_grey) - (old_grey < new_grey) for i in range(n): row_comb[i] += new_vector[i] * direction sign = -sign old_grey = new_grey return total @jit(nopython=True) def perm_bbfg(M): # pragma: no cover """ Returns the permanent of a matrix using the bbfg formula in Gray ordering The code is a re-implementation from a Python 2 code found in `Permanent code golf <https://codegolf.stackexchange.com/questions/97060/calculate-the-permanent-as-quickly-as-possible>`_ using Numba. Args: M (array) : a square array. Returns: float or complex: the permanent of a matrix ``M`` """ n = len(M) if n == 0: return M.dtype.type(1.0) row_comb = np.sum(M, 0) total = 0 old_gray = 0 sign = +1 binary_power_dict = [2**i for i in range(n)] num_loops = 2 ** (n - 1) for bin_index in range(1, num_loops + 1): reduced = np.prod(row_comb) total += sign * reduced new_gray = bin_index ^ (bin_index // 2) gray_diff = old_gray ^ new_gray gray_diff_index = binary_power_dict.index(gray_diff) new_vector = M[gray_diff_index] direction = 2 * ((old_gray > new_gray) - (old_gray < new_gray)) for i in range(n): row_comb[i] += new_vector[i] * direction sign = -sign old_gray = new_gray return total / num_loops
[docs] def permanent_repeated(A, rpt): r"""Calculates the permanent of matrix :math:`A`, where the ith row/column of :math:`A` is repeated :math:`rpt_i` times. This function constructs the matrix .. math:: B = \begin{bmatrix} 0 & A\\ A^T & 0 \end{bmatrix}, and then calculates :math:`perm(A)=haf(B)`, by calling >>> hafnian_repeated(B, rpt*2, loop=False) Args: A (array): matrix of size [N, N] rpt (Sequence): sequence of N positive integers indicating the corresponding rows/columns of A to be repeated. Returns: int or float or complex: the permanent of matrix ``A`` """ n = A.shape[0] O = np.zeros([n, n]) B = np.vstack([np.hstack([O, A]), np.hstack([A.T, O])]) return hafnian_repeated(B, rpt * 2, loop=False)
[docs] @jit(nopython=True, parallel=True) def brs(A, E): # pragma: no cover r""" Calculates the Bristolian, a matrix function introduced for calculating the threshold detector statistics on measurements of Fock states interfering in linear optical interferometers. See the paper 'Threshold detector statistics of Bosonic states' for more detail (to be published soon) Args: A (array): matrix of size [m, n] E (array): matrix of size [r, n] Returns: int or float or complex: the Bristol of matrices A and E """ m = A.shape[0] steps = 2**m ones = np.ones(m, dtype=np.int8) total = 0 for j in prange(steps): kept_rows = np.where(find_kept_edges(j, ones) != 0)[0] Ay = A[kept_rows, :] plusminus = (-1) ** ((m - len(kept_rows)) % 2) total += plusminus * perm_bbfg(Ay.conj().T @ Ay + E) return total
[docs] @jit(nopython=True, parallel=True) def ubrs(A): # pragma: no cover r""" Calculates the Unitary Bristolian, a matrix function introduced for calculating the threshold detector statistics on measurements of Fock states interfering in lossless linear optical interferometers. See the paper 'Threshold detector statistics of Bosonic states' for more detail (to be published soon) Args: A (array): matrix of size [m, n] Returns: int or float or complex: the Unitary Bristol of matrix A """ m = A.shape[0] steps = 2**m ones = np.ones(m, dtype=np.int8) total = 0 for j in prange(1, steps): kept_rows = np.where(find_kept_edges(j, ones) != 0)[0] Az = A[kept_rows, :] plusminus = (-1) ** ((m - len(kept_rows)) % 2) total += plusminus * perm_bbfg(Az.conj().T @ Az) return total
def fock_prob(n, m, U): r""" Calculates the probability of a an input Fock state, n, scattering to an output Fock state, m, through an interferometer described by matrix U. The matrix U does not need to be unitary, but the total photon number at the input and the output must be equal. Args: n (sequence[int]): length-M list giving the input Fock state occupancy of each mode m (sequence[int]): length-M list giving the output Fock state occupancy of each mode U (array): M x M matrix describing the a linear optical transformation Returns: float: probability of Fock state, n, scattering to m, through an interferometer, U """ if sum(n) != sum(m): raise ValueError("number of input photons must equal number of output photons") in_modes = np.array(list(chain(*[[i] * j for i, j in enumerate(n) if j > 0]))) out_modes = np.array(list(chain(*[[i] * j for i, j in enumerate(m) if j > 0]))) Umn = U[np.ix_(out_modes, in_modes)] n = np.array(n) m = np.array(m) return abs(perm(Umn)) ** 2 / ( np.prod(factorial(n), dtype=np.float64) * np.prod(factorial(m), dtype=np.float64) ) def fock_threshold_prob(n, d, T): r""" Calculates the probability of a an M_in mode input Fock state, n, scattering through an interferometer described by T, being detected by M_out threshold detectors, with outcome given by M_out-length list, d. T is an M_out x M_in matrix. It does not need to be unitary but M_out <= M_in. Args: n (sequence[int]): length-M_in list giving the input Fock state occupancy of each mode d (sequence[int]): length-M_out list giving the outputs of threshold detectors T (array): M_out x M_in matrix describing the a linear optical transformation, M_out <= M_in Returns: float: probability of Fock state, n, scattering through an interferometer, T, to give threshold detector outcome, d """ n = np.array(n) d = np.array(d) if len(n) != T.shape[1]: raise ValueError("length of n must matrix number of input modes of T") if len(d) != T.shape[0]: raise ValueError("length of d must match number of output modes of T") if T.shape[0] > T.shape[1]: raise ValueError("number of output modes cannot be larger than number of input modes") fac_prod = np.prod(factorial(n), dtype=np.float64) in_modes = np.array(list(chain(*[[i] * j for i, j in enumerate(n) if j > 0]))) C = np.where(d > 0)[0] A = T[np.ix_(C, in_modes)] # if matrix is unitary, use the Unitary Bristolian E = np.eye(T.shape[1]) - T.conj().T @ T if np.allclose(E, np.zeros((T.shape[1], T.shape[1]))): U_dn = T[np.ix_(C, in_modes)] return ubrs(U_dn).real / fac_prod E_n = E[np.ix_(in_modes, in_modes)] return brs(A, E_n).real / fac_prod

Contents

Getting started

Background

The Walrus API