How to use sf.fock.probe in the Quantum Optics Neural Networks?

I need some help!
Does Pennylane have a function like sf.fock.probe in StrawberryField? I want to know the probabilities of some certain fock states(but not the probabilities of the computational basis by qml.probs) from the output of the [Quantum Optics Neural Networks] and train the neural network based on the probabilities. In the previous work of [Quantum Optics Neural Networks], they use the exact output(e.g. [0, 1, 1, 0]) to train the QONN, but what if the output is the superposition of some certain fock states(e.g. the output is 0.70|0110>+0.70|1010> )? Could you please give me some advice? Thanks very much.

Hi @Hamitonian!

As you correctly point out, you can use qml.probs to extract all computational basis probabilities up to a certain cutoff.

If instead you would like the probabilities of a single Fock state, you can instead use qml.FockStateProjector:

import pennylane as qml
from pennylane import numpy as np

dev = qml.device("strawberryfields.fock", cutoff_dim=5, wires=2)

@qml.qnode(dev)
def circuit(x, state):
    qml.Squeezing(x[0], 0, wires=0)
    qml.Squeezing(x[1], 0, wires=1)
    qml.Beamsplitter(x[2], x[3], wires=[0, 1])
    return qml.expval(qml.FockStateProjector(state))

x = np.array([0.1, 0.2, 0.3, 0.4], requires_grad=True)
states = [
    np.array([0, 2], requires_grad=False),
    np.array([2, 0], requires_grad=False)
]

Executing this QNode:

>>> circuit(x, states[0])
0.01692903

Your advice does help a lot! Thanks.
I modified the code in the tutorial Optimizing a Quantum Optics Neural Network, but it seems that the code has something wrong.
I want the label of training set to be the probabilities of some certain Fock states, but it says that Each wire in the quantum circuit can only be measured once.
However I need to know the probabilities of the states |1010>,|0101>,|1001>,|0110> in the same circuits. My codes are as follows, could you please help me to modify it?

import pennylane as qml
from pennylane import numpy as np
from strawberryfields.ops import *
import nlopt
dev = qml.device("strawberryfields.fock", wires=4, cutoff_dim=4)

def layer(theta, phi, wires):
    M = len(wires)
    phi_nonlinear = np.pi / 2

    qml.templates.Interferometer(
        theta, phi, np.zeros(M), wires=wires, mesh="rectangular", beamsplitter="pennylane"
    )

    for i in wires:
        qml.Kerr(phi_nonlinear, wires=i)

@qml.qnode(dev)
def quantum_neural_net(var, x, states = [
    np.array([1, 0, 1, 0], requires_grad=False),
    np.array([0, 1, 0, 1], requires_grad=False),
    np.array([1, 0, 0, 1], requires_grad=False),
    np.array([0, 1, 1, 0], requires_grad=False)
]):
    wires = list(range(len(x)))

    # Encode input x into a sequence of quantum fock states
    for i in wires:
        qml.FockState(x[i], wires=i)

    # "layer" subcircuits
    for i, v in enumerate(var):
        layer(v[: len(v) // 2], v[len(v) // 2 :], wires)

    return [qml.expval(qml.FockStateProjector(w)) for w in states]

def square_loss(labels, predictions):
    term = 0
    for l, p in zip(labels, predictions):
        lnorm = l / np.linalg.norm(l)
        pnorm = p / np.linalg.norm(p)

        term = term + np.abs(np.dot(lnorm, pnorm.T)) ** 2

    return 1 - term / len(labels)

def cost(var, data_input, labels):
    predictions = np.array([quantum_neural_net(var, x) for x in data_input])
    sl = square_loss(labels, predictions)

    return sl


X = np.array([[1, 0, 1, 0]], requires_grad=False)

Y = np.array([[0.1, 0.5, 0.2, 0.2]], requires_grad=False)

num_layers = 4
M = len(X[0])
num_variables_per_layer = M * (M - 1)

var_init = (4 * np.random.rand(num_layers, num_variables_per_layer) - 2) * np.pi
print(var_init)

cost_grad = qml.grad(cost)

print_every = 1

evals = 0
def cost_wrapper(var, grad=[ ]):
    global evals
    evals += 1

    if grad.size > 0:
        # Get the gradient for `var` (idx 0) by first "unflattening" it
        var_grad = cost_grad(var.reshape((num_layers, num_variables_per_layer)), X, Y)[0]
        grad[:] = var_grad.flatten()
    cost_val = cost(var.reshape((num_layers, num_variables_per_layer)), X, Y)

    if evals % print_every == 0:
        print(f"Iter: {evals:4d}    Cost: {cost_val:.4e}")

    return float(cost_val)


opt_algorithm = nlopt.LD_LBFGS  # Gradient-based
#opt_algorithm = nlopt.LN_BOBYQA  # Gradient-free

opt = nlopt.opt(opt_algorithm, num_layers*num_variables_per_layer)

opt.set_min_objective(cost_wrapper)

opt.set_lower_bounds(-2*np.pi * np.ones(num_layers*num_variables_per_layer))
opt.set_upper_bounds(2*np.pi * np.ones(num_layers*num_variables_per_layer))

var = opt.optimize(var_init.flatten())
var = var.reshape(var_init.shape)

Another problem is that I can only get the probability of one fock state, because of the QuantumFunctionError: Each wire in the quantum circuit can only be measured once. I want to know the probabilities of some certain fock states simultaneously, for example the states:

states = [
    np.array([1, 0, 1, 0], requires_grad=False),
    np.array([0, 1, 0, 1], requires_grad=False),
    np.array([1, 0, 0, 1], requires_grad=False),
    np.array([0, 1, 1, 0], requires_grad=False)
]

in the same quantum circuits and wires. Thanks.

Hi @Hamitonian!

Yes this is a good catch , we currently do not support measurements on more than one wire.

However, this is something we can work to support with the Strawberry Fields devices – I’ll forward this on to the PennyLane team :slight_smile:

Thanks a lot:grinning: