submodule of graphix to run measurement pattern on discrete-variable optical QPU
from .converter import to_perceval

graphix_perceval/__init__.py
TeamGraphix/graphix-perceval
import numpy as np
import perceval.components as comp
import sympy as sp

# perceval representation of Clifford gates.
# see graphix.clifford module for the definitions and details of Clifford operatos for each index.
CLIFFORD_TO_PERCEVAL_BS = [
    [comp.BS(theta=0.0)],
    [comp.BS(theta=sp.pi, phi_bl=-sp.pi / 2, phi_br=-sp.pi / 2)],
    [comp.BS(theta=sp.pi, phi_tl=-sp.pi / 2, phi_bl=sp.pi / 2, phi_tr=sp.pi / 2, phi_br=sp.pi / 2)],
    [comp.BS(theta=0.0, phi_bl=sp.pi / 2, phi_br=sp.pi / 2)],
    [comp.BS(theta=0.0, phi_br=sp.pi / 2)],
    [comp.BS(theta=0.0, phi_br=-sp.pi / 2)],
    [comp.BS(theta=sp.pi / 2, phi_bl=3 * sp.pi / 2, phi_br=3 * sp.pi / 2)],
    [comp.BS(theta=sp.pi / 2, phi_tl=sp.pi / 2, phi_bl=-sp.pi / 2, phi_tr=-sp.pi / 2, phi_br=sp.pi / 2)],
    [comp.BS(theta=sp.pi / 2, phi_bl=sp.pi / 2, phi_br=3 * sp.pi / 2)],
    [comp.BS(theta=sp.pi, phi_tl=3 * sp.pi / 4, phi_tr=-3 * sp.pi / 4)],
    [comp.BS(theta=sp.pi, phi_tl=-3 * sp.pi / 4, phi_tr=3 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_bl=sp.pi / 2, phi_br=sp.pi / 2)],
    [comp.BS(theta=sp.pi / 2, phi_tl=3 * sp.pi / 4, phi_bl=sp.pi / 4, phi_tr=sp.pi / 4, phi_br=3 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_tr=sp.pi / 2, phi_br=3 * sp.pi / 2)],
    [comp.BS(theta=sp.pi / 2, phi_tl=sp.pi / 2, phi_bl=3 * sp.pi / 2)],
    [comp.BS(theta=sp.pi / 2, phi_tr=sp.pi, phi_br=3 * sp.pi)],
    [comp.BS(theta=sp.pi / 2, phi_bl=sp.pi, phi_tr=3 * sp.pi / 4, phi_br=sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_tr=3 * sp.pi / 4, phi_br=5 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_bl=sp.pi, phi_tr=sp.pi / 4, phi_br=3 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_tr=5 * sp.pi / 4, phi_br=3 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_tl=5 * sp.pi / 4, phi_bl=3 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_tl=sp.pi / 2, phi_tr=sp.pi / 4, phi_br=5 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_bl=sp.pi / 2, phi_tr=sp.pi / 4, phi_br=5 * sp.pi / 4)],
    [comp.BS(theta=sp.pi / 2, phi_tl=3 * sp.pi / 4, phi_bl=5 * sp.pi / 4)],
]

