final first commit

qibojit-benchmarks/benchmarks/libraries/__init__.py

@@ -0,0 +1,74 @@
+def parse(options):
+    """Parse options from string.
+
+    Args:
+        options (str): String with options.
+                       It should have the form 'arg1=value1,arg2=value2,...'.
+
+    Returns:
+        dict: {'arg1': value1, 'arg2': value2, ...}
+
+    """
+    kwargs = {}
+    if options is not None:
+        for parameter in options.split(","):
+            if "=" in parameter:
+                k, v = parameter.split("=")
+                kwargs[k] = v
+            else:
+                raise ValueError(f"Cannot parse parameter {parameter}.")
+    return kwargs
+
+
+def get(backend_name, options=None):
+    options = parse(options)
+    if backend_name == "qibo":
+        from benchmarks.libraries.qibo import Qibo
+        return Qibo(**options)
+
+    elif backend_name == "qiskit":
+        from benchmarks.libraries.qiskit import Qiskit
+        return Qiskit(**options)
+    elif backend_name == "qiskit-gpu":
+        from benchmarks.libraries.qiskit import QiskitGpu
+        return QiskitGpu(**options)
+
+    elif backend_name == "cirq":
+        from benchmarks.libraries.cirq import Cirq
+        return Cirq()
+    elif backend_name == "qsim":
+        from benchmarks.libraries.cirq import QSim
+        return QSim(**options)
+    elif backend_name == "qsim-gpu":
+        from benchmarks.libraries.cirq import QSimGpu
+        return QSimGpu(**options)
+    elif backend_name == "qsim-cuquantum":
+        from benchmarks.libraries.cirq import QSimCuQuantum
+        return QSimCuQuantum(**options)
+    elif backend_name == "tfq":
+        from benchmarks.libraries.cirq import TensorflowQuantum
+        return TensorflowQuantum()
+
+    elif backend_name == "qulacs":
+        from benchmarks.libraries.qulacs import Qulacs
+        return Qulacs()
+    elif backend_name == "qulacs-gpu":
+        from benchmarks.libraries.qulacs import QulacsGpu
+        return QulacsGpu()
+
+    elif backend_name == "qcgpu":
+        from benchmarks.libraries.qcgpu import QCGPU
+        return QCGPU()
+
+    elif backend_name == "projectq":
+        from benchmarks.libraries.projectq import ProjectQ
+        return ProjectQ(**options)
+
+    elif backend_name == "hybridq":
+        from benchmarks.libraries.hybridq import HybridQ
+        return HybridQ(**options)
+    elif backend_name == "hybridq-gpu":
+        from benchmarks.libraries.hybridq import HybridQGPU
+        return HybridQGPU(**options)
+
+    raise KeyError(f"Unknown simulation library {backend_name}.")

qibojit-benchmarks/benchmarks/libraries/abstract.py

