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

def test_interferometers(self, tol):
        """Test that the compilation correctly decomposes the interferometer using
        the rectangular_symmetric mesh"""
        prog = sf.Program(8)
        U = random_interferometer(4)

        with prog.context as q:
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[0], q[4])
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[1], q[5])
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[2], q[6])
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[3], q[7])
            ops.Interferometer(U) | (q[0], q[1], q[2], q[3])
            ops.Interferometer(U) | (q[4], q[5], q[6], q[7])
            ops.MeasureFock() | q

        res = prog.compile("chip2")

        expected = sf.Program(8)

        with expected.context as q:
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[0], q[4])
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[1], q[5])
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[2], q[6])
            ops.S2gate(SQ_AMPLITUDE, 0) | (q[3], q[7])
            ops.Interferometer(U, mesh="rectangular_symmetric", drop_identity=False) | (q[0], q[1], q[2], q[3])
            ops.Interferometer(U, mesh="rectangular_symmetric", drop_identity=False) | (q[4], q[5], q[6], q[7])
            ops.MeasureFock() | q

        expected = expected.compile(DummyCircuit())

def test_gaussian_program(depth, width):
    """Tests that a circuit and its compiled version produce the same Gaussian state"""
    eng = sf.LocalEngine(backend="gaussian")
    eng1 = sf.LocalEngine(backend="gaussian")
    circuit = sf.Program(width)
    with circuit.context as q:
        for _ in range(depth):
            U, s, V, alphas = random_params(width, 2.0 / depth, 1.0)
            ops.Interferometer(U) | q
            for i in range(width):
                ops.Sgate(s[i]) | q[i]
            ops.Interferometer(V) | q
            for i in range(width):
                ops.Dgate(alphas[i]) | q[i]
    compiled_circuit = circuit.compile("gaussian_unitary")
    cv = eng.run(circuit).state.cov()
    mean = eng.run(circuit).state.means()

    cv1 = eng1.run(compiled_circuit).state.cov()
    mean1 = eng1.run(compiled_circuit).state.means()
    assert np.allclose(cv, cv1)
    assert np.allclose(mean, mean1)

def test_merge(self, tol):
        """Test that two interferometers merge: U = U1 @ U2"""
        n = 3
        U1 = random_interferometer(n)
        U2 = random_interferometer(n)

        int1 = ops.Interferometer(U1)
        int1inv = ops.Interferometer(U1.conj().T)
        int2 = ops.Interferometer(U2)

        # an interferometer merged with its inverse is identity
        assert int1.merge(int1inv) is None

        # two merged unitaries are the same as their product
        assert np.allclose(int1.merge(int2).p[0].x, U2 @ U1, atol=tol, rtol=0)

def test_no_s2gates(self, tol):
        """Test identity S2gates are inserted when no S2gates
        are provided."""
        prog = sf.Program(8)
        U = random_interferometer(4)

        with prog.context as q:
            ops.Interferometer(U) | (q[0], q[1], q[2], q[3])
            ops.Interferometer(U) | (q[4], q[5], q[6], q[7])
            ops.MeasureFock() | q

        expected = sf.Program(8)

        with expected.context as q:
            ops.S2gate(0) | (q[0], q[4])
            ops.S2gate(0) | (q[1], q[5])
            ops.S2gate(0) | (q[2], q[6])
            ops.S2gate(0) | (q[3], q[7])
            ops.Interferometer(U) | (q[0], q[1], q[2], q[3])
            ops.Interferometer(U) | (q[4], q[5], q[6], q[7])
            ops.MeasureFock() | q

        res = prog.compile("chip2")
        expected = expected.compile("chip2")

def test_gaussian_program(depth, width):
    """Tests that a circuit and its compiled version produce the same Gaussian state"""
    eng = sf.LocalEngine(backend="gaussian")
    eng1 = sf.LocalEngine(backend="gaussian")
    circuit = sf.Program(width)
    with circuit.context as q:
        for _ in range(depth):
            U, s, V, alphas = random_params(width, 2.0 / depth, 1.0)
            ops.Interferometer(U) | q
            for i in range(width):
                ops.Sgate(s[i]) | q[i]
            ops.Interferometer(V) | q
            for i in range(width):
                ops.Dgate(alphas[i]) | q[i]
    compiled_circuit = circuit.compile("gaussian_unitary")
    cv = eng.run(circuit).state.cov()
    mean = eng.run(circuit).state.means()

    cv1 = eng1.run(compiled_circuit).state.cov()
    mean1 = eng1.run(compiled_circuit).state.means()
    assert np.allclose(cv, cv1)
    assert np.allclose(mean, mean1)

