microsoft/qdk
Publicmirrored from https://github.com/microsoft/qdkAvailable
source/pip/tests/test_cpu_simulator.py
428lines · modecode
| 1 | # Copyright (c) Microsoft Corporation. |
| 2 | # Licensed under the MIT License. |
| 3 | |
| 4 | from collections import Counter |
| 5 | from pathlib import Path |
| 6 | from typing import Sequence, cast |
| 7 | import math |
| 8 | import random |
| 9 | |
| 10 | import pytest |
| 11 | |
| 12 | from qsharp._native import Result |
| 13 | |
| 14 | import qsharp |
| 15 | from qsharp import TargetProfile |
| 16 | from qsharp import openqasm |
| 17 | |
| 18 | from qsharp._simulation import run_qir_cpu, NoiseConfig |
| 19 | |
| 20 | current_file_path = Path(__file__) |
| 21 | # Get the directory of the current file |
| 22 | current_dir = current_file_path.parent |
| 23 | |
| 24 | |
| 25 | def read_file(file_name: str) -> str: |
| 26 | return Path(file_name).read_text(encoding="utf-8") |
| 27 | |
| 28 | |
| 29 | def read_file_relative(file_name: str) -> str: |
| 30 | return Path(current_dir / file_name).read_text(encoding="utf-8") |
| 31 | |
| 32 | |
| 33 | def result_array_to_string(results: Sequence[Result]) -> str: |
| 34 | chars = [] |
| 35 | for value in results: |
| 36 | if value == Result.Zero: |
| 37 | chars.append("0") |
| 38 | elif value == Result.One: |
| 39 | chars.append("1") |
| 40 | else: |
| 41 | chars.append("-") |
| 42 | return "".join(chars) |
| 43 | |
| 44 | |
| 45 | def test_cpu_seeding_no_noise(): |
| 46 | qsharp.init(target_profile=TargetProfile.Base) |
| 47 | qsharp.eval( |
| 48 | """ |
| 49 | operation BellTest() : Result[] { |
| 50 | use qs = Qubit[2]; |
| 51 | H(qs[0]); |
| 52 | CNOT(qs[0], qs[1]); |
| 53 | MResetEachZ(qs) |
| 54 | } |
| 55 | """ |
| 56 | ) |
| 57 | |
| 58 | qir = str(qsharp.compile("BellTest()")) |
| 59 | |
| 60 | results = [run_qir_cpu(qir, 1, None, seed)[0] for seed in range(100)] |
| 61 | print(results) |
| 62 | |
| 63 | # Results will be an array of 100 lists [Result, Result] |
| 64 | # Each result should be [Zero, Zero] or [One, One] |
| 65 | # As evident from a manual experiment running with the seeds of 0..99 |
| 66 | # gives 41:59 results split. Experiment should be repeatable for fixed seeds. |
| 67 | |
| 68 | # Verify we have 6 of each result |
| 69 | count_00 = sum(1 for r in results if r == [Result.Zero, Result.Zero]) |
| 70 | count_11 = sum(1 for r in results if r == [Result.One, Result.One]) |
| 71 | assert count_00 == 41 |
| 72 | assert count_11 == 100 - 41 |
| 73 | # TODO: count_00 is always suspiciously lower than count_11 for MANY ranges of seeds. |
| 74 | # Investigate if there's some bias in the simulator. Technically this isn't indication of a fault: |
| 75 | # we need roughly equal counts for shots, not for seeds. |
| 76 | |
| 77 | |
| 78 | def test_cpu_no_noise(): |
| 79 | """Simple test that CPU simulator works without noise.""" |
| 80 | qsharp.init(target_profile=TargetProfile.Base) |
| 81 | qsharp.eval(read_file_relative("CliffordIsing.qs")) |
| 82 | |
| 83 | input = qsharp.compile( |
| 84 | "IsingModel2DEvolution(4, 4, PI() / 2.0, PI() / 2.0, 10.0, 10)" |
| 85 | ) |
| 86 | |
| 87 | output = run_qir_cpu(str(input)) |
| 88 | print(output) |
| 89 | # Expecting deterministic output, no randomization seed needed. |
| 90 | assert output == [[Result.Zero] * 16], "Expected result of 0s with pi/2 angles." |
| 91 | |
| 92 | |
| 93 | def test_cpu_bitflip_noise(): |
| 94 | """Bitflip noise for CPU simulator.""" |
| 95 | qsharp.init(target_profile=TargetProfile.Base) |
| 96 | qsharp.eval(read_file_relative("CliffordIsing.qs")) |
| 97 | |
| 98 | input = qsharp.compile( |
| 99 | "IsingModel2DEvolution(4, 4, PI() / 2.0, PI() / 2.0, 10.0, 10)" |
| 100 | ) |
| 101 | |
| 102 | p_noise = 0.005 |
| 103 | noise = NoiseConfig() |
| 104 | noise.rx.set_bitflip(p_noise) |
| 105 | noise.rzz.set_pauli_noise("XX", p_noise) |
| 106 | noise.mresetz.set_bitflip(p_noise) |
| 107 | |
| 108 | output = run_qir_cpu(str(input), shots=3, noise=noise, seed=17) |
| 109 | result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] |
| 110 | print(result) |
| 111 | # Reasonable results obtained from manual run |
| 112 | assert result == ["1000010001000001", "0000000000000000", "0001000001100000"] |
| 113 | |
| 114 | |
| 115 | def test_cpu_mixed_noise(): |
| 116 | qsharp.init(target_profile=TargetProfile.Base) |
| 117 | qsharp.eval(read_file_relative("CliffordIsing.qs")) |
| 118 | |
| 119 | input = qsharp.compile( |
| 120 | "IsingModel2DEvolution(4, 4, PI() / 2.0, PI() / 2.0, 4.0, 4)" |
| 121 | ) |
| 122 | |
| 123 | noise = NoiseConfig() |
| 124 | noise.rz.set_bitflip(0.008) |
| 125 | noise.rz.loss = 0.005 |
| 126 | noise.rzz.set_depolarizing(0.008) |
| 127 | noise.rzz.loss = 0.005 |
| 128 | |
| 129 | output = run_qir_cpu(str(input), shots=3, noise=noise, seed=53) |
| 130 | result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] |
| 131 | print(result) |
| 132 | # Reasonable results obtained from manual run |
| 133 | assert result == ["000000000--00000", "0010000000000000", "0000000000100000"] |
| 134 | |
| 135 | |
| 136 | def test_cpu_isolated_loss(): |
| 137 | qsharp.init(target_profile=TargetProfile.Base) |
| 138 | program = """ |
| 139 | import Std.Math.PI; |
| 140 | operation Main() : Result[] { |
| 141 | use qs = Qubit[3]; |
| 142 | X(qs[0]); |
| 143 | X(qs[1]); |
| 144 | CNOT(qs[0], qs[1]); |
| 145 | // When loss is configured for X gate, qubit 2 should be unaffected. |
| 146 | Rx(PI() / 2.0, qs[2]); |
| 147 | Rx(PI() / 2.0, qs[2]); |
| 148 | MeasureEachZ(qs) |
| 149 | } |
| 150 | """ |
| 151 | qsharp.eval(program) |
| 152 | |
| 153 | input = qsharp.compile( |
| 154 | "Main()" |
| 155 | ) |
| 156 | |
| 157 | noise = NoiseConfig() |
| 158 | noise.x.loss = 0.1 |
| 159 | |
| 160 | output = run_qir_cpu(str(input), shots=1000, noise=noise) |
| 161 | result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] |
| 162 | histogram = Counter(result) |
| 163 | total = sum(histogram.values()) |
| 164 | allowed_percent = { |
| 165 | "101": 0.81, |
| 166 | "1-1": 0.09, |
| 167 | "-11": 0.09, |
| 168 | "--1": 0.01, |
| 169 | } |
| 170 | tolerance = 0.2 * total |
| 171 | for bitstring, actual_count in histogram.items(): |
| 172 | assert bitstring in allowed_percent, f"Unexpected measurement string: '{bitstring}'." |
| 173 | expected_count = allowed_percent[bitstring] * total |
| 174 | assert abs(actual_count - expected_count) <= tolerance, ( |
| 175 | f"Count for {bitstring} outside 20% tolerance. " |
| 176 | f"Actual={actual_count}, Expected≈{expected_count:.0f}, Shots={total}." |
| 177 | ) |
| 178 | # We don't check for missing strings, as low-probability strings may not appear in finite shots. |
| 179 | |
| 180 | |
| 181 | def test_cpu_isolated_loss_and_noise(): |
| 182 | qsharp.init(target_profile=TargetProfile.Base) |
| 183 | program = """ |
| 184 | import Std.Math.PI; |
| 185 | operation Main() : Result[] { |
| 186 | use qs = Qubit[5]; |
| 187 | for _ in 1..100 { |
| 188 | X(qs[0]); |
| 189 | X(qs[1]); |
| 190 | CNOT(qs[0], qs[1]); |
| 191 | } |
| 192 | Rx(PI() / 2.0, qs[4]); |
| 193 | Rx(PI() / 2.0, qs[4]); |
| 194 | MeasureEachZ(qs) |
| 195 | } |
| 196 | """ |
| 197 | qsharp.eval(program) |
| 198 | |
| 199 | input = qsharp.compile( |
| 200 | "Main()" |
| 201 | ) |
| 202 | |
| 203 | noise = NoiseConfig() |
| 204 | noise.x.set_bitflip(0.001) |
| 205 | noise.x.loss = 0.001 |
| 206 | |
| 207 | output = run_qir_cpu(str(input), shots=1000, noise=noise) |
| 208 | result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] |
| 209 | histogram = Counter(result) |
| 210 | total = sum(histogram.values()) |
| 211 | assert total > 0, "No measurement results recorded." |
| 212 | for bitstring in histogram: |
| 213 | assert bitstring.endswith("001"), f"Unexpected suffix in '{bitstring}'." |
| 214 | probability_00001 = histogram.get("00001", 0) / total |
| 215 | assert 0.5 < probability_00001 < 0.8, ( |
| 216 | f"Probability of 00001 outside expected range. " |
| 217 | f"Actual={probability_00001:.2%}, Shots={total}." |
| 218 | ) |
| 219 | |
| 220 | |
| 221 | def build_x_chain_qir(n_instances: int, n_x: int) -> str: |
| 222 | # Construct multiple instances of x gate chains |
| 223 | prefix = f""" |
| 224 | OPENQASM 3.0; |
| 225 | include "stdgates.inc"; |
| 226 | bit[{n_instances}] c; |
| 227 | qubit[{n_instances}] q; |
| 228 | """ |
| 229 | |
| 230 | infix = """ |
| 231 | x q; |
| 232 | """ |
| 233 | |
| 234 | suffix = """ |
| 235 | c = measure q; |
| 236 | """ |
| 237 | |
| 238 | src_parallel = prefix + infix * n_x + suffix |
| 239 | |
| 240 | # Compile resulting program |
| 241 | qsharp.init(target_profile=TargetProfile.Base) |
| 242 | qir_parallel = openqasm.compile(src_parallel) |
| 243 | return str(qir_parallel) |
| 244 | |
| 245 | |
| 246 | @pytest.mark.parametrize( |
| 247 | "p_noise, n_x, n_instances, n_shots, max_percent", |
| 248 | [ |
| 249 | (0.001, 200, 6, 4096, 2.0), |
| 250 | (0.01, 200, 6, 4096, 2.0), |
| 251 | (0.001, 50, 12, 1024, 4.0), # 50 shots is low, so higher error tolerated |
| 252 | ], |
| 253 | ) |
| 254 | def test_cpu_x_chain( |
| 255 | p_noise: float, n_x: int, n_instances: int, n_shots: int, max_percent: float |
| 256 | ): |
| 257 | """ |
| 258 | Simulate multi-instance X-chain with bitflip noise many times |
| 259 | Compare result frequencies with analytically computed probabilities |
| 260 | """ |
| 261 | # Use the CPU simulator with noise |
| 262 | noise = NoiseConfig() |
| 263 | noise.x.set_bitflip(p_noise) |
| 264 | |
| 265 | qir = build_x_chain_qir(n_instances, n_x) |
| 266 | output = run_qir_cpu(qir, shots=n_shots, noise=noise, seed=18) |
| 267 | histogram = [0 for _ in range(n_instances + 1)] |
| 268 | for shot in output: |
| 269 | shot_results = cast(Sequence[Result], shot) |
| 270 | count_1 = shot_results.count(Result.One) |
| 271 | histogram[count_1] += 1 |
| 272 | |
| 273 | # Probability of obtaining 0 and 1 at the end of the X chain. |
| 274 | p_0 = ((2.0 * p_noise - 1.0) ** n_x + 1.0) / 2.0 |
| 275 | p_1 = 1.0 - p_0 |
| 276 | |
| 277 | # Number of results with k ones that should be there. |
| 278 | p_N = [ |
| 279 | p_0 ** ((n_instances - k)) * (p_1**k) * math.comb(n_instances, k) * n_shots |
| 280 | for k in range(n_instances + 1) |
| 281 | ] |
| 282 | |
| 283 | # Error % for deviation from analytical value |
| 284 | error_percent = [abs(a - b) * 100.0 / n_shots for (a, b) in zip(histogram, p_N)] |
| 285 | print(", ".join(f"{a} (Δ≈{b:.1f}%)" for (a, b) in zip(histogram, error_percent))) |
| 286 | |
| 287 | # We tolerate configured percentage error. |
| 288 | assert all( |
| 289 | err < max_percent for err in error_percent |
| 290 | ), f"Error percent too high: {error_percent}" |
| 291 | |
| 292 | |
| 293 | def generate_op_sequence( |
| 294 | n_qubits: int, n_ops: int, n_rand: int |
| 295 | ) -> list[tuple[int, int]]: |
| 296 | """Return operation tuples and randomly swap neighboring pairs n_rand times.""" |
| 297 | if n_qubits < 0 or n_ops < 0 or n_rand < 0: |
| 298 | raise ValueError("Tuple bounds must be non-negative") |
| 299 | |
| 300 | ops = [(q, op) for op in range(n_ops) for q in range(n_qubits)] |
| 301 | |
| 302 | if len(ops) < 2 or n_rand == 0: |
| 303 | return ops |
| 304 | |
| 305 | max_index = len(ops) - 1 |
| 306 | for _ in range(n_rand): |
| 307 | idx = random.randrange(max_index) |
| 308 | left, right = ops[idx], ops[idx + 1] |
| 309 | if left[0] != right[0]: |
| 310 | ops[idx], ops[idx + 1] = right, left |
| 311 | |
| 312 | return ops |
| 313 | |
| 314 | |
| 315 | @pytest.mark.parametrize("noisy_gate, noise_number", [(0, 2), (1, 1), (2, 2), (3, 2)]) |
| 316 | def test_cpu_permuted_rotations(noisy_gate: int, noise_number: int): |
| 317 | qsharp.init(target_profile=TargetProfile.Base) |
| 318 | |
| 319 | n_shots = 2000 |
| 320 | n_qubits = 11 |
| 321 | seed = 2026 |
| 322 | p_loss = 0.1 |
| 323 | tolerance_percent = 2.0 |
| 324 | assert n_qubits >= 2, "Need at least two qubits" |
| 325 | |
| 326 | random.seed(seed) |
| 327 | i1, i2 = random.sample(range(n_qubits), 2) |
| 328 | prefix = f""" |
| 329 | operation tiny_coeffs() : Result[] {{ |
| 330 | use q = Qubit[{n_qubits}]; |
| 331 | let i1 = {i1}; |
| 332 | let i2 = {i2}; |
| 333 | """ |
| 334 | |
| 335 | # The following sequence of rotations is equivalent to identity: |
| 336 | # 0. H <- could be any rotation |
| 337 | # 1. Rx(1.123456789) |
| 338 | # 2. Ry(1.212121212) |
| 339 | # 3. Rz(1.14856940153986) |
| 340 | # 4. Ry(-1.41836046203971) |
| 341 | # 5. Rz(-0.325946593598928) |
| 342 | # 6. H <- adjoint to step 0 |
| 343 | # We will perform these rotations on every qubit, but randomly intermix sequences for different qubits. |
| 344 | # This should still result in identity on all qubits as gates on different qubits commute. |
| 345 | # noise_number = how many times noisy gate appears in sequence. |
| 346 | |
| 347 | n_ops = 7 |
| 348 | ops = generate_op_sequence(n_qubits, n_ops, n_qubits * n_ops * 100) |
| 349 | infix = "" |
| 350 | for qubit, op in ops: |
| 351 | match op: |
| 352 | case 0 | 6: |
| 353 | infix += f" H(q[{qubit}]);\n" |
| 354 | case 1: |
| 355 | infix += f" Rx(1.123456789, q[{qubit}]);\n" |
| 356 | case 2: |
| 357 | infix += f" Ry(1.212121212, q[{qubit}]);\n" |
| 358 | case 3: |
| 359 | infix += f" Rz(1.14856940153986, q[{qubit}]);\n" |
| 360 | case 4: |
| 361 | infix += f" Ry(-1.41836046203971, q[{qubit}]);\n" |
| 362 | case 5: |
| 363 | infix += f" Rz(-0.325946593598928, q[{qubit}]);\n" |
| 364 | |
| 365 | suffix = """ |
| 366 | let m1 = M(q[i1]); |
| 367 | let m2 = M(q[i2]); |
| 368 | ResetAll(q); |
| 369 | return [m1, m2]; |
| 370 | } |
| 371 | """ |
| 372 | |
| 373 | program = prefix + infix + suffix |
| 374 | qsharp.eval(program) |
| 375 | input = qsharp.compile("tiny_coeffs()") |
| 376 | |
| 377 | noise = NoiseConfig() |
| 378 | p_combined_loss = 1.0 - ((1.0 - p_loss) ** noise_number) |
| 379 | match noisy_gate: |
| 380 | case 0: |
| 381 | noise.h.loss = p_loss |
| 382 | case 1: |
| 383 | noise.rx.loss = p_loss |
| 384 | case 2: |
| 385 | noise.ry.loss = p_loss |
| 386 | case 3: |
| 387 | noise.rz.loss = p_loss |
| 388 | case _: |
| 389 | raise ValueError("Invalid noisy_gate value") |
| 390 | |
| 391 | output = run_qir_cpu(str(input), shots=n_shots, noise=noise, seed=seed) |
| 392 | result_strings = [ |
| 393 | result_array_to_string(cast(Sequence[Result], shot)) for shot in output |
| 394 | ] |
| 395 | assert ( |
| 396 | len(result_strings) == n_shots |
| 397 | ), f"Shot count mismatch. Actual={len(result_strings)}, Expected={n_shots}" |
| 398 | |
| 399 | p_minus = p_combined_loss |
| 400 | p_0 = 1.0 - p_minus |
| 401 | allowed = [ |
| 402 | ("00", n_shots * p_0 * p_0), |
| 403 | ("0-", n_shots * p_0 * p_minus), |
| 404 | ("-0", n_shots * p_minus * p_0), |
| 405 | ("--", n_shots * p_minus * p_minus), |
| 406 | ] |
| 407 | |
| 408 | counts = {pattern: 0 for pattern, _ in allowed} |
| 409 | for entry in result_strings: |
| 410 | assert ( |
| 411 | entry in counts |
| 412 | ), f"Unexpected measurement string: '{entry}'. Program={program}." |
| 413 | counts[entry] += 1 |
| 414 | |
| 415 | tolerance = tolerance_percent / 100.0 * n_shots |
| 416 | print( |
| 417 | f"Permuted rotations test: n_qubits={n_qubits}, n_shots={n_shots}, seed={seed}, noise#{noise_number}, Δ<={tolerance:.0f} i1={i1}, i2={i2}" |
| 418 | ) |
| 419 | summary_msg = ", ".join( |
| 420 | f"'{pattern}': {counts[pattern]} (Δ={abs(counts[pattern] - expected_count):.0f})" |
| 421 | for pattern, expected_count in allowed |
| 422 | ) |
| 423 | print(summary_msg) |
| 424 | for pattern, expected_count in allowed: |
| 425 | actual_count = counts[pattern] |
| 426 | assert ( |
| 427 | abs(actual_count - expected_count) <= tolerance |
| 428 | ), f"Count for {pattern} off by more than {tolerance_percent:.1f}% of shots. Actual={actual_count}, Expected={expected_count:.0f}, noise#{noise_number}, Program={program}." |
| 429 | |