@@ -0,0 +1,173 @@
+from abc import ABC, abstractmethod
+import numpy as np
+
+
+class AbstractBackend(ABC):
+
+    def __init__(self):
+        self.name = None
+        self.__version__ = None
+
+    @abstractmethod
+    def from_qasm(self, qasm):
+        raise NotImplementedError
+
+    @abstractmethod
+    def __call__(self, circuit):
+        raise NotImplementedError
+
+    def transpose_state(self, x):
+        """Switch order of qubits in state vector to be compatible to Qibo."""
+        shape = tuple(x.shape)
+        nqubits = int(np.log2(shape[0]))
+        x = np.reshape(x, nqubits * (2,))
+        x = np.transpose(x, range(nqubits - 1, -1, -1))
+        return np.reshape(x, shape)
+
+    @abstractmethod
+    def get_precision(self):
+        raise NotImplementedError
+
+    def set_precision(self, precision):
+        raise NotImplementedError(f"Cannot set precision for {self.name} backend.")
+
+    @abstractmethod
+    def get_device(self):
+        raise NotImplementedError
+
+
+class ParserBackend(AbstractBackend):
+
+    QASM_GATES = {"h": "H", "x": "X", "y": "Y", "z": "Z",
+                  "rx": "RX", "ry": "RY", "rz": "RZ",
+                  "u1": "U1", "u2": "U2", "u3": "U3",
+                  "cx": "CNOT", "swap": "SWAP", "cz": "CZ",
+                  "crx": "CRX", "cry": "CRY", "crz": "CRZ",
+                  "cu1": "CU1", "cu3": "CU3", "rzz": "RZZ",
+                  "ccx": "TOFFOLI", "id": "I"}
+    PARAMETRIZED_GATES = {"rx", "ry", "rz", "u1", "u2", "u3",
+                          "crx", "cry", "crz", "cu1", "cu3", "rzz"}
+
+    def parse(self, qasm_code):
+        """Extracts circuit information from QASM script.
+
+        Args:
+            qasm_code: String with the QASM code to parse.
+
+        Returns:
+            nqubits: The total number of qubits in the circuit.
+            gate_list: List that specifies the gates of the circuit.
+                Contains tuples of the form
+                (Qibo gate name, qubit IDs, optional additional parameter).
+                The additional parameter is the ``register_name`` for
+                measurement gates or ``theta`` for parametrized gates.
+        """
+        import re
+        def read_args(args):
+            _args = iter(re.split(r"[\[\],]", args))
+            for name in _args:
+                if name:
+                    index = next(_args)
+                    if not index.isdigit():
+                        raise ValueError("Invalid QASM qubit arguments: {}".format(args))
+                    yield name, int(index)
+
+        # Remove comment lines
+        lines = "".join(line for line in qasm_code.split("\n")
+                        if line and line[:2] != "//")
+        lines = (line for line in lines.split(";") if line)
+
+        if next(lines) != "OPENQASM 2.0":
+            raise ValueError("QASM code should start with 'OPENQASM 2.0'.")
+
+        qubits = {} # Dict[Tuple[str, int], int]: map from qubit tuple to qubit id
+        cregs_size = {} # Dict[str, int]: map from `creg` name to its size
+        registers = {} # Dict[str, List[int]]: map from register names to target qubit ids
+        gate_list = [] # List[Tuple[str, List[int]]]: List of (gate name, list of target qubit ids)
+        for line in lines:
+            command, args = line.split(None, 1)
+            # remove spaces
+            command = command.replace(" ", "")
+            args = args.replace(" ", "")
+
+            if command == "include":
+                pass
+
+            elif command == "qreg":
+                for name, nqubits in read_args(args):
+                    for i in range(nqubits):
+                        qubits[(name, i)] = len(qubits)
+
+            elif command == "creg":
+                for name, nqubits in read_args(args):
+                    cregs_size[name] = nqubits
+
+            elif command == "measure":
+                args = args.split("->")
+                if len(args) != 2:
+                    raise ValueError("Invalid QASM measurement: {}".format(line))
+                qubit = next(read_args(args[0]))
+                if qubit not in qubits:
+                    raise ValueError("Qubit {} is not defined in QASM code."
+                                     "".format(qubit))
+
+                register, idx = next(read_args(args[1]))
+                if register not in cregs_size:
+                    raise ValueError("Classical register name {} is not defined "
+                                     "in QASM code.".format(register))
+                if idx >= cregs_size[register]:
+                    raise ValueError("Cannot access index {} of register {} "
+                                     "with {} qubits."
+                                     "".format(idx, register, cregs_size[register]))
+                if register in registers:
+                    if idx in registers[register]:
+                        raise KeyError("Key {} of register {} has already "
+                                       "been used.".format(idx, register))
+                    registers[register][idx] = qubits[qubit]
+                else:
+                    registers[register] = {idx: qubits[qubit]}
+                    gate_list.append(("M", register))
+
+            else:
+                pieces = [x for x in re.split("[()]", command) if x]
+                if len(pieces) == 1:
+                    gatename, params = pieces[0], None
+                    if gatename not in self.QASM_GATES:
+                        raise ValueError("QASM command {} is not recognized."
+                                         "".format(command))
+                    if gatename in self.PARAMETRIZED_GATES:
+                        raise ValueError("Missing parameters for QASM "
+                                         "gate {}.".format(gatename))
+
+                elif len(pieces) == 2:
+                    gatename, params = pieces
+                    if gatename not in self.PARAMETRIZED_GATES:
+                        raise ValueError("Invalid QASM command {}."
+                                         "".format(command))
+                    params = params.replace(" ", "").split(",")
+                    try:
+                        for i, p in enumerate(params):
+                            if 'pi' in p:
+                                import math
+                                from operator import mul
+                                from functools import reduce
+                                s = p.replace('pi', str(math.pi)).split('*')
+                                p = reduce(mul, [float(j) for j in s], 1)
+                            params[i] = float(p)
+                    except ValueError:
+                        raise ValueError("Invalid value {} for gate parameters."
+                                         "".format(params))
+
+                else:
+                    raise ValueError("QASM command {} is not recognized."
+                                     "".format(command))
+
+                # Add gate to gate list
+                qubit_list = []
+                for qubit in read_args(args):
+                    if qubit not in qubits:
+                        raise ValueError("Qubit {} is not defined in QASM "
+                                         "code.".format(qubit))
+                    qubit_list.append(qubits[qubit])
+                gate_list.append((self.QASM_GATES[gatename], list(qubit_list), params))
+        return len(qubits), gate_list

qibojit-benchmarks/benchmarks/libraries/cirq.py