# unitary representation on modes [0, 1] and [2, 3].
        U01 = multi_dot(U_list01)
        U23 = multi_dot(U_list23)

        # check unitaries are equal
        if not np.allclose(U01, U23):
            raise CircuitError(
                "Interferometer on modes [0, 1] must be identical to interferometer on modes [2, 3]."
            )

        U = block_diag(U01, U23)

        # replace B with an interferometer
        B = [
            Command(ops.Interferometer(U01), registers[:2]),
            Command(ops.Interferometer(U23), registers[2:]),
        ]

        # decompose the interferometer, using Mach-Zehnder interferometers
        B = self.decompose(B)

        # Do a final circuit topology check
        # ---------------------------------
        seq = super().compile(A + B + C, registers)
        return seq

def _decompose(self, reg, **kwargs):
        cmds = []
        mesh = kwargs.get("mesh", "rectangular")

        if self.active:
            if not self.vacuum:
                cmds = [Command(Interferometer(self.U2), reg)]

            for n, expr in enumerate(self.Sq):
                if np.abs(expr - 1) >= _decomposition_tol:
                    r = np.abs(np.log(expr))
                    phi = np.angle(np.log(expr))
                    cmds.append(Command(Sgate(-r, phi), reg[n]))

            cmds.append(Command(Interferometer(self.U1, mesh=mesh), reg))
        else:
            if not self.vacuum:
                cmds = [Command(Interferometer(self.U1, mesh=mesh), reg)]

        return cmds

_operation_map = {
        'CoherentState': Coherent,
        'DisplacedSqueezedState': DisplacedSqueezed,
        'SqueezedState': Squeezed,
        'ThermalState': Thermal,
        'GaussianState': Gaussian,
        'Beamsplitter': BSgate,
        'ControlledAddition': CXgate,
        'ControlledPhase': CZgate,
        'Displacement': Dgate,
        'QuadraticPhase': Pgate,
        'Rotation': Rgate,
        'TwoModeSqueezing': S2gate,
        'Squeezing': Sgate,
        'Interferometer': Interferometer
    }

    _observable_map = {
        'NumberOperator': mean_photon,
        'TensorN': number_expectation,
        'X': homodyne(0),
        'P': homodyne(np.pi/2),
        'QuadOperator': homodyne(),
        'PolyXP': poly_xp,
        'FockStateProjector': fock_state,
        'Identity': identity
    }

    _circuits = {}

    def __init__(self, wires, *, analytic=True, cutoff_dim=10, shots=1000, hbar=2):

# multiply all unitaries together, to get the final
        # unitary representation on modes [0, 1] and [2, 3].
        U01 = multi_dot(U_list01)
        U23 = multi_dot(U_list23)

        # check unitaries are equal
        if not np.allclose(U01, U23):
            raise CircuitError(
                "Interferometer on modes [0, 1] must be identical to interferometer on modes [2, 3]."
            )

        U = block_diag(U01, U23)

        # replace B with an interferometer
        B = [
            Command(ops.Interferometer(U01), registers[:2]),
            Command(ops.Interferometer(U23), registers[2:]),
        ]

        # decompose the interferometer, using Mach-Zehnder interferometers
        B = self.decompose(B)

        # Do a final circuit topology check
        # ---------------------------------
        seq = super().compile(A + B + C, registers)
        return seq

# unitary representation on modes [0, 1] and [2, 3].
        U0 = multi_dot(U_list0)
        U4 = multi_dot(U_list4)

        # check unitaries are equal
        if not np.allclose(U0, U4):
            raise CircuitError(
                "Interferometer on modes [0, 1, 2, 3] must be identical to interferometer on modes [4, 5, 6, 7]."
            )

        U = block_diag(U0, U4)

        # replace B with an interferometer
        B = [
            Command(ops.Interferometer(U0), registers[:4]),
            Command(ops.Interferometer(U4), registers[4:]),
        ]

        # decompose the interferometer, using Mach-Zehnder interferometers
        B = self.decompose(B)

        # Do a final circuit topology check
        # ---------------------------------
        seq = super().compile(A + B + C, registers)
        return seq