Skip to content

Commit

Permalink
Gaussian state preparation + debugging (#522)
Browse files Browse the repository at this point in the history 
* Gaussian state preparation + debugging

- Bosonic backend now supports Thermal, Vac, Coherent, Squeezed, DispSqueezed, Gaussian state prep
- Minor debugging in homodyne/heterodyne, thermal_loss, mbsqueezing

* Apply suggestions from code review

Deleting qmat and Amat since they are only relevant for gaussian backend.

* Update strawberryfields/backends/bosonicbackend/bosoniccircuit.py

Deleting two lines from Amat that were missed

* Ran black
  • Loading branch information
@elib20
elib20 committed Jan 19, 2021
1 parent 4b9d828 commit 106372f
Show file tree
Hide file tree
Showing 2 changed files with 87 additions and 140 deletions.
139 changes: 70 additions & 69 deletions strawberryfields/backends/bosonicbackend/backend.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,7 +265,7 @@ def begin_circuit(self, num_subsystems, **kwargs):
self.circuit = BosonicModes(num_subsystems)

def add_mode(self, n=1):
self.circuit.add_mode(n)
self.circuit.add_mode([n])

def del_mode(self, modes):
self.circuit.del_mode(modes)
Expand All @@ -277,28 +277,23 @@ def reset(self, pure=True, **kwargs):
self.circuit.reset(self._init_modes)

def prepare_thermal_state(self, nbar, mode):
# self.circuit.init_thermal(nbar, mode)
pass
self.circuit.init_thermal(nbar, mode)

def prepare_vacuum_state(self, mode):
# self.circuit.loss(0.0, mode)
pass
self.circuit.loss(0.0, mode)

def prepare_coherent_state(self, r, phi, mode):
# self.circuit.loss(0.0, mode)
# self.circuit.displace(r, phi, mode)
pass
self.circuit.loss(0.0, mode)
self.circuit.displace(r, phi, mode)

def prepare_squeezed_state(self, r, phi, mode):
# self.circuit.loss(0.0, mode)
# self.circuit.squeeze(r, phi, mode)
pass
self.circuit.loss(0.0, mode)
self.circuit.squeeze(r, phi, mode)

def prepare_displaced_squeezed_state(self, r_d, phi_d, r_s, phi_s, mode):
# self.circuit.loss(0.0, mode)
# self.circuit.squeeze(r_s, phi_s, mode)
# self.circuit.displace(r_d, phi_d, mode)
pass
self.circuit.loss(0.0, mode)
self.circuit.squeeze(r_s, phi_s, mode)
self.circuit.displace(r_d, phi_d, mode)

def prepare_cat(self, alpha, phi, cutoff, desc, D):
""" Prepares the arrays of weights, means and covs for a cat state"""
Expand Down Expand Up @@ -598,13 +593,13 @@ def measure_homodyne(self, phi, mode, shots=1, select=None, **kwargs):
self.circuit.phase_shift(-phi, mode)

if select is None:
qs = self.circuit.homodyne(mode, **kwargs)[0, 0]
val = self.circuit.homodyne(mode, **kwargs)[0, 0]
else:
val = select * 2 / np.sqrt(2 * self.circuit.hbar)
qs = self.circuit.post_select_homodyne(mode, val, **kwargs)
self.circuit.post_select_homodyne(mode, val, **kwargs)

# `qs` will always be a single value since multiple shots is not supported
return np.array([[qs * np.sqrt(2 * self.circuit.hbar) / 2]])
return np.array([[val * np.sqrt(2 * self.circuit.hbar) / 2]])

def measure_heterodyne(self, mode, shots=1, select=None):

Expand Down Expand Up @@ -662,59 +657,65 @@ def thermal_loss(self, T, nbar, mode):
self.circuit.thermal_loss(T, nbar, mode)

def measure_fock(self, modes, shots=1, select=None, **kwargs):
if select is not None:
raise NotImplementedError(
"Gaussian backend currently does not support " "postselection"
)
if shots != 1:
warnings.warn(
"Cannot simulate non-Gaussian states. "
"Conditional state after Fock measurement has not been updated."
)

mu = self.circuit.mean
mean = self.circuit.smeanxp()
cov = self.circuit.scovmatxp()

x_idxs = np.array(modes)
p_idxs = x_idxs + len(mu)
modes_idxs = np.concatenate([x_idxs, p_idxs])
reduced_cov = cov[np.ix_(modes_idxs, modes_idxs)]
reduced_mean = mean[modes_idxs]

# check we are sampling from a gaussian state with zero mean
if np.allclose(mu, np.zeros_like(mu)):
samples = hafnian_sample_state(reduced_cov, shots)
else:
samples = hafnian_sample_state(reduced_cov, shots, mean=reduced_mean)

return samples
# if select is not None:
# raise NotImplementedError(
# "Gaussian backend currently does not support " "postselection"
# )
# if shots != 1:
# warnings.warn(
# "Cannot simulate non-Gaussian states. "
# "Conditional state after Fock measurement has not been updated."
# )

# mu = self.circuit.mean
# mean = self.circuit.smeanxp()
# cov = self.circuit.scovmatxp()

# x_idxs = np.array(modes)
# p_idxs = x_idxs + len(mu)
# modes_idxs = np.concatenate([x_idxs, p_idxs])
# reduced_cov = cov[np.ix_(modes_idxs, modes_idxs)]
# reduced_mean = mean[modes_idxs]

# # check we are sampling from a gaussian state with zero mean
# if np.allclose(mu, np.zeros_like(mu)):
# samples = hafnian_sample_state(reduced_cov, shots)
# else:
# samples = hafnian_sample_state(reduced_cov, shots, mean=reduced_mean)

# return samples

# TODO
pass

def measure_threshold(self, modes, shots=1, select=None, **kwargs):
if shots != 1:
if select is not None:
raise NotImplementedError(
"Gaussian backend currently does not support " "postselection"
)
warnings.warn(
"Cannot simulate non-Gaussian states. "
"Conditional state after Threshold measurement has not been updated."
)

mu = self.circuit.mean
cov = self.circuit.scovmatxp()
# check we are sampling from a gaussian state with zero mean
if not np.allclose(mu, np.zeros_like(mu)):
raise NotImplementedError(
"Threshold measurement is only supported for " "Gaussian states with zero mean"
)
x_idxs = np.array(modes)
p_idxs = x_idxs + len(mu)
modes_idxs = np.concatenate([x_idxs, p_idxs])
reduced_cov = cov[np.ix_(modes_idxs, modes_idxs)]
samples = torontonian_sample_state(reduced_cov, shots)

return samples
# if shots != 1:
# if select is not None:
# raise NotImplementedError(
# "Gaussian backend currently does not support " "postselection"
# )
# warnings.warn(
# "Cannot simulate non-Gaussian states. "
# "Conditional state after Threshold measurement has not been updated."
# )

# mu = self.circuit.mean
# cov = self.circuit.scovmatxp()
# # check we are sampling from a gaussian state with zero mean
# if not np.allclose(mu, np.zeros_like(mu)):
# raise NotImplementedError(
# "Threshold measurement is only supported for " "Gaussian states with zero mean"
# )
# x_idxs = np.array(modes)
# p_idxs = x_idxs + len(mu)
# modes_idxs = np.concatenate([x_idxs, p_idxs])
# reduced_cov = cov[np.ix_(modes_idxs, modes_idxs)]
# samples = torontonian_sample_state(reduced_cov, shots)

# return samples

# TODO
pass

def state(self, modes=None, peaks=None, **kwargs):
"""Returns the state of the quantum simulation.
Expand Down
88 changes: 17 additions & 71 deletions strawberryfields/backends/bosonicbackend/bosoniccircuit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,10 +177,11 @@ def mbsqueeze(self, k, r, phi, r_anc, eta_anc, avg):
if self.active[k] is None:
raise ValueError("Cannot squeeze mode, mode does not exist")

phi = phi + (1 - np.sign(r)) * np.pi / 2
if r < 0:
phi += np.pi
r = np.abs(r)
theta = np.arccos(np.exp(-r))
self.phase_shift(phi / 2, k)
self.phase_shift(-phi / 2, k)

if avg:
X = np.diag([np.cos(theta), 1 / np.cos(theta)])
Expand Down Expand Up @@ -212,7 +213,7 @@ def mbsqueeze(self, k, r, phi, r_anc, eta_anc, avg):
self.active = self.active[:new_mode]
prefac = -np.tan(theta) / np.sqrt(2 * self.hbar * eta_anc)
self.displace(prefac * val[0][0], np.pi / 2, k)
self.phase_shift(-phi / 2, k)
self.phase_shift(phi / 2, k)

if not avg:
return val
Expand Down Expand Up @@ -281,12 +282,11 @@ def fromsmean(self, r, modes=None):
r (array): vector of means in :math:`(x_1,p_1,x_2,p_2,\dots)` ordering
modes (Sequence): sequence of modes corresponding to the vector of means
"""
mode_list = modes
if modes is None:
mode_list = range(self.nlen)
modes = range(self.nlen)

for idx, mode in enumerate(mode_list):
self.mean[mode] = 0.5 * (r[2 * idx] + 1j * r[2 * idx + 1])
mode_ind = np.sort(np.append(2 * np.array(modes), 2 * np.array(modes) + 1))
self.means[:, mode_ind] = r

def fromscovmat(self, V, modes=None):
r"""Updates the circuit's state when a standard covariance matrix is provided.
Expand All @@ -309,50 +309,12 @@ def fromscovmat(self, V, modes=None):
if n > self.nlen:
raise ValueError("Covariance matrix is larger than the number of subsystems.")

# convert to xp ordering
rotmat = changebasis(n)
VV = np.dot(np.dot(np.transpose(rotmat), V), rotmat)

A = VV[0:n, 0:n]
B = VV[0:n, n : 2 * n]
C = VV[n : 2 * n, n : 2 * n]
Bt = np.transpose(B)

if n < self.nlen:
# reset modes to be prepared back to the vacuum state
for mode in modes:
self.loss(0.0, mode)

rows = modes.reshape(-1, 1)
cols = modes.reshape(1, -1)
self.nmat[rows, cols] = 0.25 * (A + C + 1j * (B - Bt) - 2 * np.identity(n))
self.mmat[rows, cols] = 0.25 * (A - C + 1j * (B + Bt))

def qmat(self, modes=None):
""" Construct the covariance matrix for the Q function"""
if modes is None:
modes = list(range(self.nlen))

rows = np.reshape(modes, [-1, 1])
cols = np.reshape(modes, [1, -1])

sigmaq = (
np.concatenate(
(
np.concatenate(
(self.nmat[rows, cols], np.conjugate(self.mmat[rows, cols])),
axis=1,
),
np.concatenate(
(self.mmat[rows, cols], np.conjugate(self.nmat[rows, cols])),
axis=1,
),
),
axis=0,
)
+ np.identity(2 * len(modes))
)
return sigmaq
mode_ind = np.sort(np.append(2 * np.array(modes), 2 * np.array(modes) + 1))
# Delete current entries of covariance and any correlations to other modes
self.covs[:, mode_ind, :] = 0
self.covs[:, :, mode_ind] = 0
# Set new covariance elements
self.covs[np.ix_(np.arange(self.covs.shape[0], dtype=int), mode_ind, mode_ind)] = V

def fidelity_coherent(self, alpha, modes=None):
""" Returns a function that evaluates the Q function of the given state """
Expand Down Expand Up @@ -406,21 +368,6 @@ def parity_val(self, modes=None):
parity = np.sum(weighted_exp)
return parity

def Amat(self):
""" Constructs the A matrix from Hamilton's paper"""
######### this needs to be conjugated
sigmaq = (
np.concatenate(
(
np.concatenate((np.transpose(self.nmat), self.mmat), axis=1),
np.concatenate((np.transpose(np.conjugate(self.mmat)), self.nmat), axis=1),
),
axis=0,
)
+ np.identity(2 * self.nlen)
)
return np.dot(Xmat(self.nlen), np.identity(2 * self.nlen) - np.linalg.inv(sigmaq))

def loss(self, T, k):
r"""Implements a loss channel in mode k by amplitude loss amount \sqrt{T}
(energy loss amount T)"""
Expand All @@ -440,14 +387,13 @@ def thermal_loss(self, T, nbar, k):
if self.active[k] is None:
raise ValueError("Cannot apply loss channel, mode does not exist")
X = np.sqrt(T) * np.identity(2)
Y = self.hbar * (1 - T) * nbar * np.identity(2) / 2
Y = self.hbar * (1 - T) * (2 * nbar + 1) * np.identity(2) / 2
X2, Y2 = self.expandXY([k], X, Y)
self.apply_channel(X2, Y2)

def init_thermal(self, population, mode):
def init_thermal(self, nbar, mode):
""" Initializes a state of mode in a thermal state with the given population"""
# self.loss(0.0, mode)
# self.nmat[mode][mode] = population
self.thermal_loss(0.0, nbar, mode)

def is_vacuum(self, tol=0.0):
""" Checks if the state is vacuum by calculating its fidelity with vacuum """
Expand Down Expand Up @@ -632,7 +578,7 @@ def post_select_heterodyne(self, n, alpha_val):

covmat = self.hbar * np.identity(2) / 2
indices = [n]
vals = np.array(alpha_val.real, alpha_val.imag)
vals = np.array([alpha_val.real, alpha_val.imag])
self.post_select_generaldyne(covmat, indices, vals)
return

Expand Down

0 comments on commit 106372f

Please sign in to comment.