@@ -0,0 +1,167 @@
+import numpy as np
+from benchmarks.libraries import abstract
+
+
+class Cirq(abstract.ParserBackend):
+
+    def __init__(self):
+        import cirq
+        self.name = "cirq"
+        self.__version__ = cirq.__version__
+        self.cirq = cirq
+        self.precision = "double"
+        self.simulator = cirq.Simulator(dtype=np.complex128)
+
+    def RX(self, theta):
+        return self.cirq.rx(theta)
+
+    def RY(self, theta):
+        return self.cirq.ry(theta)
+
+    def RZ(self, theta):
+        return self.cirq.rz(theta)
+
+    def CU1(self, theta):
+        return self.cirq.CZPowGate(exponent=theta / np.pi)
+
+    def CU3(self, theta, phi, lam):
+        gate = self.cirq.circuits.qasm_output.QasmUGate(theta / np.pi, phi / np.pi, lam / np.pi)
+        return gate.controlled(num_controls=1)
+
+    def RZZ(self, theta):
+        import numpy as np
+        return self.cirq.ZZPowGate(exponent=theta / np.pi, global_shift=-0.5)
+
+    def __getattr__(self, x):
+        return getattr(self.cirq, x)
+
+    def __getitem__(self, x):
+        return getattr(self.cirq, x)
+
+    def from_qasm(self, qasm):
+        from cirq.contrib.qasm_import import circuit_from_qasm, exception
+        try:
+            return circuit_from_qasm(qasm)
+        except exception.QasmException:
+            nqubits, gatelist = self.parse(qasm)
+            qubits = [self.cirq.GridQubit(i, 0) for i in range(nqubits)]
+            circuit = self.cirq.Circuit()
+            for gatename, qid, params in gatelist:
+                if params is not None:
+                    gate = getattr(self, gatename)(*params)
+                else:
+                    gate = getattr(self, gatename)
+                circuit.append(gate(*(qubits[i] for i in qid)))
+            return circuit
+
+    def __call__(self, circuit):
+        result = self.simulator.simulate(circuit)
+        return result.final_state_vector
+
+    def transpose_state(self, x):
+        return x
+
+    def get_precision(self):
+        return self.precision
+
+    def set_precision(self, precision):
+        import numpy as np
+        self.precision = precision
+        if precision == "single":
+            self.simulator = self.cirq.Simulator(dtype=np.complex64)
+        else:
+            self.simulator = self.cirq.Simulator(dtype=np.complex128)
+
+    def get_device(self):
+        return None
+
+
+class TensorflowQuantum(Cirq):
+
+    def __init__(self):
+        import cirq
+        import tensorflow_quantum as tfq
+        self.name = "tfq"
+        self.cirq = cirq
+        self.precision = "single"
+        self.__version__ = tfq.__version__
+        self.state_layer = tfq.layers.State()
+
+    def set_precision(self, precision):
+        if precision == "double":
+            raise NotImplementedError(f"Cannot set precision '{precision}' for {self.name} backend.")
+
+    def from_qasm(self, qasm):
+        circuit = super().from_qasm(qasm)
+        # change `NamedQubit`s to `GridQubit`s as TFQ understands only `GridQubit`
+        qubit_map = {}
+        for q in circuit.all_qubits():
+            if isinstance(q, self.cirq.NamedQubit):
+                i = int(str(q).split("_")[-1])
+                qubit_map[q] = self.cirq.GridQubit(i, 0)
+        if qubit_map:
+            return circuit.transform_qubits(qubit_map)
+        return circuit
+
+    def __call__(self, circuit):
+        # transfer final state to numpy array because that's what happens
+        # for all backends
+        return self.state_layer(circuit)[0].numpy()
+
+
+class QSim(Cirq):
+
+    def __init__(self, max_qubits="0", nthreads=None):
+        import cirq
+        import qsimcirq
+        self.name = "qsim"
+        self.cirq = cirq
+        self.qsimcirq = qsimcirq
+        self.precision = "single"
+        self.__version__ = qsimcirq.__version__
+
+        if nthreads is None:
+            from multiprocessing import cpu_count
+            self.nthreads = cpu_count()
+        else:
+            self.nthreads = int(nthreads)
+        self.max_qubits = int(max_qubits)
+
+        self.simulator = self.get_simulator()
+
+    def get_simulator(self):
+        return self.qsimcirq.QSimSimulator({'t': self.nthreads, 'f': self.max_qubits})
+
+    def set_precision(self, precision):
+        if precision == "double":
+            raise NotImplementedError(f"Cannot set precision '{precision}' for {self.name} backend.")
+
+
+class QSimGpu(QSim):
+
+    def __init__(self, max_qubits="0"):
+        super().__init__(max_qubits)
+        self.name = "qsim-gpu"
+
+    def get_simulator(self):
+        qsim_options = self.qsimcirq.QSimOptions(
+                use_gpu=True,
+                gpu_mode=0,
+                max_fused_gate_size=self.max_qubits
+            )
+        return self.qsimcirq.QSimSimulator(qsim_options)
+
+
+class QSimCuQuantum(QSim):
+
+    def __init__(self, max_qubits="0"):
+        super().__init__(max_qubits)
+        self.name = "qsim-cuquantum"
+
+    def get_simulator(self):
+        qsim_options = self.qsimcirq.QSimOptions(
+                use_gpu=True,
+                gpu_mode=1,
+                max_fused_gate_size=self.max_qubits
+            )
+        return self.qsimcirq.QSimSimulator(qsim_options)

qibojit-benchmarks/benchmarks/libraries/hybridq.py

