Tidy up the docs, fix typos and formatting issues

docs/code/libhaf.rst

@@ -77,17 +77,18 @@ The following functions are intended as the main interface to the C++ Hafnian li
 
 .. rst-class:: longtable docutils
 
-======================================    ==============================================
-:cpp:func:`hafnian::hafnian_recursive`    Returns the hafnian of a matrix using the recursive algorithm described in *Counting perfect matchings as fast as Ryser* :cite:`bjorklund2012counting`.
-:cpp:func:`hafnian::hafnian`              Returns the hafnian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__.
-:cpp:func:`hafnian::loop_hafnian`         Returns the loop hafnian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__.
-:cpp:func:`hafnian::hafnian_rpt`          Returns the hafnian of a matrix with repeated rows and columns using the algorithm described in *From moments of sum to moments of product*, `doi:10.1016/j.jmva.2007.01.013 <https://dx.doi.org/10.1016/j.jmva.2007.01.013>`__.
-:cpp:func:`hafnian::hafnian_approx`       Returns the approximate hafnian of a matrix with non-negative entries by sampling over determinants. The higher the number of samples, the better the accuracy.
-:cpp:func:`hafnian::torontonian`          Returns the Torontonian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__.
-:cpp:func:`hafnian::torontonian_fsum`     Returns the torontonian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__, with increased accuracy via the ``fsum`` summation algorithm.
-:cpp:func:`hafnian::permanent`            Returns the permanent of a matrix using Ryser's algorithm with Gray code ordering.
-:cpp:func:`hafnian::perm_fsum`            Returns the permanent of a matrix using Ryser's algorithm with Gray code ordering, with increased accuracy via the ``fsum`` summation algorithm.
-======================================    ==============================================
+=================================================            ==============================================
+:cpp:func:`hafnian::hafnian_recursive`                       Returns the hafnian of a matrix using the recursive algorithm described in *Counting perfect matchings as fast as Ryser* :cite:`bjorklund2012counting`.
+:cpp:func:`hafnian::hafnian`                                 Returns the hafnian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__.
+:cpp:func:`hafnian::loop_hafnian`                            Returns the loop hafnian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__.
+:cpp:func:`hafnian::hafnian_rpt`                             Returns the hafnian of a matrix with repeated rows and columns using the algorithm described in *From moments of sum to moments of product*, `doi:10.1016/j.jmva.2007.01.013 <https://dx.doi.org/10.1016/j.jmva.2007.01.013>`__.
+:cpp:func:`hafnian::hafnian_approx`                          Returns the approximate hafnian of a matrix with non-negative entries by sampling over determinants. The higher the number of samples, the better the accuracy.
+:cpp:func:`hafnian::torontonian`                             Returns the Torontonian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__.
+:cpp:func:`hafnian::torontonian_fsum`                        Returns the torontonian of a matrix using the algorithm described in *A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer*, `arxiv:1805.12498 <https://arxiv.org/abs/1805.12498>`__, with increased accuracy via the ``fsum`` summation algorithm.
+:cpp:func:`hafnian::permanent`                               Returns the permanent of a matrix using Ryser's algorithm with Gray code ordering.
+:cpp:func:`hafnian::perm_fsum`                               Returns the permanent of a matrix using Ryser's algorithm with Gray code ordering, with increased accuracy via the ``fsum`` summation algorithm.
+:cpp:func:`hafnian::hermite_multidimensional_cpp`            Returns photon number statistics of a Gaussian state for a given covariance matrix as described in *Multidimensional Hermite polynomials and photon distribution for polymode mixed light* `arxiv:9308033 <https://arxiv.org/abs/hep-th/9308033>`__.
+=================================================            ==============================================
 
 
 API

docs/conf.py

@@ -94,7 +94,7 @@ class TypeMock(type):
     # TIP: if using the sphinx-bootstrap-theme, you need
     # "treeViewIsBootstrap": True,
     "exhaleExecutesDoxygen": True,
-    "exhaleDoxygenStdin":    "INPUT = ../src/stdafx.h ../src/hafnian.hpp ../src/version.hpp ../src/eigenvalue_hafnian.hpp ../src/hafnian_approx.hpp ../src/recursive_hafnian.hpp ../src/repeated_hafnian.hpp  ../src/torontonian.hpp ../src/permanent.hpp",
+    "exhaleDoxygenStdin":    "INPUT = ../src/stdafx.h ../src/hafnian.hpp ../src/version.hpp ../src/eigenvalue_hafnian.hpp ../src/hafnian_approx.hpp ../src/recursive_hafnian.hpp ../src/repeated_hafnian.hpp  ../src/torontonian.hpp ../src/permanent.hpp ../src/hermite_multidimensional.hpp",
     # "exhaleUseDoxyfile": True
 }
 