CLIFFORD_TO_PERCEVAL_POLAR = [
    [comp.WP(delta=0.0, xsi=0.0)],
    [comp.WP(delta=sp.pi / 2, xsi=sp.pi / 4), comp.PS(-sp.pi / 2)],
    [comp.WP(delta=sp.pi / 2, xsi=0.0), comp.WP(delta=sp.pi / 2, xsi=sp.pi / 4), comp.PS(-sp.pi / 2)],
    [comp.WP(delta=sp.pi / 2, xsi=0.0), comp.PS(-sp.pi / 2)],
    [comp.WP(delta=-sp.pi / 4, xsi=0.0), comp.PS(sp.pi / 4)],
    [comp.WP(delta=sp.pi / 4, xsi=0.0), comp.PS(7 * sp.pi / 4)],
    [comp.WP(delta=sp.pi / 2, xsi=np.pi / 8), comp.PS(3 * sp.pi / 2)],
    [comp.WP(delta=3 * sp.pi / 4, xsi=np.pi / 4), comp.PS(sp.pi)],
    [comp.WP(delta=sp.pi / 2, xsi=sp.pi / 8), comp.WP(delta=sp.pi / 2, xsi=sp.pi / 4), comp.PS(sp.pi)],
    [comp.WP(delta=-sp.pi / 4, xsi=sp.pi), comp.WP(delta=sp.pi / 2, xsi=sp.pi / 4), comp.PS(sp.pi)],
    [comp.WP(delta=sp.pi / 4, xsi=sp.pi), comp.WP(delta=sp.pi / 2, xsi=sp.pi / 4), comp.PS(-sp.pi)],
    [comp.WP(delta=sp.pi / 2, xsi=3 * sp.pi / 8), comp.PS(sp.pi / 2)],
    [comp.WP(delta=sp.pi / 2, xsi=3 * sp.pi / 8), comp.WP(delta=sp.pi / 2, xsi=sp.pi / 4)],
    [comp.WP(delta=sp.pi / 4, xsi=sp.pi / 4), comp.WP(delta=sp.pi / 2, xsi=0.0)],
    [comp.WP(delta=sp.pi / 4, xsi=3 * sp.pi / 4), comp.WP(delta=sp.pi / 2, xsi=0.0)],
    [comp.WP(delta=sp.pi / 4, xsi=sp.pi / 4), comp.PS(sp.pi)],
    [comp.WP(delta=sp.pi / 4, xsi=0.0), comp.WP(delta=sp.pi / 2, xsi=sp.pi / 8)],
    [comp.WP(delta=sp.pi / 4, xsi=0.0), comp.WP(delta=sp.pi / 2, xsi=3 * sp.pi / 8), comp.PS(sp.pi)],
    [comp.WP(delta=sp.pi / 4, xsi=sp.pi / 2), comp.WP(delta=sp.pi / 2, xsi=3 * sp.pi / 8), comp.PS(sp.pi)],
    [comp.WP(delta=sp.pi / 4, xsi=sp.pi / 2), comp.WP(delta=sp.pi / 2, xsi=5 * sp.pi / 8)],
    [comp.WP(delta=sp.pi / 4, xsi=0.0), comp.WP(delta=sp.pi / 4, xsi=sp.pi / 4), comp.PS(sp.pi)],
    [comp.WP(delta=sp.pi / 4, xsi=sp.pi / 4), comp.WP(delta=sp.pi / 2, xsi=sp.pi / 8)],
    [comp.WP(delta=sp.pi / 4, xsi=0.0), comp.WP(delta=sp.pi / 4, xsi=3 * sp.pi / 4)],
    [comp.WP(delta=sp.pi / 4, xsi=sp.pi / 2), comp.WP(delta=sp.pi / 4, xsi=sp.pi / 4), comp.PS(sp.pi)],
]

graphix_perceval/clifford.py
TeamGraphix/graphix-perceval
from __future__ import annotations

import graphix
import perceval as pcvl
import sympy as sp
from graphix.extraction import ResourceGraph, ResourceType, get_fusion_network_from_graph
from perceval import components as comp

from graphix_perceval.clifford import CLIFFORD_TO_PERCEVAL_POLAR
from graphix_perceval.experiment import PercevalExperiment, Photon, PhotonType


def pattern2graphstate(
    pattern: graphix.Pattern,
) -> tuple[graphix.GraphState, dict[int, float], list[int]]:
    """Create a graph state from a MBQC pattern.

    Parameters
    ----------
    pattern : :class:`graphix.Pattern` object
        MBQC pattern to be run on the device

    Returns
    -------
    graph_state : :class:`graphix.GraphState` object
        Graph state corresponding to the pattern.
    phasedict : dict
        Dictionary of phases for each node.
    output_nodes : list
        List of output nodes.
    """
    nodes, edges = pattern.get_graph()
    vop_init = pattern.get_vops()
    graph_state = graphix.GraphState(nodes=nodes, edges=edges, vops=vop_init)
    phasedict = {}
    for command in pattern.get_measurement_commands():
        phasedict[command[1]] = command[3]

    output_nodes = pattern.output_nodes
    return graph_state, phasedict, output_nodes


def to_perceval(pattern: graphix.Pattern) -> PercevalExperiment:
    """Convert a graphix.Pattern to a perceval.Circuit.

    Parameters
    ----------
    pattern : graphix.Pattern
        GraphState to be converted to a perceval.Circuit

    Returns
    -------
    experiment : PercevalExperiment
        :class:`graphix_perceval.experiment.PercevalExperiment` object
    """
    if not isinstance(pattern, graphix.Pattern):
        raise TypeError("pattern must be a graphix.Pattern object")
    graph_state, phasedict, output_nodes = pattern2graphstate(pattern)
    ResourceGraphs = get_fusion_network_from_graph(graph_state)
    vops = pattern.get_vops()

    pcc = PercevalCircuitConstructor()
    for ResourceGraph in ResourceGraphs:
        pcc.add_ResourceGraph(ResourceGraph, phasedict, output_nodes)

    pcc.add_fusions()

    pcc.apply_local_clifford(vops)

    perceval_circuit = pcc.setup_perceval_circuit()

    exp = PercevalExperiment(perceval_circuit, pcc.photons)

    return exp


