Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rectangular gkp #668

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 14 additions & 0 deletions .github/CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,18 @@

<h3>New features since last release</h3>

* Rectangular GKP states are now supported.
[(#668)](https://github.com/XanaduAI/strawberryfields/pull/668)

For example,

```python
prog = sf.Program(1)

with prog.context as q:
sf.ops.GKP(shape="rectangular", alpha=2) | q[0]
```

<h3>Breaking Changes</h3>

<h3>Bug fixes</h3>
Expand All @@ -12,6 +24,8 @@

This release contains contributions from (in alphabetical order):

J. Eli Bourassa

# Release 0.21.0 (current release)

<h3>New features since last release</h3>
Expand Down
27 changes: 21 additions & 6 deletions strawberryfields/backends/bosonicbackend/backend.py
Original file line number Diff line number Diff line change
Expand Up @@ -540,7 +540,9 @@ def prepare_cat_real_rep(self, a, theta, p, ampl_cutoff, D):

return weights, means, cov

def prepare_gkp(self, state, epsilon, ampl_cutoff, representation="real", shape="square"):
def prepare_gkp(
self, state, epsilon, ampl_cutoff, representation="real", shape="square", alpha=1
):
r"""Prepares the arrays of weights, means and covs for a finite energy GKP state.

GKP states are qubits, with the qubit state defined by:
Expand All @@ -555,6 +557,7 @@ def prepare_gkp(self, state, epsilon, ampl_cutoff, representation="real", shape=
ampl_cutoff (float): this determines how many terms to keep
representation (str): ``'real'`` or ``'complex'`` reprsentation
shape (str): shape of the lattice; default 'square'
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)

Returns:
tuple: arrays of the weights, means and covariances for the state
Expand All @@ -566,16 +569,24 @@ def prepare_gkp(self, state, epsilon, ampl_cutoff, representation="real", shape=
if representation == "complex":
raise NotImplementedError("The complex description of GKP is not implemented")

if shape != "square":
raise NotImplementedError("Only square GKP are implemented for now")
if shape not in ["square", "rectangular"]:
ilan-tz marked this conversation as resolved.
Show resolved Hide resolved
ilan-tz marked this conversation as resolved.
Show resolved Hide resolved
raise NotImplementedError("Only square and rectangular GKP are implemented.")

if shape == "square":
if alpha != 1:
raise ValueError(
"For square GKPs, alpha must be 1. For alpha not equal to "
+ "1, use shape='rectangular'."
elib20 marked this conversation as resolved.
Show resolved Hide resolved
)

theta, phi = state[0], state[1]

def coeff(peak_loc):
def coeff(peak_loc, alpha):
elib20 marked this conversation as resolved.
Show resolved Hide resolved
"""Returns the value of the weight for a given peak.

Args:
elib20 marked this conversation as resolved.
Show resolved Hide resolved
peak_loc (array): location of the ideal peak in phase space
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)

Returns:
float: weight of the peak
Expand Down Expand Up @@ -610,7 +621,7 @@ def coeff(peak_loc):
prefactor = np.exp(
-np.pi
* 0.25
* (l ** 2 + m ** 2)
* ((l * np.sqrt(alpha)) ** 2 + (m / np.sqrt(alpha)) ** 2)
* (1 - np.exp(-2 * epsilon))
/ (1 + np.exp(-2 * epsilon))
)
Expand Down Expand Up @@ -649,15 +660,19 @@ def coeff(peak_loc):
)

# Calculate the weights for each peak
weights = coeff(means)
weights = coeff(means, alpha)
filt = abs(weights) > ampl_cutoff
weights = weights[filt]

weights /= np.sum(weights)
# Apply finite energy effect to means
means = means[filt]

means[:, 0] *= np.sqrt(alpha)
means[:, 1] /= np.sqrt(alpha)

means *= 0.5 * damping * np.sqrt(np.pi * self.circuit.hbar)