@@ -0,0 +1,142 @@
+import os
+import numpy as np
+from benchmarks.libraries import abstract
+
+
+class HybridQ(abstract.ParserBackend):
+
+    def __init__(self, max_qubits="0", simplify="False"):
+        from hybridq.gate import Gate, MatrixGate
+        self.name = "hybridq"
+        self.__version__ = "0.7.7.post2"
+        self.Gate = Gate
+        self.MatrixGate = MatrixGate
+        self.max_qubits = int(max_qubits)
+        if simplify in ("true", "True"):
+            self.simplify = True
+        else:
+            self.simplify = False
+        self.complex_type = "complex128"
+        self.max_qubits = int(max_qubits)
+
+    def H(self, q):
+        return self.Gate('H', qubits=(q,))
+
+    def X(self, q):
+        return self.Gate('X', qubits=(q,))
+
+    def Y(self, q):
+        return self.Gate('Y', qubits=(q,))
+
+    def Z(self, q):
+        return self.Gate('Z', qubits=(q,))
+
+    def RX(self, q, theta):
+        return self.Gate('RX', params=[theta], qubits=(q,))
+
+    def RY(self, q, theta):
+        return self.Gate('RY', params=[theta], qubits=(q,))
+
+    def RZ(self, q, theta):
+        return self.Gate('RZ', params=[theta], qubits=(q,))
+
+    def U1(self, q, theta):
+        phase = np.exp(1j * theta)
+        matrix = np.diag([1, phase])
+        return self.MatrixGate(U=matrix, qubits=(q,))
+
+    def U2(self, q, phi, lam):
+        plus = np.exp(0.5j * (phi + lam))
+        minus = np.exp(0.5j * (phi - lam))
+        matrix = np.array([[np.conj(plus), -np.conj(minus)], [minus, plus]]) / np.sqrt(2)
+        return self.MatrixGate(U=matrix, qubits=(q,))
+
+    def U3(self, q, theta, phi, lam):
+        return self.Gate('U3', params=[theta, phi, lam], qubits=(q,))
+
+    def CNOT(self, q1, q2):
+        return self.Gate('CNOT', qubits=(q1, q2))
+
+    def SWAP(self, q1, q2):
+        return self.Gate('SWAP', qubits=(q1, q2))
+
+    def CZ(self, q1, q2):
+        return self.Gate('CZ', qubits=(q1, q2))
+
+    def CU1(self, q1, q2, theta):
+        return self.Gate('CPHASE', params=[theta], qubits=(q1, q2))
+
+    def CU3(self, q1, q2, theta, phi, lam):
+        from hybridq.gate import Control
+        cost, sint = np.cos(theta / 2.0), np.sin(theta / 2.0)
+        pplus, pminus = np.exp(0.5j * (phi + lam)), np.exp(0.5j * (phi - lam))
+        matrix = np.array([[np.conj(pplus) * cost, -np.conj(pminus) * sint],
+                           [pminus * sint, pplus * cost]])
+        gate = self.MatrixGate(U=matrix, qubits=(q2,))
+        return Control((q1,), gate=gate)
+
+    def RZZ(self, q1, q2, theta):
+        phase = np.exp(0.5j * theta)
+        phasec = np.conj(phase)
+        matrix = np.diag([phasec, phase, phase, phasec])
+        return self.MatrixGate(U=matrix, qubits=(q1, q2))
+
+    def from_qasm(self, qasm):
+        from hybridq.circuit import Circuit
+        nqubits, gatelist = self.parse(qasm)
+        circuit = Circuit()
+        for gatename, qubits, params in gatelist:
+            args = list(qubits)
+            if params:
+                args.extend(params)
+            gate = getattr(self, gatename)(*args)
+            circuit.append(gate)
+        return circuit
+
+    def __call__(self, circuit):
+        from hybridq.circuit.simulation import simulate
+        initial_state = len(circuit.all_qubits()) * '0'
+        final_state = simulate(circuit, optimize="evolution",
+                               initial_state=initial_state,
+                               complex_type=self.complex_type,
+                               simplify=self.simplify,
+                               compress=self.max_qubits,
+                               max_largest_intermediate=2**40)
+        return final_state.ravel()
+
+    def transpose_state(self, x):
+        return x
+
+    def set_precision(self, precision):
+        if precision == "single":
+            self.complex_type = "complex64"
+        else:
+            self.complex_type = "complex128"
+
+    def get_precision(self):
+        if self.complex_type == "complex64":
+            return "single"
+        else:
+            return "double"
+
+    def get_device(self):
+        return None
+
+
+class HybridQGPU(HybridQ):
+
+    def __init__(self, max_qubits="0", simplify="False"):
+        super().__init__(max_qubits=max_qubits, simplify=simplify)
+        self.name = "hybridq-gpu"
+
+    def __call__(self, circuit):
+        from hybridq.circuit.simulation import simulate
+        initial_state = len(circuit.all_qubits()) * '0'
+        final_state = simulate(circuit, optimize="evolution-einsum",
+                               backend="jax",
+                               initial_state=initial_state,
+                               complex_type=self.complex_type,
+                               simplify=self.simplify,
+                               compress=self.max_qubits,
+                               max_largest_intermediate=2**40)
+        return final_state.ravel()

qibojit-benchmarks/benchmarks/libraries/projectq.py