class PercevalCircuitConstructor:
    def __init__(self):
        self.num_photons = 0
        self.ResourceGraphs: list[ResourceGraph] = []
        self.fusion_pairs: list[tuple[Photon, Photon]] = []
        self.photons: list[Photon] = []
        self.node_id2photon_ids: dict[int, list[int]] = {}
        self.clifford_comps: list[comp.BS] = []
        self._is_fused: bool = False
        self._clifford_applied: bool = False

    def add_ResourceGraph(self, ResourceGraph: ResourceGraph, phasedict: dict[int, float], readouts: list) -> None:
        if self._is_fused:
            raise RuntimeError("Cannot add ResourceGraph after fusion")
        for node_id in ResourceGraph.graph.nodes:
            if node_id in readouts:
                ph = Photon(
                    exp_id=self.num_photons, type=PhotonType.READOUT, node_id=node_id, angle=phasedict.get(node_id)
                )
            else:
                ph = Photon(
                    exp_id=self.num_photons, type=PhotonType.COMPUTE, node_id=node_id, angle=phasedict.get(node_id)
                )
            if self.node_id2photon_ids.get(node_id) is not None:
                self.node_id2photon_ids[node_id].append(self.num_photons)
            else:
                self.node_id2photon_ids[node_id] = [self.num_photons]
            self.num_photons += 1
            self.photons.append(ph)

        if ResourceGraph.type in (ResourceType.GHZ, ResourceType.LINEAR):
            self.ResourceGraphs.append(ResourceGraph)
        else:
            raise TypeError(f"ResourceType {ResourceGraph.type} is not supported")

    def get_readouts(self) -> list[Photon]:
        return [ph for ph in self.photons if ph.type == PhotonType.READOUT]

    def get_computes(self) -> list[Photon]:
        return [ph for ph in self.photons if ph.type == PhotonType.COMPUTE]

    def get_witnesses(self) -> list[Photon]:
        return [ph for ph in self.photons if ph.type == PhotonType.WITNESS]

    def add_fusions(self) -> None:
        """Find edges that connects two ResourceGraphs.
        If the two ResourceGraphs share a same node id, then they are fused.
        Type-1 fusion.
        """
        for _, photon_ids in self.node_id2photon_ids.items():
            fusing_photons = sorted(photon_ids)
            for idx in range(len(fusing_photons) - 1):
                self.fusion_pairs.append((self.photons[fusing_photons[idx]], self.photons[fusing_photons[idx + 1]]))
                # Note that the photon with the larger index is the witness
                self.photons[fusing_photons[idx + 1]].type = PhotonType.WITNESS

        self._is_fused = True

    def get_all_ResourceGraphs(self) -> list[ResourceGraph]:
        return self.ghz_ResourceGraphs | self.linear_ResourceGraphs

    def setup_perceval_circuit(self, name: str | None = None, merge: bool = False) -> pcvl.Circuit:
        if not self._is_fused:
            raise RuntimeError("Must fuse before setting up perceval circuit")
        if not self._clifford_applied:
            raise RuntimeError("Must apply local clifford before setting up perceval circuit")
        circ = pcvl.Circuit(self.num_photons * 2, name=name)
        # Create circuits for all the ResourceGraphs
        photon_idx = 0
        for cl in self.ResourceGraphs:
            if cl.type == ResourceType.GHZ:
                circ.add(
                    [idx for idx in range(photon_idx, photon_idx + len(cl.graph.nodes))],
                    ghz_circuit(len(cl.graph.nodes)),
                    merge,
                )
            elif cl.type == ResourceType.LINEAR:
                circ.add(
                    [idx for idx in range(photon_idx, photon_idx + len(cl.graph.nodes))],
                    linear_circuit(len(cl.graph.nodes)),
                    merge,
                )
            photon_idx += len(cl.graph.nodes)

        circ.add(0, comp.PERM(list(range(self.num_photons))))  # work as a barrier

        # Create circuits for all the Fusions
        for ph1, ph2 in self.fusion_pairs:
            circ.add(list(range(ph1.id, ph2.id + 1)), fusion_circuit(ph1, ph2), merge)

        circ.add(0, comp.PERM(list(range(self.num_photons))))  # work as a barrier

        # Add local clifford
        for photon_id, clifford_comp in self.clifford_comps:
            circ.add(photon_id, clifford_comp)

        circ.add(0, comp.PERM(list(range(self.num_photons * 2))))  # work as a barrier

        # Convert measurement basis
        for ph in self.photons:
            if ph.type == PhotonType.COMPUTE:
                circ.add(ph.id, comp.QWP(ph.angle[0]))
                circ.add(ph.id, comp.HWP(ph.angle[1]))

        # Currently, Perceval does not support measurement in polarization, so we need to convert it to dual-rail encoding.
        # |{P:H},0> -> |0,1> = |0> (this is the opposite of the definition in perceval)
        # |{P:V},0> -> |1,0> = |1> (this is the opposite of the definition in perceval)
        circ.add(
            0,
            comp.PERM(
                sum(
                    list([2 * i] for i in range(0, self.num_photons))
                    + list([2 * i + 1] for i in range(0, self.num_photons)),
                    [],
                )
            ),
        )
        for i in range(0, self.num_photons):
            circ.add(i * 2, comp.PBS())

        return circ

    def apply_local_clifford(self, vops: dict[int, int]) -> None:
        if not self._is_fused:
            raise RuntimeError("Must fuse before applying local clifford")
        for node_id, cid in vops.items():
            for ph_id in self.node_id2photon_ids[node_id]:
                if self.photons[ph_id].type != PhotonType.WITNESS:
                    self.clifford_comps.append((ph_id, local_clifford_circuit(cid)))

        self._clifford_applied = True


