diff --git a/.gitignore b/.gitignore index c14e9a3a747b0d74dfbc5f565a1dae400ec0c3a7..87560ae1ef61b0a4d4b2df573d9e229698d51e01 100644 --- a/.gitignore +++ b/.gitignore @@ -140,4 +140,6 @@ cython_debug/ .idea .vscode .DS_Store -.dylibs \ No newline at end of file +.dylibs +c_fusion/ +wheelhouse/ \ No newline at end of file diff --git a/cqlib/simulator/mergy.py b/cqlib/simulator/mergy.py index 5cf11fd4210ac27cb9f9b7e4d4fb85a2ea6e13d9..6775babc394b00f19ade20bb250721dcce66e183 100644 --- a/cqlib/simulator/mergy.py +++ b/cqlib/simulator/mergy.py @@ -55,13 +55,22 @@ class Gate: else: self.mat = np.asarray(mat) + def __repr__(self) -> str: + return (f"\nname: {self.name}\nqubits: {self.qubits}\n" + f"theta: {self.theta}\nmat: {self.mat}\n") + class FusionGate: """ Represents a fused gate """ - def __init__(self, name: str, qubits: list[int] = None, idx: list = None): + def __init__( + self, + name: str, + qubits: list[int] = None, + idx: list = None + ): """ Initializes a FusionGate object. @@ -72,14 +81,16 @@ class FusionGate: """ self.name = name if qubits is None: - self.qubits = [] - else: - self.qubits = qubits + raise TypeError("Qubits is needed.") + self.qubits = qubits if idx is None: self.idx = [] else: self.idx = idx + def __repr__(self) -> str: + return f"\nname: {self.name}\nqubits: {self.qubits}\nidx: {self.idx}\n" + def to_dag(gates: Sequence[Gate]) -> nx.DiGraph: """ @@ -93,36 +104,38 @@ def to_dag(gates: Sequence[Gate]) -> nx.DiGraph: """ dg = nx.DiGraph() pre_nodes = defaultdict(lambda: -1) + for idx, gate in enumerate(gates): - if len(gate.qubits) == 1: - qubit = gate.qubits[0] - node = pre_nodes[qubit] - if node == -1: - dg.add_node(idx, gate=FusionGate(gate.name, gate.qubits.copy(), [idx])) - pre_nodes[qubit] = idx - else: - dg.nodes[node]["gate"].idx.append(idx) - dg.nodes[node]["gate"].name = "fgate" + pre = [pre_nodes[q] for q in gate.qubits] + dg.add_node(idx, gate=FusionGate(gate.name, gate.qubits.copy(), [idx])) + for p in pre: + if p != -1: + dg.add_edge(p, idx) + predecessors = list(dg.predecessors(idx)) + + if len(predecessors) == 1: + gate1 = dg.nodes[idx]["gate"] + gate2 = dg.nodes[predecessors[0]]["gate"] + if set(gate1.qubits) & set(gate2.qubits) == set(gate1.qubits): + dg.nodes[predecessors[0]]["gate"].idx += gate1.idx + dg.nodes[predecessors[0]]["gate"].name = "fgate" + dg.remove_node(idx) else: - c, t = gate.qubits - if pre_nodes[c] == pre_nodes[t] and pre_nodes[c] != -1: - # dg.nodes[pre_nodes[c]]['gate'].mat = gate.mat @ dg.nodes[pre_nodes[c]]['gate'].mat - dg.nodes[pre_nodes[c]]["gate"].idx.append(idx) - dg.nodes[pre_nodes[c]]["gate"].name = "fgate" - else: - dg.add_node(idx, gate=FusionGate(gate.name, gate.qubits.copy(), [idx])) - for predecessor in [pre_nodes[c], pre_nodes[t]]: - if predecessor == -1: - continue - gate2 = dg.nodes[predecessor]["gate"] - if len(gate2.qubits) == 1: + gate1 = dg.nodes[idx]["gate"] + for predecessor in predecessors: + gate2 = dg.nodes[predecessor]["gate"] + if set(gate1.qubits) & set(gate2.qubits) == set(gate2.qubits): + if all(pre_nodes[q] == predecessor for q in gate2.qubits): dg.nodes[idx]["gate"].idx += gate2.idx dg.nodes[idx]["gate"].name = "fgate" + for j in dg.predecessors(predecessor): + dg.add_edge(j, idx) dg.remove_node(predecessor) - else: - dg.add_edge(predecessor, idx) - pre_nodes[c] = idx - pre_nodes[t] = idx + + if dg.has_node(idx): + for q in gate.qubits: + pre_nodes[q] = idx + return dg @@ -191,15 +204,11 @@ class CostBasedFusion: list: A list of fused gates. """ fusions = [] - costs = [] - fusion_to = [] - fusion_to.append(fusion_start) - costs.append( - np.power( - self.cost_factor, - max(len(fusion_gates[fusion_start].qubits) - 1, 1), - ) - ) + fusion_to = [fusion_start] + costs = [np.power( + self.cost_factor, + max(len(fusion_gates[fusion_start].qubits) - 1, 1), + )] for i in range(fusion_start + 1, fusion_end): fusion_to.append(i) @@ -235,7 +244,8 @@ class CostBasedFusion: # return fusions return self.convert(fusions) - def convert(self, fusions: list): + @staticmethod + def convert(fusions: list): dg = nx.DiGraph() pre_nodes = defaultdict(lambda: -1) for idx, gate in enumerate(fusions): @@ -269,62 +279,59 @@ class CostBasedFusion: for q in gate.qubits: pre_nodes[q] = idx - fusionsgate = [] + fusions_gate = [] for node in dg.nodes: - fusionsgate.append(dg.nodes[node]["gate"]) - fusionsgate[-1].idx.sort() - return fusionsgate + fusions_gate.append(dg.nodes[node]["gate"]) + fusions_gate[-1].idx.sort() + return fusions_gate -def get_mat(gate: Gate, qnum: int, qmap: dict) -> np.ndarray: +def get_mat( + gates: list[Gate], + fusion: FusionGate, + qnum: int, + qmap: dict +) -> np.ndarray: """ Generates the matrix for a given gate and maps it to the appropriate qubits. Args: - gate (Gate): The quantum gate. + gates (list[Gate]): The quantum gate list. qnum (int): The total number of qubits. qmap (dict): A dictionary mapping the qubit indices to the correct position. Returns: np.ndarray: The matrix representation of the gate applied to the correct qubits. """ - qubits = [qmap[q] for q in gate.qubits] - if len(qubits) == 1: - q = qubits[0] - mat = np.eye(2 ** q) - mat = np.kron(mat, gate.mat) - mat = np.kron(mat, np.eye(2 ** (qnum - 1 - q))) - return mat - else: - c, t = qubits - mat1 = np.eye(2 ** c) - mat1 = np.kron(mat1, np.diag([1, 0])) - mat1 = np.kron(mat1, np.eye(2 ** (qnum - 1 - c))) - - if c < t: - mat2 = np.eye(2 ** c) - mat2 = np.kron(mat2, np.diag([0, 1])) - mat2 = np.kron(mat2, np.eye(2 ** (t - 1 - c))) - if gate.name in ["CNOT", "CX"]: - mat2 = np.kron(mat2, x_mat) - if gate.name == "CY": - mat2 = np.kron(mat2, y_mat) - if gate.name == "CZ": - mat2 = np.kron(mat2, z_mat) - mat2 = np.kron(mat2, np.eye(2 ** (qnum - 1 - t))) - return mat1 + mat2 - else: - mat2 = np.eye(2 ** t) - if gate.name in ["CNOT", "CX"]: - mat2 = np.kron(mat2, x_mat) - if gate.name == "CY": - mat2 = np.kron(mat2, y_mat) - if gate.name == "CZ": - mat2 = np.kron(mat2, z_mat) - mat2 = np.kron(mat2, np.eye(2 ** (c - 1 - t))) - mat2 = np.kron(mat2, np.diag([0, 1])) - mat2 = np.kron(mat2, np.eye(2 ** (qnum - 1 - c))) - return mat1 + mat2 + shape = [2 ** qnum] + [2 for _ in range(qnum)] + mat = np.eye(2 ** qnum).reshape(shape) + for i in fusion.idx: + op = gates[i] + qubits = [qmap[q] for q in op.qubits] + subscripts = ( + [0] + + [qubit + 1 for qubit in qubits] + + [ + i + for i in range(1, qnum + 1) + if i not in {qubit + 1 for qubit in qubits} + ] + ) + mat = np.reshape(np.transpose(mat, subscripts), (2 ** qnum, 2 ** len(qubits), -1)) + mat = np.matmul(np.reshape(op.mat, (2 ** len(qubits), -1)), mat) + mat = np.transpose( + np.reshape(mat, shape), + _inv_subscripts(subscripts), + ) + mat = mat.reshape(2 ** qnum, -1).T.copy() + return mat + + +def _inv_subscripts(subscripts: list[int]) -> list[int]: + inv = [0] * len(subscripts) + for i, sub in enumerate(subscripts): + inv[sub] = i + return inv def fusion_to_gate(gates: list, fusion: FusionGate) -> Gate: @@ -339,7 +346,6 @@ def fusion_to_gate(gates: list, fusion: FusionGate) -> Gate: Returns: Gate: A new Gate object representing the fused operation. """ - if fusion.name == "fgate": qubits = list(fusion.qubits) qnum = len(qubits) @@ -347,12 +353,9 @@ def fusion_to_gate(gates: list, fusion: FusionGate) -> Gate: qmap = {} for i, v in enumerate(qubits): qmap[v] = i - mat = np.eye(2 ** qnum) - for i in fusion.idx: - mat = get_mat(gates[i], qnum, qmap) @ mat + mat = get_mat(gates, fusion, qnum, qmap) return Gate(fusion.name, qubits, theta=[], mat=mat) - else: - return gates[fusion.idx[0]] + return gates[fusion.idx[0]] def merge_gate(ori_gates: list[Gate], max_qubit: int) -> list[Gate]: @@ -370,17 +373,17 @@ def merge_gate(ori_gates: list[Gate], max_qubit: int) -> list[Gate]: if max_qubit == 1: return ori_gates dg = to_dag(ori_gates) - fusiongates = [] + fusion_gates = [] for node in dg.nodes: - fusiongates.append(dg.nodes[node]["gate"]) + fusion_gates.append(dg.nodes[node]["gate"]) fusion_method = CostBasedFusion() - fusiongates = fusion_method.aggregate_gates( - fusiongates, 0, len(fusiongates), max_qubit + fusion_gates = fusion_method.aggregate_gates( + fusion_gates, 0, len(fusion_gates), max_qubit ) new_gates = [] - for fusion in fusiongates: + for fusion in fusion_gates: new_gates.append(fusion_to_gate(ori_gates, fusion)) return new_gates diff --git a/cqlib/simulator/statevector_simulator.py b/cqlib/simulator/statevector_simulator.py index 8cc12e3e75abd811b219b0199180ead83b007783..8ee99183e19076b1397cb8b4b6663dc95f47e553 100644 --- a/cqlib/simulator/statevector_simulator.py +++ b/cqlib/simulator/statevector_simulator.py @@ -33,14 +33,17 @@ The results will be returned in the order of Q1, Q0. """ import ctypes +import time from collections import Counter import numpy as np from cqlib.circuits.circuit import Circuit from cqlib.circuits.parameter import Parameter +from cqlib.exceptions import CqlibError from cqlib.simulator.mergy import Gate, merge_gate -from cqlib.simulator.statevector_simulator_c import get_measure, get_probs, get_sample, get_state +from cqlib.simulator.statevector_simulator_c import get_measure, get_probs, \ + get_sample, get_state, free_memory gate_name_map = { "H": 72, @@ -61,12 +64,19 @@ gate_name_map = { "RX": 8288, "RY": 8289, "RZ": 8290, + "U3": 8551, "U": 8551, "fgate": 102, "XY": 8889, "XY2M": 88895077, "XY2P": 88895080, "RXY": 828889, + "CCX": 102, + "SWAP": 102, + "CRX": 102, + "CRY": 102, + "CRZ": 102, + "custom_gate": 102, } @@ -136,6 +146,10 @@ class StatevectorSimulator: self.measure_qubits = [] self.qubit_mapping = {q: i for i, q in enumerate(circuit.qubits)} self._parse_circuit() + self.state_ptr_capsule = None + self.probs_ptr_capsule = None + self.measure_ptr_capsule = None + self.samples_ptr_capsule = None def _parse_circuit(self): """ @@ -156,10 +170,10 @@ class StatevectorSimulator: param = float(param.value(params=self.circuit.parameters_value)) ps.append(param) name = item.instruction.name - if name in ['CRX', 'CRY', 'CRZ', 'CCX', 'SWAP']: - name = 'fgate' - if self.is_fusion: - raise ValueError(f"gate {item.instruction.name} not support fusion.") + # if name in ['CRX', 'CRY', 'CRZ', 'CCX', 'SWAP']: + # name = 'fgate' + # if self.is_fusion: + # raise ValueError(f"gate {item.instruction.name} not support fusion.") self.gates.append(Gate( name=name, @@ -210,7 +224,8 @@ class StatevectorSimulator: """ self._check_fusion() gates_list = self.get_gates() - state_array = get_state(self.nq, gates_list, self.omp_threads) + state_array, state_ptr_capsule = get_state(self.nq, 2 ** self.nq, gates_list, self.omp_threads) + self.state_ptr_capsule = state_ptr_capsule state = {np.binary_repr(i, width=self.nq): val for i, val in enumerate(state_array)} return state @@ -223,8 +238,10 @@ class StatevectorSimulator: np.ndarray: A list of probabilities for each possible outcome. """ self._check_fusion() - gates_list = self.get_gates() - probs_array = get_probs(self.nq, gates_list, self.omp_threads) + if self.state_ptr_capsule is None: + self.statevector() + probs_array, probs_ptr_capsule = get_probs(self.nq, 2 ** self.nq, self.omp_threads, self.state_ptr_capsule) + self.probs_ptr_capsule = probs_ptr_capsule return {np.binary_repr(i, width=self.nq): val for i, val in enumerate(probs_array)} def measure(self) -> np.array: @@ -235,16 +252,23 @@ class StatevectorSimulator: Returns: np.ndarray: A probability distribution of the measurement results. """ + self._check_fusion() + if self.probs_ptr_capsule is None: + self.probs() + if not self.measure_qubits: # all measured - return self.probs() - mq_len = len(self.measure_qubits) - assert mq_len == len(set(self.measure_qubits)) - if self.measure_qubits == sorted(self.measure_qubits.copy()) and mq_len == self.nq: - return self.probs() + raise CqlibError("Measured qubits are empty, please add measurement gates first.") - self._check_fusion() - gates_list = self.get_gates() - measure_array = get_measure(self.nq, gates_list, self.omp_threads, self.measure_qubits) + mq_ordered = sorted(self.measure_qubits.copy()) + # Call get_measure from the C extension module + measure_array, measure_ptr_capsule = get_measure( + self.nq, + self.omp_threads, + self.measure_qubits, + mq_ordered, + self.probs_ptr_capsule + ) + self.measure_ptr_capsule = measure_ptr_capsule return {np.binary_repr(i, width=len(self.measure_qubits)): val for i, val in enumerate(measure_array)} def sample( @@ -253,6 +277,7 @@ class StatevectorSimulator: is_sorted: bool = False, sample_block_th: int = 10, is_raw_data: bool = False, + rng_seed: int = None ) -> np.ndarray | dict[str, int]: """ Samples the quantum circuit multiple times, returning either the raw sampled data or a @@ -263,6 +288,7 @@ class StatevectorSimulator: is_sorted (bool): Whether to return the results sorted by state (default: False). sample_block_th (int): Block threshold for sampling optimization (default: 10). is_raw_data (bool): If True, returns raw sample data instead of a frequency dictionary (default: False). + rng_seed (int): Seed for the random number generator (default: None). Returns: np.ndarray | dict[str, int]: The sampled results, either as raw data or as a frequency distribution. @@ -271,21 +297,24 @@ class StatevectorSimulator: self.measure_qubits = list(range(self.nq)) assert shots < 4294967296 # uint32_t - measure_dict = self.measure() # Should return a dictionary {bitstring: probability} - mq_len = len(self.measure_qubits) - # Convert measure_dict to a list of probabilities - measure_probs = np.zeros(2 ** mq_len, dtype=np.float64) - for bitstring, prob in measure_dict.items(): - idx = int(bitstring, 2) - measure_probs[idx] = prob + if self.measure_ptr_capsule is None: + # If measure_ptr_capsule doesn't exist, perform measurement first + self.measure() + + if rng_seed is None: + rng_seed = int(time.time()) - samples_array = get_sample( + # Call get_sample from the C extension module + samples_array, samples_ptr_capsule = get_sample( shots, - measure_probs, + self.measure_ptr_capsule, self.omp_threads, self.nq, - sample_block_th + sample_block_th, + rng_seed ) + # Store samples_ptr_capsule to keep the samples_ptr alive + self.samples_ptr_capsule = samples_ptr_capsule # Check if samples_array is None if samples_array is None: @@ -293,19 +322,24 @@ class StatevectorSimulator: if is_raw_data: return samples_array + + samples_array = samples_array.astype(np.uint64) mq_len = len(self.measure_qubits) - powers_of_two = 2 ** np.arange(mq_len)[::-1] - samples_int = samples_array.dot(powers_of_two).astype(np.uint64) - counts = Counter(samples_int) + # powers_of_two = 2 ** np.arange(mq_len)[::-1] + # samples_int = samples_array.dot(powers_of_two).astype(np.uint64) + counts = Counter(samples_array) # Convert counts to binary strings result = { - np.binary_repr(k, width=mq_len): v for k, v in counts.items() + np.binary_repr(int(k), width=mq_len): v for k, v in counts.items() } if is_sorted: result = dict(sorted(result.items())) return result - def get_gates(self): + def get_gates(self) -> list[dict]: + """ + Retrieves detailed information about all gates in the current quantum circuit. + """ gates_list = [] for gate in self.gates: gate_dict = { diff --git a/setup.py b/setup.py index 8fafce3ec0cd7383e50f90e1b53949f6051a9f53..67b4d5883e2e340429d551ea3e4faea19e3d4ca9 100644 --- a/setup.py +++ b/setup.py @@ -64,10 +64,15 @@ class CustomBuildExt(build_ext_orig): ext_modules = [ Extension( 'cqlib.simulator.statevector_simulator_c', - include_dirs=['cqlib/simulator/c_fusion/include', python_paths['include'], numpy.get_include()], + include_dirs=[ + 'cqlib/simulator/c_fusion/include', + python_paths['include'], + numpy.get_include() + ], library_dirs=[python_paths['platlib']], # 链接库文件所在目录 sources=[ 'cqlib/simulator/c_fusion/src/c_fusion.c', + 'cqlib/simulator/c_fusion/src/state_vector.c', ], define_macros=[('Py_LIMITED_API', '0x030A0000')], py_limited_api=True, diff --git a/tests/simulator/test_statevector_sim.py b/tests/simulator/test_statevector_sim.py index 6de611ed74f9393b0690a92d5d719f76c86922c9..17eec02cd51610c678117f2d5d475a733ea0819d 100644 --- a/tests/simulator/test_statevector_sim.py +++ b/tests/simulator/test_statevector_sim.py @@ -586,7 +586,6 @@ def test_qcis_gate_circuit_2(): def test_swap_circuit(): """ - todo: bug """ c = Circuit(2) c.h(0)