@@ -0,0 +1,129 @@
+import numpy as np
+from benchmarks.libraries import abstract
+
+class ProjectQ(abstract.ParserBackend):
+
+    def __init__(self, max_qubits="0", local_optimizer="0"):
+        """Initialize data members.
+
+        Args:
+            max_qubits (str): if "0", gate fusion is disabled, otherwise it's enabled.
+                              Note that it's not possible to set the maximum fused gate size.
+            local_optimizer (str): if "0", local optimization of circuits is disabled,
+                              otherwise it's enabled.
+        """
+        import projectq
+        self.name = "projectq"
+        self.projectq = projectq
+        self.__version__ = None
+        self.gate_fusion = int(max_qubits) > 0
+        self.local_optimizer = bool(int(local_optimizer))
+
+    def RX(self, theta):
+        return self.projectq.ops.Rx(theta)
+
+    def RY(self, theta):
+        return self.projectq.ops.Ry(theta)
+
+    def RZ(self, theta):
+        return self.projectq.ops.Rz(theta)
+
+    def U1(self, theta):
+        return self.projectq.ops.R(theta)
+
+    def U2(self, phi, lam):
+        pplus, pminus = np.exp(0.5j * (phi + lam)), np.exp(0.5j * (phi - lam))
+        matrix = np.array([[np.conj(pplus), -np.conj(pminus)],
+                           [pminus, pplus]])
+        matrix /= np.sqrt(2)
+        return self.projectq.ops.MatrixGate(matrix)
+
+    def U3(self, theta, phi, lam):
+        cost, sint = np.cos(theta / 2.0), np.sin(theta / 2.0)
+        pplus, pminus = np.exp(0.5j * (phi + lam)), np.exp(0.5j * (phi - lam))
+        matrix = np.array([[np.conj(pplus) * cost, -np.conj(pminus) * sint],
+                           [pminus * sint, pplus * cost]])
+        return self.projectq.ops.MatrixGate(matrix)
+
+
+    def SWAP(self):
+        return self.projectq.ops.Swap
+    
+    def CRX(self, theta):
+        return self.projectq.ops.C(self.RX(theta))
+
+    def CRY(self, theta):
+        return self.projectq.ops.C(self.RY(theta))
+
+    def CRZ(self, theta):
+        return self.projectq.ops.CRz(theta)
+    
+    def CU1(self, theta):
+        U1 = self.projectq.ops.R(theta)
+        return self.projectq.ops.C(U1, n_qubits=1)
+
+    def CU3(self, theta):
+        raise NotImplementedError
+
+    def RZZ(self, theta):
+        return self.projectq.ops.Rzz(theta)
+
+    def __getattr__(self, x):
+        return getattr(self.projectq.ops, x)
+
+    def __item__(self, x):
+        return getattr(self.projectq.ops, x)
+
+    def from_qasm(self, qasm):
+        nqubits, gatelist = self.parse(qasm)
+        backend = self.projectq.backends.Simulator(gate_fusion=self.gate_fusion)
+        if self.local_optimizer:
+            self.eng = self.projectq.MainEngine(
+                backend=backend, engine_list=[self.projectq.cengines.LocalOptimizer()]
+            )
+        else:
+            self.eng = self.projectq.MainEngine(backend=backend)
+        qureg = self.eng.allocate_qureg(nqubits)
+        for gatename, qubits, params in gatelist:
+            gate = getattr(self, gatename)
+            if params is not None:
+                parameters = list(params)
+                if len(qubits) > 1:
+                    gate(*parameters) | tuple(qureg[i] for i in qubits)
+                else:
+                    gate(*parameters) | qureg[qubits[0]]
+            elif len(qubits) > 1:
+                if gatename == "SWAP":
+                    gate() | tuple(qureg[i] for i in qubits)
+                else:
+                    gate | tuple(qureg[i] for i in qubits)
+            else:
+                gate | qureg[qubits[0]]
+
+        return qureg 
+
+    def __call__(self, qureg):
+        self.eng.flush()
+        self.qubit_id, wave = self.eng.backend.cheat()
+        # measure everything to avoid error when running
+        self.projectq.ops.All(self.projectq.ops.Measure) | qureg
+        return np.array(wave)
+
+    def transpose_state(self, x):
+        shape = tuple(x.shape)
+        nqubits = int(np.log2(shape[0]))
+        x = np.reshape(x, nqubits * (2,))
+        x = np.transpose(x, range(nqubits - 1, -1, -1))
+        x = np.transpose(x, tuple(self.qubit_id[key] for key in self.qubit_id))
+        x = np.reshape(x, shape)
+        return x
+
+    def set_precision(self, precision):
+        if precision != "double":
+            raise NotImplementedError(f"Cannot set {precision} precision for {self.name} backend.")
+
+    def get_precision(self):
+        return "double"
+
+    def get_device(self):
+        return None

qibojit-benchmarks/benchmarks/libraries/qcgpu.py

@@ -0,0 +1,85 @@
+import numpy as np
+from benchmarks.libraries import abstract
+
+
+class QCGPU(abstract.ParserBackend):
+
+    def __init__(self):
+        import os
+        os.environ["PYOPENCL_CTX"] = "0"
+        import qcgpu
+        self.name = "qcgpu"
+        self.qcgpu = qcgpu
+        self.__version__ = None
+
+    def RX(self, target, theta):
+        cost, sint = np.cos(theta / 2.0), np.sin(theta / 2.0)
+        matrix = np.array([[cost, -1j * sint], [-1j * sint, cost]])
+        gate = self.qcgpu.Gate(matrix)
+        return ("apply_gate", (gate, target))
+
+    def RY(self, target, theta):
+        cost, sint = np.cos(theta / 2.0), np.sin(theta / 2.0)
+        matrix = np.array([[cost, -sint], [sint, cost]])
+        gate = self.qcgpu.Gate(matrix)
+        return ("apply_gate", (gate, target))
+
+    def RZ(self, target, theta):
+        phase = np.exp(0.5j * theta)
+        matrix = np.diag([np.conj(phase), phase])
+        gate = self.qcgpu.Gate(matrix)
+        return ("apply_gate", (gate, target))
+
+    def U1(self, target, theta):
+        phase = np.exp(1j * theta)
+        matrix = np.diag([1, phase])
+        gate = self.qcgpu.Gate(matrix)
+        return ("apply_gate", (gate, target))
+
+    def CU1(self, control, target, theta):
+        phase = np.exp(1j * theta)
+        matrix = np.diag([1, phase])
+        gate = self.qcgpu.Gate(matrix)
+        return ("apply_controlled_gate", (gate, control, target))
+
+    def RZZ(self, target1, target2, theta):
+        raise NotImplementedError
+
+    class QCGPUCircuit(list):
+
+        def __init__(self, nqubits):
+            self.nqubits = nqubits
+
+    def from_qasm(self, qasm):
+        nqubits, gatelist = self.parse(qasm)
+        circuit = self.QCGPUCircuit(nqubits)
+        for gate, qubits, params in gatelist:
+            args = list(qubits)
+            if params is not None:
+                args.extend(params)
+            if gate == "SWAP":
+                target1, target2 = qubits
+                circuit.append(("cx", (target1, target2)))
+                circuit.append(("cx", (target2, target1)))
+                circuit.append(("cx", (target1, target2)))
+            elif gate in {"RX", "RY", "RZ", "U1", "CU1"}:
+                circuit.append(getattr(self, gate)(*args))
+            else:
+                circuit.append((gate.lower(), args))
+        return circuit
+
+    def __call__(self, circuit):
+        state = self.qcgpu.State(circuit.nqubits)
+        for gate, args in circuit:
+            getattr(state, gate)(*args)
+        return state.amplitudes()
+
+    def set_precision(self, precision):
+        if precision != "single":
+            raise NotImplementedError(f"Cannot set {precision} precision for {self.name} backend.")
+
+    def get_precision(self):
+        return "single"
+
+    def get_device(self):
+        return None