def local_clifford_circuit(clifford_id: int) -> pcvl.Circuit:
    """Create a Perceval Circuit for a local clifford.

    Parameters
    ----------
    mode_id : int
        Mode id.
    clifford_id : int
        Clifford id.

    Returns
    -------
    perceval.Circuit
        Perceval Circuit for a local clifford.
    """
    if not 0 <= clifford_id <= 23:
        raise ValueError("clifford_id must be in [0, 23]")
    circ = pcvl.Circuit(m=1, name="LOCAL CLIFFORD ID:" + str(clifford_id))
    for comps in CLIFFORD_TO_PERCEVAL_POLAR[clifford_id]:
        circ.add(0, comps)
    return circ


def fusion_circuit(ph1: Photon, ph2: Photon) -> pcvl.Circuit:
    """Create a Perceval Circuit for fusing two photons.

    Parameters
    ----------
    ph1 : Photon
        First photon.
    ph2 : Photon
        Second photon.

    Returns
    -------
    perceval.Circuit
        Perceval Circuit for fusing two photons.
    """
    if not isinstance(ph1, Photon) or not isinstance(ph2, Photon):
        raise TypeError("ph1 and ph2 must be Photon objects")
    if ph1.type == PhotonType.WITNESS and ph1.id > ph2.id:
        ph1, ph2 = ph2, ph1
    if ph2.type != PhotonType.WITNESS:
        raise ValueError("The second photon must be a witness")
    l = ph2.id - ph1.id
    circ = pcvl.Circuit(m=l + 1, name="FUSE " + str(ph1.id) + "-" + str(ph2.id))
    # If the photons are not neighbors, we swap the ph2 and the photon next to ph1,
    # do the fusion and swap back.
    if l > 1:
        a, *b, c = list(range(0, l))
        perm = comp.PERM([c, *b, a])
        circ.add(1, perm)
    circ.add((0, 1), comp.PBS())
    circ.add(1, comp.HWP(sp.pi / 8))
    if l > 1:
        circ.add(1, perm)
    return circ


def linear_circuit(num_photons: int, name: str = "") -> pcvl.Circuit:
    """Create a Perceval Circuit for a linear ResourceGraph.

    Parameters
    ----------
    num_photons : int
        Number of photons.

    Returns
    -------
    perceval.Circuit
        Perceval Circuit for a linear ResourceGraph.
    """
    if not isinstance(num_photons, int):
        raise TypeError("num_photons must be an integer")
    circ = pcvl.Circuit(m=num_photons, name="LINEAR " + name)
    for i in range(num_photons):
        circ.add(i, comp.HWP(sp.pi / 8))

    for i in range(num_photons - 1):
        circ.add((i, i + 1), comp.PBS())
        if i >= 1 and i != num_photons - 2:
            circ.add(i + 1, comp.HWP(sp.pi / 8))
    circ.add(0, comp.PERM(list(range(num_photons))))  # work as a barrier
    circ.add(0, comp.HWP(sp.pi / 8))
    circ.add(num_photons - 1, comp.HWP(sp.pi / 8))

    return circ


def ghz_circuit(num_photons: int, name: str = "") -> pcvl.Circuit:
    """Create a Perceval Circuit for a GHZ ResourceGraph.

    Parameters
    ----------
    num_photons : int
        Number of photons.

    Returns
    -------
    perceval.Circuit
        Perceval Circuit for a GHZ ResourceGraph.
    """
    if not isinstance(num_photons, int):
        raise TypeError("num_photons must be an integer")
    circ = pcvl.Circuit(m=num_photons, name="GHZ " + name)
    for i in range(num_photons):
        circ.add(i, comp.HWP(sp.pi / 8))

    for i in range(num_photons - 1):
        circ.add((i, i + 1), comp.PBS())

    for i in range(1, num_photons):
        circ.add(i, comp.HWP(sp.pi / 8))  # Hadamard

    return circ

graphix_perceval/converter.py
TeamGraphix/graphix-perceval
from __future__ import annotations

import collections
import itertools
import sys
import warnings
from enum import Enum

import perceval as pcvl
import sympy as sp
from _collections_abc import dict_items
from perceval.algorithm import Sampler
from perceval.utils import PostSelect
from tabulate import tabulate

IS_NOTEBOOK = "ipykernel" in sys.modules
if IS_NOTEBOOK:
    from IPython.display import HTML, display  # type: ignore


class PhotonType(Enum):
    READOUT = "READOUT"
    COMPUTE = "COMPUTE"
    WITNESS = "WITNESS"
    LOSS = "LOSS"
    NONE = None

    def __str__(self) -> str:
        return self.name


