Adds functionality to calculate click cumulants

thewalrus/_torontonian.py

@@ -16,8 +16,70 @@
 """
 import numpy as np
 import numba
-from .quantum import Qmat, Xmat, Amat
-from . import reduction
+#from .quantum import Xmat, Amat
+#from thewalrus import reduction
+
+def Xmat(N):
+    r"""Returns the matrix :math:`X_n = \begin{bmatrix}0 & I_n\\ I_n & 0\end{bmatrix}`
+
+    Args:
+        N (int): positive integer
+
+    Returns:
+        array: :math:`2N\times 2N` array
+    """
+    I = np.identity(N)
+    O = np.zeros_like(I)
+    X = np.block([[O, I], [I, O]])
+    return X
+
+
+def Amat(cov, hbar=2, cov_is_qmat=False):
+    r"""Returns the :math:`A` matrix of the Gaussian state whose hafnian gives the photon number probabilities.
+
+    Args:
+        cov (array): :math:`2N\times 2N` covariance matrix
+        hbar (float): the value of :math:`\hbar` in the commutation
+            relation :math:`[\x,\p]=i\hbar`.
+        cov_is_qmat (bool): if ``True``, it is assumed that ``cov`` is in fact the Q matrix.
+
+    Returns:
+        array: the :math:`A` matrix.
+    """
+    # number of modes
+    N = len(cov) // 2
+    X = Xmat(N)
+
+    # inverse Q matrix
+    if cov_is_qmat:
+        Q = cov
+    else:
+        Q = Qmat_numba(cov, hbar=hbar)
+
+    Qinv = np.linalg.inv(Q)
+
+    # calculate Hamilton's A matrix: A = X.(I-Q^{-1})*
+    A = X @ (np.identity(2 * N) - Qinv).conj()
+    return A
+
+def reduction(A, rpt):
+    r"""Calculates the reduction of an array by a vector of indices.
+
+    This is equivalent to repeating the ith row/column of :math:`A`, :math:`rpt_i` times.
+
+    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:
+        array: the reduction of A by the index vector rpt
+    """
+    rows = [i for sublist in [[idx] * j for idx, j in enumerate(rpt)] for i in sublist]
+
+    if A.ndim == 1:
+        return A[rows]
+
+    return A[:, rows][rows]
 
 
 def tor(A):
@@ -247,7 +309,7 @@ def threshold_detection_prob(mu, cov, det_pattern, hbar=2, atol=1e-10, rtol=1e-1
     if np.allclose(mu, 0, atol=atol, rtol=rtol):
         # no displacement
         n_modes = cov.shape[0] // 2
-        Q = Qmat(cov, hbar)
+        Q = Qmat_numba(cov, hbar)
         O = Xmat(n_modes) @ Amat(cov, hbar=hbar)
         rpt2 = np.concatenate((det_pattern, det_pattern))
         Os = reduction(O, rpt2)

thewalrus/quantum/__init__.py

@@ -107,6 +107,7 @@
     mean_clicks
     variance_clicks
     photon_number_cumulant
+    click_cumulant
 
 Photon number distributions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -175,6 +176,7 @@
     mean_clicks,
     variance_clicks,
     s_ordered_expectation,
+    click_cumulant,
 )
 
 from .photon_number_distributions import (

thewalrus/quantum/means_and_variances.py

@@ -23,6 +23,7 @@
 import numpy as np
 
 from .._hafnian import hafnian, reduction
+from .._torontonian import threshold_detection_prob_displacement
 
 from .conversions import (
     reduced_gaussian,
@@ -408,7 +409,7 @@ def _list_to_freq_dict(words):
     return {i : words.count(i) for i in set(words)}
 
 def photon_number_cumulant(mu, cov, modes, hbar=2):
-    r"""Calculates the cumulant of the modes in the Gaussian state.
+    r"""Calculates the photon-number cumulant of the modes in the Gaussian state.
 
     Args:
         mu (array): length-:math:`2N` means vector in xp-ordering.
@@ -432,3 +433,33 @@ def photon_number_cumulant(mu, cov, modes, hbar=2):
         kappa += term
 
     return kappa
+
+
+def click_cumulant(mu, cov, modes, hbar=2):
+    r"""Calculates the cumulant of the modes in the Gaussian state.
+
+    Args:
+        mu (array): length-:math:`2N` means vector in xp-ordering.
+        cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering.
+        modes (list or array): list of modes.
+        hbar (float): value of hbar in the uncertainty relation.
+
+    Returns:
+        (float): the cumulant
+    """
+
+    modes = list(modes)  # turns modes from array to list if passed in as array
+    num_modes = len(mu) // 2
+    kappa = 0
+    for pi in partition(modes):
+        size = len(pi)
+        term = factorial(size - 1) * (-1) ** (size - 1)
+        for B in pi:
+            B = list(set(B))  # remove repetitions
+            pattern = np.ones([len(B)])
+            mu_red, cov_red = reduced_gaussian(mu, cov, B)
+            summand = threshold_detection_prob_displacement(mu_red, cov_red, pattern, hbar=hbar)
+            term *= summand
+        kappa += term
+
+    return kappa

thewalrus/tests/test_quantum.py

@@ -31,7 +31,7 @@
 )
 
 from thewalrus.random import random_covariance, random_interferometer
-
+from thewalrus._torontonian import threshold_detection_prob_displacement
 from thewalrus.quantum.fock_tensors import _prefactor
 from thewalrus.quantum.photon_number_distributions import _squeezed_state_distribution
 from thewalrus.quantum.adjacency_matrices import _mean_clicks_adj
@@ -75,6 +75,7 @@
     characteristic_function,
     photon_number_moment,
     n_body_marginals,
+    click_cumulant,
 )
 
 @pytest.mark.parametrize("n", [0, 1, 2])
@@ -1813,3 +1814,32 @@ def test_n_body_marginals_too_high_correlation():
         ValueError, match="The order of the correlations is higher than the number of modes"
     ):
         n_body_marginals(mu, cov, 4, 4)
+
+
+
+def test_click_cumulants():
+    """Tests the correctness of the click cumulant function"""
+    M = 3
+    cov = random_covariance(M)
+    mu = np.zeros([2 * M])
+    mu0, cov0 = reduced_gaussian(mu, cov, [0])
+    mu1, cov1 = reduced_gaussian(mu, cov, [1])
+    mu2, cov2 = reduced_gaussian(mu, cov, [2])
+    mu01, cov01 = reduced_gaussian(mu, cov, [0, 1])
+    mu02, cov02 = reduced_gaussian(mu, cov, [0, 2])
+    mu12, cov12 = reduced_gaussian(mu, cov, [1, 2])
+    expected = (
+        threshold_detection_prob_displacement(mu, cov, [1, 1, 1])
+        - threshold_detection_prob_displacement(mu01, cov01, [1, 1])
+        * threshold_detection_prob_displacement(mu2, cov2, [1])
+        - threshold_detection_prob_displacement(mu02, cov02, [1, 1])
+        * threshold_detection_prob_displacement(mu1, cov1, [1])
+        - threshold_detection_prob_displacement(mu12, cov12, [1, 1])
+        * threshold_detection_prob_displacement(mu0, cov0, [1])
+        + 2
+        * threshold_detection_prob_displacement(mu2, cov2, [1])
+        * threshold_detection_prob_displacement(mu1, cov1, [1])
+        * threshold_detection_prob_displacement(mu0, cov0, [1])
+    )
+    obtained = click_cumulant(mu, cov, [0, 1, 2])
+    assert np.allclose(expected, obtained)