# Covariances all the same
covs = (
0.5
Expand Down
16 changes: 12 additions & 4 deletions strawberryfields/backends/fockbackend/backend.py
Original file line number Diff line number Diff line change
Expand Up @@ -295,7 +295,7 @@ def measure_fock(self, modes, shots=1, select=None, **kwargs):
return self.circuit.measure_fock(self._remap_modes(modes), select=select)

def prepare_gkp(
self, state, epsilon, ampl_cutoff, representation="real", shape="square", mode=None
self, state, epsilon, ampl_cutoff, representation="real", shape="square", alpha=1, mode=None
):
r"""Prepares the Fock representation of a finite energy GKP state.

Expand All @@ -311,6 +311,7 @@ def prepare_gkp(
amplcutoff (float): this determines how many terms to keep
representation (str): ``'real'`` or ``'complex'`` reprsentation
shape (str): shape of the lattice; default 'square'
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)

Returns:
tuple: arrays of the weights, means and covariances for the state
Expand All @@ -322,8 +323,15 @@ def prepare_gkp(
if representation == "complex":
raise NotImplementedError("The complex description of GKP is not implemented")

if shape != "square":
raise NotImplementedError("Only square GKP are implemented for now")
if shape not in ["square", "rectangular"]:
raise NotImplementedError("Only square and rectangular GKP are implemented.")

if shape == "square":
if alpha != 1:
raise ValueError(
"For square GKPs, alpha must be 1. For alpha not equal to "
+ "1, use shape='rectangular'."
elib20 marked this conversation as resolved.
Show resolved Hide resolved
)

theta, phi = state[0], state[1]
self.circuit.prepare_gkp(theta, phi, epsilon, ampl_cutoff, self._remap_modes(mode))
self.circuit.prepare_gkp(theta, phi, epsilon, ampl_cutoff, alpha, self._remap_modes(mode))
elib20 marked this conversation as resolved.
Show resolved Hide resolved
8 changes: 5 additions & 3 deletions strawberryfields/backends/fockbackend/circuit.py
Original file line number Diff line number Diff line change
Expand Up @@ -799,13 +799,15 @@ def measure_homodyne(self, phi, mode, select=None, **kwargs):
# `homodyne_sample` will always be a single value since multiple shots is not supported
return np.array([[homodyne_sample]])

def prepare_gkp(self, theta, phi, epsilon, ampl_cutoff, mode):
def prepare_gkp(self, theta, phi, epsilon, ampl_cutoff, alpha, mode):
"""
Prepares a mode in a GKP state.
"""

if self._pure:
self.prepare(ops.square_gkp_state(theta, phi, epsilon, ampl_cutoff, self._trunc), mode)
self.prepare(
ops.rect_gkp_state(theta, phi, epsilon, ampl_cutoff, alpha, self._trunc), mode
)
else:
st = ops.square_gkp_state(theta, phi, epsilon, ampl_cutoff, self._trunc)
st = ops.rect_gkp_state(theta, phi, epsilon, ampl_cutoff, alpha, self._trunc)
self.prepare(np.outer(st, st.conjugate()), mode)
30 changes: 17 additions & 13 deletions strawberryfields/backends/fockbackend/ops.py
Original file line number Diff line number Diff line change
Expand Up @@ -515,7 +515,7 @@ def hermiteVals(q_mag, num_bins, m_omega_over_hbar, trunc):
return q_tensor, Hvals


def gkp_displacements(t, k, epsilon):
def gkp_displacements(t, k, epsilon, alpha):
"""
Helper function to generate the displacements parameters associated with the teeth of
GKP computational basis state k.
Expand All @@ -524,14 +524,15 @@ def gkp_displacements(t, k, epsilon):
t (array): the teeth of GKP computational basis
k (int): a computational basis state label, can be either 0 or 1
epsilon (float): finite energy parameter of the state
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)

Returns:
array: the displacements
"""
return np.sqrt(0.5 * np.pi) * (2 * t + k) / np.cosh(epsilon)
return np.sqrt(0.5 * np.pi * alpha) * (2 * t + k) / np.cosh(epsilon)


def gkp_coeffs(t, k, epsilon):
def gkp_coeffs(t, k, epsilon, alpha):
"""
Helper function to generate the coefficient parameters associated with the teeth of
GKP computational basis state k.
Expand All @@ -540,39 +541,41 @@ def gkp_coeffs(t, k, epsilon):
t (array): the teeth of GKP computational basis
k (int): a computational basis state label, can be either 0 or 1
epsilon (float): finite energy parameter of the state
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)

Returns:
array: the coefficients
"""
return np.exp(-0.5 * np.pi * np.tanh(epsilon) * (k + 2 * t) ** 2)
return np.exp(-0.5 * np.pi * alpha * np.tanh(epsilon) * (k + 2 * t) ** 2)


@functools.lru_cache()
def square_gkp_basis_state(i, epsilon, ampl_cutoff, cutoff):
def rect_gkp_basis_state(i, epsilon, ampl_cutoff, alpha, cutoff):
elib20 marked this conversation as resolved.
Show resolved Hide resolved
"""
Generate the Fock expansion of a (subnormalized) computational GKP basis state. Normalization occurs in the ``square_gkp_state`` method.
Generate the Fock expansion of a (subnormalized) computational GKP basis state. Normalization occurs in the ``rect_gkp_state`` method.

