From 7eee9c65c0d5b4f413d12999ab664c941c266fb7 Mon Sep 17 00:00:00 2001 From: rachelchadwick <33873960+rachelchadwick@users.noreply.github.com> Date: Wed, 10 Aug 2022 14:19:15 -0400 Subject: [PATCH 1/4] Update generate_hafnian_sample --- thewalrus/samples.py | 50 +++++++++++++++++++------------------------- 1 file changed, 22 insertions(+), 28 deletions(-) diff --git a/thewalrus/samples.py b/thewalrus/samples.py index addb0c1ff..45416706f 100644 --- a/thewalrus/samples.py +++ b/thewalrus/samples.py @@ -112,56 +112,50 @@ def generate_hafnian_sample( N = len(cov) // 2 result = [] prev_prob = 1.0 - nmodes = N if mean is None: local_mu = np.zeros(2 * N) else: local_mu = mean - for k in range(nmodes): - probs1 = np.zeros([cutoff + 1], dtype=np.float64) + for k in range(N): + total_photons = int(np.sum(result)) + prob_cumulative = 0 kk = np.arange(k + 1) mu_red, V_red = reduced_gaussian(local_mu, cov, kk) + sample = np.random.random() if approx: Q = Qmat(V_red, hbar=hbar) A = Amat(Q, hbar=hbar, cov_is_qmat=True) + A = np.real_if_close(A) + A[np.where((-1e-8 < A) & (A < 0))] = 0 + normalisation = np.sqrt(np.linalg.det(Q).real) - for i in range(cutoff): + for i in range(min(max_photons + 1 - total_photons, cutoff) + 1): + if i == min(max_photons + 1 - total_photons, cutoff) + 0: + return -1 indices = result + [i] ind2 = indices + indices if approx: factpref = np.prod(fac(indices)) mat = reduction(A, ind2) - probs1[i] = ( - hafnian(np.abs(mat.real), approx=True, num_samples=approx_samples) / factpref + prob_i = ( + hafnian(mat, approx=True, num_samples=int(approx_samples)) + / factpref + / normalisation ) else: - probs1[i] = density_matrix_element( + prob_i = density_matrix_element( mu_red, V_red, indices, indices, include_prefactor=True, hbar=hbar ).real - if approx: - probs1 = probs1 / np.sqrt(np.linalg.det(Q).real) - - probs2 = probs1 / prev_prob - probs3 = np.maximum( - probs2, np.zeros_like(probs2) - ) # pylint: disable=assignment-from-no-return - ssum = np.sum(probs3) - if ssum < 1.0: - probs3[-1] = 1.0 - ssum - - # The following normalization of probabilities is needed to prevent np.random.choice error - if ssum > 1.0: - probs3 = probs3 / ssum - - result.append(np.random.choice(a=range(len(probs3)), p=probs3)) - if result[-1] == cutoff: - return -1 - if np.sum(result) > max_photons: - return -1 - prev_prob = probs1[result[-1]] + prob_cumulative += prob_i / prev_prob + if prob_cumulative >= sample: + result.append(i) + prev_prob = prob_i + break + + prev_prob = prob_i return result From b55377a59ad33436fb23630c4727f10bcfc643d0 Mon Sep 17 00:00:00 2001 From: rachelchadwick <33873960+rachelchadwick@users.noreply.github.com> Date: Wed, 10 Aug 2022 15:26:23 -0400 Subject: [PATCH 2/4] Run black --- thewalrus/samples.py | 30 ++++++++++++++++++++++++------ 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/thewalrus/samples.py b/thewalrus/samples.py index 45416706f..25da8bf4e 100644 --- a/thewalrus/samples.py +++ b/thewalrus/samples.py @@ -269,7 +269,9 @@ def hafnian_sample_state( np.array[int]: photon number samples from the Gaussian state """ if parallel: - params = [[cov, 1, mean, hbar, cutoff, max_photons, approx, approx_samples]] * samples + params = [ + [cov, 1, mean, hbar, cutoff, max_photons, approx, approx_samples] + ] * samples compute_list = [] for p in params: compute_list.append(dask.delayed(_hafnian_sample)(p)) @@ -283,7 +285,14 @@ def hafnian_sample_state( def hafnian_sample_graph( - A, n_mean, samples=1, cutoff=5, max_photons=30, approx=False, approx_samples=1e5, parallel=False + A, + n_mean, + samples=1, + cutoff=5, + max_photons=30, + approx=False, + approx_samples=1e5, + parallel=False, ): r"""Returns samples from the Gaussian state specified by the adjacency matrix :math:`A` and with total mean photon number :math:`n_{mean}` @@ -427,7 +436,9 @@ def _torontonian_sample(args): j = 0 while j < samples: - result = generate_torontonian_sample(cov, mu, hbar=hbar, max_photons=max_photons) + result = generate_torontonian_sample( + cov, mu, hbar=hbar, max_photons=max_photons + ) if result != -1: samples_array.append(result) j = j + 1 @@ -435,7 +446,9 @@ def _torontonian_sample(args): return np.vstack(samples_array) -def torontonian_sample_state(cov, samples, mu=None, hbar=2, max_photons=30, parallel=False): +def torontonian_sample_state( + cov, samples, mu=None, hbar=2, max_photons=30, parallel=False +): r"""Returns samples from the Torontonian of a Gaussian state Args: @@ -550,7 +563,10 @@ def torontonian_sample_classical_state(cov, samples, mean=None, hbar=2, atol=1e- np.array[int]: threshold samples from the Gaussian state with covariance cov and vector means mean. """ return np.where( - hafnian_sample_classical_state(cov, samples, mean=mean, hbar=hbar, atol=atol) > 0, 1, 0 + hafnian_sample_classical_state(cov, samples, mean=mean, hbar=hbar, atol=atol) + > 0, + 1, + 0, ) @@ -574,7 +590,9 @@ def photon_number_sampler(probabilities, num_samples, out_of_bounds=False): probabilities = probabilities.flatten() / sum_p vals = np.arange(cutoff**num_modes, dtype=int) return [ - np.unravel_index(np.random.choice(vals, p=probabilities), [cutoff] * num_modes) + np.unravel_index( + np.random.choice(vals, p=probabilities), [cutoff] * num_modes + ) for _ in range(num_samples) ] From 96341bf08329932609a9a48af0340f88fff22d70 Mon Sep 17 00:00:00 2001 From: rachelchadwick <33873960+rachelchadwick@users.noreply.github.com> Date: Wed, 10 Aug 2022 18:55:31 -0400 Subject: [PATCH 3/4] Change kwarg values --- thewalrus/samples.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/thewalrus/samples.py b/thewalrus/samples.py index 25da8bf4e..a9db5f94d 100644 --- a/thewalrus/samples.py +++ b/thewalrus/samples.py @@ -88,7 +88,13 @@ # pylint: disable=too-many-branches def generate_hafnian_sample( - cov, mean=None, hbar=2, cutoff=6, max_photons=30, approx=False, approx_samples=1e5 + cov, + mean=None, + hbar=2, + cutoff=None, + max_photons=30, + approx=False, + approx_samples=1e5, ): # pylint: disable=too-many-branches r"""Returns a single sample from the Hafnian of a Gaussian state. @@ -112,6 +118,8 @@ def generate_hafnian_sample( N = len(cov) // 2 result = [] prev_prob = 1.0 + if cutoff is None: + cutoff = max_photons + 1 if mean is None: local_mu = np.zeros(2 * N) else: From 0e6ea9daa48f964d6711ec0597ab62ac47a2c75d Mon Sep 17 00:00:00 2001 From: rachelchadwick <33873960+rachelchadwick@users.noreply.github.com> Date: Wed, 10 Aug 2022 20:51:08 -0400 Subject: [PATCH 4/4] Pure state final mode --- thewalrus/samples.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/thewalrus/samples.py b/thewalrus/samples.py index a9db5f94d..95b455ef7 100644 --- a/thewalrus/samples.py +++ b/thewalrus/samples.py @@ -66,6 +66,8 @@ is_classical_cov, reduced_gaussian, density_matrix_element, + is_pure_cov, + pure_state_amplitude, ) __all__ = [ @@ -124,6 +126,7 @@ def generate_hafnian_sample( local_mu = np.zeros(2 * N) else: local_mu = mean + pure_cov = is_pure_cov(cov, hbar=hbar) for k in range(N): total_photons = int(np.sum(result)) @@ -153,9 +156,20 @@ def generate_hafnian_sample( / normalisation ) else: - prob_i = density_matrix_element( - mu_red, V_red, indices, indices, include_prefactor=True, hbar=hbar - ).real + if (k == N - 1) and (pure_cov == True): + amp_i = pure_state_amplitude( + mu_red, cov, indices, hbar=hbar, check_purity=False + ) + prob_i = np.abs(amp_i) ** 2 + else: + prob_i = density_matrix_element( + mu_red, + V_red, + indices, + indices, + include_prefactor=True, + hbar=hbar, + ).real prob_cumulative += prob_i / prev_prob if prob_cumulative >= sample: