Adding batchhafnian code #21

bgupt · 2019-06-28T11:52:58Z

Description of the Change: Code to compute photon number statistics for given Gaussian state with a covariance matrix.

Benefits: Computes probabilities of all photon statistics upto a given cutoff

Possible Drawbacks: None

Related GitHub Issues: None

…googletest and python test modules. Reformattes with black and astyle.

nquesada · 2019-06-30T21:12:30Z

@bgupt could the following functionality be built into the C++ code?

Generating an example:

from hafnian import batchhafnian
import numpy as np

V =  np.array([[ 1.5 , 1.41421356,  0. , 0.],
     [ 1.41421356, 1.5  ,0.  ,0.        ],
     [ 0. ,0. ,1.5 , -1.41421356],
     [ 0. ,0. ,-1.41421356,  1.5       ]])
rs = np.zeros(2*n)

res = 4
indv = res*np.ones(2*n,dtype=int)
A = batchhafnian(V,rs,res)

This is the part that would be great to move into the C++ code

If I understand your code correctly you encode a tensor $T_{i_1,\ldots,i_n}$ of $n$ indices each with local dimensions base into an a linear array $A$ of size $M = \text{base}^n$.

The function find_rep below allows you to find, given the position of an element in the linear array $A$ to which element of the tensor $T_{i_1,\ldots,i_n}$ it corresponds.

Note that $n = 2 \times \ell$ where $\ell$ is the number of modes

def find_rep(val, base, n):
    x = np.empty([n]) # In C++ this would be a vector of size n
    local_val = val
    x[0] = 1
    for i in range(1,n):
        x[i] = x[i-1]*base
    digits = np.zeros([n], dtype = int)
    for i in range(n):
        digits[i] = local_val // x[-i-1] # Careful here this is integer division, not float division, although I think C++ does this by default using / 
        local_val = local_val - digits[i]*x[-i-1]

    return digits

find_rep(85,4,4) # the 85th element of the array with 2 modes (n = 2*2) and with local base size base = 4 corrsponds to the 1,1,1,1 element of the tensor

array([1, 1, 1, 1])

The function below does the following: We know that each element of $A$ is connected to an element of $T_{i_1,\ldots,i_n}$. The following function transform $A$ so that it now corresponds to the element
$T_{i_1,\ldots,i_n} \to S_{i_1,\ldots,i_n} = T_{i_1,\ldots,i_n}/\sqrt{i_!i_2! \ldots i_n!}$

def renorm(tn, m, res):
    facts = np.ones([res])
    for i in range(1,res):
        facts[i] = i*facts[i-1]
    for i in range(1,res):
        facts[i] = 1.0/np.sqrt(facts[i])
    for i in range(res**m):
        digits = find_rep(i, res, m)
        pref = 1.0
        for j in digits:
            pref = pref*facts[j]
        tn[i] = tn[i]*pref 
    return tn

B = np.copy(A)
renorm(B, 2*n, res);


For example, the last element of $A$ which corresponds to the 3,3,3,3 element of T is


```python
A[-1]

(4.499999999999993+0j)

B[-1]

(0.12499999999999988+0j)

A[-1]/B[-1] # which is also sqrt(3! 3! 3! 3!)

(35.999999999999986+0j)

Would it be possible to make the function batch hafnian in C++ take care of the division by all the factorials, i.e. have it return the array B constructed Above directly?

nquesada · 2019-06-30T21:21:16Z

One final thing is that it would be great if the inclusion of this factorial parameters could be made optional, i.e., add an extra argument to the C++ function...

hafnian/__init__.py

codecov · 2019-07-01T11:56:40Z

Codecov Report

Merging #21 into master will decrease coverage by 0.01%.
The diff coverage is 96.9%.

@@            Coverage Diff            @@
##           master     #21      +/-   ##
=========================================
- Coverage   95.21%   95.2%   -0.02%     
=========================================
  Files           8       9       +1     
  Lines         564     625      +61     
=========================================
+ Hits          537     595      +58     
- Misses         27      30       +3

Impacted Files	Coverage Δ
hafnian/_hafnian.py	`100% <ø> (ø)`	⬆️
hafnian/__init__.py	`100% <100%> (ø)`	⬆️
hafnian/quantum.py	`97.82% <95.65%> (-0.82%)`	⬇️
hafnian/_hermite_multidimensional.py	`98% <98%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0bed206...ee2ca7b. Read the comment docs.

…etest accordingly.