Args:
i (int): a computational basis state label, can be either 0 or 1
epsilon (float): finite energy parameter of the state
ampl_cutoff (float): this determines how many terms to keep in the Hilbert space expansion
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
cutoff (int): Fock space truncation

Returns
(array): the expansion of the ith computational basis state in the Fock basis

"""
z_max = int(np.ceil(np.sqrt(-0.25 / np.pi * np.log(ampl_cutoff) / np.tanh(epsilon))))
coeffs = [gkp_coeffs(t, i, epsilon) for t in range(-z_max, z_max + 1)]
coeffs = [gkp_coeffs(t, i, epsilon, alpha) for t in range(-z_max, z_max + 1)]
r = -0.5 * np.log(np.tanh(epsilon))
alphas = [gkp_displacements(t, i, epsilon) for t in range(-z_max, z_max + 1)]
num_kets = len(alphas)
ket = [coeffs[j] * displacedSqueezed(alphas[j], 0, r, 0, cutoff) for j in range(num_kets)]
disps = [gkp_displacements(t, i, epsilon, alpha) for t in range(-z_max, z_max + 1)]
num_kets = len(disps)
ket = [coeffs[j] * displacedSqueezed(disps[j], 0, r, 0, cutoff) for j in range(num_kets)]
return sum(ket)


@functools.lru_cache()
def square_gkp_state(theta, phi, epsilon, ampl_cutoff, cutoff):
def rect_gkp_state(theta, phi, epsilon, ampl_cutoff, alpha, cutoff):
elib20 marked this conversation as resolved.
Show resolved Hide resolved
r"""
Generate the Fock expansion of an abitrary GKP state parametrized as
:math:`|\psi\rangle = \cos{\tfrac{\theta}{2}} \vert 0 \rangle_{\rm gkp} + e^{-i \phi} \sin{\tfrac{\theta}{2}} \vert 1 \rangle_{\rm gkp}`.
Expand All @@ -582,15 +585,16 @@ def square_gkp_state(theta, phi, epsilon, ampl_cutoff, cutoff):
phi (float): the longitude with respect to the x-axis in the Bloch sphere
epsilon (float): finite energy parameter of the state
ampl_cutoff (float): this determines how many terms to keep
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
cutoff (int): Fock space truncation

Returns:
tuple: arrays of the weights, means and covariances for the state
"""
qubit_coeff0 = np.cos(theta / 2)
qubit_coeff1 = np.sin(theta / 2) * np.exp(-1j * phi)
ket0 = square_gkp_basis_state(0, epsilon, ampl_cutoff, cutoff)
ket1 = square_gkp_basis_state(1, epsilon, ampl_cutoff, cutoff)
ket0 = rect_gkp_basis_state(0, epsilon, ampl_cutoff, alpha, cutoff)
ket1 = rect_gkp_basis_state(1, epsilon, ampl_cutoff, alpha, cutoff)
ket = qubit_coeff0 * ket0 + qubit_coeff1 * ket1
ket /= np.linalg.norm(ket)
return ket
11 changes: 9 additions & 2 deletions strawberryfields/ops.py
Original file line number Diff line number Diff line change
Expand Up @@ -949,14 +949,21 @@ class GKP(Preparation):
cutoff (float): this determines how many terms to keep
representation (str): ``'real'`` or ``'complex'`` reprsentation
shape (str): shape of the lattice; default ``'square'``
alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
"""

def __init__(
self, state=None, epsilon=0.2, ampl_cutoff=1e-12, representation="real", shape="square"
self,
state=None,
epsilon=0.2,
ampl_cutoff=1e-12,
representation="real",
shape="square",
alpha=1,
):
if state is None:
state = [0, 0]
super().__init__([state, epsilon, ampl_cutoff, representation, shape])
super().__init__([state, epsilon, ampl_cutoff, representation, shape, alpha])

def _apply(self, reg, backend, **kwargs):
backend.prepare_gkp(*self.p, mode=reg[0])
Expand Down
57 changes: 57 additions & 0 deletions tests/backend/test_bosonic_backend.py
Original file line number Diff line number Diff line change
Expand Up @@ -474,6 +474,35 @@ def test_gkp_complex(self):
with pytest.raises(NotImplementedError):
backend.run_prog(prog)

@pytest.mark.parametrize("eps", EPS_VALS)
ilan-tz marked this conversation as resolved.
Show resolved Hide resolved
def test_sensor_gkp(self, eps):
r"""Checks that the GKP sensor state is invariant under Fourier gate."""
x = np.linspace(-3 * np.sqrt(np.pi), 3 * np.sqrt(np.pi), 40)
p = np.copy(x)

# Prepare GKP 0 and apply Hadamard
prog_sensor = sf.Program(1)
with prog_sensor.context as q:
sf.ops.GKP(state=[np.pi / 2, 0], epsilon=eps, shape="rectangular", alpha=2) | q[0]

