How to use the strawberryfields.ops.BSgate function in StrawberryFields

def setup_two_mode_circuit(self, setup_eng, cutoff):
        """Create the circuit for following tests"""
        eng_ref, p0 = setup_eng(2)

        S = ops.Sgate(2)
        B = ops.BSgate(2.234, -1.165)

        initial_state = np.complex64(
            np.random.rand(*[cutoff] * 4) + 1j * np.random.rand(*[cutoff] * 4)
        )

        with p0.context as q:
            ops.DensityMatrix(initial_state) | q

        prog = sf.Program(p0)
        with prog.context as q:
            S | q[0]
            B | q
            S | q[1]
            B | q

        rho = eng_ref.run([p0, prog]).state.dm()

def single_input_circuit(x):
                eng.reset()
                with eng:
                    Dgate(x[0], 0.) | q[0]
                    Dgate(x[1], 0.) | q[1]
                    BSgate(phi=params[0]) | (q[0], q[1])
                    BSgate() | (q[0], q[1])
                state = eng.run('fock', cutoff_dim=10, eval=True)
                p0 = state.fock_prob([0, 2])
                p1 = state.fock_prob([2, 0])
                normalisation = p0 + p1 + 1e-10
                outp = p1 / normalisation
                return outp

def test_is_unitary_with_channel(self, prog):
        """test that the is_unitary function returns False if channels are present"""
        with prog.context as q:
            Sgate(0.4) | q[0]
            LossChannel(0.4) | q[0]
            BSgate(0.4) | q
        assert not utils.is_unitary(prog)

def test_S2gate_decomp_equal(self, setup_eng, r, tol):
        """Tests that the S2gate gives the same transformation as its decomposition."""
        eng, prog = setup_eng(2)

        phi = 0.273
        BS = ops.BSgate(np.pi / 4, 0)
        with prog.context as q:
            ops.S2gate(r, phi) | q
            # run decomposition with reversed arguments
            BS | q
            ops.Sgate(-r, phi) | q[0]
            ops.Sgate(r, phi) | q[1]
            BS.H | q

        eng.run(prog)
        assert np.all(eng.backend.is_vacuum(tol))

def dummy_func(r, theta, phi, q):
            Sgate(r) | q[0]
            BSgate(theta, phi) | (q[0], q[1])

def single_input_circuit(x):

        eng.reset()
        with eng:
            Dgate(x[0], 0.) | q[0]
            Dgate(x[1], 0.) | q[1]
            BSgate(phi=params[0]) | (q[0], q[1])
            BSgate() | (q[0], q[1])
        state = eng.run('fock', cutoff_dim=10, eval=True)

        p0 = state.fock_prob([0, 2])
        p1 = state.fock_prob([2, 0])
        normalization = p0 + p1 + 1e-10
        output = p1 / normalization
        return output

# 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]))

    # output state
    rhoC = rho[m1, m1, m2, m2]/prob

    #fidelity with the target
    fidelity = np.abs(np.trace(np.einsum('ij,jk->ik', rhoC, ONdm)))
    return (fidelity, prob, normIn, normOut, rhoC)

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]))

    # output state
    rhoC = rho[m1, m1, m2, m2]/prob

    #fidelity with the target
    fidelity = np.abs(np.trace(np.einsum('ij,jk->ik', rhoC, ONdm)))
    return (fidelity, prob, normIn, normOut, rhoC)

for n, m, theta, phi, _ in BS1:
                theta = theta if np.abs(theta) >= _decomposition_tol else 0
                phi = phi if np.abs(phi) >= _decomposition_tol else 0

                if "symmetric" in mesh:
                    # Mach-Zehnder interferometers
                    cmds.append(Command(MZgate(np.mod(phi, 2*np.pi), np.mod(theta, 2*np.pi)), (reg[n], reg[m])))

                else:
                    # Clements style beamsplitters
                    if not (drop_identity and phi == 0):
                        cmds.append(Command(Rgate(phi), reg[n]))

                    if not (drop_identity and theta == 0):
                        cmds.append(Command(BSgate(theta, 0), (reg[n], reg[m])))

            for n, expphi in enumerate(R):
                # local phase shifts
                q = np.log(expphi).imag if np.abs(expphi - 1) >= _decomposition_tol else 0
                if not (drop_identity and q == 0):
                    cmds.append(Command(Rgate(q), reg[n]))

            if BS2 is not None:
                # Clements style beamsplitters

                for n, m, theta, phi, _ in reversed(BS2):
                    theta = theta if np.abs(theta) >= _decomposition_tol else 0
                    phi = phi if np.abs(phi) >= _decomposition_tol else 0

                    if not (drop_identity and theta == 0):
                        cmds.append(Command(BSgate(-theta, 0), (reg[n], reg[m])))

def decompose(self, reg):
        # make BS gate
        theta = self.layer['BS'][0]
        phi = self.layer['BS'][1]
        BS = BSgate(theta, phi)

        # make cross-Kerr gate
        CK = None
        param = self.layer.get('CK', [0])[0]
        if param != 0:
            CK = CKgate(param)

        # make Kerr gate
        K = None
        param = self.layer.get('K', [0])[0]
        if param != 0:
            K = Kgate(param)

        # make rotation gate
        R = None
        param = self.layer.get('R', [0])[0]