…test when renorm=False. Test for True needs to be written.

nquesada · 2019-07-09T17:04:51Z

src/batchhafnian.hpp

@@ -117,6 +136,10 @@ inline std::vector<T> hermite_multidimensional(std::vector<T> &R_mat, std::vecto
                std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
                ullint prevCoordinate = vec2index(tmpjump, resolution);
                H[nextCoordinate] = H[nextCoordinate] - (static_cast<T>(jumpFrom[ii]-1))*static_cast<T>(R(k,ii))*H[prevCoordinate];
+
+				if (renorm) {


the renormalization need to occur outside the loop for (ullint jj = 0; jj < Hdim-1; jj++) otherwise you break the recursion relation. So you need a second loop for (ullint jj = 0; jj < Hdim-1; jj++) in which do the renormalization after the first loop is completed.

Computed the renormalization factor in a separate array in the same loop and used a separate loop to renormalize if renorm = 1.

nquesada

Hi @bgupt . I left a few comments that should be addressed. In particular because you renomarlize inside the recursion loops you get wrong answers. Also, how hard would it be allow for these calculation to also be done in quad precision?

nquesada · 2019-07-09T17:26:34Z

src/batchhafnian.hpp

+
+				if (renorm) {
+					H[nextCoordinate] = H[nextCoordinate]/renorm_factor(nextPos);
+				}


Also, you should not be calling a function to calculate factorials everytime. Instead just save the square root factorials all the way to resolution and just multiply those that you need for each element of the tensor. For this you will need to write a function that allows you to find the indices correspoding to a given element of the tensor inside the vector H. Such function in python is:

def find_rep(val, base, n): x = np.empty([n]) # In C++ this would be a vector of size n local_val = val x[0] = 1 for i in range(1,n): x[i] = x[i-1]*base digits = np.zeros([n], dtype = int) for i in range(n): digits[i] = local_val // x[-i-1] # Careful here this is integer division, not float division, although I think C++ does this by default using / local_val = local_val - digits[i]*x[-i-1] return digits

Thanks for the suggestion. Implemented that.

nquesada

I am not sure if it is a good idea to use integers to represent factorials.

nquesada · 2019-07-09T17:32:50Z

src/batchhafnian.hpp

+
+}
+
+ullint factorial(int n)


What is the largest factorial that you can calculate in integer precision ullint

I am not too sure about that. But factorial of 21 exceeds the range of unsigned long long int so I suppose the max would be 20.

But I changed the factorial function from ullint to long double. So, it should be better now.

bgupt · 2019-07-10T17:36:53Z

I am not sure if it is a good idea to use integers to represent factorials.

Changed it to long double.

… Removed gradhaf mentions. Changed factorial data type from unsigned long long int to long double. More efficient computation of the renormalization factors.

…on tests.

…to batchhafnian

hafnian/_hafnian.py

hafnian/__init__.py

josh146 · 2019-07-13T03:57:46Z

hafnian/_hafnian.py

-        raise ValueError("Input matrix must be symmetric.")
+    if check_symmetry:
+        if np.linalg.norm(A - A.T) >= tol:
+            raise ValueError("Input matrix must be symmetric.")


How comes this check was added? Is there a use case for attempting to calculate the hafnian of non-symmetric matrices?

That is just to make sure that the input matrix is symmetric. It raises error if not.

Oh, I meant the check_symmetry flag. Is there a reason to turn off this check?

There was because I thought multidimensional hermite polynomials are defined for any matrix. In reality they only make sense for symmetric matrices. Will remove this option

hafnian/_hermite_multidimensional.py

hafnian/tests/test_hermite_multidimensional.py

setup.py

src/hermite_multidimensional.hpp

src/tests/hafnian_unittest.cpp

…es in source code.

…for the exceptions

josh146 · 2019-07-16T08:03:24Z

hafnian/__init__.py

+    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
+


isort probably moved this down, but this will have to be moved back to before the line

from ._hafnian import (

to make sure it imports correctly on windows :)

@bgupt can you look into this?

josh146 · 2019-07-16T08:04:40Z

hafnian/_hafnian.py

+    r"""
+    THIS NEEDS TO BE WRITTEN
+
+    Returns the hafnian of matrix with repeated rows/columns.


Is this a copy/paste mistake 🙂 Probably should undo this change.

It is. Thanks for catching it up!

hafnian/_hermite_multidimensional.py

josh146 · 2019-07-16T08:05:55Z

hafnian/_hermite_multidimensional.py

+        C (array): An array
+        index (array): A set of indices
+    Returns:
+        complex: The product of the array elements determines by index


Suggested change

complex: The product of the array elements determines by index

complex: The product of the array elements determined by index

hafnian/_hermite_multidimensional.py

josh146 · 2019-07-16T08:07:17Z

hafnian/_hermite_multidimensional.py

+def expansion_coeff(alpha, resolution, renorm=True):
+    """
+    Returns the (quasi) geometric series as a vector with components alpha^i/sqrt(i!) for 0 <= i < resolution.
+    If renorm is false it omits the division by factorials


Can probably delete this line since it's already discussed in the args: list.

Suggested change

If renorm is false it omits the division by factorials

hafnian/_hermite_multidimensional.py

josh146 · 2019-07-16T08:08:58Z

hafnian/_hermite_multidimensional.py

+        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
+    Returns:
+        (array): The multidimensional Hermite polynomials.


Note that you'll need blank lines in three places in order for this to render properly in the docs.

After the first line of the docstring. So after "Returns the multidimensional Hermite polynomials :math:H_k^{(R)}(y).", there should be an empty line.

Before "Args:"

Before "Returns"

josh146 · 2019-07-16T08:09:40Z

hafnian/_hermite_multidimensional.py

+        make_tensor: If False returns a flattened one dimensional array instead of a tensor with the values of the polynomial
+    Returns:
+        (array): The values of the hafnians.
+    """


same comment as above re: line breaks. You can build the docs locally (make doc) to ensure the docstrings render correctly :)

josh146 · 2019-07-16T08:10:44Z

hafnian/hafnian.pyx

+    for i in range(n):
+        y_mat.push_back(d[i])
+
+    return hermite_multidimensional_cpp(R_mat, y_mat, resolution, renorm)


I'm a bit confused; if this is the 'batch_hafnian', should we just call the function 'batch_hafnian' instead of hermite_multidimensional?

No. As far as I know this connection has not been made and I think is better to keep them separated. The multidimensional hermite polynomials are a thing that has existed for centuries. The fact that they can be used to calculate batches of hafnians is a new finding...

Co-Authored-By: Josh Izaac <josh146@gmail.com>

…urns and also removed the option check_symmetry from input validation

…to batchhafnian

Adding batchhafnian serial code to the repo along with corresponidng …

11e3cd5

…googletest and python test modules. Reformattes with black and astyle.

bgupt requested review from josh146 and nquesada June 28, 2019 11:53

bgupt added 4 commits June 28, 2019 21:49

Adding batchhafnian files along with googletests and python tests.

aa7ae01

Adding python test for batchhafnian.

f3dcdae

Adding permanent file.

24420a7

Increasing the tolerance for hafnian_approx test.

96971a6

Merge branch 'master' into batchhafnian

bd67897

josh146 reviewed Jul 1, 2019

View reviewed changes

hafnian/__init__.py Outdated Show resolved Hide resolved

josh146 reviewed Jul 1, 2019

View reviewed changes

hafnian/__init__.py Outdated Show resolved Hide resolved

josh146 added 2 commits July 1, 2019 19:40

Apply suggestions from code review

5918566

fix unittest

a758609

bgupt added 2 commits July 8, 2019 22:47

Changing input matrices in batchhafnians to R and y. Modify the googl…

c013920

…etest accordingly.

Adding renormalization to batchhafnian functions. Current tests only …

89abd16

…test when renorm=False. Test for True needs to be written.

nquesada reviewed Jul 9, 2019

View reviewed changes

nquesada requested changes Jul 9, 2019

View reviewed changes

nquesada reviewed Jul 9, 2019

View reviewed changes

bgupt and others added 9 commits July 10, 2019 23:15

Fixed the renormalization in batch hafnian as per Nico's suggestions.…

63059ce

… Removed gradhaf mentions. Changed factorial data type from unsigned long long int to long double. More efficient computation of the renormalization factors.

Minor changes in hafnian.hpp

4c471ea

Adding renormalization to batchhafnian along with googletest and Pyth…

986080d

…on tests.

Merge branch 'batchhafnian' of https://github.com/xanaduai/hafnian in…

82d216a

…to batchhafnian

Renaming functions to streamline python library.

ba966b6

Removes white space

8de965b

Fixed merge

515c4d5

Adding _hermite_multidimensional.py

364ba0c

Merge branch 'batchhafnian' of https://github.com/xanaduai/hafnian in…

72b4e41

…to batchhafnian

bgupt and others added 4 commits July 12, 2019 22:04

Renamed src/batchhafnian.hpp -> src/hermite_multidimensional.hpp

020aaeb

Changes to hafnian.cpp

7bedc7b

Merge branch 'batchhafnian' of https://github.com/xanaduai/hafnian in…

5392b78

…to batchhafnian

Added new functions hafnian_batched and multiple tests

3dcdb70

josh146 requested changes Jul 13, 2019

View reviewed changes

bgupt and others added 5 commits July 15, 2019 11:53

Added assert statements in python tests, reformatting and minor chang…

b89adb2

…es in source code.

The test_coherent_squeezed is still not passing

86006cd

Test pass now linting cleaning

f3a1f5f

Linted and blacked. Also the test go through all the new code except …

76f8884

…for the exceptions

Removed lingering printing commands

355fe38

josh146 reviewed Jul 16, 2019

View reviewed changes

nquesada and others added 15 commits July 16, 2019 09:34

Update hafnian/_hermite_multidimensional.py

4d2d9d7

Co-Authored-By: Josh Izaac <josh146@gmail.com>

Update hafnian/_hermite_multidimensional.py

865af95

Co-Authored-By: Josh Izaac <josh146@gmail.com>

Update hafnian/_hermite_multidimensional.py

a2395e4

Co-Authored-By: Josh Izaac <josh146@gmail.com>

Applied code review for docstrings adding returns before args and ret…

3447dd7

…urns and also removed the option check_symmetry from input validation

Added suggestions from codefactor

6073a88

Added suggestions from codefactor

0e11876

removes superflous import

8f9ca82

Adding -fopenmp to makefile for googletest

cfd9a39

Removing -lomp from googletest Makefile.

c50ed27

Merge branch 'master' into batchhafnian

98af9e7

increasing rtol to 0.2 in sample_test_hafnian.

206b7af

Changing the random seed.

d7ae9f1

Merge branch 'batchhafnian' of https://github.com/XanaduAI/hafnian in…

3238b65

…to batchhafnian

incresing tolerance on tests.

0b20f9f

Changes test probability vacuum

ee2ca7b

nquesada merged commit fa17e0b into master Jul 17, 2019

nquesada deleted the batchhafnian branch July 17, 2019 14:11

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Adding batchhafnian code #21

Adding batchhafnian code #21

bgupt commented Jun 28, 2019

nquesada commented Jun 30, 2019

nquesada commented Jun 30, 2019

codecov bot commented Jul 1, 2019 •

edited

Loading

nquesada Jul 9, 2019

bgupt Jul 10, 2019 •

edited

Loading

nquesada left a comment

nquesada Jul 9, 2019

bgupt Jul 10, 2019

nquesada left a comment

nquesada Jul 9, 2019

bgupt Jul 10, 2019

bgupt Jul 10, 2019

bgupt commented Jul 10, 2019

josh146 Jul 13, 2019

bgupt Jul 15, 2019 •

edited

Loading

josh146 Jul 16, 2019

nquesada Jul 16, 2019

josh146 Jul 16, 2019

nquesada Jul 16, 2019

josh146 Jul 16, 2019

nquesada Jul 16, 2019

josh146 Jul 16, 2019

josh146 Jul 16, 2019

josh146 Jul 16, 2019

josh146 Jul 16, 2019

josh146 Jul 16, 2019

nquesada Jul 16, 2019

	complex: The product of the array elements determines by index
	complex: The product of the array elements determined by index

Adding batchhafnian code #21

Adding batchhafnian code #21

Conversation

bgupt commented Jun 28, 2019

nquesada commented Jun 30, 2019

This is the part that would be great to move into the C++ code

Would it be possible to make the function batch hafnian in C++ take care of the division by all the factorials, i.e. have it return the array B constructed Above directly?

nquesada commented Jun 30, 2019

codecov bot commented Jul 1, 2019 • edited Loading

Codecov Report

Choose a reason for hiding this comment

bgupt Jul 10, 2019 • edited Loading

Choose a reason for hiding this comment

nquesada left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nquesada left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bgupt commented Jul 10, 2019

Choose a reason for hiding this comment

bgupt Jul 15, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Jul 1, 2019 •

edited

Loading

bgupt Jul 10, 2019 •

edited

Loading

bgupt Jul 15, 2019 •

edited

Loading