diff --git a/wai/wai-quantum-service/Cargo.toml b/wai/wai-quantum-service/Cargo.toml index d09f06dfe..18b304ea2 100644 --- a/wai/wai-quantum-service/Cargo.toml +++ b/wai/wai-quantum-service/Cargo.toml @@ -17,7 +17,7 @@ serde = { version = "1", features = ["derive"] } serde_json = { version = "1", features = ["float_roundtrip"] } blake3 = "1" getrandom = "0.4" -wai-quantum = { path = "../wai-quantum", default-features = false, features = ["quantum", "quantum_toolchain", "quantum_receipt", "quantum_sv", "quantum_tn", "quantum_tableau", "quantum_mpsim", "quantum_spd", "quantum_noisy", "quantum_sv_metal"] } +wai-quantum = { path = "../wai-quantum", default-features = false, features = ["quantum", "quantum_toolchain", "quantum_receipt", "quantum_sv", "quantum_tn", "quantum_tableau", "quantum_cutrank", "quantum_mpsim", "quantum_spd", "quantum_noisy", "quantum_sv_metal"] } # proc_pid_rusage and waitid, for each job's energy (src/energy.rs). [target.'cfg(target_os = "macos")'.dependencies] diff --git a/wai/wai-quantum-service/README.md b/wai/wai-quantum-service/README.md index cf6393608..053a6c406 100644 --- a/wai/wai-quantum-service/README.md +++ b/wai/wai-quantum-service/README.md @@ -41,12 +41,19 @@ The program may use the standard gates, user `gate` definitions, and the control | `pauli` | any circuit, up to 1024 qubits | expectations | bounded | | `mps` | any circuit, up to 4096 qubits, at bond dimension up to 256; best when entanglement stays low | counts, expectations, amplitudes | exact when nothing is truncated, otherwise heuristic | | `tensor_network` | any circuit | amplitudes | exact | +| `cutrank` | a Clifford circuit, then one layer of one-qubit gates at any angles, up to 4096 qubits, as wide as the state's cut-rank allows | counts, amplitudes | exact (amplitudes up to one shared phase); single when only a single-precision table fits | | `density_matrix` | noisy circuits, up to 16 qubits within a job's memory | counts, expectations | exact | | `trajectories` | noisy circuits, up to 32 qubits | counts, expectations | estimate | How the methods run: - **statevector** groups consecutive gates into passes over the state. Each pass applies its gates tile by tile, in a buffer that fits a core's cache. The grouping never changes a bit of the result. +- **cutrank** takes a Clifford circuit followed, on each qubit, by one-qubit gates at any angles (a qubit's gates from its first non-Clifford one on act on that qubit alone). + - It writes the Clifford part's state as a graph state with local Cliffords, and sums the amplitude one qubit at a time along an order of the qubits, in a table of `2^cut` entries. Here `cut` is the state's cut-rank, its entanglement in bits, across the cut the sum has reached. + - The planner searches for an order with narrower cuts. The cost is set by the widest cut, `max_cut` in the result, not by the number of qubits. + - Counts are exact samples, drawn qubit by qubit from the table's conditional probabilities. + - Amplitudes are exact up to one phase shared by all of the job's amplitudes, because a stabilizer tableau does not fix the state's global phase; the result's `global_phase` says so. + - The table is held in double precision when it fits the job's memory, and in single precision otherwise (accuracy class `single`). - **pauli** reports each value with its truncation `bound`. It holds at most as many Pauli terms as the job's memory allows; a sum that would outgrow it raises its threshold until it fits (`threshold_reached`), and the bound covers everything that dropped. Unless `ladder` is `false`, it also reports a `ladder`: the value and bound at up to three thresholds a decade above the one asked for, so you can see the value settle as the threshold drops. - **stabilizer**, with noise, samples by Pauli frames against one noiseless outcome of the tableau, 64 shots to a machine word. Its counts have the law of independent noisy shots. Its expectation values are estimates with standard errors, from at least 10,000 shots unless `trajectories` sets the number. Amplitude damping is not a Pauli channel and is refused. - **trajectories** treats every shot as an independent noisy trajectory in law. Shots that took the same errors share one state, so the error-free part of a circuit, and the error-free majority of the shots, run once. `histories` in the result counts the distinct error histories that ran. @@ -71,7 +78,7 @@ curl -s https://wai.transaction.science/quantum-jobs/v1/run -H 'content-type: ap "observables": ["Z2", "X2"], "shots": 1000 }' ``` -`auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it. +`auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `cutrank`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it. The accuracy classes: diff --git a/wai/wai-quantum-service/smoke.sh b/wai/wai-quantum-service/smoke.sh index 39611dc5d..dc2c74b77 100644 --- a/wai/wai-quantum-service/smoke.sh +++ b/wai/wai-quantum-service/smoke.sh @@ -22,7 +22,7 @@ py 'import json,sys d=json.load(open(sys.argv[1])) assert d["service"]=="wai-quantum-service", d ms={m["method"] for m in d["methods"]} -assert {"exact","statevector","stabilizer","pauli","tensor_network","density_matrix","trajectories"} <= ms, ms +assert {"exact","statevector","stabilizer","pauli","tensor_network","cutrank","density_matrix","trajectories"} <= ms, ms assert len(d["receipt_key"])==64 print("describe: ok,", len(ms), "methods, receipt key", d["receipt_key"][:16]+"…")' "$tmp/describe.json" @@ -236,6 +236,21 @@ assert d["status"]=="done" and abs(e[0]-math.cos(0.7))<1e-12 and abs(e[1]-math.s assert r["dynamic"]["branches_weighed"]==4 and sum(r["counts"].values())==2000, r print("teleport (mid-circuit measurement, if): ok,", r["method"], e)' "$tmp/tele.out" +# A Clifford circuit on 60 qubits, then T and H on each: past the state +# vector, and not Clifford, so auto samples it exactly by cut-rank. Its graph +# (a ring with chords of length 5) keeps every cut narrow. +py 'import json +n=60; q="OPENQASM 3.0; qubit[%d] q; bit[%d] c; h q;" % (n, n) +q+="".join(" cz q[%d], q[%d]; cz q[%d], q[%d];" % (i, (i+1)%n, i, (i+5)%n) for i in range(n)) +q+="".join(" t q[%d]; h q[%d];" % (i, i) for i in range(n)) + " c = measure q;" +print(json.dumps({"qasm": q, "shots": 100, "seed": 3}))' > "$tmp/cut.json" +post /v1/run "$tmp/cut.json" > "$tmp/cut.out" +py 'import json,sys +d=json.load(open(sys.argv[1])); r=d["result"] +assert d["status"]=="done" and r["method"]=="cutrank" and sum(r["counts"].values())==100, d +assert r["accuracy"]["class"]=="exact" and r["accuracy"]["max_cut"]<=12, r["accuracy"] +print("cutrank (60 qubits, Clifford then T and H): ok, max cut", r["accuracy"]["max_cut"], "bits,", r["distinct_outcomes"], "distinct outcomes")' "$tmp/cut.out" + # A refusal says why. py 'import json; print(json.dumps({"qasm": "qubit[1] q; bit c; c = measure q; while (c) { x q; c = measure q; }"}))' > "$tmp/bad.json" code=$(curl -sS -o "$tmp/bad.out" -w '%{http_code}' -X POST "$U/v1/jobs" -H 'content-type: application/json' --data-binary "@$tmp/bad.json") diff --git a/wai/wai-quantum-service/src/catalog.rs b/wai/wai-quantum-service/src/catalog.rs index eb3711fc8..b1ea0cc62 100644 --- a/wai/wai-quantum-service/src/catalog.rs +++ b/wai/wai-quantum-service/src/catalog.rs @@ -14,12 +14,13 @@ //! | `pauli` | any, ≤ 1024 qubits | expectations | bounded | //! | `mps` | any, ≤ 4096 qubits, best on low-entanglement circuits | counts, expectations, amplitudes | exact, or heuristic with its fidelity estimate | //! | `tensor_network` | any | amplitudes | exact | +//! | `cutrank` | Clifford, then one layer of one-qubit gates at any angles, ≤ 4096 qubits, as wide as its cut-rank allows | counts, amplitudes | exact (amplitudes up to one shared phase), or single | //! | `density_matrix` | any, noisy, ≤ 16 qubits in a job's memory | counts, expectations | exact | //! | `trajectories` | any, noisy, ≤ 32 qubits | counts, expectations | estimate | //! //! `auto` picks, in order: with noise, `density_matrix`, `stabilizer`, then `trajectories`; //! for counts, `stabilizer`, `exact`, `statevector`, `statevector_single`, -//! `mps`; for expectations only, `stabilizer`, `exact`, `statevector`, +//! `cutrank`, `mps`; for expectations only, `stabilizer`, `exact`, `statevector`, //! `statevector_single`, `pauli`, `mps`; for amplitudes only, `exact`, //! `statevector`, `statevector_single`, `tensor_network`, `mps`. A method that cannot take the circuit is //! skipped, and the plan says why the one chosen was chosen. @@ -53,6 +54,18 @@ pub const MPS_MAX_WORK: u64 = 400_000_000_000; /// shots): about a minute. pub const STABILIZER_MAX_WORK: u64 = 500_000_000_000; +/// The most qubits `cutrank` takes: its graph form is an `n × n` bit matrix +/// reduced in `O(n³/64)` word operations. +pub const CUTRANK_MAX_QUBITS: u32 = 4096; +/// The widest cut a `cutrank` table may span, whatever the memory: `2³⁴` +/// entries. +pub const CUTRANK_MAX_CUT: usize = 34; +/// The table entries a `cutrank` job's passes may touch: about ten minutes. +pub const CUTRANK_MAX_WORK: u64 = 4_000_000_000_000; +/// The word operations the planner spends searching for an order of the +/// qubits with narrower cuts: about a second. +const CUTRANK_ANNEAL_OPS: u64 = 2_000_000_000; + const MIB: u64 = 1 << 20; const GIB: u64 = 1 << 30; @@ -90,6 +103,7 @@ pub fn describe(job_bytes: u64) -> Value { { "method": "pauli", "accuracy": "bounded", "circuits": "any", "max_qubits": PAULI_MAX_QUBITS, "answers": ["expectations"], "options": { "threshold": "coefficients below it are dropped (default 1e-5); the bound covers what was dropped", "max_weight": "Pauli strings heavier than this are dropped", "ladder": "also report the value at up to three thresholds a decade apart above it, to show it settling (default true)", "max_terms": "the most Pauli terms to hold (default: what the job's memory holds); past it the threshold rises, and the bound covers what drops" } }, { "method": "mps", "accuracy": "exact when nothing is truncated, otherwise heuristic", "circuits": "any; best on low-entanglement circuits", "max_qubits": MPS_MAX_QUBITS, "answers": ["counts", "expectations", "amplitudes"], "options": { "max_bond": "the largest bond dimension kept (default 64, at most 256)", "cutoff": "singular values whose squared share is below it are dropped (default 1e-12)" } }, { "method": "tensor_network", "accuracy": "exact", "circuits": "any", "max_qubits": qasm::MAX_QUBITS, "answers": ["amplitudes"] }, + { "method": "cutrank", "accuracy": "exact (amplitudes up to one phase shared by the job's amplitudes); single when only a single-precision table fits", "circuits": "a Clifford circuit, then one layer of one-qubit gates at any angles: each qubit's gates from its first non-Clifford one on act on that qubit alone; the cost is set by the state's cut-rank along an order of the qubits, not by their number", "max_qubits": CUTRANK_MAX_QUBITS, "answers": ["counts", "amplitudes"] }, { "method": "density_matrix", "accuracy": "exact", "circuits": "any, with noise", "max_qubits": density_max_qubits(job_bytes), "answers": ["counts", "expectations"] }, { "method": "trajectories", "accuracy": "estimate", "circuits": "any, with noise", "max_qubits": TRAJECTORY_MAX_QUBITS.min(statevector_max_qubits(job_bytes)), "answers": ["counts", "expectations"], "options": { "trajectories": "how many (default: the shots, or 1000 for expectations alone); shots that take the same errors share one state" } }, ]) @@ -260,6 +274,39 @@ pub fn mps_work(r: &Request, chi: u64) -> u64 { updates.saturating_mul((2 * chi).pow(3)).saturating_add(r.shots.saturating_mul(r.program.n_qubits as u64 * chi * chi * 2)) } +/// A `cutrank` job's state: the program's Clifford part as a graph form with +/// its closing one-qubit gates, summed along an order found by a seeded +/// search. The planner and the run build it alike, so they agree on the order. +pub fn product_measured(r: &Request) -> Result { + use wai_quantum::quantum_cutrank::ProductMeasured; + use wai_quantum::quantum_tableau::Tableau; + let n = r.program.n_qubits; + if n > CUTRANK_MAX_QUBITS { + return Err(format!("cutrank holds at most {CUTRANK_MAX_QUBITS} qubits")); + } + let (bulk, last) = r.program.product_measured()?; + let t = Tableau::run(n, &bulk).map_err(|e| e.to_string())?; + let pm = ProductMeasured::new(&t, &last).map_err(|e| e.to_string())?; + // One evaluation of an order rebuilds a basis at each cut: about + // k · cut_k · words operations at qubit k. + let words = (n as u64).div_ceil(64); + let per_order: u64 = pm.cut_ranks().iter().enumerate().map(|(k, &c)| (k as u64 + 1) * (c as u64 + 1) * words).sum(); + let moves = (CUTRANK_ANNEAL_OPS / per_order.max(1)).min(20_000) as usize; + if moves == 0 { + return Ok(pm); + } + let order = pm.anneal_order(moves, 1); + pm.with_order(&order).map_err(|e| e.to_string()) +} + +/// The memory a `cutrank` job takes: its largest table with a margin, the +/// graph form, and room to run. +pub fn cutrank_bytes(pm: &wai_quantum::quantum_cutrank::ProductMeasured, precision: wai_quantum::quantum_cutrank::Precision) -> u64 { + let n = pm.qubits() as u64; + let table = pm.table_bytes(precision); + table + table / 16 + n * n.div_ceil(64) * 8 * 4 + 256 * MIB +} + /// Whether `method` can take the job, and the plan if so; otherwise why not. /// A program that measures in the middle, resets after a qubit's first gate, /// or acts on measured bits ([`crate::dynamic`]): the dense state vector and @@ -412,6 +459,27 @@ fn try_method(method: &str, r: &Request, dyadic: &Result { + use wai_quantum::quantum_cutrank::Precision; + if wants_expectations { + return Err("cutrank answers counts and amplitudes".into()); + } + let pm = product_measured(r)?; + let cut = pm.max_cut(); + let precision = if cutrank_bytes(&pm, Precision::Double) <= cap { Precision::Double } else { Precision::Single }; + let memory = cutrank_bytes(&pm, precision); + if cut > CUTRANK_MAX_CUT || memory > cap { + return Err(format!("the widest cut along the best order found is {cut} bits: a table of 2^{cut} entries, more than a job's memory")); + } + // A sample is a read of the table and a pass for each qubit. + let per = pm.work().ceil() as u64; + let work = per.saturating_mul(r.amplitudes.len() as u64 + 3 * r.shots); + if work > CUTRANK_MAX_WORK { + return Err(format!("{work} table entries touched is more than one job may take ({CUTRANK_MAX_WORK}); ask for fewer shots or amplitudes")); + } + let held = if precision == Precision::Double { "double" } else { "single" }; + Ok(plan(memory, work, format!("Clifford circuit then one layer of one-qubit gates on {n} qubits: summed exactly by cut-rank, the widest cut {cut} bits, the table in {held} precision"))) + } "tensor_network" => { if !wants_amplitudes || counts || wants_expectations { return Err("tensor_network answers amplitudes only".into()); @@ -498,7 +566,7 @@ pub fn plan(v: &Value, job_bytes: u64) -> Result { let order: &[&str] = if r.noise.is_some() { &["density_matrix", "stabilizer", "trajectories"] } else if r.shots > 0 { - &["stabilizer", "exact", "statevector", "statevector_single", "mps"] + &["stabilizer", "exact", "statevector", "statevector_single", "cutrank", "mps"] } else if !r.observables.is_empty() && r.amplitudes.is_empty() { &["stabilizer", "exact", "statevector", "statevector_single", "pauli", "mps"] } else { diff --git a/wai/wai-quantum-service/src/discovery.rs b/wai/wai-quantum-service/src/discovery.rs index 01ba21291..bb40633c3 100644 --- a/wai/wai-quantum-service/src/discovery.rs +++ b/wai/wai-quantum-service/src/discovery.rs @@ -108,7 +108,7 @@ pub fn openapi(server: &str) -> Value { "required": ["qasm"], "properties": { "qasm": { "type": "string", "description": "The program: OpenQASM 3, or OpenQASM 2. Registers, standard gates (h, x, y, z, s, sdg, t, tdg, sx, rx, ry, rz, p, u, cx, cy, cz, ch, swap, iswap, cp, crx, cry, crz, ccx, ccz, cswap, rxx, ryy, rzz), the control modifiers ctrl @, ctrl(n) @ and negctrl @ on x y z h s sdg t tdg p and the controlled names, user gate definitions, measurement anywhere, reset anywhere, and `if (c[0]) …` / `if (c == 3) { … } else { … }` on measured bits (OpenQASM 2's `if(c==n) gate;` too). A program that measures in the middle, resets after a gate or acts on measured bits runs on statevector or stabilizer: counts split their shots at each measurement, expectation values are the exact mixture over the branches, of the state before the final measurements. Loops, inv @ and pow @ are refused, with the line." }, - "method": { "type": "string", "enum": ["auto", "exact", "statevector", "statevector_single", "stabilizer", "pauli", "mps", "tensor_network", "density_matrix", "trajectories"], "default": "auto" }, + "method": { "type": "string", "enum": ["auto", "exact", "statevector", "statevector_single", "stabilizer", "pauli", "mps", "tensor_network", "cutrank", "density_matrix", "trajectories"], "default": "auto" }, "shots": { "type": "integer", "minimum": 0, "maximum": 1000000, "description": "Measurement shots (default 1024 when nothing else is asked for)." }, "seed": { "type": "integer", "default": 1, "description": "The sampling seed: the same job and seed give the same counts." }, "observables": { "type": "array", "items": pauli }, diff --git a/wai/wai-quantum-service/src/mcp.rs b/wai/wai-quantum-service/src/mcp.rs index 81b266551..71976a925 100644 --- a/wai/wai-quantum-service/src/mcp.rs +++ b/wai/wai-quantum-service/src/mcp.rs @@ -31,7 +31,7 @@ fn tools() -> Value { "noise": { "type": "object", "description": "Optional noise, each a probability: depolarizing_1q, depolarizing_2q, amplitude_damping, dephasing (after each gate on its qubits), readout (each measured bit).", "properties": { "depolarizing_1q": { "type": "number" }, "depolarizing_2q": { "type": "number" }, "amplitude_damping": { "type": "number" }, "dephasing": { "type": "number" }, "readout": { "type": "number" } } }, - "method": { "type": "string", "enum": ["auto", "exact", "statevector", "statevector_single", "stabilizer", "pauli", "mps", "tensor_network", "density_matrix", "trajectories"], "description": "Leave as auto unless you need a particular method." }, + "method": { "type": "string", "enum": ["auto", "exact", "statevector", "statevector_single", "stabilizer", "pauli", "mps", "tensor_network", "cutrank", "density_matrix", "trajectories"], "description": "Leave as auto unless you need a particular method." }, "seed": { "type": "integer", "description": "Sampling seed (default 1); the same job and seed give the same counts." }, "options": { "type": "object", "description": "Advanced, per method: max_bond (mps), threshold and max_terms (pauli), fusion (statevector), trajectories." } }); diff --git a/wai/wai-quantum-service/src/qasm.rs b/wai/wai-quantum-service/src/qasm.rs index 1a4cbfb9d..b084e6ea6 100644 --- a/wai/wai-quantum-service/src/qasm.rs +++ b/wai/wai-quantum-service/src/qasm.rs @@ -1732,12 +1732,54 @@ impl Program { .collect()) } + /// The program as a Clifford circuit, then one layer of one-qubit gates: + /// each qubit's gates from its first that is not Clifford on, multiplied + /// into one matrix (the identity for a qubit with none). The input of the + /// `cutrank` method. Refused when a gate on two or more qubits comes after + /// a qubit's first non-Clifford gate, or is not Clifford itself. + pub fn product_measured(&self) -> Result<(Vec, Vec<[C; 4]>), String> { + let one = complex(1.0, 0.0); + let zero = complex(0.0, 0.0); + let mul = |a: &[C], b: &[C; 4]| -> [C; 4] { + let m = |x: C, y: C| complex(x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re); + let add = |x: C, y: C| complex(x.re + y.re, x.im + y.im); + [add(m(a[0], b[0]), m(a[1], b[2])), add(m(a[0], b[1]), m(a[1], b[3])), add(m(a[2], b[0]), m(a[3], b[2])), add(m(a[2], b[1]), m(a[3], b[3]))] + }; + let mut bulk = Vec::new(); + let mut tail: Vec> = vec![None; self.n_qubits as usize]; + for op in &self.ops { + if let [q] = op.qubits[..] { + let q = q as usize; + if let Some(t) = tail[q] { + tail[q] = Some(mul(&op.matrix().matrix, &t)); + } else { + match Self::lower_clifford(std::slice::from_ref(op)) { + Ok((gates, _)) => bulk.extend(gates), + Err(_) => tail[q] = Some(mul(&op.matrix().matrix, &[one, zero, zero, one])), + } + } + continue; + } + if let Some(&q) = op.qubits.iter().find(|&&q| tail[q as usize].is_some()) { + return Err(format!("{:?} on qubits {:?} comes after a gate on qubit {q} that is not Clifford; cutrank takes a Clifford circuit, then one layer of one-qubit gates", op.kind, op.qubits)); + } + let (gates, _) = Self::lower_clifford(std::slice::from_ref(op)).map_err(|e| format!("{e}; cutrank takes a Clifford circuit, then one layer of one-qubit gates"))?; + bulk.extend(gates); + } + Ok((bulk, tail.into_iter().map(|t| t.unwrap_or([one, zero, zero, one])).collect())) + } + /// The Clifford lowering, and where each op's gates end in it. fn clifford_bounds(&self) -> Result<(Vec, Vec), String> { + Self::lower_clifford(&self.ops) + } + + /// The Clifford lowering of `ops`, and where each op's gates end in it. + fn lower_clifford(ops: &[Op]) -> Result<(Vec, Vec), String> { use Kind::*; let mut out = Vec::new(); - let mut ends = Vec::with_capacity(self.ops.len()); - for op in &self.ops { + let mut ends = Vec::with_capacity(ops.len()); + for op in ops { let q = &op.qubits; let p = &op.params; let not = || format!("{:?}{} is not a Clifford gate", op.kind, if p.is_empty() { String::new() } else { format!("({})", p.iter().map(|x| format!("{x}")).collect::>().join(", ")) }); diff --git a/wai/wai-quantum-service/src/run.rs b/wai/wai-quantum-service/src/run.rs index 9bacd9cc4..d71579490 100644 --- a/wai/wai-quantum-service/src/run.rs +++ b/wai/wai-quantum-service/src/run.rs @@ -361,6 +361,36 @@ pub fn run(request: &Value, plan: &Plan) -> Result { out["amplitudes"] = Value::Array(entries); } } + "cutrank" => { + use wai_quantum::quantum_cutrank::Precision; + // The same order the planner found, and the precision its memory allows. + let pm = catalog::product_measured(&r)?; + let precision = if plan.memory_bytes >= catalog::cutrank_bytes(&pm, Precision::Double) { Precision::Double } else { Precision::Single }; + let cut = pm.max_cut(); + out["accuracy"] = json!({ + "class": if precision == Precision::Double { "exact" } else { "single" }, + "max_cut": cut, + "order_work": pm.work(), + "note": "the state's sum taken qubit by qubit along an order of the qubits, in a table of 2^max_cut entries (single class: the table stored in single precision, its arithmetic in double); counts are exact samples, seeded; amplitudes are exact up to one phase shared by all of the job's amplitudes, since a stabilizer tableau does not fix the state's global phase" + }); + if r.shots > 0 { + let mut map: HashMap = HashMap::new(); + for (bits, _) in pm.sample_in(r.shots as usize, r.seed, cut, precision).map_err(|e| e.to_string())? { + *map.entry(r.program.bitstring(|q| bits[q as usize])).or_insert(0) += 1; + } + merge(&mut out, counts(map)); + } + if !r.amplitudes.is_empty() { + let mut entries = Vec::new(); + for b in &r.amplitudes { + let x: Vec = (0..n).map(|q| b.as_bytes()[(n - 1 - q) as usize] == b'1').collect(); + let a = pm.amplitude_in(&x, cut, precision).map_err(|e| e.to_string())?; + entries.push(json!({ "bits": b, "re": a.re, "im": a.im })); + } + out["amplitudes"] = Value::Array(entries); + out["global_phase"] = json!("not fixed: the amplitudes share one unknown phase; their magnitudes and relative phases are exact"); + } + } "tensor_network" => { use wai_quantum::quantum_tn::{Network, PathOptions, contract_plan, find_path}; let gates = r.program.matrices(); @@ -885,6 +915,113 @@ mod tests { assert_eq!(counts.values().map(|v| v.as_u64().unwrap()).sum::(), 1000); } + /// A Clifford circuit on `n` qubits (Hadamards, phases, controlled-Z and + /// CX between pseudo-random pairs), then one layer of one-qubit gates at + /// generic angles on every qubit. + fn clifford_then_layer(n: u32, pairs: usize, seed: u64) -> String { + let mut s = seed; + let mut next = |m: u32| { + s = s.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1_442_695_040_888_963_407); + ((s >> 33) % m as u64) as u32 + }; + let mut q = format!("OPENQASM 3.0; qubit[{n}] q; h q;"); + for _ in 0..pairs { + let a = next(n); + let b = (a + 1 + next(n - 1)) % n; + q += &match next(4) { + 0 => format!(" cz q[{a}], q[{b}];"), + 1 => format!(" cx q[{a}], q[{b}]; s q[{b}];"), + 2 => format!(" h q[{a}]; cz q[{a}], q[{b}];"), + _ => format!(" sdg q[{a}]; cx q[{b}], q[{a}];"), + }; + } + for i in 0..n { + q += &format!(" t q[{i}]; ry({}) q[{i}]; rz({}) q[{i}];", 0.3 + 0.1 * i as f64, 1.1 - 0.07 * i as f64); + } + q + } + + #[test] + fn cutrank_amplitudes_are_the_state_vectors_up_to_one_phase() { + for (n, seed) in [(6u32, 1u64), (9, 2), (11, 3)] { + let qasm = clifford_then_layer(n, 3 * n as usize, seed); + let bits: Vec = (0..24u64).map(|k| (0..n).map(|q| if (k.wrapping_mul(0x9e37_79b9) >> q) & 1 == 1 { '1' } else { '0' }).collect()).collect(); + let cr = go(json!({ "qasm": qasm, "method": "cutrank", "amplitudes": bits })); + let sv = go(json!({ "qasm": qasm, "method": "statevector", "amplitudes": bits })); + assert_eq!(cr["method"], "cutrank", "{cr}"); + assert_eq!(cr["accuracy"]["class"], "exact", "{cr}"); + assert!(cr["global_phase"].as_str().is_some_and(|g| g.starts_with("not fixed")), "{cr}"); + let amp = |v: &Value, i: usize| (v["amplitudes"][i]["re"].as_f64().unwrap(), v["amplitudes"][i]["im"].as_f64().unwrap()); + // The shared phase, from the largest amplitude. + let k = (0..bits.len()).max_by(|&a, &b| { + let (x, y) = (amp(&sv, a), amp(&sv, b)); + (x.0 * x.0 + x.1 * x.1).total_cmp(&(y.0 * y.0 + y.1 * y.1)) + }).unwrap(); + let ((cr_k_re, cr_k_im), (sv_k_re, sv_k_im)) = (amp(&cr, k), amp(&sv, k)); + let norm = sv_k_re * sv_k_re + sv_k_im * sv_k_im; + // phase = cr_k / sv_k + let (pr, pi) = ((cr_k_re * sv_k_re + cr_k_im * sv_k_im) / norm, (cr_k_im * sv_k_re - cr_k_re * sv_k_im) / norm); + assert!(((pr * pr + pi * pi) - 1.0).abs() < 1e-9, "n {n}: |phase| {}", (pr * pr + pi * pi).sqrt()); + for (i, b) in bits.iter().enumerate() { + let ((cr_re, cr_im), (sv_re, sv_im)) = (amp(&cr, i), amp(&sv, i)); + let (want_re, want_im) = (pr * sv_re - pi * sv_im, pr * sv_im + pi * sv_re); + assert!((cr_re - want_re).abs() < 1e-10 && (cr_im - want_im).abs() < 1e-10, "n {n}, {b}: cutrank {cr_re} {cr_im}, statevector × phase {want_re} {want_im}"); + } + } + } + + #[test] + fn cutrank_counts_follow_the_exact_distribution() { + let n = 5u32; + let qasm = clifford_then_layer(n, 9, 7); + let all: Vec = (0..1u32 << n).map(|x| format!("{x:05b}")).collect(); + let sv = go(json!({ "qasm": qasm, "method": "statevector", "amplitudes": all })); + let shots = 40_000u64; + let cr = go(json!({ "qasm": qasm, "method": "cutrank", "shots": shots, "seed": 5 })); + let counts = cr["counts"].as_object().unwrap(); + assert_eq!(counts.values().map(|v| v.as_u64().unwrap()).sum::(), shots); + for (i, b) in all.iter().enumerate() { + let (re, im) = (sv["amplitudes"][i]["re"].as_f64().unwrap(), sv["amplitudes"][i]["im"].as_f64().unwrap()); + let p = re * re + im * im; + let got = counts.get(b).and_then(Value::as_u64).unwrap_or(0) as f64; + let sd = (shots as f64 * p * (1.0 - p)).sqrt().max(1.0); + assert!((got - shots as f64 * p).abs() < 5.0 * sd, "{b}: {got} of {shots}, expected {:.1}", shots as f64 * p); + } + // Seeded. + let again = go(json!({ "qasm": qasm, "method": "cutrank", "shots": shots, "seed": 5 })); + assert_eq!(again["counts"], cr["counts"]); + } + + #[test] + fn auto_samples_a_wide_clifford_then_layer_circuit_by_cutrank() { + // 40 qubits: past the dense state vector in a test job's memory, and not Clifford. + let n = 40u32; + let mut qasm = format!("OPENQASM 3.0; qubit[{n}] q; bit[{n}] c; h q;"); + for i in 0..n { + qasm += &format!(" cz q[{i}], q[{}];", (i + 1) % n); + qasm += &format!(" cz q[{i}], q[{}];", (i + 5) % n); + } + for i in 0..n { + qasm += &format!(" t q[{i}]; h q[{i}];"); + } + qasm += " c = measure q;"; + let r = go(json!({ "qasm": qasm, "shots": 200, "seed": 3 })); + assert_eq!(r["method"], "cutrank", "{r}"); + assert_eq!(r["counts"].as_object().unwrap().values().map(|v| v.as_u64().unwrap()).sum::(), 200); + assert!(r["accuracy"]["max_cut"].as_u64().unwrap() <= 12, "{r}"); + } + + #[test] + fn cutrank_refuses_what_it_cannot_sum() { + let plan_of = |j: Value| plan(&j, CAP); + let e = plan_of(json!({ "qasm": "qubit[2] q; rx(0.3) q[0]; cx q[0], q[1];", "method": "cutrank", "shots": 10 })).unwrap_err(); + assert!(e.contains("comes after a gate on qubit 0 that is not Clifford"), "{e}"); + let e = plan_of(json!({ "qasm": "qubit[3] q; ccx q[0], q[1], q[2];", "method": "cutrank", "shots": 10 })).unwrap_err(); + assert!(e.contains("not a Clifford gate"), "{e}"); + let e = plan_of(json!({ "qasm": "qubit[2] q; h q[0]; t q[1];", "method": "cutrank", "observables": ["Z0"] })).unwrap_err(); + assert!(e.contains("counts and amplitudes"), "{e}"); + } + #[test] fn dynamic_runs_repeat_with_their_seed() { let qasm = TELEPORT.replace("THETA", "1.1");