class Photon:
    def __init__(self, exp_id: int, type: PhotonType, node_id: int = 0, angle: float | None = None):
        self.id = exp_id
        self.type = type
        self.node_id = node_id
        # angle for QWP and HWP
        if angle:
            self.angle = [sp.pi / 4, (sp.pi - (2 * angle * sp.pi)) / 8]
        else:
            self.angle = [sp.pi / 4, sp.pi / 8]  # X-basis measurement

    def __str__(self) -> str:
        return f"Photon(ID:{str(self.id)} Node:{str(self.node_id)} ({str(self.type)}))"

    def __repr__(self) -> str:
        return self.__str__()


class PercevalExperiment:
    """PercevalExperiment class for running MBQC patterns on Perceval simulators and Quandela devices.

    Attributes
    ----------
    pattern: :class:`graphix.Pattern` object
        MBQC pattern to be run on the device
    circ: :class:`perceval.Circuit` object
        Perceval circuit corresponding to the pattern.
    backend : str
        Name of a Perceval simulator or Quandela device
    """

    def __init__(self, circuit: pcvl.Circuit, photons: list[Photon]):
        """

        Parameters
        ----------
        pattern: :class:`graphix.Pattern` object
            MBQC pattern to be run on the Quandela device or Perceval simulator.
        """
        self.circ = circuit
        self.photons = photons
        self.processor = None
        self.input_state = None
        self.output_states: dict[str, str] | None = None

    def set_local_processor(self, backend: str, source: pcvl.Source = pcvl.Source(), name: str = None):
        """Set the local computing backend.

        Parameters
        ----------
        backend : str
            Name of a local backend.
        source : :class:`perceval.Source` object, optional
            Setting of single-photon source.
        name : str, optional
            Name for the processor.
        """
        if self.circ is None:
            warnings.warn("The circuit has not been converted to Perceval circuit. It will be converted automatically.")
            self.to_perceval()
        if self.processor is not None:
            warnings.warn("The processor has already been set. The previous processor will be overwritten.")
        self.processor = pcvl.Processor(backend=backend, m_circuit=self.circ, source=source, name=name)
        self.backend = backend

        self.set_input_state()
        self.set_output_states()

    def set_remote_processor(self, backend: str, token: str):
        """Set the remote computing backend.

        Parameters
        ----------
        backend : str
            Name of a remote backend.
        token : str
            Token for the remote processor.
        """
        if self.circ is None:
            warnings.warn("The circuit has not been converted to Perceval circuit. It will be converted automatically.")
            self.to_perceval()
        if self.processor is not None:
            warnings.warn("The processor has already been set. The previous processor will be overwritten.")
        self.processor = pcvl.RemoteProcessor(name=backend, token=token)
        self.processor.set_circuit(self.circ)
        self.backend = backend

        self.set_input_state()
        self.set_output_states()

    def set_input_state(self):
        """Set the input states for the processor.
        The default input state is |{P:H}> for each photon and |0> for each ancillary mode.
        """
        if self.processor is None:
            raise Exception(
                "No processor has been set. Please set a processor by `set_local_processor` or `set_remote_procesor` before running the experiment."
            )

        input_state = "|"
        input_state = input_state + ",".join([r"{P:H}" for _ in range(len(self.photons))])
        input_state = input_state + "," + ",".join(["0"] * len(self.photons))
        input_state = input_state + ">"

        self.input_state = pcvl.BasicState(input_state)
        self.processor.with_polarized_input(self.input_state)  # not with_input (it will not work for polarized input)

    def set_output_states(self):
        r"""Set the output states.
        Currently, Perceval does not support feed-forward opetations,
        so we postselect the output states where

        - The witness photons are in |{P:H}> and translated to |0,1>
        - The computing photons are in |{P:H}> and translated to |0,1>
        - The readout photons are in |{P:H}> or |{P:V}>
        """
        if self.processor is None:
            raise Exception(
                "No processor has been set. Please set a processor by `set_local_processor` or `set_remote_procesor` before running the experiment."
            )
        (readouts, witnesses, comps) = (
            self.get_readout_photons(),
            self.get_witness_photons(),
            self.get_compute_photons(),
        )
        out_states = {}
        x = 0
        (zero, one) = ([0, 1], [1, 0])
        for st in itertools.product([zero, one], repeat=len(readouts)):
            basic_out_state = [[]] * len(self.photons)
            for w in witnesses:
                basic_out_state[w.id] = zero
            for c in comps:
                basic_out_state[c.id] = zero
            for i in range(len(readouts)):
                basic_out_state[readouts[i].id] = st[i]
            out_states[
                str(pcvl.BasicState(list(itertools.chain.from_iterable(basic_out_state))))
            ] = f"|{x:0{len(readouts)}b}>"
            x = x + 1
        self.output_states = out_states

    def get_probability_distribution(self, format_result=True, postselection=True) -> PhotonDistribution:
        r"""Get the probability distribution of the measurement results.

        Parameters
        ----------
        format_result : bool, optional
            whether to format the result so that only the result corresponding to the output qubit is taken out.
        postselection : bool, optional
            whether to postselect the results.

        Returns
        -------
        result : PhotonDistribution
            Probability distribution of the measurement results.
        """
        if self.processor is None:
            raise Exception(
                "No processor has been set. Please set a processor by `set_local_processor` or `set_remote_procesor` before running the experiment."
            )
        if postselection:
            self.set_postselection()

        sampler = Sampler(self.processor)
        probs = PhotonDistribution(sampler.probs()["results"])

        if format_result:
            probs.replace_keys(self.output_states)

        return probs

    def sample(self, num_samples=1024, format_result=True, postselection=True) -> PhotonCount:
        """Run the MBQC pattern on IBMQ devices

        Parameters
        ----------
        num_samples : int, optional
            Number of samples.
        format_result : bool, optional
            whether to format the result so that only the result corresponding to the output qubit is taken out.
        postselection : bool, optional
            whether to postselect the results.

        Returns
        -------
        result : PhotonCount
            Measurement result.
        """
        if self.processor is None:
            raise Exception(
                "No processor has been set. Please set a processor by `set_local_processor` or `set_remote_procesor` before running the experiment."
            )
        if postselection:
            self.set_postselection()

        sampler = Sampler(self.processor)
        sample_result = PhotonCount(collections.Counter(sampler.samples(num_samples)["results"]))

        if format_result:
            sample_result.replace_keys(self.output_states)

        return sample_result

    def set_postselection(self):
        """Postselect the results according to the pattern."""
        ps = PostSelect()
        for ph in self.get_readout_photons():
            ps.eq([2 * ph.id, 2 * ph.id + 1], 1)
        for ph in self.get_compute_photons():
            ps.eq([2 * ph.id], 0).eq([2 * ph.id + 1], 1)
        for ph in self.get_witness_photons():
            ps.eq([2 * ph.id], 0).eq([2 * ph.id + 1], 1)

        self.processor.set_postselection(ps)

    def get_readout_photons(self):
        return [ph for ph in self.photons if ph.type == PhotonType.READOUT]

    def get_compute_photons(self):
        return [ph for ph in self.photons if ph.type == PhotonType.COMPUTE]

    def get_witness_photons(self):
        return [ph for ph in self.photons if ph.type == PhotonType.WITNESS]


