How to use the strawberryfields.Engine function in StrawberryFields

def SF_gate_reference(sf_op, wires, *args):
    """SF reference circuit for gate tests"""
    eng = sf.Engine("gaussian")
    prog = sf.Program(2)
    with prog.context as q:
        sf.ops.S2gate(0.1) | q
        sf_op(*args) | [q[i] for i in wires]

    state = eng.run(prog).state
    return state.mean_photon(0)[0], state.mean_photon(1)[0]

def SF_expectation_reference(sf_expectation, wires, num_wires, *args):
    """SF reference circuit for expectation tests"""
    eng = sf.Engine("gaussian")

    # Allows returning the variance of tensor number for 3 modes
    prog = sf.Program(num_wires)
    with prog.context as q:
        sf.ops.Dgate(0.1) | q[0]
        sf.ops.S2gate(0.1) | (q[0], q[1])

    state = eng.run(prog).state
    return sf_expectation(state, wires, args)[0]

def build_circuit(self):
        params_counter = 0
        sgates = []
        dgates = []
        kgates = []
        for gate_structure in self.gates_structure:
            if gate_structure[0] is Sgate:
                sgates.append(ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2]), make_param(**gate_structure[3])]))
            if gate_structure[0] is Dgate:
                dgates.append(ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2]), make_param(**gate_structure[3])]))
            if gate_structure[0] is Kgate:
                kgates.append(ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2])]))

        eng, q = sf.Engine(self.n_qumodes)

        rl, U = takagi(self.adj_matrix)
        initial_squeezings = np.tanh(rl)

        with eng:
            for i ,squeeze_value in enumerate(initial_squeezings):
                Sgate(squeeze_value) | i

            Interferometer(U) | q

            for gate in sgates:
                gate.gate(gate.params[0], gate.params[1]) | gate.qumodes

            Interferometer(self.interferometer_matrix) | q

            for gate in dgates:

current_layer = int(gate_structure[2]['name'].split('_')[-1][0])
            if gate_structure[0] is Sgate:
                current_gate = ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2]), make_param(**gate_structure[3])])
                all_sgates[current_layer].append(current_gate)
            if gate_structure[0] is Dgate:
                current_gate = ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2]), make_param(**gate_structure[3])])
                all_dgates[current_layer].append(current_gate)
            if gate_structure[0] is Kgate:
                current_gate = ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2])])
                all_kgates[current_layer].append(current_gate)
            if gate_structure[0] is Vgate:
                current_gate = ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2])])
                all_vgates[current_layer].append(current_gate)


        eng, q = sf.Engine(self.n_qumodes)
        rl, U = takagi(adj_matrix)
        initial_squeezings = np.arctanh(rl)

        with eng:
            for i ,squeeze_value in enumerate(initial_squeezings):
                Sgate(squeeze_value) | i

            Interferometer(U) | q
            for layer in range(number_of_layers):
                sgates = all_sgates[layer]
                dgates = all_dgates[layer]
                kgates = all_kgates[layer]
                vgates = all_vgates[layer]

                if len(sgates) != 0:
                    Interferometer(self.interferometer_matrix) | q

def build_circuit(self):
        params_counter = 0
        sgates = []
        dgates = []
        kgates = []
        for gate_structure in self.gates_structure:
            if gate_structure[0] is Sgate:
                sgates.append(ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2]), make_param(**gate_structure[3])]))
            if gate_structure[0] is Dgate:
                dgates.append(ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2]), make_param(**gate_structure[3])]))
            if gate_structure[0] is Kgate:
                kgates.append(ParametrizedGate(gate_structure[0], gate_structure[1], [make_param(**gate_structure[2])]))

        eng, q = sf.Engine(self.n_qumodes)

        rl, U = takagi(self.adj_matrix)
        initial_squeezings = np.arctanh(rl)

        with eng:
            for i ,squeeze_value in enumerate(initial_squeezings):
                Sgate(squeeze_value) | i

            Interferometer(U) | q

            for gate in sgates:
                gate.gate(gate.params[0], gate.params[1]) | gate.qumodes

            Interferometer(self.interferometer_matrix) | q

            for gate in dgates:

# displacement magnitudes (assume displacement is real for now)
    d_r = params[6:9]
    # beamsplitter theta
    bs_theta1, bs_theta2, bs_theta3 = params[9:12]
    # beamsplitter phi
    bs_phi1, bs_phi2, bs_phi3 = params[12:]

    # quantum circuit prior to entering the beamsplitter
    prog = sf.Program(3)

    with prog.context as q:
        for k in range(3):
            Sgate(sq_r[k], sq_phi[k]) | q[k]
            Dgate(d_r[k]) | q[k]

    eng = sf.Engine("fock", backend_options={"cutoff_dim": cutoff})
    stateIn = eng.run(prog).state
    normIn = np.abs(stateIn.trace())

    # norm of output state and probability
    prog_BS = sf.Program(3)
    with prog_BS.context as q:
        BSgate(bs_theta1, bs_phi1) | (q[0], q[1])
        BSgate(bs_theta2, bs_phi2) | (q[1], q[2])
        BSgate(bs_theta3, bs_phi3) | (q[0], q[1])

    stateOut = eng.run(prog_BS).state
    normOut = np.abs(stateOut.trace())
    rho = stateOut.dm()

    # probability of meausring m1 and m2
    prob = np.abs(np.trace(rho[m1, m1, m2, m2]))

phi = make_param(name='phi', stdev=0.2, shape=[depth], regularize=False)
    theta = make_param(name='theta', stdev=0.2, shape=[depth], regularize=False)
    a = make_param(name='a', stdev=0.2, shape=[depth], regularize=True, monitor=True)
    rtheta = make_param(name='rtheta', stdev=0.2, shape=[depth], regularize=False, monitor=True)
    r = make_param(name='r', stdev=0.2, shape=[depth], regularize=True, monitor=True)
    kappa = make_param(name='kappa', stdev=0.2, shape=[depth], regularize=True, monitor=True)

    # Define a single layer of gates
    def layer(l):
        BSgate(phi[l], theta[l]) | (q[0], q[1])
        Dgate(a[l]) | q[0]
        Rgate(rtheta[l]) | q[0]
        Sgate(r[l]) | q[0]
        Kgate(kappa[l]) | q[0]

    eng, q = sf.Engine(2)

    with eng:
        # Make depth layers
        for d in range(depth):
            layer(d)

    state = eng.run('tf', cutoff_dim=7, eval=False)
    circuit_output = state.all_fock_probs()

    return circuit_output

#!/usr/bin/env python3
import strawberryfields as sf
from strawberryfields.ops import *
from numpy import pi

# initialize engine and program objects
eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 5})
ham_simulation = sf.Program(2)

# set the Hamiltonian parameters
J = 1           # hopping transition
U = 1.5         # on-site interaction
k = 20          # Lie product decomposition terms
t = 1.086       # timestep
theta = -J*t/k  
r = -U*t/(2*k)

with ham_simulation.context as q:
    # prepare the initial state
    Fock(2) | q[0]

    # Two node tight-binding
    # Hamiltonian simulation