Switching from np.empty to np.zeros (#309)

**Description of the Change:**
`np.empty` returns an array of given shape and type, without
re-initializing its values.
This means that the array's entries are not zeros (they are equal to
whatever was in memory before that memory slot became part of the
returned array).
This creates instabilities in some of the hermite renormalized methods.

For this reason, we are switching to `np.zeros`

.github/CHANGELOG.md

@@ -69,6 +69,8 @@ which uses the old Numba code. When setting to a higher value, the new Julia cod
 [(#294)](https://github.com/XanaduAI/MrMustard/pull/294)
 * Fixed the documentations for loss_XYd and amp_XYd functions for Gaussian channels.
 [(#305)](https://github.com/XanaduAI/MrMustard/pull/305)
+* Replaced all instances of `np.empty` with `np.zeros` to fix instabilities.
+[(#309)](https://github.com/XanaduAI/MrMustard/pull/309)
 
 ### Documentation
 

mrmustard/math/lattice/strategies/compactFock/diagonal_amps.py

@@ -36,9 +36,9 @@ def use_offDiag_pivot(
     K_l = SQRT[pivot]
     K_i = SQRT[pivot + 1]
     if B.ndim == 1:
-        G_in = np.empty(2 * M, dtype=np.complex128)
+        G_in = np.zeros(2 * M, dtype=np.complex128)
     elif B.ndim == 2:
-        G_in = np.empty((2 * M, B.shape[1]), dtype=np.complex128)
+        G_in = np.zeros((2 * M, B.shape[1]), dtype=np.complex128)
 
     ########## READ ##########
     GB = arr1[(2 * d,) + params] * B
@@ -100,9 +100,9 @@ def use_diag_pivot(A, B, M, cutoffs, params, arr0, arr1):  # pragma: no cover
     K_l = SQRT[pivot]
     K_i = SQRT[pivot + 1]
     if B.ndim == 1:
-        G_in = np.empty(2 * M, dtype=np.complex128)
+        G_in = np.zeros(2 * M, dtype=np.complex128)
     elif B.ndim == 2:
-        G_in = np.empty((2 * M, B.shape[1]), dtype=np.complex128)
+        G_in = np.zeros((2 * M, B.shape[1]), dtype=np.complex128)
 
     ########## READ ##########
     GB = arr0[params] * B
@@ -180,27 +180,27 @@ def fock_representation_diagonal_amps(A, B, G0, M, cutoffs):
     list_type = numba.types.ListType(tuple_type)
 
     if B.ndim == 1:
-        arr0 = np.empty(cutoffs, dtype=np.complex128)
-        arr2 = np.empty((M,) + cutoffs, dtype=np.complex128)
-        arr1 = np.empty((2 * M,) + cutoffs, dtype=np.complex128)
+        arr0 = np.zeros(cutoffs, dtype=np.complex128)
+        arr2 = np.zeros((M,) + cutoffs, dtype=np.complex128)
+        arr1 = np.zeros((2 * M,) + cutoffs, dtype=np.complex128)
         if M == 1:
-            arr1010 = np.empty((1, 1, 1), dtype=np.complex128)
-            arr1001 = np.empty((1, 1, 1), dtype=np.complex128)
+            arr1010 = np.zeros((1, 1, 1), dtype=np.complex128)
+            arr1001 = np.zeros((1, 1, 1), dtype=np.complex128)
         else:
-            arr1010 = np.empty((M, M - 1) + cutoffs, dtype=np.complex128)
-            arr1001 = np.empty((M, M - 1) + cutoffs, dtype=np.complex128)
+            arr1010 = np.zeros((M, M - 1) + cutoffs, dtype=np.complex128)
+            arr1001 = np.zeros((M, M - 1) + cutoffs, dtype=np.complex128)
 
     elif B.ndim == 2:
         batch_length = B.shape[1]
-        arr0 = np.empty(cutoffs + (batch_length,), dtype=np.complex128)
-        arr2 = np.empty((M,) + cutoffs + (batch_length,), dtype=np.complex128)
-        arr1 = np.empty((2 * M,) + cutoffs + (batch_length,), dtype=np.complex128)
+        arr0 = np.zeros(cutoffs + (batch_length,), dtype=np.complex128)
+        arr2 = np.zeros((M,) + cutoffs + (batch_length,), dtype=np.complex128)
+        arr1 = np.zeros((2 * M,) + cutoffs + (batch_length,), dtype=np.complex128)
         if M == 1:
-            arr1010 = np.empty((1, 1, 1) + (batch_length,), dtype=np.complex128)
-            arr1001 = np.empty((1, 1, 1) + (batch_length,), dtype=np.complex128)
+            arr1010 = np.zeros((1, 1, 1) + (batch_length,), dtype=np.complex128)
+            arr1001 = np.zeros((1, 1, 1) + (batch_length,), dtype=np.complex128)
         else:
-            arr1010 = np.empty((M, M - 1) + cutoffs + (batch_length,), dtype=np.complex128)
-            arr1001 = np.empty((M, M - 1) + cutoffs + (batch_length,), dtype=np.complex128)
+            arr1010 = np.zeros((M, M - 1) + cutoffs + (batch_length,), dtype=np.complex128)
+            arr1001 = np.zeros((M, M - 1) + cutoffs + (batch_length,), dtype=np.complex128)
 
     arr0[(0,) * M] = G0
     return fock_representation_diagonal_amps_NUMBA(

mrmustard/math/lattice/strategies/compactFock/helperFunctions.py

@@ -20,7 +20,7 @@ def repeat_twice(params):
     Returns:
         (1D array): [a,a,b,b,c,c,...]
     """
-    pivot = np.empty(2 * len(params), dtype=np.int64)
+    pivot = np.zeros(2 * len(params), dtype=np.int64)
     for i, val in enumerate(params):
         pivot[2 * i] = val
         pivot[2 * i + 1] = val

mrmustard/math/lattice/strategies/compactFock/singleLeftoverMode_amps.py

@@ -107,7 +107,7 @@ def use_offDiag_pivot(
 
     ########## READ ##########
     read_GB = (2 * d,) + params
-    GB = np.empty((cutoff_leftoverMode, cutoff_leftoverMode, len(B)), dtype=np.complex128)
+    GB = np.zeros((cutoff_leftoverMode, cutoff_leftoverMode, len(B)), dtype=np.complex128)
     for m in range(cutoff_leftoverMode):
         for n in range(cutoff_leftoverMode):
             GB[m, n] = arr1[(m, n) + read_GB] * B
@@ -180,7 +180,7 @@ def use_diag_pivot(A, B, M, cutoff_leftoverMode, cutoffs_tail, params, arr0, arr
 
     ########## READ ##########
     read_GB = params
-    GB = np.empty((cutoff_leftoverMode, cutoff_leftoverMode, len(B)), dtype=np.complex128)
+    GB = np.zeros((cutoff_leftoverMode, cutoff_leftoverMode, len(B)), dtype=np.complex128)
     for m in range(cutoff_leftoverMode):
         for n in range(cutoff_leftoverMode):
             GB[m, n] = arr0[(m, n) + read_GB] * B
@@ -300,26 +300,24 @@ def fock_representation_1leftoverMode_amps(A, B, G0, M, cutoffs):
     list_type = numba.types.ListType(tuple_type)
     zero_tuple = (0,) * (M - 1)
 
-    arr0 = np.zeros(
-        (cutoff_leftoverMode, cutoff_leftoverMode) + cutoffs_tail, dtype=np.complex128
-    )  # doesn't work with np.empty
+    arr0 = np.zeros((cutoff_leftoverMode, cutoff_leftoverMode) + cutoffs_tail, dtype=np.complex128)
     arr0[(0,) * (M + 1)] = G0
-    arr2 = np.empty(
+    arr2 = np.zeros(
         (cutoff_leftoverMode, cutoff_leftoverMode) + (M - 1,) + cutoffs_tail, dtype=np.complex128
     )
-    arr1 = np.empty(
+    arr1 = np.zeros(
         (cutoff_leftoverMode, cutoff_leftoverMode) + (2 * (M - 1),) + cutoffs_tail,
         dtype=np.complex128,
     )
     if M == 2:
-        arr1010 = np.empty((1, 1, 1, 1, 1), dtype=np.complex128)
-        arr1001 = np.empty((1, 1, 1, 1, 1), dtype=np.complex128)
+        arr1010 = np.zeros((1, 1, 1, 1, 1), dtype=np.complex128)
+        arr1001 = np.zeros((1, 1, 1, 1, 1), dtype=np.complex128)
     else:
-        arr1010 = np.empty(
+        arr1010 = np.zeros(
             (cutoff_leftoverMode, cutoff_leftoverMode) + (M - 1, M - 2) + cutoffs_tail,
             dtype=np.complex128,
         )
-        arr1001 = np.empty(
+        arr1001 = np.zeros(
             (cutoff_leftoverMode, cutoff_leftoverMode) + (M - 1, M - 2) + cutoffs_tail,
             dtype=np.complex128,
         )

tests/test_lab/test_detectors.py

@@ -309,7 +309,7 @@ def test_sampling_mean_and_var(
             state = State(dm=state.dm(cutoffs=[40]))
         detector = Homodyne(0.0)
 
-        results = np.empty((self.N_MEAS, 2))
+        results = np.zeros((self.N_MEAS, 2))
         for i in range(self.N_MEAS):
             _ = state << detector
             results[i] = math.asnumpy(detector.outcome)