class PhotonCount(dict):
    """PhotonCount class for storing the counts of the measurement results.

    perceval.BSCount does not seem to show fock state with one qubit properly."""

    def __init__(self, counts: dict[str, int] = {}):
        if not isinstance(counts, dict):
            raise TypeError("counts must be a dictionary.")
        super().__init__()
        self.counts = dict(counts)

    def __str__(self) -> str:
        return str(self.counts)

    def __getitem__(self, key: str) -> int:
        if not isinstance(key, str):
            raise TypeError("key must be a string.")
        return self.counts[key]

    def __setitem__(self, key: str, value: int):
        if not isinstance(key, str):
            raise TypeError("key must be a string.")
        if not (isinstance(value, int) and value >= 0):
            raise TypeError("value must be a positive integer.")
        self.counts[key] = value

    def items(self) -> dict_items:
        return self.counts.items()

    def draw(self, sort: bool = True):
        """Draw the counts result in a table.
        If the code is run in a Jupyter notebook, the table will be displayed in HTML format.
        If the code is run in a terminal, the table will be displayed in ASCII format.

        Parameters
        ----------
        sort : bool, optional
            Whether to sort the counts by the key.
        """
        headers = ["state", "counts"]
        d = []
        for key, value in self.counts.items():
            d.append([str(key), value])
        if sort:
            d.sort()
        if IS_NOTEBOOK:
            table = tabulate(d, headers=headers, tablefmt="html")
            display(HTML(table))
        else:
            table = tabulate(d, headers=headers, tablefmt="pretty")
            print(table)

    def replace_keys(self, replace_dict: dict[str, str]):
        """Replace the keys of the counts.

        Parameters
        ----------
        replace_dict : dict
            Dictionary of the replacement.
        """
        replaced = {}
        # Iterate over original measurement results
        for key, value in self.counts.items():
            if str(key) not in replace_dict:
                continue
            replaced[replace_dict[str(key)]] = value
        self.counts = replaced


