diff --git a/cspell.json b/cspell.json index c428002..ee4bf33 100644 --- a/cspell.json +++ b/cspell.json @@ -18,6 +18,10 @@ "splinalg", "spsparse", "eigvals", + "eigenvals", + "cumtime", + "csvfile", + "yscale", "eigvecs", "eigh", "finfo", diff --git a/qpiai_quantum/quantum_info/density_matrix.py b/qpiai_quantum/quantum_info/density_matrix.py index a7b744d..85c2979 100644 --- a/qpiai_quantum/quantum_info/density_matrix.py +++ b/qpiai_quantum/quantum_info/density_matrix.py @@ -49,29 +49,23 @@ def _get_num_qubits(self) -> int: def _init_from_array(self, data: list | np.ndarray): """Initialize from array.""" - # Handle nested list format - if isinstance(data, list) and len(data) > 0: - if isinstance(data[0], list): - # Could be density matrix or nested statevector - if isinstance(data[0], list): - # Check for dict-based density matrix - if ( - isinstance(data[0][0], dict) - and "real" in data[0][0] - and "imag" in data[0][0] - ): - self._init_from_dict_density_matrix(data) - return - - # 2D numeric density matrix - if isinstance(data[0][0], list): - self.data = np.array(data, dtype=complex) - else: - # Nested statevector [[c1], [c2], ...] + if isinstance(data, np.ndarray): + self.data = data.astype(complex) + elif isinstance(data, list) and len(data) > 0: + first = data[0] + if isinstance(first, list): + if len(first) > 0 and isinstance(first[0], dict) and "real" in first[0]: + self._init_from_dict_density_matrix(data) + return + elif len(first) == 1 and not isinstance(first[0], (list, dict)): + # Nested column statevector [[c1], [c2], ...] statevector = np.array([item[0] for item in data], dtype=complex) self._init_from_statevector(statevector) return - elif isinstance(data[0], dict) and "real" in data[0] and "imag" in data[0]: + else: + # 2D numeric / complex density matrix + self.data = np.array(data, dtype=complex) + elif isinstance(first, dict) and "real" in first and "imag" in first: # Dict format statevector statevector = np.array( [complex(item["real"], item["imag"]) for item in data], diff --git a/qpiai_quantum/simulator/__init__.py b/qpiai_quantum/simulator/__init__.py index a0e273b..a1bf33f 100644 --- a/qpiai_quantum/simulator/__init__.py +++ b/qpiai_quantum/simulator/__init__.py @@ -7,9 +7,14 @@ from .base_simulator import BaseSimulator from .statevector import StatevectorSimulator from .result import QasmSimulatorResult +from .state import QuantumState +from .density_matrix import DensityMatrixState, DensityMatrixSimulator __all__ = [ "BaseSimulator", "StatevectorSimulator", "QasmSimulatorResult", + "QuantumState", + "DensityMatrixState", + "DensityMatrixSimulator", ] diff --git a/qpiai_quantum/simulator/benchmark_statevector_speed.py b/qpiai_quantum/simulator/benchmark_statevector_speed.py new file mode 100644 index 0000000..64a7636 --- /dev/null +++ b/qpiai_quantum/simulator/benchmark_statevector_speed.py @@ -0,0 +1,761 @@ +import cProfile +import csv +import io +import json +import pstats +import statistics +import sys +import time +import types +from dataclasses import dataclass +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.backends.backend_pdf import PdfPages + + +ROOT = Path(__file__).resolve().parents[2] +if "qpiai_quantum" not in sys.modules: + qpiai_quantum = types.ModuleType("qpiai_quantum") + qpiai_quantum.__path__ = [str(ROOT / "qpiai_quantum")] + sys.modules["qpiai_quantum"] = qpiai_quantum + +from qpiai_quantum.circuit.circuit import Circuit # noqa: E402 +from qpiai_quantum.icr.circuitoperation import OperationType # noqa: E402 +from qpiai_quantum.simulator.gates import ( # noqa: E402 + DECOMPOSED_GATES, + decompose, + gate_spec, +) +from qpiai_quantum.simulator.result import QasmSimulatorResult # noqa: E402 +from qpiai_quantum.simulator.statevector import StatevectorSimulator # noqa: E402 + + +OPTIMIZED_GATES = ("h", "x", "y", "z", "s", "rx", "ry", "rz") +PARAMETRIC_GATES = {"rx", "ry", "rz"} +ITERATIONS = 30 +WARMUPS = 5 +RNG_SEED = 20260711 + + +@dataclass(frozen=True) +class BenchmarkCase: + name: str + workload: str + n_qubits: int + specs: tuple[tuple[str, tuple[int, ...], tuple[float, ...]], ...] + + +class OriginalStatevectorSimulator(StatevectorSimulator): + def run( + self, + circuit: Circuit, + shots: int = 1024, + seed: int | None = None, + name: str | None = None, + ) -> QasmSimulatorResult: + start_time = time.perf_counter() + n_qubits = circuit.num_qubits + n_cbits = circuit.num_clbits + + if n_qubits == 0: + raise ValueError("Cannot simulate a circuit with 0 qubits.") + + state = np.zeros(2**n_qubits, dtype=complex) + state[0] = 1.0 + measure_map: dict[int, int] = {} + + def _apply_gate(gate_name: str, params: list[float], qubits: list[int]) -> None: + nonlocal state + gate_name_lower = gate_name.lower() + + if gate_name_lower in DECOMPOSED_GATES: + for sub_name, sub_params, sub_qubits in decompose( + gate_name_lower, qubits + ): + _apply_gate(sub_name, sub_params, sub_qubits) + return + + _, unitary = gate_spec(gate_name_lower, params, num_qubits=len(qubits)) + state = self._apply_unitary(state, n_qubits, qubits, unitary) + + def _apply_op(op) -> None: + nonlocal state + if op.operation_type == OperationType.BARRIER: + return + + if op.operation_type == OperationType.MEASURE: + if op.qubits and op.clbits: + for q, c in zip(op.qubits, op.clbits): + measure_map[q] = c + return + + if op.operation_type in ( + OperationType.N_QUBIT_NON_PARAMETRIC, + OperationType.N_QUBIT_PARAMETRIC, + OperationType.SWAP, + ): + _apply_gate(op.gate_name, op.params or [], op.qubits) + return + + if op.operation_type == OperationType.OPERATION: + if hasattr(op, "order") and op.order is not None: + for sub_op in op.order: + _apply_op(sub_op) + else: + _apply_gate(op.gate_name, op.params or [], op.qubits) + return + + if op.gate_name.lower() == "reset" and op.qubits: + for q in op.qubits: + state = self._apply_reset(state, n_qubits, q) + else: + raise ValueError(f"Unsupported operation type: {op.operation_type}") + + for op in circuit.icr.evolve: + _apply_op(op) + + counts = ( + self._sample_counts(state, n_qubits, n_cbits, measure_map, shots, seed) + if n_cbits > 0 and measure_map + else {} + ) + + return QasmSimulatorResult( + name=name or circuit.name, + counts=counts, + statevector=state.tolist(), + shots=shots, + executionTime=time.perf_counter() - start_time, + method="statevector-original", + job_status="completed", + n_qubits=n_qubits, + n_cbits=n_cbits, + ) + + +def apply_spec(circuit: Circuit, spec): + gate, qubits, params = spec + if gate in PARAMETRIC_GATES: + getattr(circuit, gate)(qubits[0], params[0]) + elif gate == "cx": + circuit.cx(qubits[0], qubits[1]) + elif gate == "cz": + circuit.cz(qubits[0], qubits[1]) + elif gate == "swap": + circuit.swap(qubits[0], qubits[1]) + else: + getattr(circuit, gate)(qubits[0]) + + +def build_circuit(n_qubits: int, specs) -> Circuit: + circuit = Circuit(n_qubits) + for spec in specs: + apply_spec(circuit, spec) + return circuit + + +def statevector(simulator, circuit: Circuit) -> np.ndarray: + return np.asarray(simulator.run(circuit, shots=1).statevector, dtype=complex) + + +def describe_specs(specs) -> str: + lines = [] + for index, (gate, qubits, params) in enumerate(specs, start=1): + suffix = f", params={list(params)}" if params else "" + lines.append(f"{index}. {gate.upper()} qubits={list(qubits)}{suffix}") + return "\n".join(lines) + + +def first_divergence(n_qubits: int, specs) -> str: + original = OriginalStatevectorSimulator() + optimized = StatevectorSimulator() + for index in range(1, len(specs) + 1): + prefix = specs[:index] + circuit = build_circuit(n_qubits, prefix) + original_state = statevector(original, circuit) + optimized_state = statevector(optimized, circuit) + if not np.allclose(original_state, optimized_state): + gate, qubits, params = specs[index - 1] + return ( + f"first divergent gate {index}: " + f"{gate.upper()} qubits={list(qubits)} params={list(params)}" + ) + return "no divergent prefix found" + + +def verify_case(case: BenchmarkCase): + circuit = build_circuit(case.n_qubits, case.specs) + original_state = statevector(OriginalStatevectorSimulator(), circuit) + optimized_state = statevector(StatevectorSimulator(), circuit) + + if not np.allclose(original_state, optimized_state): + detail = first_divergence(case.n_qubits, case.specs) + raise RuntimeError( + f"Verification failed for {case.name}\n" + f"Circuit:\n{describe_specs(case.specs)}\n" + f"Original statevector:\n{original_state}\n" + f"Optimized statevector:\n{optimized_state}\n" + f"{detail}\n" + "Reason: the optimized pair-wise kernel is not mathematically " + "equivalent to the original matrix path for this prefix." + ) + + if case.workload == "Bell state": + expected = np.zeros(4, dtype=complex) + expected[0] = 1 / np.sqrt(2) + expected[3] = 1 / np.sqrt(2) + if not np.allclose(original_state, expected) or not np.allclose( + optimized_state, expected + ): + raise RuntimeError( + "Bell-state expected-value verification failed.\n" + f"Original statevector: {original_state}\n" + f"Optimized statevector: {optimized_state}\n" + f"Expected statevector: {expected}" + ) + + return circuit + + +def summarize(samples: list[float]) -> dict[str, float]: + return { + "mean": statistics.fmean(samples), + "median": statistics.median(samples), + "min": min(samples), + "max": max(samples), + "stddev": statistics.stdev(samples), + } + + +def benchmark_case(case: BenchmarkCase) -> dict[str, object]: + circuit = verify_case(case) + original = OriginalStatevectorSimulator() + optimized = StatevectorSimulator() + + for _ in range(WARMUPS): + original.run(circuit, shots=1) + optimized.run(circuit, shots=1) + + original_times = [] + optimized_times = [] + for _ in range(ITERATIONS): + start = time.perf_counter() + original.run(circuit, shots=1) + original_times.append(time.perf_counter() - start) + + start = time.perf_counter() + optimized.run(circuit, shots=1) + optimized_times.append(time.perf_counter() - start) + + original_stats = summarize(original_times) + optimized_stats = summarize(optimized_times) + return { + "name": case.name, + "workload": case.workload, + "n_qubits": case.n_qubits, + "depth": len(case.specs), + "iterations": ITERATIONS, + "warmups": WARMUPS, + "original": original_stats, + "optimized": optimized_stats, + "speedup": original_stats["mean"] / optimized_stats["mean"], + } + + +def repeated_gate_specs(gate: str, n_qubits: int, depth: int, seed: int): + rng = np.random.default_rng(seed) + specs = [] + for index in range(depth): + if n_qubits >= 15: + target = n_qubits // 2 + (index % max(1, n_qubits - n_qubits // 2)) + else: + target = index % n_qubits + params = ( + (float(rng.uniform(-2 * np.pi, 2 * np.pi)),) + if gate in PARAMETRIC_GATES + else () + ) + specs.append((gate, (target,), params)) + return tuple(specs) + + +def random_single_qubit_specs(n_qubits: int, depth: int, seed: int): + rng = np.random.default_rng(seed) + low_target = n_qubits // 2 if n_qubits >= 15 else 0 + specs = [] + for _ in range(depth): + gate = str(rng.choice(OPTIMIZED_GATES)) + target = int(rng.integers(low_target, n_qubits)) + params = ( + (float(rng.uniform(-2 * np.pi, 2 * np.pi)),) + if gate in PARAMETRIC_GATES + else () + ) + specs.append((gate, (target,), params)) + return tuple(specs) + + +def mixed_specs(n_qubits: int, depth: int, seed: int): + rng = np.random.default_rng(seed) + low_target = n_qubits // 2 if n_qubits >= 15 else 0 + specs = [] + for index in range(depth): + if n_qubits > 1 and index % 9 == 8: + first = int(rng.integers(0, n_qubits)) + second = int((first + rng.integers(1, n_qubits)) % n_qubits) + specs.append(("swap", (first, second), ())) + elif n_qubits > 1 and index % 5 == 4: + control = int(rng.integers(0, n_qubits)) + target = int((control + rng.integers(1, n_qubits)) % n_qubits) + specs.append(("cx", (control, target), ())) + elif n_qubits > 1 and index % 7 == 6: + control = int(rng.integers(0, n_qubits)) + target = int((control + rng.integers(1, n_qubits)) % n_qubits) + specs.append(("cz", (control, target), ())) + else: + gate = str(rng.choice(OPTIMIZED_GATES)) + target = int(rng.integers(low_target, n_qubits)) + params = ( + (float(rng.uniform(-2 * np.pi, 2 * np.pi)),) + if gate in PARAMETRIC_GATES + else () + ) + specs.append((gate, (target,), params)) + return tuple(specs) + + +def build_cases() -> list[BenchmarkCase]: + cases: list[BenchmarkCase] = [] + stress_depths = {5: 120, 10: 60, 15: 8, 20: 1} + random_depths = {5: 140, 10: 70, 15: 10, 20: 1} + mixed_depths = {5: 100, 10: 50, 15: 10, 20: 1} + + for n_qubits, depth in stress_depths.items(): + for gate in OPTIMIZED_GATES: + cases.append( + BenchmarkCase( + name=f"{gate.upper()} repeated, {n_qubits}q", + workload="Single-gate stress", + n_qubits=n_qubits, + specs=repeated_gate_specs( + gate, n_qubits, depth, RNG_SEED + n_qubits + ), + ) + ) + + for n_qubits, depth in random_depths.items(): + cases.append( + BenchmarkCase( + name=f"Random optimized gates, {n_qubits}q", + workload="Random single-qubit", + n_qubits=n_qubits, + specs=random_single_qubit_specs( + n_qubits, depth, RNG_SEED + 100 + n_qubits + ), + ) + ) + + for n_qubits, depth in mixed_depths.items(): + cases.append( + BenchmarkCase( + name=f"Mixed 1q + CX/CZ/SWAP, {n_qubits}q", + workload="Mixed", + n_qubits=n_qubits, + specs=mixed_specs(n_qubits, depth, RNG_SEED + 200 + n_qubits), + ) + ) + + cases.extend( + [ + BenchmarkCase( + name="Bell state H(0), CX(0, 1)", + workload="Bell state", + n_qubits=2, + specs=(("h", (0,), ()), ("cx", (0, 1), ())), + ), + BenchmarkCase( + name="Large random optimized gates, 12q depth 240", + workload="Large stress", + n_qubits=12, + specs=random_single_qubit_specs(12, 240, RNG_SEED + 300), + ), + BenchmarkCase( + name="Large random optimized gates, 15q depth 16", + workload="Large stress", + n_qubits=15, + specs=random_single_qubit_specs(15, 16, RNG_SEED + 301), + ), + BenchmarkCase( + name="Large mixed circuit, 12q depth 160", + workload="Large stress", + n_qubits=12, + specs=mixed_specs(12, 160, RNG_SEED + 302), + ), + ] + ) + return cases + + +def profile_optimized() -> dict[str, object]: + case = BenchmarkCase( + name="Profile random optimized gates, 12q depth 240", + workload="Profiling", + n_qubits=12, + specs=random_single_qubit_specs(12, 240, RNG_SEED + 400), + ) + circuit = verify_case(case) + simulator = StatevectorSimulator() + for _ in range(3): + simulator.run(circuit, shots=1) + + profiler = cProfile.Profile() + profiler.enable() + simulator.run(circuit, shots=1) + profiler.disable() + + stream = io.StringIO() + stats = pstats.Stats(profiler, stream=stream).strip_dirs().sort_stats("cumulative") + stats.print_stats(12) + entries = [] + total = stats.total_tt + for key, value in stats.stats.items(): + filename, line, func_name = key + call_count, primitive, total_time, cumulative_time, callers = value + entries.append( + { + "function": f"{func_name} ({filename}:{line})", + "calls": call_count, + "total_s": total_time, + "cumulative_s": cumulative_time, + "percent_total": (total_time / total * 100) if total else 0.0, + } + ) + entries.sort(key=lambda item: item["cumulative_s"], reverse=True) + return { + "case": case.name, + "total_profiled_s": total, + "top_entries": entries[:12], + "raw": stream.getvalue(), + } + + +def machine_info() -> dict[str, str]: + import platform + + cpu = platform.processor() or platform.machine() or "not reported by OS" + return { + "python": platform.python_version(), + "numpy": np.__version__, + "platform": platform.platform(), + "cpu": cpu, + } + + +def fmt_ms(value: float) -> str: + return f"{value * 1000:.3f}" + + +def write_csv(results: list[dict[str, object]], output: Path): + with output.open("w", newline="", encoding="utf-8") as csv_file: + writer = csv.writer(csv_file) + writer.writerow( + [ + "workload", + "name", + "n_qubits", + "depth", + "iterations", + "warmups", + "original_mean_s", + "optimized_mean_s", + "original_median_s", + "optimized_median_s", + "original_min_s", + "optimized_min_s", + "original_max_s", + "optimized_max_s", + "original_stddev_s", + "optimized_stddev_s", + "speedup", + ] + ) + for row in results: + writer.writerow( + [ + row["workload"], + row["name"], + row["n_qubits"], + row["depth"], + row["iterations"], + row["warmups"], + row["original"]["mean"], + row["optimized"]["mean"], + row["original"]["median"], + row["optimized"]["median"], + row["original"]["min"], + row["optimized"]["min"], + row["original"]["max"], + row["optimized"]["max"], + row["original"]["stddev"], + row["optimized"]["stddev"], + row["speedup"], + ] + ) + + +def draw_text_page(pdf, title: str, lines: list[str]): + fig = plt.figure(figsize=(8.27, 11.69)) + fig.patch.set_facecolor("white") + plt.axis("off") + fig.text(0.07, 0.96, title, fontsize=15, fontweight="bold", va="top") + y = 0.91 + for line in lines: + if y < 0.05: + pdf.savefig(fig, bbox_inches="tight") + plt.close(fig) + fig = plt.figure(figsize=(8.27, 11.69)) + fig.patch.set_facecolor("white") + plt.axis("off") + y = 0.96 + is_heading = line.endswith(":") + fig.text( + 0.08, + y, + line, + fontsize=10.2, + fontweight="bold" if is_heading else "normal", + va="top", + ) + y -= 0.034 if line else 0.018 + pdf.savefig(fig, bbox_inches="tight") + plt.close(fig) + + +def draw_table_page(pdf, title: str, rows: list[dict[str, object]]): + columns = [ + "Workload", + "Case", + "Qubits", + "Depth", + "Orig mean ms", + "Opt mean ms", + "Speedup", + "Opt std ms", + ] + per_page = 17 + for start in range(0, len(rows), per_page): + chunk = rows[start : start + per_page] + fig, ax = plt.subplots(figsize=(11.69, 8.27)) + ax.axis("off") + ax.set_title(title, fontsize=15, fontweight="bold", pad=18) + data = [ + [ + row["workload"], + row["name"], + row["n_qubits"], + row["depth"], + fmt_ms(row["original"]["mean"]), + fmt_ms(row["optimized"]["mean"]), + f"{row['speedup']:.3f}x", + fmt_ms(row["optimized"]["stddev"]), + ] + for row in chunk + ] + table = ax.table( + cellText=data, + colLabels=columns, + loc="center", + cellLoc="left", + colLoc="left", + ) + table.auto_set_font_size(False) + table.set_fontsize(7.4) + table.scale(1, 1.48) + pdf.savefig(fig, bbox_inches="tight") + plt.close(fig) + + +def draw_runtime_plot(pdf, results: list[dict[str, object]]): + representative = [ + row + for row in results + if row["workload"] in ("Random single-qubit", "Mixed", "Large stress") + ] + labels = [row["name"].replace(" gates, ", "\n") for row in representative] + x = np.arange(len(labels)) + original = [row["original"]["mean"] * 1000 for row in representative] + optimized = [row["optimized"]["mean"] * 1000 for row in representative] + + fig, ax = plt.subplots(figsize=(11.69, 8.27)) + width = 0.38 + ax.bar(x - width / 2, original, width, label="Original") + ax.bar(x + width / 2, optimized, width, label="Optimized pair-wise") + ax.set_yscale("log") # cspell:disable-line + ax.set_ylabel("Mean runtime (ms, log scale)") + ax.set_title("Runtime comparison for representative workloads") + ax.set_xticks(x) + ax.set_xticklabels(labels, rotation=35, ha="right") + ax.legend() + ax.grid(axis="y", alpha=0.25) + pdf.savefig(fig, bbox_inches="tight") + plt.close(fig) + + +def draw_speedup_plot(pdf, results: list[dict[str, object]]): + fig, ax = plt.subplots(figsize=(11.69, 8.27)) + for workload in ("Random single-qubit", "Mixed"): + rows = [row for row in results if row["workload"] == workload] + rows.sort(key=lambda row: row["n_qubits"]) + ax.plot( + [row["n_qubits"] for row in rows], + [row["speedup"] for row in rows], + marker="o", + label=workload, + ) + ax.axhline(1.0, color="black", linewidth=1, linestyle="--", label="break-even") + ax.set_xlabel("Qubit count") + ax.set_ylabel("Speedup (Original / Optimized)") + ax.set_title("Speedup versus qubit count") + ax.legend() + ax.grid(alpha=0.25) + pdf.savefig(fig, bbox_inches="tight") + plt.close(fig) + + +def write_report( + results: list[dict[str, object]], + profile: dict[str, object], + output: Path, +): + info = machine_info() + speedups = [float(row["speedup"]) for row in results] + faster = [row for row in results if float(row["speedup"]) > 1.0] + materially_faster = [row for row in results if float(row["speedup"]) >= 1.10] + random_rows = [row for row in results if row["workload"] == "Random single-qubit"] + mixed_rows = [row for row in results if row["workload"] == "Mixed"] + large_rows = [row for row in results if row["workload"] == "Large stress"] + recommendation = ( + "Submit as a PR only if the goal is code clarity or a foundation for future " + "vectorized kernels; based on these measurements, the performance case alone " + "is not strong enough." + if len(materially_faster) <= len(results) // 2 + else "The optimization shows enough broad performance improvement to justify a PR." + ) + + with PdfPages(output) as pdf: + intro = [ + "Benchmark methodology:", + "Original path: every gate uses gate_spec(...) followed by _apply_unitary(...).", + "Optimized path: H, X, Y, Z, S, RX, RY, and RZ use the current pair-wise traversal kernel.", + "Each benchmark uses the same Circuit object for both simulators.", + "Before timing, each case verifies np.allclose(original_statevector, optimized_statevector).", + "The Bell benchmark also checks the expected (|00> + |11>) / sqrt(2) state.", + f"Warm-up runs per case: {WARMUPS}.", + f"Measured iterations per case: {ITERATIONS}.", + "Timer: time.perf_counter().", + "Random seeds are fixed for reproducibility.", + "For 15- and 20-qubit workloads, single-qubit targets are limited to middle/high indices so the pair-wise implementation is computationally feasible under 30 iterations.", + "This is still a fair original-vs-optimized comparison because both simulators execute the same circuits.", + "", + "Environment:", + f"Python version: {info['python']}", + f"NumPy version: {info['numpy']}", + f"Platform: {info['platform']}", + f"CPU: {info['cpu']}", + "", + "Correctness verification summary:", + f"Verified benchmark cases: {len(results)}.", + "All benchmark cases passed original-vs-optimized statevector equality.", + "No benchmark timings were recorded for an incorrect case.", + ] + draw_text_page( + pdf, + "QpiAI Quantum SDK: Pair-Wise Single-Qubit Gate Benchmark", + intro, + ) + + summary = [ + "Summary statistics:", + f"Cases faster than original: {len(faster)} / {len(results)}.", + f"Cases at least 10 percent faster: {len(materially_faster)} / {len(results)}.", + f"Median speedup: {statistics.median(speedups):.3f}x.", + f"Mean speedup: {statistics.fmean(speedups):.3f}x.", + f"Best speedup: {max(speedups):.3f}x.", + f"Worst speedup: {min(speedups):.3f}x.", + "", + "Interpretation:", + "Speedup above 1.0 means the pair-wise optimized implementation is faster.", + "Speedup below 1.0 means the original NumPy matrix path is faster.", + "The pair-wise traversal removes redundant pair visits and avoids per-amplitude bit/flipped calculations.", + "However, this implementation still allocates a result array per gate and creates offset index arrays inside each stride block.", + "For low target qubits, many small np.arange/indexing operations create noticeable overhead.", + "For high target qubits, fewer but larger indexed slices are used, shifting cost toward array allocation and vectorized arithmetic.", + "The original _apply_unitary path delegates most arithmetic and reshaping work to optimized NumPy internals.", + "", + "Final recommendation:", + recommendation, + "The improvement is statistically meaningful only for cases where the speedup is consistently above 1.0 and materially beyond run-to-run noise.", + "The added code complexity is justified only if follow-up work removes allocation/index creation overhead or specializes common targets further.", + ] + draw_text_page(pdf, "Analysis and Recommendation", summary) + + draw_table_page(pdf, "Runtime and speedup table", results) + draw_runtime_plot(pdf, results) + draw_speedup_plot(pdf, results) + + profile_lines = [ + "Profiling setup:", + f"Profiled case: {profile['case']}.", + f"Profiled total Python time: {profile['total_profiled_s']:.6f} seconds.", + "", + "Top cumulative-time functions:", + ] + for entry in profile["top_entries"]: + profile_lines.append( + f"{entry['function']}: calls={entry['calls']}, " + f"total={entry['total_s']:.6f}s, cumulative={entry['cumulative_s']:.6f}s" + ) + profile_lines.extend( + [ + "", + "Profiling interpretation:", + "The dominant cost is the optimized gate dispatcher/kernel and the NumPy operations it performs inside each gate.", + "The main bottlenecks are repeated result allocation, offset index construction, advanced indexing gathers/scatters, and arithmetic on temporary arrays.", + "Dispatcher overhead is present but smaller than per-gate array work for larger statevectors.", + ] + ) + draw_text_page(pdf, "Profiling Summary", profile_lines) + + draw_table_page(pdf, "Random-circuit workloads", random_rows) + draw_table_page(pdf, "Mixed workloads", mixed_rows) + draw_table_page(pdf, "Large stress workloads", large_rows) + + +def main(): + output = Path(__file__).with_name("result_speed.pdf") + cases = build_cases() + results = [] + for index, case in enumerate(cases, start=1): + print(f"[{index}/{len(cases)}] {case.name}") + results.append(benchmark_case(case)) + + print("[profile] optimized implementation") + profile = profile_optimized() + + payload = { + "machine": machine_info(), + "iterations": ITERATIONS, + "warmups": WARMUPS, + "results": results, + "profile": profile, + } + output.with_suffix(".json").write_text(json.dumps(payload, indent=2), encoding="utf-8") + write_csv(results, output.with_suffix(".csv")) + write_report(results, profile, output) + print(output) + + +if __name__ == "__main__": + main() diff --git a/qpiai_quantum/simulator/density_matrix.py b/qpiai_quantum/simulator/density_matrix.py new file mode 100644 index 0000000..c705d1c --- /dev/null +++ b/qpiai_quantum/simulator/density_matrix.py @@ -0,0 +1,490 @@ +""" +Density Matrix State & Local Simulator +======================================= +Implements the DensityMatrixState representation and the DensityMatrixSimulator +backend for executing quantum circuits using density matrix formalism. +""" + +from __future__ import annotations + +import time +from typing import TYPE_CHECKING +import numpy as np + +from qpiai_quantum.icr.circuitoperation import OperationType +from qpiai_quantum.simulator.base_simulator import BaseSimulator +from qpiai_quantum.simulator.gates import ( + ALL_KNOWN_GATES, + DECOMPOSED_GATES, + H, + S, + SDG, + SWAP, + T, + TDG, + X, + Y, + Z, + controlled, + decompose, + gate_spec, + rx_matrix, + ry_matrix, + rz_matrix, +) +from qpiai_quantum.simulator.result import QasmSimulatorResult +from qpiai_quantum.simulator.state import QuantumState + +if TYPE_CHECKING: + from qpiai_quantum.circuit import Circuit + + +class DensityMatrixState(QuantumState): + """ + Density matrix representation of a quantum state (2^n x 2^n complex matrix). + + Supports pure and mixed quantum states, single- and multi-qubit unitary + evolutions, Kraus noise channels, partial traces, and measurement statistics. + """ + + def __init__(self, num_qubits: int, data: np.ndarray | None = None) -> None: + """ + Initialize a DensityMatrixState for n qubits. + + Args: + num_qubits: Number of qubits. + data: Optional 2D complex numpy array of shape (2^n, 2^n). + Defaults to the pure zero state |0...0><0...0|. + """ + if num_qubits <= 0: + raise ValueError(f"num_qubits must be positive, got {num_qubits}") + + self._num_qubits = num_qubits + dim = 1 << num_qubits + + if data is None: + self._data = np.zeros((dim, dim), dtype=complex) + self._data[0, 0] = 1.0 + 0.0j + else: + arr = np.asarray(data, dtype=complex) + if arr.ndim == 1: + # Convert statevector to density matrix + if arr.shape[0] != dim: + raise ValueError( + f"Statevector size {arr.shape[0]} does not match 2^{num_qubits} = {dim}" + ) + norm = np.linalg.norm(arr) + if norm > 0: + arr = arr / norm + self._data = np.outer(arr, np.conj(arr)) + elif arr.ndim == 2: + if arr.shape != (dim, dim): + raise ValueError( + f"Density matrix shape {arr.shape} does not match ({dim}, {dim})" + ) + self._data = arr.copy() + else: + raise ValueError( + f"Expected 1D statevector or 2D density matrix, got {arr.ndim}D array" + ) + + @property + def num_qubits(self) -> int: + """Number of qubits represented in the density matrix.""" + return self._num_qubits + + @property + def data(self) -> np.ndarray: + """Underlying 2D numpy array representing the density matrix.""" + return self._data + + @data.setter + def data(self, new_data: np.ndarray) -> None: + dim = 1 << self._num_qubits + arr = np.asarray(new_data, dtype=complex) + if arr.shape != (dim, dim): + raise ValueError(f"Shape must be ({dim}, {dim}), got {arr.shape}") + self._data = arr + + # ---------------------------------------------------------------------- + # Factory Constructors + # ---------------------------------------------------------------------- + + @classmethod + def pure_state(cls, statevector: np.ndarray | list[complex]) -> DensityMatrixState: + """ + Construct a DensityMatrixState from a pure statevector. + + Args: + statevector: 1D array of amplitudes of length 2^n. + """ + sv = np.asarray(statevector, dtype=complex) + dim = sv.shape[0] + num_qubits = int(np.log2(dim)) + if 1 << num_qubits != dim: + raise ValueError(f"Statevector length {dim} must be a power of 2.") + return cls(num_qubits=num_qubits, data=sv) + + @classmethod + def maximally_mixed(cls, num_qubits: int) -> DensityMatrixState: + """ + Construct a maximally mixed state ρ = I / 2^n. + + Args: + num_qubits: Number of qubits. + """ + dim = 1 << num_qubits + matrix = np.eye(dim, dtype=complex) / dim + return cls(num_qubits=num_qubits, data=matrix) + + # ---------------------------------------------------------------------- + # Single-Qubit Gates + # ---------------------------------------------------------------------- + + def apply_h(self, target: int) -> None: + """Apply Hadamard gate to target qubit.""" + self.apply_unitary(H, [target]) + + def apply_x(self, target: int) -> None: + """Apply Pauli-X gate to target qubit.""" + self.apply_unitary(X, [target]) + + def apply_y(self, target: int) -> None: + """Apply Pauli-Y gate to target qubit.""" + self.apply_unitary(Y, [target]) + + def apply_z(self, target: int) -> None: + """Apply Pauli-Z gate to target qubit.""" + self.apply_unitary(Z, [target]) + + def apply_s(self, target: int) -> None: + """Apply Phase (S) gate to target qubit.""" + self.apply_unitary(S, [target]) + + def apply_sdg(self, target: int) -> None: + """Apply S† gate to target qubit.""" + self.apply_unitary(SDG, [target]) + + def apply_t(self, target: int) -> None: + """Apply T gate to target qubit.""" + self.apply_unitary(T, [target]) + + def apply_tdg(self, target: int) -> None: + """Apply T† gate to target qubit.""" + self.apply_unitary(TDG, [target]) + + def apply_rx(self, target: int, theta: float) -> None: + """Apply RX rotation gate to target qubit.""" + self.apply_unitary(rx_matrix(theta), [target]) + + def apply_ry(self, target: int, theta: float) -> None: + """Apply RY rotation gate to target qubit.""" + self.apply_unitary(ry_matrix(theta), [target]) + + def apply_rz(self, target: int, theta: float) -> None: + """Apply RZ rotation gate to target qubit.""" + self.apply_unitary(rz_matrix(theta), [target]) + + # ---------------------------------------------------------------------- + # Multi-Qubit Gates + # ---------------------------------------------------------------------- + + def apply_cnot(self, control: int, target: int) -> None: + """Apply CNOT (CX) gate between control and target qubits.""" + self.apply_unitary(controlled(X, 1), [control, target]) + + def apply_cz(self, control: int, target: int) -> None: + """Apply Controlled-Z (CZ) gate between control and target qubits.""" + self.apply_unitary(controlled(Z, 1), [control, target]) + + def apply_swap(self, target1: int, target2: int) -> None: + """Apply SWAP gate between target1 and target2 qubits.""" + self.apply_unitary(SWAP, [target1, target2]) + + def apply_cswap(self, control: int, target1: int, target2: int) -> None: + """Apply Fredkin (CSWAP) gate.""" + self.apply_unitary(controlled(SWAP, 1), [control, target1, target2]) + + # ---------------------------------------------------------------------- + # Generic Unitary & Kraus Transformations + # ---------------------------------------------------------------------- + + def apply_unitary(self, matrix: np.ndarray, qubits: list[int]) -> None: + """ + Apply a unitary transformation U to specified qubits: ρ -> U ρ U†. + + Args: + matrix: 2^k x 2^k complex matrix where k = len(qubits). + qubits: List of qubit indices. + """ + n = self._num_qubits + k = len(qubits) + if matrix.shape != (2**k, 2**k): + raise ValueError( + f"Matrix shape {matrix.shape} does not match 2^{k} x 2^{k}" + ) + for q in qubits: + if not (0 <= q < n): + raise ValueError(f"Qubit index {q} out of bounds for n={n}") + + self._data = self._transform_density_matrix(self._data, n, qubits, matrix) + + def apply_kraus(self, kraus_ops: list[np.ndarray], qubits: list[int]) -> None: + """ + Apply a quantum channel represented by Kraus operators: ρ -> Σ_i K_i ρ K_i†. + + Args: + kraus_ops: List of Kraus matrices K_i of shape (2^k, 2^k). + qubits: List of target qubits. + """ + n = self._num_qubits + dim = 1 << n + new_rho = np.zeros((dim, dim), dtype=complex) + for k_op in kraus_ops: + new_rho += self._transform_density_matrix(self._data, n, qubits, k_op) + self._data = new_rho + + @staticmethod + def _transform_density_matrix( + rho: np.ndarray, n: int, qubits: list[int], M: np.ndarray + ) -> np.ndarray: + """ + Transform density matrix by operator M: ρ -> M ρ M† on target qubits. + """ + k = len(qubits) + # Reshape into 2n tensor: n ket axes followed by n bra axes + tensor = rho.reshape([2] * (2 * n)) + + # Little-endian convention: qubit q corresponds to axis (n - 1 - q) + row_axes = [n - 1 - q for q in qubits] + col_axes = [n + (n - 1 - q) for q in qubits] + + # 1. Apply M on row (ket) indices: (M @ flat_row) + tensor = np.moveaxis(tensor, row_axes, list(range(k))) + shape_row = tensor.shape + flat_row = tensor.reshape(2**k, -1) + flat_row = M @ flat_row + tensor = flat_row.reshape(shape_row) + tensor = np.moveaxis(tensor, list(range(k)), row_axes) + + # 2. Apply M* on col (bra) indices: equivalent to ρ' @ M† + tensor = np.moveaxis(tensor, col_axes, list(range(k))) + shape_col = tensor.shape + flat_col = tensor.reshape(2**k, -1) + flat_col = np.conj(M) @ flat_col + tensor = flat_col.reshape(shape_col) + tensor = np.moveaxis(tensor, list(range(k)), col_axes) + + dim = 1 << n + return tensor.reshape((dim, dim)) + + # ---------------------------------------------------------------------- + # Quantum Information Diagnostics & Operations + # ---------------------------------------------------------------------- + + def partial_trace(self, qubits_to_trace: list[int]) -> DensityMatrixState: + """ + Compute the reduced density matrix by tracing out the specified qubits. + + Args: + qubits_to_trace: Qubits to trace out. + + Returns: + A new DensityMatrixState for the remaining subsystem. + """ + n = self._num_qubits + for q in qubits_to_trace: + if not (0 <= q < n): + raise ValueError(f"Qubit index {q} out of bounds for n={n}") + + remaining_qubits = [q for q in range(n) if q not in qubits_to_trace] + if not remaining_qubits: + raise ValueError("Cannot trace out all qubits in partial_trace.") + + tensor = self._data.reshape([2] * (2 * n)) + + # For each qubit to trace, contract its row axis with its col axis + traced_tensor = tensor + # We trace out one by one from highest qubit index down to keep indices consistent + for q in sorted(qubits_to_trace, reverse=True): + current_n = traced_tensor.ndim // 2 + # Find axis of q in current subsystem + row_axis = current_n - 1 - q + col_axis = current_n + row_axis + traced_tensor = np.trace(traced_tensor, axis1=row_axis, axis2=col_axis) + + rem_n = len(remaining_qubits) + rem_dim = 1 << rem_n + reduced_matrix = traced_tensor.reshape((rem_dim, rem_dim)) + return DensityMatrixState(num_qubits=rem_n, data=reduced_matrix) + + def purity(self) -> float: + """Return purity γ = Tr(ρ²).""" + return float(np.real(np.trace(self._data @ self._data))) + + def von_neumann_entropy(self) -> float: + """Return von Neumann entropy S(ρ) = -Tr(ρ log₂ ρ).""" + eigenvalues = np.linalg.eigvalsh(self._data) + eigenvalues = eigenvalues[eigenvalues > 1e-15] + return float(-np.sum(eigenvalues * np.log2(eigenvalues))) + + def get_probabilities(self) -> np.ndarray: + """Return measurement probabilities for all 2^n computational basis states.""" + probs = np.real(np.diag(self._data)) + # Clamp negative rounding errors + probs = np.maximum(probs, 0.0) + total = probs.sum() + if total > 0: + probs = probs / total + return probs + + def reset(self) -> None: + """Reset the density matrix state back to |0...0><0...0|.""" + dim = 1 << self._num_qubits + self._data = np.zeros((dim, dim), dtype=complex) + self._data[0, 0] = 1.0 + 0.0j + + +class DensityMatrixSimulator(BaseSimulator): + """ + Local Density Matrix Quantum Simulator. + + Simulates quantum circuits using the density matrix formalism, supporting + both pure states and mixed states/noise channels. + """ + + @property + def name(self) -> str: + """Name of the simulator backend.""" + return "QpiAI-QDM-Local" + + def run( + self, + circuit: Circuit, + shots: int = 1024, + seed: int | None = None, + name: str | None = None, + initial_state: np.ndarray | DensityMatrixState | None = None, + ) -> QasmSimulatorResult: + """ + Execute a quantum circuit using density matrix simulation. + + Args: + circuit: The quantum circuit to execute. + shots: Number of measurement shots to sample. + seed: Optional RNG seed. + name: Optional name for the result object. + initial_state: Optional initial density matrix or statevector. + + Returns: + QasmSimulatorResult containing density matrix and measurement counts. + """ + start_time = time.perf_counter() + n_qubits = circuit.num_qubits + n_cbits = circuit.num_clbits + + if n_qubits <= 0: + raise ValueError( + f"Circuit must have at least 1 qubit, got {n_qubits} qubits" + ) + + # Initialize density matrix state + if initial_state is None: + state = DensityMatrixState(n_qubits) + elif isinstance(initial_state, DensityMatrixState): + if initial_state.num_qubits != n_qubits: + raise ValueError( + f"initial_state qubits {initial_state.num_qubits} != circuit qubits {n_qubits}" + ) + state = DensityMatrixState(n_qubits, data=initial_state.data) + else: + state = DensityMatrixState(n_qubits, data=initial_state) + + measure_map: dict[int, int] = {} + + def _apply_gate(gate_name: str, params: list[float], qubits: list[int]) -> None: + norm_name = gate_name.lower() + if norm_name in DECOMPOSED_GATES: + for sub_name, sub_p, sub_q in decompose(norm_name, qubits): + _apply_gate(sub_name, sub_p, sub_q) + return + + if norm_name in ALL_KNOWN_GATES: + _, matrix = gate_spec(norm_name, params, len(qubits)) + state.apply_unitary(matrix, qubits) + elif norm_name == "barrier": + pass + else: + raise ValueError(f"Unknown gate '{gate_name}'") + + def _apply_op(op) -> None: + if op.operation_type == OperationType.MEASURE: + if op.qubits and op.clbits: + for q, c in zip(op.qubits, op.clbits): + measure_map[q] = c + elif op.operation_type in ( + OperationType.N_QUBIT_NON_PARAMETRIC, + OperationType.N_QUBIT_PARAMETRIC, + OperationType.SWAP, + ): + _apply_gate(op.gate_name, op.params or [], op.qubits) + elif op.operation_type == OperationType.OPERATION: + if hasattr(op, "order") and op.order is not None: + for sub_op in op.order: + _apply_op(sub_op) + else: + _apply_gate(op.gate_name, op.params or [], op.qubits) + else: + raise ValueError(f"Unsupported operation type: {op.operation_type}") + + # Execute circuit evolution + for op in circuit.icr.evolve: + _apply_op(op) + + # Sample measurement counts if measurements were registered + probs = state.get_probabilities() + if n_cbits > 0 and measure_map: + counts = self._sample_counts( + probs, n_qubits, n_cbits, measure_map, shots, seed + ) + else: + counts = {} + + elapsed_time = time.perf_counter() - start_time + + return QasmSimulatorResult( + name=name or circuit.name, + counts=counts, + statevector=None, + density_matrix=state.data.tolist(), + shots=shots, + execution_time=elapsed_time, + method="density_matrix", + job_status="completed", + n_qubits=n_qubits, + n_cbits=n_cbits, + ) + + @staticmethod + def _sample_counts( + probs: np.ndarray, + n_qubits: int, + n_cbits: int, + measure_map: dict[int, int], + shots: int, + seed: int | None = None, + ) -> dict[str, int]: + """Sample measurement outcomes from probability distribution.""" + rng = np.random.default_rng(seed) + outcomes = rng.choice(len(probs), size=shots, p=probs) + + counts: dict[str, int] = {} + for outcome in outcomes: + cbits = ["0"] * n_cbits + for qubit, cbit in measure_map.items(): + bit = (outcome >> qubit) & 1 + cbits[cbit] = str(bit) + bitstring = "".join(reversed(cbits)) if n_cbits > 0 else "" + counts[bitstring] = counts.get(bitstring, 0) + 1 + + return counts diff --git a/qpiai_quantum/simulator/state.py b/qpiai_quantum/simulator/state.py new file mode 100644 index 0000000..7769c85 --- /dev/null +++ b/qpiai_quantum/simulator/state.py @@ -0,0 +1,133 @@ +""" +Quantum State Interface +======================= +Abstract base class defining the standard interface for quantum state +representations (StateVector, DensityMatrix, Stabilizer, etc.) in the +QpiAI Quantum SDK. +""" + +from abc import ABC, abstractmethod +import numpy as np + + +class QuantumState(ABC): + """Abstract interface for quantum state representations.""" + + @property + @abstractmethod + def num_qubits(self) -> int: + """Return the number of qubits in the quantum state.""" + pass + + @property + @abstractmethod + def data(self) -> np.ndarray: + """Return the underlying state array (e.g. 1D vector or 2D matrix).""" + pass + + # ---------------------------------------------------------------------- + # Single-Qubit Gates + # ---------------------------------------------------------------------- + + @abstractmethod + def apply_h(self, target: int) -> None: + """Apply Hadamard gate to target qubit.""" + pass + + @abstractmethod + def apply_x(self, target: int) -> None: + """Apply Pauli-X gate to target qubit.""" + pass + + @abstractmethod + def apply_y(self, target: int) -> None: + """Apply Pauli-Y gate to target qubit.""" + pass + + @abstractmethod + def apply_z(self, target: int) -> None: + """Apply Pauli-Z gate to target qubit.""" + pass + + @abstractmethod + def apply_s(self, target: int) -> None: + """Apply Phase (S) gate to target qubit.""" + pass + + @abstractmethod + def apply_sdg(self, target: int) -> None: + """Apply S† gate to target qubit.""" + pass + + @abstractmethod + def apply_t(self, target: int) -> None: + """Apply T (π/8) gate to target qubit.""" + pass + + @abstractmethod + def apply_tdg(self, target: int) -> None: + """Apply T† gate to target qubit.""" + pass + + @abstractmethod + def apply_rx(self, target: int, theta: float) -> None: + """Apply RX rotation gate by theta radians to target qubit.""" + pass + + @abstractmethod + def apply_ry(self, target: int, theta: float) -> None: + """Apply RY rotation gate by theta radians to target qubit.""" + pass + + @abstractmethod + def apply_rz(self, target: int, theta: float) -> None: + """Apply RZ rotation gate by theta radians to target qubit.""" + pass + + # ---------------------------------------------------------------------- + # Multi-Qubit Gates + # ---------------------------------------------------------------------- + + @abstractmethod + def apply_cnot(self, control: int, target: int) -> None: + """Apply CNOT (CX) gate between control and target qubits.""" + pass + + @abstractmethod + def apply_cz(self, control: int, target: int) -> None: + """Apply Controlled-Z (CZ) gate between control and target qubits.""" + pass + + @abstractmethod + def apply_swap(self, target1: int, target2: int) -> None: + """Apply SWAP gate between target1 and target2 qubits.""" + pass + + @abstractmethod + def apply_cswap(self, control: int, target1: int, target2: int) -> None: + """Apply Fredkin (CSWAP) gate.""" + pass + + # ---------------------------------------------------------------------- + # Generic Unitary & Channel Evolution + # ---------------------------------------------------------------------- + + @abstractmethod + def apply_unitary(self, matrix: np.ndarray, qubits: list[int]) -> None: + """Apply an arbitrary k-qubit unitary matrix to the specified qubits.""" + pass + + # ---------------------------------------------------------------------- + # Diagnostics & Measurement + # ---------------------------------------------------------------------- + + @abstractmethod + def get_probabilities(self) -> np.ndarray: + """Return measurement probabilities for all 2^n computational basis states.""" + pass + + @abstractmethod + def reset(self) -> None: + """Reset the quantum state back to |0...0>.""" + pass + diff --git a/tests/test_density_matrix_simulator.py b/tests/test_density_matrix_simulator.py new file mode 100644 index 0000000..9319df1 --- /dev/null +++ b/tests/test_density_matrix_simulator.py @@ -0,0 +1,183 @@ +""" +Unit tests for DensityMatrixState and DensityMatrixSimulator. +""" + +import math +import numpy as np + +from qpiai_quantum.circuit import Circuit +from qpiai_quantum.simulator.density_matrix import ( + DensityMatrixSimulator, + DensityMatrixState, +) + + +class TestDensityMatrixState: + """Tests for the DensityMatrixState representation.""" + + def test_initial_zero_state(self): + dm = DensityMatrixState(num_qubits=2) + assert dm.num_qubits == 2 + assert dm.data.shape == (4, 4) + expected = np.zeros((4, 4), dtype=complex) + expected[0, 0] = 1.0 + np.testing.assert_allclose(dm.data, expected) + assert math.isclose(dm.purity(), 1.0, abs_tol=1e-7) + assert math.isclose(dm.von_neumann_entropy(), 0.0, abs_tol=1e-7) + + def test_maximally_mixed_state(self): + dm = DensityMatrixState.maximally_mixed(num_qubits=2) + expected = 0.25 * np.eye(4, dtype=complex) + np.testing.assert_allclose(dm.data, expected) + assert math.isclose(dm.purity(), 0.25, abs_tol=1e-7) + assert math.isclose(dm.von_neumann_entropy(), 2.0, abs_tol=1e-7) + + def test_pure_state_constructor(self): + # |+> state + plus_sv = np.array([1 / np.sqrt(2), 1 / np.sqrt(2)], dtype=complex) + dm = DensityMatrixState.pure_state(plus_sv) + expected = 0.5 * np.ones((2, 2), dtype=complex) + np.testing.assert_allclose(dm.data, expected) + assert math.isclose(dm.purity(), 1.0, abs_tol=1e-7) + + def test_single_qubit_gates(self): + # Pauli X on |0> -> |1><1| + dm = DensityMatrixState(num_qubits=1) + dm.apply_x(0) + np.testing.assert_allclose(dm.data, [[0, 0], [0, 1]]) + + # Hadamard on |0> -> |+><+| + dm.reset() + dm.apply_h(0) + np.testing.assert_allclose(dm.data, [[0.5, 0.5], [0.5, 0.5]]) + + # Phase (S) on |+> + dm.apply_s(0) + np.testing.assert_allclose(dm.data, [[0.5, -0.5j], [0.5j, 0.5]]) + + # S† undoes S + dm.apply_sdg(0) + np.testing.assert_allclose(dm.data, [[0.5, 0.5], [0.5, 0.5]]) + + # Rotations + dm.reset() + dm.apply_rx(0, np.pi) # RX(pi) ~ -i X, so |0><0| -> |1><1| + np.testing.assert_allclose(dm.data, [[0, 0], [0, 1]], atol=1e-12) + + def test_two_qubit_bell_state_and_partial_trace(self): + dm = DensityMatrixState(num_qubits=2) + # Create Bell state (|00> + |11>) / sqrt(2) + dm.apply_h(0) + dm.apply_cnot(0, 1) + + # Expected: 0.5 * (|00><00| + |00><11| + |11><00| + |11><11|) + expected = np.zeros((4, 4), dtype=complex) + expected[0, 0] = 0.5 + expected[0, 3] = 0.5 + expected[3, 0] = 0.5 + expected[3, 3] = 0.5 + np.testing.assert_allclose(dm.data, expected, atol=1e-12) + assert math.isclose(dm.purity(), 1.0, abs_tol=1e-7) + + # Tracing out qubit 1 yields maximally mixed state on qubit 0 + reduced_0 = dm.partial_trace([1]) + assert reduced_0.num_qubits == 1 + np.testing.assert_allclose(reduced_0.data, 0.5 * np.eye(2), atol=1e-12) + assert math.isclose(reduced_0.purity(), 0.5, abs_tol=1e-7) + assert math.isclose(reduced_0.von_neumann_entropy(), 1.0, abs_tol=1e-7) + + def test_swap_and_cz_gates(self): + # Prepare |10> = X on qubit 1 + dm = DensityMatrixState(num_qubits=2) + dm.apply_x(1) # state is |10> (basis index 2 in little-endian) + assert math.isclose(dm.get_probabilities()[2], 1.0, abs_tol=1e-7) + + # SWAP(0, 1) -> |01> (basis index 1 in little-endian) + dm.apply_swap(0, 1) + assert math.isclose(dm.get_probabilities()[1], 1.0, abs_tol=1e-7) + + # Prepare |11>, apply CZ -> |11> with phase -1 (in DM, (-1)*(-1) = +1) + dm.apply_x(1) # now |11> + dm.apply_cz(0, 1) + assert math.isclose(dm.get_probabilities()[3], 1.0, abs_tol=1e-7) + + def test_kraus_channel(self): + # Amplitude damping on |1><1| + dm = DensityMatrixState(num_qubits=1) + dm.apply_x(0) # |1><1| + + gamma = 0.3 + k0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]], dtype=complex) + k1 = np.array([[0, np.sqrt(gamma)], [0, 0]], dtype=complex) + + dm.apply_kraus([k0, k1], [0]) + np.testing.assert_allclose( + dm.data, [[gamma, 0], [0, 1 - gamma]], atol=1e-12 + ) + assert math.isclose(np.trace(dm.data).real, 1.0, abs_tol=1e-7) + + +class TestDensityMatrixSimulator: + """Tests for the DensityMatrixSimulator backend.""" + + def test_simulator_bell_state_run(self): + qc = Circuit(2, 2) + qc.h(0) + qc.cx(0, 1) + qc.measure(0, 0) + qc.measure(1, 1) + + sim = DensityMatrixSimulator() + result = sim.run(qc, shots=500, seed=42) + + assert result.method == "density_matrix" + assert result.job_status == "completed" + assert result.density_matrix is not None + + dm_arr = np.array(result.density_matrix) + assert dm_arr.shape == (4, 4) + + # Measurement counts should only contain '00' and '11' + assert set(result.counts.keys()).issubset({"00", "11"}) + assert sum(result.counts.values()) == 500 + + def test_simulator_ghz_3qubit(self): + qc = Circuit(3, 3) + qc.h(0) + qc.cx(0, 1) + qc.cx(1, 2) + qc.measure_all() + + sim = DensityMatrixSimulator() + result = sim.run(qc, shots=1000, seed=123) + + assert set(result.counts.keys()).issubset({"000", "111"}) + assert result.counts.get("000", 0) > 400 + assert result.counts.get("111", 0) > 400 + + dm_arr = np.array(result.density_matrix) + assert dm_arr.shape == (8, 8) + assert math.isclose(dm_arr[0, 0].real, 0.5, abs_tol=1e-7) + assert math.isclose(dm_arr[7, 7].real, 0.5, abs_tol=1e-7) + assert math.isclose(dm_arr[0, 7].real, 0.5, abs_tol=1e-7) + assert math.isclose(dm_arr[7, 0].real, 0.5, abs_tol=1e-7) + + def test_simulator_result_interoperability_with_quantum_info_density_matrix(self): + from qpiai_quantum.quantum_info.density_matrix import DensityMatrix + + qc = Circuit(2) + qc.h(0) + qc.cx(0, 1) + + sim = DensityMatrixSimulator() + result = sim.run(qc) + + # Interoperability: loading nested list into quantum_info.DensityMatrix + dm = DensityMatrix(result.density_matrix) + assert dm.num_qubits == 2 + assert dm.is_pure() is True + assert math.isclose(dm.purity(), 1.0, abs_tol=1e-7) + probs = dm.probabilities_dict() + assert math.isclose(probs.get("00", 0.0), 0.5, abs_tol=1e-7) + assert math.isclose(probs.get("11", 0.0), 0.5, abs_tol=1e-7) + diff --git a/tests/test_statevector_single_qubit_optimization.py b/tests/test_statevector_single_qubit_optimization.py new file mode 100644 index 0000000..56ad44c --- /dev/null +++ b/tests/test_statevector_single_qubit_optimization.py @@ -0,0 +1,264 @@ +import time +import sys +import types +import unittest +from pathlib import Path + +import numpy as np + +# Avoid executing qpiai_quantum.__init__, which imports optional algorithm +# dependencies that are unrelated to this simulator equivalence test. +if "qpiai_quantum" not in sys.modules: + qpiai_quantum = types.ModuleType("qpiai_quantum") + qpiai_quantum.__path__ = [str(Path(__file__).resolve().parents[1] / "qpiai_quantum")] + sys.modules["qpiai_quantum"] = qpiai_quantum + +from qpiai_quantum.circuit.circuit import Circuit +from qpiai_quantum.icr.circuitoperation import OperationType +from qpiai_quantum.simulator.gates import DECOMPOSED_GATES, decompose, gate_spec +from qpiai_quantum.simulator.result import QasmSimulatorResult +from qpiai_quantum.simulator.statevector import StatevectorSimulator + + +OPTIMIZED_GATES = ("h", "x", "y", "z", "s", "rx", "ry", "rz") +PARAMETRIC_GATES = {"rx", "ry", "rz"} + + +class OriginalStatevectorSimulator(StatevectorSimulator): + def run( + self, + circuit: Circuit, + shots: int = 1024, + seed: int | None = None, + name: str | None = None, + ) -> QasmSimulatorResult: + start_time = time.perf_counter() + n_qubits = circuit.num_qubits + n_cbits = circuit.num_clbits + + if n_qubits == 0: + raise ValueError("Cannot simulate a circuit with 0 qubits.") + + dim = 2**n_qubits + state = np.zeros(dim, dtype=complex) + state[0] = 1.0 + measure_map: dict[int, int] = {} + + def _apply_gate(gate_name: str, params: list[float], qubits: list[int]) -> None: + nonlocal state + gate_name_lower = gate_name.lower() + + if gate_name_lower in DECOMPOSED_GATES: + for sub_name, sub_params, sub_qubits in decompose( + gate_name_lower, qubits + ): + _apply_gate(sub_name, sub_params, sub_qubits) + return + + _, U = gate_spec(gate_name_lower, params, num_qubits=len(qubits)) + state = self._apply_unitary(state, n_qubits, qubits, U) + + def _apply_op(op) -> None: + nonlocal state + if op.operation_type == OperationType.BARRIER: + return + + if op.operation_type == OperationType.MEASURE: + if op.qubits and op.clbits: + for q, c in zip(op.qubits, op.clbits): + measure_map[q] = c + return + + if op.operation_type in ( + OperationType.N_QUBIT_NON_PARAMETRIC, + OperationType.N_QUBIT_PARAMETRIC, + OperationType.SWAP, + ): + _apply_gate(op.gate_name, op.params or [], op.qubits) + return + + if op.operation_type == OperationType.OPERATION: + if hasattr(op, "order") and op.order is not None: + for sub_op in op.order: + _apply_op(sub_op) + else: + _apply_gate(op.gate_name, op.params or [], op.qubits) + return + + if op.gate_name.lower() == "reset" and op.qubits: + for q in op.qubits: + state = self._apply_reset(state, n_qubits, q) + else: + raise ValueError(f"Unsupported operation type: {op.operation_type}") + + for op in circuit.icr.evolve: + _apply_op(op) + + counts = ( + self._sample_counts(state, n_qubits, n_cbits, measure_map, shots, seed) + if n_cbits > 0 and measure_map + else {} + ) + + return QasmSimulatorResult( + name=name or circuit.name, + counts=counts, + statevector=state.tolist(), + shots=shots, + executionTime=time.perf_counter() - start_time, + method="statevector-original", + job_status="completed", + n_qubits=n_qubits, + n_cbits=n_cbits, + ) + + +def apply_spec(circuit: Circuit, spec: tuple[str, tuple[int, ...], tuple[float, ...]]): + gate, qubits, params = spec + if gate in PARAMETRIC_GATES: + getattr(circuit, gate)(qubits[0], params[0]) + elif gate == "cx": + circuit.cx(qubits[0], qubits[1]) + else: + getattr(circuit, gate)(qubits[0]) + + +def build_circuit(n_qubits: int, specs) -> Circuit: + circuit = Circuit(n_qubits) + for spec in specs: + apply_spec(circuit, spec) + return circuit + + +def statevector(simulator, circuit: Circuit) -> np.ndarray: + return np.asarray(simulator.run(circuit, shots=1).statevector, dtype=complex) + + +def describe_specs(specs) -> str: + lines = [] + for index, (gate, qubits, params) in enumerate(specs, start=1): + params_text = f", params={list(params)}" if params else "" + lines.append(f"{index}. {gate.upper()} qubits={list(qubits)}{params_text}") + return "\n".join(lines) + + +def first_divergence(n_qubits: int, specs) -> str: + original = OriginalStatevectorSimulator() + optimized = StatevectorSimulator() + for index in range(1, len(specs) + 1): + prefix = specs[:index] + original_state = statevector(original, build_circuit(n_qubits, prefix)) + optimized_state = statevector(optimized, build_circuit(n_qubits, prefix)) + if not np.allclose(original_state, optimized_state): + gate, qubits, params = specs[index - 1] + return ( + f"first divergence at gate {index}: " + f"{gate.upper()} qubits={list(qubits)} params={list(params)}" + ) + return "no prefix divergence found" + + +def assert_equivalent(n_qubits: int, specs, expected: np.ndarray | None = None): + circuit = build_circuit(n_qubits, specs) + original_state = statevector(OriginalStatevectorSimulator(), circuit) + optimized_state = statevector(StatevectorSimulator(), circuit) + + if expected is not None: + assert np.allclose(original_state, expected), ( + "Original simulator did not match expected state.\n" + f"Circuit:\n{describe_specs(specs)}\n" + f"Original statevector:\n{original_state}\n" + f"Expected statevector:\n{expected}" + ) + assert np.allclose(optimized_state, expected), ( + "Optimized simulator did not match expected state.\n" + f"Circuit:\n{describe_specs(specs)}\n" + f"Optimized statevector:\n{optimized_state}\n" + f"Expected statevector:\n{expected}" + ) + + assert np.allclose(original_state, optimized_state), ( + "Original and optimized statevectors differ.\n" + f"Circuit:\n{describe_specs(specs)}\n" + f"Original statevector:\n{original_state}\n" + f"Optimized statevector:\n{optimized_state}\n" + f"{first_divergence(n_qubits, specs)}" + ) + + +class TestStatevectorSingleQubitOptimization(unittest.TestCase): + def test_individual_non_parametric_gates(self): + for gate in ("h", "x", "y", "z", "s"): + for n_qubits in (1, 2, 5, 10): + with self.subTest(gate=gate, n_qubits=n_qubits): + target = min(n_qubits - 1, max(0, n_qubits // 2)) + assert_equivalent(n_qubits, [(gate, (target,), ())]) + + def test_individual_rotation_gates_multiple_angles(self): + angles = np.random.default_rng(20260710).uniform( + -4 * np.pi, 4 * np.pi, size=12 + ) + for gate in ("rx", "ry", "rz"): + for n_qubits in (1, 2, 5, 10): + for theta in angles: + with self.subTest(gate=gate, n_qubits=n_qubits, theta=theta): + target = min(n_qubits - 1, max(0, n_qubits // 2)) + assert_equivalent(n_qubits, [(gate, (target,), (float(theta),))]) + + def test_multiple_sequential_single_qubit_gates(self): + specs = [ + ("h", (0,), ()), + ("rx", (0,), (0.125,)), + ("y", (0,), ()), + ("ry", (0,), (-1.75,)), + ("s", (0,), ()), + ("z", (0,), ()), + ("rz", (0,), (2.25,)), + ("x", (0,), ()), + ] + assert_equivalent(1, specs) + + def test_single_qubit_gates_on_different_targets(self): + for n_qubits in (2, 5, 10): + with self.subTest(n_qubits=n_qubits): + specs = [ + ("h", (0,), ()), + ("x", (n_qubits - 1,), ()), + ("rx", (1,), (0.37,)), + ("ry", (n_qubits // 2,), (-0.91,)), + ("rz", (n_qubits - 1,), (1.23,)), + ("s", (0,), ()), + ("y", (n_qubits // 2,), ()), + ("z", (1,), ()), + ] + assert_equivalent(n_qubits, specs) + + def test_random_circuits_containing_only_optimized_gates(self): + for n_qubits, depth in ((1, 40), (2, 60), (5, 100), (10, 120)): + with self.subTest(n_qubits=n_qubits, depth=depth): + rng = np.random.default_rng(12345 + n_qubits) + specs = [] + for _ in range(depth): + gate = str(rng.choice(OPTIMIZED_GATES)) + target = int(rng.integers(0, n_qubits)) + params = ( + (float(rng.uniform(-2 * np.pi, 2 * np.pi)),) + if gate in PARAMETRIC_GATES + else () + ) + specs.append((gate, (target,), params)) + assert_equivalent(n_qubits, specs) + + def test_bell_state_matches_original_and_expected_state(self): + expected = np.zeros(4, dtype=complex) + expected[0] = 1 / np.sqrt(2) + expected[3] = 1 / np.sqrt(2) + assert_equivalent(2, [("h", (0,), ()), ("cx", (0, 1), ())], expected) + + def test_bell_state_variation_matches_original(self): + specs = [("x", (1,), ()), ("h", (0,), ()), ("cx", (0, 1), ())] + assert_equivalent(2, specs) + + +if __name__ == "__main__": + unittest.main()