How to use the strawberryfields.program_utils.Command function in StrawberryFields

sq, U, V = dec.bipartite_graph_embed(B, mean_photon_per_mode=mean_photon_per_mode, atol=tol, rtol=0)

        if not self.identity or not drop_identity:
            for m, s in enumerate(sq):
                s = s if np.abs(s) >= _decomposition_tol else 0

                if not (drop_identity and s == 0):
                    cmds.append(Command(S2gate(-s), (reg[m], reg[m+N])))

            for X, _reg in ((U, reg[:N]), (V, reg[N:])):

                if np.allclose(X, np.identity(len(X)), atol=_decomposition_tol, rtol=0):
                    X = np.identity(len(X))

                if not (drop_identity and np.all(X == np.identity(len(X)))):
                    cmds.append(Command(Interferometer(X, mesh=mesh, drop_identity=drop_identity, tol=tol), _reg))

        return cmds

def _decompose(self, reg, **kwargs):
        # two opposite squeezers sandwiched between 50% beamsplitters
        S = Sgate(self.p[0], self.p[1])
        BS = BSgate(np.pi / 4, 0)
        return [
            Command(BS, reg),
            Command(S, reg[0]),
            Command(S.H, reg[1]),
            Command(BS.H, reg)
        ]

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])))
                    if not (drop_identity and phi == 0):
                        cmds.append(Command(Rgate(-phi), reg[n]))

        return cmds

Args:
        prog (Program): quantum program
        cutoff_dim (int): the Fock basis truncation

    Returns:
        Program: modified program
    """
    prog = copy.deepcopy(prog)
    prog.locked = False  # unlock the copy so we can modify it
    N = prog.init_num_subsystems
    prog._add_subsystems(N)  # pylint: disable=protected-access
    prog.init_num_subsystems = 2 * N
    I = _interleaved_identities(N, cutoff_dim)
    # prepend the circuit with the I ket preparation
    prog.circuit.insert(0, Command(Ket(I), list(prog.reg_refs.values())))
    return prog

D = np.diag(V)
        is_diag = np.all(V == np.diag(D))

        BD = changebasis(self.ns) @ V @ changebasis(self.ns).T
        BD_modes = [BD[i*2:(i+1)*2, i*2:(i+1)*2]
                    for i in range(BD.shape[0]//2)]
        is_block_diag = (not is_diag) and np.all(BD == block_diag(*BD_modes))

        if self.pure and is_diag:
            # covariance matrix consists of x/p quadrature squeezed state
            for n, expr in enumerate(D[:self.ns]):
                if np.abs(expr - 1) >= _decomposition_tol:
                    r = np.abs(np.log(expr)/2)
                    cmds.append(Command(Squeezed(r, 0), reg[n]))
                else:
                    cmds.append(Command(Vac, reg[n]))

        elif self.pure and is_block_diag:
            # covariance matrix consists of rotated squeezed states
            for n, v in enumerate(BD_modes):
                if not np.all(v - np.identity(2) < _decomposition_tol):
                    r = np.abs(np.arccosh(np.sum(np.diag(v)) / 2)) / 2
                    phi = np.arctan(2 * v[0, 1] / np.sum(np.diag(v) * [1, -1]))
                    cmds.append(Command(Squeezed(r, phi), reg[n]))
                else:
                    cmds.append(Command(Vac, reg[n]))

        elif not self.pure and is_diag and np.all(D[:self.ns] == D[self.ns:]):
            # covariance matrix consists of thermal states
            for n, nbar in enumerate(0.5 * (D[:self.ns] - 1.0)):
                if nbar >= _decomposition_tol:
                    cmds.append(Command(Thermal(nbar), reg[n]))

a = q[i]
                b = q[i+1]
                # the ops must have equal size and act on the same wires
                if a.op.ns == b.op.ns and a.reg == b.reg:
                    if a.op.ns != 1:
                        # ns > 1 is tougher. on no wire must there be anything
                        # between them, also deleting is more complicated
                        # todo treat it as a failed merge for now
                        i += 1
                        continue
                    op = a.op.merge(b.op)
                    # merge was successful, delete the old ops
                    del q[i:i+2]
                    # insert the merged op (unless it's identity)
                    if op is not None:
                        q.insert(i, Command(op, a.reg))
                    # move one spot backwards to try another merge
                    if i > 0:
                        i -= 1
                    _print_list(i, q)
                    continue
            except MergeFailure:
                pass
            i += 1  # failed at merging the ops, move forward

    # convert the circuit back into a list (via a DAG)
    DAG = grid_to_DAG(grid)
    return DAG_to_list(DAG)

# 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])))
                    if not (drop_identity and phi == 0):
                        cmds.append(Command(Rgate(-phi), reg[n]))

        return cmds

# multiply all unitaries together, to get the final
        # 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

def _decompose(self, reg, **kwargs):
        # into local phase shifts and two 50-50 beamsplitters
        return [
            Command(Rgate(self.p[0]), reg[0]),
            Command(BSgate(np.pi/4, np.pi/2), reg),
            Command(Rgate(self.p[1]), reg[0]),
            Command(BSgate(np.pi/4, np.pi/2), reg)
        ]

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