qibojit-benchmarks/benchmarks/libraries/qibo.py

@@ -0,0 +1,295 @@
+from benchmarks.libraries import abstract
+from benchmarks.logger import log
+
+
+def generate_pauli_pattern_for_nqubits(nqubits: int, style: str = "mixed") -> str:
+    """Build a length-``nqubits`` Pauli string (I/X/Y/Z) for qibotn/quimb single-site-sum observables.
+
+    Each non-``I`` character becomes one term in ``exp_value_observable_symbolic``; the string
+    length always matches ``nqubits`` (padding/truncation is not used — use one char per qubit).
+
+    Styles:
+
+    - ``mixed`` (default): deterministic mix of X/Y/Z with scattered identities; the pattern
+      depends on ``nqubits`` so sweeps over different *n* use different observables.
+    - ``dense``: repeating XYZ on every qubit (no identities).
+    - ``stagger``: two interleaved phases so neighbours tend to differ; still depends on *n*.
+    """
+    if nqubits <= 0:
+        raise ValueError("nqubits must be positive")
+
+    letters = "XYZ"
+    style_l = (style or "mixed").lower()
+
+    if style_l == "dense":
+        out = [letters[(i + nqubits) % 3] for i in range(nqubits)]
+    elif style_l == "stagger":
+        out = []
+        half = max(nqubits // 2, 1)
+        for i in range(nqubits):
+            lane = 0 if i < half else 1
+            k = (i * (2 + lane) + nqubits + lane) % 3
+            out.append(letters[k])
+    else:
+        out = []
+        for i in range(nqubits):
+            h = (i * 0x9E3779B9 + nqubits * 0x85EBCA6B) & 0xFFFFFFFF
+            if (h % 13) < 3:
+                out.append("I")
+            else:
+                rot = (i ^ (nqubits >> 1)) + ((h >> 8) % 3)
+                out.append(letters[rot % 3])
+        if all(c == "I" for c in out):
+            out[-1] = "X"
+
+    return "".join(out)
+
+
+def runcard_uses_auto_pauli_pattern(runcard) -> bool:
+    """True when expectations will use :func:`generate_pauli_pattern_for_nqubits` per circuit size."""
+    if not runcard:
+        return False
+    raw = runcard.get("pauli_pattern")
+    auto = runcard.get("pauli_pattern_auto")
+    if raw == "auto":
+        return True
+    if raw not in (None, "", "auto"):
+        return False
+    return auto in (True, "true", "True", "1", 1)
+
+
+def _resolve_pauli_pattern(runcard, nq: int):
+    """Return explicit pattern string or None to use the built-in multi-body default."""
+    if not runcard:
+        return None
+    raw = runcard.get("pauli_pattern")
+    auto = runcard.get("pauli_pattern_auto")
+    style = runcard.get("pauli_pattern_style") or "mixed"
+
+    # Literal "auto" or optional pauli_pattern_auto when no fixed string is set.
+    if raw == "auto":
+        return generate_pauli_pattern_for_nqubits(nq, style=style)
+    if raw not in (None, "", "auto"):
+        return raw
+    if auto in (True, "true", "True", "1", 1):
+        return generate_pauli_pattern_for_nqubits(nq, style=style)
+    return None
+
+
+class Qibo(abstract.AbstractBackend):
+
+    def __init__(
+        self,
+        max_qubits="0",
+        backend="qibojit",
+        platform=None,
+        accelerators="",
+        expectation=None,
+        computation_settings=None,
+    ):
+        import qibo
+
+        runcard = None
+
+        if computation_settings is not None:
+            import json
+
+            try:
+                with open(computation_settings, "r") as f:
+                    runcard = json.load(f)
+
+            except json.JSONDecodeError as e:
+                raise ValueError(f"Invalid JSON in file '{computation_settings}': {e}")
+
+            except FileNotFoundError:
+                raise FileNotFoundError(f"File not found: {computation_settings}")
+
+            if runcard["expectation_enabled"] == True:
+                expectation = True
+
+            qibo.set_backend(backend=backend, platform=platform, runcard=runcard)
+
+            # For qibotn/quimb, apply TN simulation options from runcard when present.
+            if backend == "qibotn" and platform == "quimb":
+                quimb_backend = qibo.get_backend()
+                use_mps = runcard.get("use_mps", runcard.get("MPS_enabled", True))
+                max_bond_dimension = runcard.get(
+                    "max_bond_dimension", runcard.get("max_bond", None)
+                )
+                svd_cutoff = runcard.get("svd_cutoff", 1e-10)
+                mpi_enabled = runcard.get("MPI_enabled", False)
+                quimb_backend.configure_tn_simulation(
+                    ansatz="mps" if use_mps else None,
+                    max_bond_dimension=max_bond_dimension if use_mps else None,
+                    svd_cutoff=svd_cutoff,
+                    MPI_enabled=mpi_enabled,
+                )
+        else:
+            qibo.set_backend(backend=backend, platform=platform)
+
+        from qibo import models
+
+        self.name = "qibo"
+        self.qibo = qibo
+        self.models = models
+        self.__version__ = qibo.__version__
+        self.max_qubits = int(max_qubits)
+        self.accelerators = self._parse_accelerators(accelerators)
+        self.expectation_flag = expectation
+        self.backend_name_str = backend
+        self.platform_str = platform
+        self.runcard = runcard
+
+    def from_qasm(self, qasm):
+        circuit = self.models.Circuit.from_qasm(qasm, accelerators=self.accelerators)
+        if self.max_qubits > 1:
+            if self.max_qubits > 2:
+                log.warn(
+                    "Fusion with {} qubits is not yet supported by Qibo. "
+                    "Using max_qubits=2.".format(self.max_qubits)
+                )
+            circuit = circuit.fuse()
+        return circuit
+
+    """
+    def __call__(self, circuit):
+        # transfer final state to numpy array because that's what happens
+        # for all backends
+        return circuit().state(numpy=True)
+    """
+
+    def __call__(self, circuit):
+        # transfer final state to numpy array because that's what happens
+        # for all backends
+        if self.backend_name_str == "qibojit" and self.expectation_flag is not None:
+            from qibo.symbols import X, Y, Z, I
+            from qibo.hamiltonians import SymbolicHamiltonian
+            import numpy as np
+
+            # from qibo.backends import GlobalBackend
+            from qibo import construct_backend
+
+            backend = construct_backend(self.backend_name_str)
+            # self.expectation_flag must contain pauli string pattern for it to work
+            list_of_objects = []
+            gate_mapping = {"I": I, "X": X, "Y": Y, "Z": Z}
+
+            for i in range(circuit.nqubits):
+                gate = gate_mapping[
+                    self.expectation_flag[i % len(self.expectation_flag)]
+                ]
+                list_of_objects.append(gate(i))
+            obs = np.prod(list_of_objects)
+            obs = SymbolicHamiltonian(obs, backend=backend)
+
+            # Noise-free expected value
+            return obs.expectation(circuit)
+        else:
+            if self.expectation_flag:
+                if self.backend_name_str == "qibotn" and self.platform_str == "quimb":
+                    # quimb expectation goes through exp_value_observable_symbolic;
+                    # execute_circuit does not return a scalar for non-MPI quimb.
+                    import numpy as np
+                    nq = circuit.nqubits
+
+                    # If pauli_pattern is set in the JSON config (e.g. "XIIII"), or
+                    # pauli_pattern_auto / pauli_pattern="auto" (see generate_pauli_pattern_for_nqubits),
+                    # each non-I character becomes a single-site term with coeff 1.0.
+                    # "X" on site i means X_i ⊗ I elsewhere.
+                    pauli_pattern = _resolve_pauli_pattern(self.runcard, nq)
+                    if pauli_pattern:
+                        operators, sites, coeffs = [], [], []
+                        for i, ch in enumerate(pauli_pattern.upper()):
+                            if ch != "I" and i < nq:
+                                operators.append(ch.lower())
+                                sites.append((i,))
+                                coeffs.append(1.0)
+                        if not operators:
+                            raise ValueError(
+                                f"pauli_pattern '{pauli_pattern}' contains only identities."
+                            )
+                    else:
+                        # Default observable mirrors test_mpi_quimb.py:
+                        #   z@0, x@1, zz@(2,3), yy@(3,4), xyz@(0,1,2)
+                        operators = ["z", "x"]
+                        sites     = [(0,), (min(1, nq - 1),)]
+                        coeffs    = [1.0, 0.5]
+                        if nq >= 4:
+                            operators += ["zz", "yy"]
+                            sites     += [(min(2, nq - 2), min(3, nq - 1)),
+                                          (min(3, nq - 2), min(4, nq - 1))]
+                            coeffs    += [0.8, 0.3]
+                        if nq >= 3:
+                            operators += ["xyz"]
+                            sites     += [(0, min(1, nq - 2), min(2, nq - 1))]
+                            coeffs    += [0.2]
+
+                    return np.real(
+                        self.qibo.get_backend().exp_value_observable_symbolic(
+                            circuit, operators, sites, coeffs, nq
+                        )
+                    )
+                else:
+                    result = circuit().real
+                    return result.get() if hasattr(result, "get") else result
+            else:
+                if self.backend_name_str == "qibotn":
+                    if self.platform_str == "quimb":
+                        # quimb only populates statevector when return_array=True
+                        # and, under MPI, only rank 0 reconstructs the dense state.
+                        # Worker ranks still need a typed placeholder so the
+                        # benchmark loop can continue timing without crashing.
+                        import numpy as np
+
+                        result = self.qibo.get_backend().execute_circuit(
+                            circuit, return_array=True
+                        )
+                        if result.statevector is None:
+                            return np.empty(0, dtype=self.qibo.get_dtype())
+                        return result.statevector.flatten()
+                    else:
+                        return circuit().statevector.flatten()
+                else:
+                    return circuit().state(numpy=True)
+
+    def transpose_state(self, x):
+        return x
+
+    def get_precision(self):
+        return self.qibo.get_dtype()
+
+    def set_precision(self, precision):
+        self.qibo.set_dtype(precision)
+
+    def get_device(self):
+        return self.qibo.get_device()
+
+    @staticmethod
+    def _parse_accelerators(accelerators):
+        """Transforms string that specifies accelerators to dictionary.
+
+        The string that is parsed has the following format:
+            n1device1+n2device2+n3device3,...
+        and is transformed to the dictionary:
+            {'device1': n1, 'device2': n2, 'device3': n3, ...}
+
+        Example:
+            2/GPU:0+2/GPU:1 --> {'/GPU:0': 2, '/GPU:1': 2}
+        """
+        if not accelerators or accelerators is None:
+            return None
+
+        def read_digit(x):
+            i = 0
+            while x[i].isdigit():
+                i += 1
+            return x[i:], int(x[:i])
+
+        accelerator_dict = {}
+        for entry in accelerators.split("+"):
+            device, n = read_digit(entry)
+            if device in accelerator_dict:
+                accelerator_dict[device] += n
+            else:
+                accelerator_dict[device] = n
+        return accelerator_dict

qibojit-benchmarks/benchmarks/libraries/qiskit.py

@@ -0,0 +1,64 @@
+from benchmarks.libraries import abstract
+
+
+class Qiskit(abstract.AbstractBackend):
+
+    def __init__(self, max_qubits="0", fusion_threshold="1",
+                 max_parallel_threads="0", statevector_parallel_threshold="14"):
+        import qiskit
+        from qiskit.providers.aer import StatevectorSimulator
+        self.name = "qiskit"
+        self.__version__ = qiskit.__version__
+        self.max_qubits = int(max_qubits)
+        self.sim_options = dict(
+                max_parallel_threads=int(max_parallel_threads),
+                statevector_parallel_threshold=int(statevector_parallel_threshold),
+                fusion_enable=self.max_qubits > 0,
+                fusion_max_qubit=self.max_qubits,
+                fusion_threshold=int(fusion_threshold),
+                precision="double"
+            )
+        self.simulator = StatevectorSimulator(**self.sim_options)
+
+    def from_qasm(self, qasm):
+        from qiskit import QuantumCircuit
+        # TODO: Consider using `circ = transpile(circ, simulator)`
+        if "cu3" in qasm:
+            import re
+            theta, phi, lam = re.findall(r"cu3\((.*)\)", qasm)[0].split(",")
+            gamma = - (float(phi) + float(lam)) / 2
+            qasm = re.sub(rf"cu3\((.*)\)",
+                          f"cu({theta},{phi},{lam},{gamma})",
+                          qasm)
+        return QuantumCircuit.from_qasm_str(qasm)
+
+    def __call__(self, circuit):
+        result = self.simulator.run(circuit).result()
+        return result.get_statevector(circuit)
+
+    def get_precision(self):
+        return self.sim_options.get("precision")
+
+    def set_precision(self, precision):
+        from qiskit.providers.aer import StatevectorSimulator
+        self.sim_options["precision"] = precision
+        self.simulator = StatevectorSimulator(**self.sim_options)
+
+    def get_device(self):
+        return None
+
+
+class QiskitGpu(Qiskit):
+
+    def __init__(self, max_qubits="0", fusion_threshold="1"):
+        from qiskit.providers.aer import StatevectorSimulator
+        super().__init__(max_qubits)
+        self.name = "qiskit-gpu"
+        self.sim_options = dict(
+                device="GPU",
+                fusion_enable=self.max_qubits > 0,
+                fusion_max_qubit=self.max_qubits,
+                fusion_threshold=int(fusion_threshold),
+                precision="double"
+            )
+        self.simulator = StatevectorSimulator(**self.sim_options)

qibojit-benchmarks/benchmarks/libraries/qulacs.py

@@ -0,0 +1,85 @@
+import numpy as np
+from benchmarks.libraries import abstract
+
+
+class Qulacs(abstract.ParserBackend):
+
+    def __init__(self):
+        import qulacs
+        self.name = "qulacs"
+        self.qulacs = qulacs
+        self.__version__ = None
+        self.QuantumState = self.qulacs.QuantumState
+
+    def RX(self, target, theta):
+        return self.qulacs.gate.RX(target, -theta)
+
+    def RY(self, target, theta):
+        return self.qulacs.gate.RY(target, -theta)
+
+    def RZ(self, target, theta):
+        return self.qulacs.gate.RZ(target, -theta)
+
+    def CU1(self, control, target, theta):
+        # See `https://github.com/qulacs/qulacs/issues/278` for CU1 on Qulacs
+        matrix = np.diag([1, np.exp(1j * theta)])
+        gate = self.qulacs.gate.DenseMatrix([target], matrix)
+        gate.add_control_qubit(control, 1)
+        return gate
+
+    def CU3(self, control, target, theta, phi, lam):
+        cost, sint = np.cos(theta / 2.0), np.sin(theta / 2.0)
+        pplus, pminus = np.exp(0.5j * (phi + lam)), np.exp(0.5j * (phi - lam))
+        matrix = np.array([[np.conj(pplus) * cost, -np.conj(pminus) * sint],
+                           [pminus * sint, pplus * cost]])
+        gate = self.qulacs.gate.DenseMatrix([target], matrix)
+        gate.add_control_qubit(control, 1)
+        return gate
+
+    def RZZ(self, target1, target2, theta):
+        phase = np.exp(0.5j * theta)
+        phasec = np.conj(phase)
+        matrix = np.diag([phasec, phase, phase, phasec])
+        gate = self.qulacs.gate.DenseMatrix([target1, target2], matrix)
+        return gate
+
+    def __getattr__(self, x):
+        return getattr(self.qulacs.gate, x)
+
+    def __getitem__(self, x):
+        return getattr(self.qulacs.gate, x)
+
+    def from_qasm(self, qasm):
+        nqubits, gatelist = self.parse(qasm)
+        circuit = self.qulacs.QuantumCircuit(nqubits)
+        for gatename, qubits, params in gatelist:
+            gate = getattr(self, gatename)
+            args = list(qubits)
+            if params is not None:
+                args.extend(params)
+            circuit.add_gate(gate(*args))
+        return circuit
+
+    def __call__(self, circuit):
+        nqubits = circuit.get_qubit_count()
+        state = self.QuantumState(nqubits)
+        circuit.update_quantum_state(state)
+        return state.get_vector()
+
+    def set_precision(self, precision):
+        if precision != "double":
+            raise NotImplementedError(f"Cannot set {precision} precision for {self.name} backend.")
+
+    def get_precision(self):
+        return "double"
+
+    def get_device(self):
+        return None
+
+
+class QulacsGpu(Qulacs):
+
+    def __init__(self):
+        super().__init__()
+        self.name = "qulacs-gpu"
+        self.QuantumState = self.qulacs.QuantumStateGpu