backend_sensor = bosonic.BosonicBackend()
backend_sensor.run_prog(prog_sensor)
state_sensor = backend_sensor.state()
wigner_sensor = state_sensor.wigner(0, x, p)

# Prepare GKP +
prog_sensor_F = sf.Program(1)
with prog_sensor_F.context as qF:
sf.ops.GKP(state=[np.pi / 2, 0], epsilon=eps, shape="rectangular", alpha=2) | qF[0]
sf.ops.Rgate(np.pi / 2) | qF[0]

backend_sensor_F = bosonic.BosonicBackend()
backend_sensor_F.run_prog(prog_sensor_F)
state_sensor_F = backend_sensor_F.state()
wigner_sensor_F = state_sensor_F.wigner(0, x, p)

assert np.allclose(wigner_sensor, wigner_sensor_F)
ilan-tz marked this conversation as resolved.
Show resolved Hide resolved

@pytest.mark.parametrize("fock_eps", [0.1, 0.2, 0.3])
def test_gkp_wigner_compare_backends(self, fock_eps):
"""Test correct Wigner function are generated by comparing fock and bosonic outputs.
Expand Down Expand Up @@ -501,6 +530,34 @@ def test_gkp_wigner_compare_backends(self, fock_eps):
W_bosonic = gkp_bosonic.wigner(0, xvec, xvec)
assert np.allclose(W_fock, W_bosonic)

@pytest.mark.parametrize("fock_eps", [0.1, 0.2, 0.3])
ilan-tz marked this conversation as resolved.
Show resolved Hide resolved
@pytest.mark.parametrize("rect_ratio", [0.8, 1, 1.5])
def test_rect_gkp_wigner_compare_backends(self, fock_eps, rect_ratio):
"""Test that correct Wigner function is generated by comparing fock and bosonic outputs.
Since the fock backend cannot simulate very small epsilon, we restrict to fock_eps > 0.1."""
theta = np.pi / 8
elib20 marked this conversation as resolved.
Show resolved Hide resolved
phi = np.pi / 6
ket = [theta, phi]
xvec = np.linspace(-1, 1, 200)
prog_fock = sf.Program(1)
with prog_fock.context as q:
sf.ops.GKP(epsilon=fock_eps, state=ket, shape="rectangular", alpha=rect_ratio) | q
cutoff = 100
hbar = 2
eng_fock = sf.Engine("fock", backend_options={"cutoff_dim": cutoff, "hbar": hbar})
gkp_fock = eng_fock.run(prog_fock).state
W_fock = gkp_fock.wigner(0, xvec, xvec)

prog_bosonic = sf.Program(1)

with prog_bosonic.context as q:
sf.ops.GKP(epsilon=fock_eps, state=ket, shape="rectangular", alpha=rect_ratio) | q
hbar = 2
eng_bosonic = sf.Engine("bosonic", backend_options={"hbar": hbar})
gkp_bosonic = eng_bosonic.run(prog_bosonic).state
W_bosonic = gkp_bosonic.wigner(0, xvec, xvec)
assert np.allclose(W_fock, W_bosonic)
ilan-tz marked this conversation as resolved.
Show resolved Hide resolved


class TestBosonicUserSpecifiedState:
"""Checks the Bosonic preparation method."""
Expand Down
18 changes: 14 additions & 4 deletions tests/frontend/test_ops_gate.py
Original file line number Diff line number Diff line change
Expand Up @@ -174,12 +174,22 @@ class TestGKPBasics:
@pytest.mark.parametrize("ampl", [0, 0.001])
@pytest.mark.parametrize("eps", [0, 0.001])
@pytest.mark.parametrize("r", ["real", "complex"])
@pytest.mark.parametrize("s", ["square"])
def test_gkp_str_representation(self, state, ampl, eps, r, s):
@pytest.mark.parametrize("s", ["square", "rectangular"])
@pytest.mark.parametrize("alph", [0.8, 1.1])
def test_gkp_str_representation(self, state, ampl, eps, r, s, alph):
"""Test the string representation of the GKP operation"""
assert (
str(ops.GKP(state=state, ampl_cutoff=ampl, epsilon=eps, representation=r, shape=s))
== f"GKP({str(state)}, {str(eps)}, {str(ampl)}, {r}, {s})"
str(
ops.GKP(
state=state,
ampl_cutoff=ampl,
epsilon=eps,
representation=r,
shape=s,
alpha=alph,
)
)
== f"GKP({str(state)}, {str(eps)}, {str(ampl)}, {r}, {s}, {str(alph)})"
)


Expand Down