class PhotonDistribution(dict):
    """PhotonDistribution class for storing the probability distribution of the measurement results.

    perceval.BSDistribution does not seem to show fock state with one qubit properly."""

    def __init__(self, distribution: dict[str, float] = {}):
        # TODO: use sympy.physics.secondquant.FockStateBosonBra?
        if not isinstance(distribution, dict):
            raise TypeError("distribution must be a dictionary.")
        super().__init__()
        self.distribution = dict(distribution)

    def __str__(self) -> str:
        return str(self.distribution)

    def __getitem__(self, key: str) -> float:
        if not isinstance(key, str):
            raise TypeError("key must be a string.")
        return self.distribution[key]

    def __setitem__(self, key: str, value: float):
        if not isinstance(key, str):
            raise TypeError("key must be a string.")
        if not isinstance(value, float):
            raise TypeError("value must be a float.")
        self.distribution[key] = value

    def items(self) -> dict_items:
        return self.distribution.items()

    def draw(self, sort: bool = True):
        """Draw the probability distribution in a table.
        If the code is run in a Jupyter notebook, the table will be displayed in HTML format.
        If the code is run in a terminal, the table will be displayed in ASCII format.

        Parameters
        ----------
        sort : bool, optional
            Whether to sort the distribution by the key.
        """
        headers = ["state", "probability"]
        d = []
        for key, value in self.distribution.items():
            d.append([str(key), value])
        if sort:
            d.sort()
        if IS_NOTEBOOK:
            table = tabulate(d, headers=headers, tablefmt="html")
            display(HTML(table))
        else:
            table = tabulate(d, headers=headers, tablefmt="pretty")
            print(table)

    def replace_keys(self, replace_dict: dict[str, str]):
        """Replace the keys of the distribution.

        Parameters
        ----------
        replace_dict : dict
            Dictionary of the replacement.
        """
        replaced = {}
        # Iterate over original measurement results
        for key, value in self.distribution.items():
            if str(key) not in replace_dict:
                continue
            replaced[replace_dict[str(key)]] = value
        self.distribution = replaced

graphix_perceval/experiment.py
TeamGraphix/graphix-perceval
from setuptools import find_packages, setup

with open("README.md", "r", encoding="utf-8") as fh:
    long_description = fh.read()

version = {}
with open("graphix_perceval/version.py") as fp:
    exec(fp.read(), version)

requirements = [requirement.strip() for requirement in open("requirements.txt").readlines()]

info = {
    "name": "graphix_perceval",
    "version": version["__version__"],
    "packages": find_packages(),
    "author": "Kazuki Tsuoka, Shinichi Sunami",
    "author_email": "shinichi.sunami@gmail.com",
    "maintainer": "Kazuki Tsuoka, Shinichi Sunami",
    "maintainer_email": "shinichi.sunami@gmail.com",
    "license": "Apache License 2.0",
    "description": "Perceval interface for graphix library, the MBQC compiler",
    "long_description": long_description,
    "long_description_content_type": "text/markdown",
    "url": "https://graphix.readthedocs.io",
    "project_urls": {"Bug Tracker": "https://github.com/TeamGraphix/graphix-perceval/issues"},
    "classifiers": [
        "Development Status :: 4 - Beta",
        "Environment :: Console",
        "Intended Audience :: Science/Research",
        "Programming Language :: Python :: 3.8",
        "Programming Language :: Python :: 3.9",
        "Programming Language :: Python :: 3.10",
        "License :: OSI Approved :: Apache Software License",
        "Operating System :: OS Independent",
        "Topic :: Scientific/Engineering :: Physics",
    ],
    "python_requires": ">=3.8,<3.11",
    "install_requires": requirements,
    "extras_require": {},
}

setup(**(info))

setup.py
TeamGraphix/graphix-perceval
import unittest

from graphix.clifford import CLIFFORD
from graphix_perceval.clifford import CLIFFORD_TO_PERCEVAL_BS, CLIFFORD_TO_PERCEVAL_POLAR
from sympy import matrix2numpy
import numpy as np
import perceval as pcvl


class TestConverter(unittest.TestCase):
    def test_clifford_conversion_BS(self):
        for idx in range(len(CLIFFORD_TO_PERCEVAL_BS)):
            circ = pcvl.Circuit(2)
            for component in CLIFFORD_TO_PERCEVAL_BS[idx]:
                circ.add(0, component)
            self.assertTrue(np.allclose(matrix2numpy(circ.U, dtype=complex), CLIFFORD[idx].astype(complex)))
        assert len(CLIFFORD_TO_PERCEVAL_BS) == len(CLIFFORD)

    def test_clifford_conversion_polar(self):
        for idx in range(len(CLIFFORD_TO_PERCEVAL_POLAR)):
            circ = pcvl.Circuit(1)
            for component in CLIFFORD_TO_PERCEVAL_POLAR[idx]:
                circ.add(0, component)
            self.assertTrue(np.allclose(matrix2numpy(circ.U, dtype=complex), CLIFFORD[idx].astype(complex)))
        assert len(CLIFFORD_TO_PERCEVAL_POLAR) == len(CLIFFORD)

tests/test_clifford.py
TeamGraphix/graphix-perceval
import unittest

import numpy as np
from graphix import Circuit

from graphix_perceval.converter import to_perceval
from graphix_perceval.experiment import PhotonDistribution


