code: def fock_prob(cov, mu, event, hbar=2.0):
      r"""Returns the probability of detection of a particular PNR detection event.
  
      For more details, see:
  
      * Kruse, R., Hamilton, C. S., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I.
        "A detailed study of Gaussian Boson Sampling." `arXiv:1801.07488. (2018).
        <https://arxiv.org/abs/1801.07488>`_
  
      * Hamilton, C. S., Kruse, R., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I.
        "Gaussian boson sampling." `Physical review letters, 119(17), 170501. (2017).
        <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.170501>`_
  
      Args:
          cov (array): :math:`2N\times 2N` covariance matrix
          mu (array): length-:math:`2N` means vector
          event (array): length-:math:`N` array of non-negative integers representing the
              PNR detection event of the multi-mode system.
          hbar (float): (default 2) the value of :math:`\hbar` in the commutation
              relation :math:`[\x,\p]=i\hbar`.
  
      Returns:
          float: probability of detecting the event
      """
      # number of modes
      N = len(mu) // 2
      I = np.identity(N)
  
      # mean displacement of each mode
      alpha = (mu[:N] + 1j * mu[N:]) / math.sqrt(2 * hbar)
      # the expectation values (<a_1>, <a_2>,...,<a_N>, <a^\dagger_1>, ..., <a^\dagger_N>)
      beta = np.concatenate([alpha, alpha.conj()])
  
      x = cov[:N, :N] * 2 / hbar
      xp = cov[:N, N:] * 2 / hbar
      p = cov[N:, N:] * 2 / hbar
      # the (Hermitian) matrix elements <a_i^\dagger a_j>
      aidaj = (x + p + 1j * (xp - xp.T) - 2 * I) / 4
      # the (symmetric) matrix elements <a_i a_j>
      aiaj = (x - p + 1j * (xp + xp.T)) / 4
  
      # calculate the covariance matrix sigma_Q appearing in the Q function:
      # Q(alpha) = exp[-(alpha-beta).sigma_Q^{-1}.(alpha-beta)/2]/|sigma_Q|
      Q = np.block(
          [[aidaj, aiaj.conj()], [aiaj, aidaj.conj()]]
      ) + np.identity(2 * N)
  
      # inverse Q matrix
      Qinv = np.linalg.inv(Q)
      # 1/sqrt(|Q|)
      sqrt_Qdet = 1 / math.sqrt(np.linalg.det(Q).real)
  
      prefactor = cmath.exp(-beta @ Qinv @ beta.conj() / 2)
  
      if np.all(np.array(event) == 0):
          # all PNRs detect the vacuum state
          return (prefactor * sqrt_Qdet).real / np.prod(fac(event))
  
      # the matrix X_n = [[0, I_n], [I_n, 0]]
      O = np.zeros_like(I)
      X = np.block([[O, I], [I, O]])
  
      gamma = X @ Qinv.conj() @ beta
  
      # For each mode, repeat the mode number event[i] times
      ind = [
          i
          for sublist in [[idx] * j for idx, j in enumerate(event)]
          for i in sublist
      ]
      # extend the indices for xp-ordering of the Gaussian state
      ind += [i + N for i in ind]
  
      if np.linalg.norm(beta) < tolerance:
          # state has no displacement
          part = partitions(ind, include_singles=False)
      else:
          part = partitions(ind, include_singles=True)
  
      # calculate Hamilton's A matrix: A = X.(I-Q^{-1})*
      A = X @ (np.identity(2 * N) - Qinv).conj()
      summation = np.sum(
          [
              np.prod([gamma[i[0]] if len(i) == 1 else A[i] for i in p])
              for p in part
          ]
      )
  
      return (prefactor * sqrt_Qdet * summation).real / np.prod(
          fac(event)
      )
  
imports: import math
  import cmath
  import numpy
  from scipy.special import factorial
  import pennylane
  from pennylane.ops import Identity
  from pennylane import Device
  from .._version import __version__
label_desc: Use image node_id 2 for calling a global function with example usage: fock_prob(cov, mu, event, hbar) and returns: unknown, unknown