hafnian/__init__.py

@@ -47,6 +47,10 @@
     matrices with non-negative elements. This is done by sampling determinants;
     the larger the number of samples taken, the higher the accuracy.
 
+Batched hafnian algorithm
+    An algorithm that allows to calculate the hafnians of all the reduction of
+    a given matrix up to the cutoff (resolution) provided. This algorithm uses
+    internaly the multidimensional Hermite polynomials.
 
 Python wrappers
 ---------------
@@ -55,10 +59,10 @@
     hafnian
     hafnian_repeated
     hafnian_batched
-    hermite_multidimensional
     tor
     perm
     permanent_repeated
+    hermite_multidimensional
     reduction
     version
 
@@ -69,6 +73,11 @@
 
 import numpy as np
 
+if platform.system() == "Windows":  # pragma: no cover
+    extra_dll_dir = os.path.join(os.path.dirname(__file__), ".libs")
+    if os.path.isdir(extra_dll_dir):
+        os.environ["PATH"] += os.pathsep + extra_dll_dir
+
 from ._hafnian import (
     haf_complex,
     haf_int,
@@ -84,20 +93,15 @@
 from ._torontonian import tor
 from ._version import __version__
 
-if platform.system() == "Windows":  # pragma: no cover
-    extra_dll_dir = os.path.join(os.path.dirname(__file__), ".libs")
-    if os.path.isdir(extra_dll_dir):
-        os.environ["PATH"] += os.pathsep + extra_dll_dir
-
 
 __all__ = [
     "hafnian",
     "hafnian_repeated",
+    "hafnian_batched",
     "tor",
     "perm",
     "permanent_repeated",
     "reduction",
-    "hafnian_batched",
     "hermite_multidimensional",
     "version",
 ]

hafnian/_hafnian.py

@@ -174,10 +174,7 @@ def hafnian(
 
 
 def hafnian_repeated(A, rpt, mu=None, loop=False, tol=1e-12):
-    r"""
-    THIS NEEDS TO BE WRITTEN
-
-    Returns the hafnian of matrix with repeated rows/columns.
+    r"""Returns the hafnian of matrix with repeated rows/columns.
 
     The :func:`reduction` function may be used to show the resulting matrix
     with repeated rows and columns as per ``rpt``.

hafnian/_hermite_multidimensional.py

@@ -22,61 +22,73 @@
 
 
 def return_prod(C, index):
-    r"""Given an array :math:`C_{i,j}` and an array or list of indices :math:`index = [i_1,i_2,i_3,\dots,i_n] `, returns :math:`prod_{k=1}^n C_{1,i_1}`.
+    r"""Given an array :math:`C_{i,j}` and an array or list of indices
+    :math:`index = [i_1,i_2,i_3,\dots,i_n] `, returns :math:`prod_{k=1}^n C_{k,i_k}`.
 
     Args:
         C (array): An array
         index (array): A set of indices
 
     Returns:
-        complex: The product of the array elements determines by index
+        complex: the product of the array elements determined by index
     """
     return np.prod([C[mode, val] for mode, val in enumerate(index)])
 
 
-def expansion_coeff(alpha, resolution, renorm=True):
-    r"""Returns the (quasi) geometric series as a vector with components :math:`alpha^i/sqrt(i!)` for :math:`0 \leq i < resolution`.
-
-    If renorm is false it omits the division by factorials
+def expansion_coeff(alpha, cutoff, renorm=True):
+    r"""Returns the (quasi) geometric series as a vector with components
+    :math:`\alpha^i/\sqrt{i!}` for :math:`0 \leq i < \texttt{cutoff}`.
 
     Args:
-        alpha (complex): Ratio of the geometric series
-        resoluton (int): Cutoff of the geometric series
-        renorm (bool): Decides whether to normalize by the factorials
+        alpha (complex): ratio of the geometric series
+        cutoff (int): cutoff truncation of the geometric series
+        renorm (bool): if ``False``, the components are not normalized
+            by the square-root factorials
 
     Returns:
-        array: The power series
+        array: the (quasi) geometric series
     """
-    vals = np.empty([resolution], dtype=type(alpha))
+    vals = np.empty([cutoff], dtype=type(alpha))
     vals[0] = 1.0
     if renorm:
-        for i in range(1, resolution):
+        for i in range(1, cutoff):
             vals[i] = vals[i - 1] * alpha / np.sqrt(i)
     else:
-        for i in range(1, resolution):
+        for i in range(1, cutoff):
             vals[i] = vals[i - 1] * alpha
     return vals
 
 
-def hermite_multidimensional(R, resolution, y=None, renorm=False, make_tensor=True):
+def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True):
     r"""Returns the multidimensional Hermite polynomials :math:`H_k^{(R)}(y)`.
 
-    Here :math:`R` is an :math:n \times n: square matrix,
-    :math:`y` is an :math:`n` dimensional vector. The polynomials are parametrized by the multi-index :math:`k=(k_0,k_1,\ldots,k_{n-1})
-    and are calculated for all values :math:`0 \leq k_j < \text{resolution}`, thus a tensor of dimensions :math:`\text{resolution}^n` is returned.
-    This tensor can either be flattened into a vector or returned as an actual tensor with :math:`n` indices.
-    Note that is R = np.array([[1.0+0.0j]]) then :math:`H_k^{(R)}(y)` are precisely the well known **probabilists' Hermite polynomials** :math:`He_k(y)`:,
-    and if R = 2*np.array([[1.0+0.0j]]) then :math:`H_k^{(R)}(y)` are precisely the well known **physicists' Hermite polynomials** :math:`H_k(y)`:.
+    Here :math:`R` is an :math:`n \times n` square matrix, and
+    :math:`y` is an :math:`n` dimensional vector. The polynomials are
+    parametrized by the multi-index :math:`k=(k_0,k_1,\ldots,k_{n-1})`,
+    and are calculated for all values :math:`0 \leq k_j < \text{cutoff}`,
+    thus a tensor of dimensions :math:`\text{cutoff}^n` is returned.
+
+    This tensor can either be flattened into a vector or returned as an actual
+    tensor with :math:`n` indices.
+
+    .. note::
+
+        Note that if :math:`R = (1)` then :math:`H_k^{(R)}(y)`
+        are precisely the well known **probabilists' Hermite polynomials** :math:`He_k(y)`,
+        and if :math:`R = (2)` then :math:`H_k^{(R)}(y)` are precisely the well known
+        **physicists' Hermite polynomials** :math:`H_k(y)`.
 
     Args:
-        R (array): Square matrix parametrizing the Hermite polynomial family
-        resolution (int): Maximum size of the subindices in the Hermite polynomial
-        y (array): Vector for the argument of the Hermite polynomial
-        renorm (bool): If True returns :math:`H_k^{(R)}(y)/\prod(\prod_i k_i!)`
-        make_tensor: If False returns a flattened one dimensional array instead of a tensor with the values of the polynomial
+        R (array): square matrix parametrizing the Hermite polynomial family
+        cutoff (int): maximum size of the subindices in the Hermite polynomial
+        y (array): vector argument of the Hermite polynomial
+        renorm (bool): If ``True``, normalizes the returned multidimensional Hermite
+            polynomials such that :math:`H_k^{(R)}(y)/\prod(\prod_i k_i!)`
+        make_tensor: If ``False``, returns a flattened one dimensional array
+            containing the values of the polynomial
 
     Returns:
-        (array): The multidimensional Hermite polynomials.
+        (array): the multidimensional Hermite polynomials
     """
     input_validation(R)
     n, _ = R.shape
@@ -87,59 +99,78 @@ def hermite_multidimensional(R, resolution, y=None, renorm=False, make_tensor=Tr
     if m != n:
         raise ValueError("The matrix R and vector y have incompatible dimensions")
 
-    values = np.array(hm(R, y, resolution, ren=renorm))
+    values = np.array(hm(R, y, cutoff, ren=renorm))
 
     if make_tensor:
-        shape = resolution * np.ones([n], dtype=int)
+        shape = cutoff * np.ones([n], dtype=int)
         values = np.reshape(values, shape)
 
     return values
 
 
-def hafnian_batched(A, resolution, mu=None, tol=1e-12, renorm=False, make_tensor=True):
-    r"""Returns the haf(reduction(A, k)) where k is a vector of (non-negative) integers with the same dimensions as the square matrix A
-    :math:`k = (k_0,k_1,\ldots,k_{n-1})` and where :math:`0 \leq k_j < \text{resolution}`.
+def hafnian_batched(A, cutoff, mu=None, tol=1e-12, renorm=False, make_tensor=True):
+    r"""Calculates the hafnian of :func:`reduction(A, k) <hafnian.reduction>`
+    for all possible values of vector ``k`` below the specified cutoff.
+
+    Here,
+
+    * :math:`A` is am :math:`n\times n` square matrix
+    * :math:`k` is a vector of (non-negative) integers with the same dimensions as :math:`A`,
+      i.e., :math:`k = (k_0,k_1,\ldots,k_{n-1})`, and where :math:`0 \leq k_j < \texttt{cutoff}`.
 
-    If mu is not None then it instead calculates lhaf(fill_diagonal(reduction(A, k),reduction(mu, k))),
-    this calculation can only be performed if the matrix A has an inverse.
+    The function :func:`~.hafnian_repeated` can be used to calculate the reduced hafnian
+    for a *specific* value of :math:`k`; see the documentation for more information.
+
+    .. note::
+
+        If ``mu`` is not ``None``, this function instead returns
+        ``hafnian(np.fill_diagonal(reduction(A, k), reduction(mu, k)), loop=True)``.
+        This calculation can only be performed if the matrix :math:`A` is invertible.
 
     Args:
-        R (array): Square matrix parametrizing
-        resolution (int): Maximum size of the subindices in the Hermite polynomial
-        y (array): Vector for the argument of the Hermite polynomial
-        renorm (bool): If True returns :math:`haf(reduction(A, k))/\prod(\prod_i k_i!)` or lhaf(fill_diagonal(reduction(A, k),reduction(mu, k))) is mu is not None
-        make_tensor: If False returns a flattened one dimensional array instead of a tensor with the values of the polynomial
+        A (array): a square, symmetric :math:`N\times N` array.
+        cutoff (int): maximum size of the subindices in the Hermite polynomial
+        mu (array): a vector of length :math:`N` representing the vector of means/displacement
+        renorm (bool): If ``True``, the returned hafnians are *normalized*, that is,
+            :math:`haf(reduction(A, k))/\prod(\prod_i k_i!)`
+            (or :math:`lhaf(fill\_diagonal(reduction(A, k),reduction(mu, k)))` if
+            ``mu`` is not None)
+        make_tensor: If ``False``, returns a flattened one dimensional array instead
+            of a tensor with the values of the hafnians.
 
     Returns:
-        (array): The values of the hafnians.
+        (array): the values of the hafnians for each value of :math:`k` up to the cutoff
     """
     # pylint: disable=too-many-return-statements,too-many-branches,too-many-arguments
     input_validation(A, tol=tol)
     n, _ = A.shape
+
     if not np.allclose(A, np.zeros([n, n])):
         if mu is not None:
             try:
                 yi = np.linalg.solve(A, mu)
             except np.linalg.LinAlgError:
                 raise ValueError("The matrix does not have an inverse")
             return hermite_multidimensional(
-                -A, resolution, y=-yi, renorm=renorm, make_tensor=make_tensor
+                -A, cutoff, y=-yi, renorm=renorm, make_tensor=make_tensor
             )
         yi = np.zeros([n], dtype=complex)
         return hermite_multidimensional(
-            -A, resolution, y=-yi, renorm=renorm, make_tensor=make_tensor
+            -A, cutoff, y=-yi, renorm=renorm, make_tensor=make_tensor
         )
     # Note the minus signs in the arguments. Those are intentional and are due to the fact that Dodonov et al. in PRA 50, 813 (1994) use (p,q) ordering instead of (q,p) ordering
 
     if mu is None:
-        tensor = np.zeros([resolution ** n], dtype=complex)
+        tensor = np.zeros([cutoff ** n], dtype=complex)
         tensor[0] = 1.0
     else:
-        index = resolution * np.ones([n], dtype=int)
+        index = cutoff * np.ones([n], dtype=int)
         tensor = np.empty(index, dtype=complex)
-        prim = np.array([expansion_coeff(alpha, resolution, renorm=renorm) for alpha in mu])
-        for i in product(range(resolution), repeat=n):
+        prim = np.array([expansion_coeff(alpha, cutoff, renorm=renorm) for alpha in mu])
+        for i in product(range(cutoff), repeat=n):
             tensor[i] = return_prod(prim, i)
+
     if make_tensor:
         return tensor
+
     return tensor.flatten()