class TestConverter(unittest.TestCase):
    def test_sampling_circuit_wo_postselection(self):
        circuit = Circuit(2)
        circuit.h(1)
        circuit.cnot(0, 1)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.sample(num_samples=1000, format_result=False, postselection=False)
        counts = 0
        for k, v in dist.items():
            counts += v
        self.assertEqual(counts, 1000)

    def test_sampling_circuit_w_postselection(self):
        circuit = Circuit(2)
        circuit.h(1)
        circuit.cnot(0, 1)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.sample(num_samples=1000)
        counts = 0
        for k, v in dist.items():
            self.assertTrue(k in ["|00>", "|11>"])
            counts += v
        self.assertEqual(counts, 1000)

    def test_zero_state_creation_wo_pauli_meas(self):
        circuit = Circuit(1)  # initialize with |+>
        circuit.h(0)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.get_probability_distribution()

        ans_dist = PhotonDistribution({"|0>": 1.0})
        self.assertAlmostEqual(dist["|0>"], ans_dist["|0>"])
        self.assertTrue(all(k in dist.keys() for k in ans_dist.keys()))

    def test_one_state_creation_wo_pauli_meas(self):
        circuit = Circuit(1)  # initialize with |+>
        circuit.h(0)
        circuit.x(0)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.get_probability_distribution()

        ans_dist = PhotonDistribution({"|1>": 1.0})
        self.assertAlmostEqual(dist["|1>"], ans_dist["|1>"])
        self.assertTrue(all(k in dist.keys() for k in ans_dist.keys()))

    def test_rotated_one_qubit_state_creation_wo_pauli_meas(self):
        circuit = Circuit(1)  # initialize with |+>
        circuit.h(0)
        circuit.rx(0, np.pi / 1.23)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.get_probability_distribution()

        ans_dist = PhotonDistribution({"|0>": np.cos(np.pi / 1.23 / 2) ** 2, "|1>": np.sin(np.pi / 1.23 / 2) ** 2})
        self.assertAlmostEqual(dist["|0>"], ans_dist["|0>"])
        self.assertAlmostEqual(dist["|1>"], ans_dist["|1>"])
        self.assertTrue(all(k in dist.keys() for k in ans_dist.keys()))

    def test_bell_state_phi_plus_creation_wo_pauli_meas(self):
        circuit = Circuit(2)  # initialize with |+> \otimes |+>
        circuit.h(1)
        circuit.cnot(0, 1)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.get_probability_distribution()

        ans_dist = PhotonDistribution({"|00>": 0.5, "|11>": 0.5})
        self.assertAlmostEqual(dist["|00>"], ans_dist["|00>"])
        self.assertAlmostEqual(dist["|11>"], ans_dist["|11>"])
        self.assertTrue(all(k in dist.keys() for k in ans_dist.keys()))

    def test_bell_state_phi_plus_creation_with_pauli_meas(self):
        circuit = Circuit(2)  # initialize with |+> \otimes |+>
        circuit.h(1)
        circuit.cnot(0, 1)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()
        pattern.perform_pauli_measurements()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.get_probability_distribution()

        ans_dist = PhotonDistribution({"|00>": 0.5, "|11>": 0.5})
        self.assertAlmostEqual(dist["|00>"], ans_dist["|00>"])
        self.assertAlmostEqual(dist["|11>"], ans_dist["|11>"])
        self.assertTrue(all(k in dist.keys() for k in ans_dist.keys()))

    def test_ghz_state_creation_with_pauli_meas(self):
        circuit = Circuit(3)  # initialize with |+> \otimes |+> \otimes |+>
        circuit.h(1)
        circuit.h(2)
        circuit.cnot(0, 1)
        circuit.cnot(1, 2)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()
        pattern.perform_pauli_measurements()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.get_probability_distribution()

        ans_dist = PhotonDistribution({"|000>": 0.5, "|111>": 0.5})
        self.assertAlmostEqual(dist["|000>"], ans_dist["|000>"])
        self.assertAlmostEqual(dist["|111>"], ans_dist["|111>"])
        self.assertTrue(all(k in dist.keys() for k in ans_dist.keys()))

    def test_bell_state_and_ry_with_pauli_meas(self):
        circuit = Circuit(2)  # initialize with |+> \otimes |+> \otimes |+>
        circuit.h(1)
        circuit.cnot(0, 1)
        circuit.ry(1, np.pi / 4)
        pattern = circuit.transpile()
        pattern.standardize()
        pattern.shift_signals()
        pattern.perform_pauli_measurements()

        exp = to_perceval(pattern)
        exp.set_local_processor("SLOS")
        dist = exp.get_probability_distribution()

        ans_dist = PhotonDistribution(
            {
                "|00>": (np.cos(np.pi / 8) / np.sqrt(2)) ** 2,
                "|01>": (np.sin(np.pi / 8) / np.sqrt(2)) ** 2,
                "|10>": (np.sin(np.pi / 8) / np.sqrt(2)) ** 2,
                "|11>": (np.cos(np.pi / 8) / np.sqrt(2)) ** 2,
            }
        )
        self.assertAlmostEqual(dist["|00>"], ans_dist["|00>"])
        self.assertAlmostEqual(dist["|01>"], ans_dist["|01>"])
        self.assertAlmostEqual(dist["|10>"], ans_dist["|10>"])
        self.assertAlmostEqual(dist["|11>"], ans_dist["|11>"])
        self.assertTrue(all(k in dist.keys() for k in ans_dist.keys()))

tests/test_converter.py
TeamGraphix/graphix-perceval