ruqu_core/
control_theory.rs

1//! Hybrid classical-quantum control theory engine for QEC.
2//!
3//! Models the QEC feedback loop as a discrete-time control system:
4//! `Physical qubits -> Syndrome extraction -> Classical decode -> Correction -> Repeat`
5//!
6//! If classical decoding latency exceeds the syndrome extraction period, errors
7//! accumulate faster than they are corrected (the "backlog problem").
8
9use rand::rngs::StdRng;
10use rand::{Rng, SeedableRng};
11
12#[allow(unused_imports)]
13use crate::error::{QuantumError, Result};
14
15// -- 1. Control Loop Model --------------------------------------------------
16
17/// Full QEC control loop: plant (quantum) + controller (classical) + state.
18#[derive(Debug, Clone)]
19pub struct QecControlLoop {
20    pub plant: QuantumPlant,
21    pub controller: ClassicalController,
22    pub state: ControlState,
23}
24
25/// Physical parameters of the quantum error-correction code.
26#[derive(Debug, Clone)]
27pub struct QuantumPlant {
28    pub code_distance: u32,
29    pub physical_error_rate: f64,
30    pub num_data_qubits: u32,
31    pub coherence_time_ns: u64,
32}
33
34/// Classical decoder performance characteristics.
35#[derive(Debug, Clone)]
36pub struct ClassicalController {
37    pub decode_latency_ns: u64,
38    pub decode_throughput: f64,
39    pub accuracy: f64,
40}
41
42/// Evolving state of the control loop during execution.
43#[derive(Debug, Clone)]
44pub struct ControlState {
45    pub logical_error_rate: f64,
46    pub error_backlog: f64,
47    pub rounds_decoded: u64,
48    pub total_latency_ns: u64,
49}
50
51impl ControlState {
52    pub fn new() -> Self {
53        Self { logical_error_rate: 0.0, error_backlog: 0.0, rounds_decoded: 0, total_latency_ns: 0 }
54    }
55}
56
57impl Default for ControlState {
58    fn default() -> Self { Self::new() }
59}
60
61// -- 2. Stability Analysis ---------------------------------------------------
62
63/// Result of analyzing the control loop's stability.
64#[derive(Debug, Clone)]
65pub struct StabilityCondition {
66    pub is_stable: bool,
67    pub margin: f64,
68    pub critical_latency_ns: u64,
69    pub critical_error_rate: f64,
70    pub convergence_rate: f64,
71}
72
73/// Syndrome extraction period (ns) for distance-d surface code.
74/// 6 gate layers per cycle, ~20 ns per gate layer.
75fn syndrome_period_ns(distance: u32) -> u64 {
76    6 * (distance as u64) * 20
77}
78
79/// Analyze stability: the loop is stable when `decode_latency < syndrome_period`.
80pub fn analyze_stability(config: &QecControlLoop) -> StabilityCondition {
81    let d = config.plant.code_distance;
82    let p = config.plant.physical_error_rate;
83    let t_decode = config.controller.decode_latency_ns;
84    let acc = config.controller.accuracy;
85    let t_syndrome = syndrome_period_ns(d);
86
87    let margin = if t_decode == 0 { f64::INFINITY }
88                 else { (t_syndrome as f64 / t_decode as f64) - 1.0 };
89    let is_stable = t_decode < t_syndrome;
90    let critical_latency_ns = t_syndrome;
91    let critical_error_rate = 0.01 * acc;
92    let error_injection = p * (d as f64);
93    let convergence_rate = if t_syndrome > 0 {
94        1.0 - (t_decode as f64 / t_syndrome as f64) - error_injection
95    } else { -1.0 };
96
97    StabilityCondition { is_stable, margin, critical_latency_ns, critical_error_rate, convergence_rate }
98}
99
100/// Maximum code distance stable for a given controller and physical error rate.
101/// Iterates odd distances 3, 5, 7, ... until latency exceeds syndrome period.
102pub fn max_stable_distance(controller: &ClassicalController, error_rate: f64) -> u32 {
103    let mut best = 3u32;
104    for d in (3..=201).step_by(2) {
105        if controller.decode_latency_ns >= syndrome_period_ns(d) { break; }
106        if error_rate >= 0.01 * controller.accuracy { break; }
107        best = d;
108    }
109    best
110}
111
112/// Minimum decoder throughput (syndromes/sec) to keep up with the plant.
113pub fn min_throughput(plant: &QuantumPlant) -> f64 {
114    let t_ns = syndrome_period_ns(plant.code_distance);
115    if t_ns == 0 { return f64::INFINITY; }
116    1e9 / t_ns as f64
117}
118
119// -- 3. Resource Optimization ------------------------------------------------
120
121/// Available hardware resources.
122#[derive(Debug, Clone)]
123pub struct ResourceBudget {
124    pub total_physical_qubits: u32,
125    pub classical_cores: u32,
126    pub classical_clock_ghz: f64,
127    pub total_time_budget_us: u64,
128}
129
130/// A candidate allocation on the Pareto frontier.
131#[derive(Debug, Clone)]
132pub struct OptimalAllocation {
133    pub code_distance: u32,
134    pub logical_qubits: u32,
135    pub decode_threads: u32,
136    pub expected_logical_error_rate: f64,
137    pub pareto_score: f64,
138}
139
140/// Enumerate Pareto-optimal resource allocations sorted by descending score.
141pub fn optimize_allocation(
142    budget: &ResourceBudget, error_rate: f64, min_logical: u32,
143) -> Vec<OptimalAllocation> {
144    let mut candidates = Vec::new();
145    for d in (3u32..=99).step_by(2) {
146        let qpl = 2 * d * d - 2 * d + 1;
147        if qpl == 0 { continue; }
148        let max_logical = budget.total_physical_qubits / qpl;
149        if max_logical < min_logical { continue; }
150
151        let decode_ns = if budget.classical_cores > 0 && budget.classical_clock_ghz > 0.0 {
152            ((d as f64).powi(3) / (budget.classical_cores as f64 * budget.classical_clock_ghz)) as u64
153        } else { u64::MAX };
154        let decode_threads = budget.classical_cores.min(max_logical);
155
156        let p_th = 0.01_f64;
157        let ratio = error_rate / p_th;
158        let exp = (d as f64 + 1.0) / 2.0;
159        let p_logical = if ratio < 1.0 { 0.1 * ratio.powf(exp) }
160                        else { 1.0_f64.min(ratio.powf(exp)) };
161
162        let t_syn = syndrome_period_ns(d);
163        let round_time = t_syn.max(decode_ns);
164        let budget_ns = budget.total_time_budget_us * 1000;
165        if round_time == 0 || budget_ns / round_time == 0 { continue; }
166
167        let score = if p_logical > 0.0 && max_logical > 0 {
168            (max_logical as f64).log2() - p_logical.log10()
169        } else if max_logical > 0 { (max_logical as f64).log2() + 15.0 }
170          else { 0.0 };
171
172        candidates.push(OptimalAllocation {
173            code_distance: d, logical_qubits: max_logical, decode_threads,
174            expected_logical_error_rate: p_logical, pareto_score: score,
175        });
176    }
177    candidates.sort_by(|a, b| b.pareto_score.partial_cmp(&a.pareto_score).unwrap_or(std::cmp::Ordering::Equal));
178    candidates
179}
180
181// -- 4. Latency Budget Planner -----------------------------------------------
182
183/// Breakdown of time budgets for a single QEC round.
184#[derive(Debug, Clone)]
185pub struct LatencyBudget {
186    pub syndrome_extraction_ns: u64,
187    pub decode_ns: u64,
188    pub correction_ns: u64,
189    pub total_round_ns: u64,
190    pub slack_ns: i64,
191}
192
193/// Plan the latency budget for one QEC round at the given distance and decode time.
194pub fn plan_latency_budget(distance: u32, decode_ns_per_syndrome: u64) -> LatencyBudget {
195    let extraction_ns = syndrome_period_ns(distance);
196    let correction_ns: u64 = 20;
197    let total_round_ns = extraction_ns + decode_ns_per_syndrome + correction_ns;
198    let slack_ns = extraction_ns as i64 - (decode_ns_per_syndrome as i64 + correction_ns as i64);
199    LatencyBudget { syndrome_extraction_ns: extraction_ns, decode_ns: decode_ns_per_syndrome,
200                    correction_ns, total_round_ns, slack_ns }
201}
202
203// -- 5. Backlog Simulator ----------------------------------------------------
204
205/// Full trace of a simulated control loop execution.
206#[derive(Debug, Clone)]
207pub struct SimulationTrace {
208    pub rounds: Vec<RoundSnapshot>,
209    pub converged: bool,
210    pub final_logical_error_rate: f64,
211    pub max_backlog: f64,
212}
213
214/// Snapshot of a single simulation round.
215#[derive(Debug, Clone)]
216pub struct RoundSnapshot {
217    pub round: u64,
218    pub errors_this_round: u32,
219    pub errors_corrected: u32,
220    pub backlog: f64,
221    pub decode_latency_ns: u64,
222}
223
224/// Monte Carlo simulation of the QEC control loop with seeded RNG.
225pub fn simulate_control_loop(
226    config: &QecControlLoop, num_rounds: u64, seed: u64,
227) -> SimulationTrace {
228    let mut rng = StdRng::seed_from_u64(seed);
229    let d = config.plant.code_distance;
230    let p = config.plant.physical_error_rate;
231    let n_q = config.plant.num_data_qubits;
232    let t_decode = config.controller.decode_latency_ns;
233    let acc = config.controller.accuracy;
234    let t_syn = syndrome_period_ns(d);
235
236    let mut rounds = Vec::with_capacity(num_rounds as usize);
237    let (mut backlog, mut max_backlog) = (0.0_f64, 0.0_f64);
238    let mut logical_errors = 0u64;
239
240    for r in 0..num_rounds {
241        let mut errs: u32 = 0;
242        for _ in 0..n_q { if rng.gen::<f64>() < p { errs += 1; } }
243
244        let jitter = 0.8 + 0.4 * rng.gen::<f64>();
245        let actual_lat = (t_decode as f64 * jitter) as u64;
246        let in_time = actual_lat < t_syn;
247
248        let corrected = if in_time {
249            let mut c = 0u32;
250            for _ in 0..errs { if rng.gen::<f64>() < acc { c += 1; } }
251            c
252        } else { 0 };
253
254        let uncorrected = errs.saturating_sub(corrected);
255        backlog += uncorrected as f64;
256        if in_time && backlog > 0.0 { backlog -= (backlog * acc).min(backlog); }
257        if backlog > max_backlog { max_backlog = backlog; }
258        if uncorrected > (d.saturating_sub(1)) / 2 { logical_errors += 1; }
259
260        rounds.push(RoundSnapshot {
261            round: r, errors_this_round: errs, errors_corrected: corrected,
262            backlog, decode_latency_ns: actual_lat,
263        });
264    }
265
266    let final_logical_error_rate = if num_rounds > 0 { logical_errors as f64 / num_rounds as f64 } else { 0.0 };
267    SimulationTrace { rounds, converged: backlog < 1.0, final_logical_error_rate, max_backlog }
268}
269
270// -- 6. Scaling Laws ---------------------------------------------------------
271
272/// A power-law scaling relation: `y = prefactor * x^exponent`.
273#[derive(Debug, Clone)]
274pub struct ScalingLaw {
275    pub name: String,
276    pub exponent: f64,
277    pub prefactor: f64,
278}
279
280/// Classical overhead scaling for a named decoder.
281/// Known: `"union_find"` O(n), `"mwpm"` O(n^3), `"neural"` O(n). Default: O(n^2).
282pub fn classical_overhead_scaling(decoder_name: &str) -> ScalingLaw {
283    match decoder_name {
284        "union_find" => ScalingLaw { name: "Union-Find decoder".into(), exponent: 1.0, prefactor: 1.0 },
285        "mwpm"       => ScalingLaw { name: "Minimum Weight Perfect Matching".into(), exponent: 3.0, prefactor: 0.5 },
286        "neural"     => ScalingLaw { name: "Neural network decoder".into(), exponent: 1.0, prefactor: 10.0 },
287        _            => ScalingLaw { name: format!("Generic decoder ({})", decoder_name), exponent: 2.0, prefactor: 1.0 },
288    }
289}
290
291/// Logical error rate scaling: p_L ~ prefactor * (p/p_th)^exponent per distance step.
292/// Below threshold the exponent is the suppression factor lambda = -ln(p/p_th).
293pub fn logical_error_scaling(physical_rate: f64, threshold: f64) -> ScalingLaw {
294    if threshold <= 0.0 || physical_rate <= 0.0 {
295        return ScalingLaw { name: "Logical error scaling (degenerate)".into(), exponent: 0.0, prefactor: 1.0 };
296    }
297    if physical_rate >= threshold {
298        return ScalingLaw { name: "Logical error scaling (above threshold)".into(), exponent: 0.0, prefactor: 1.0 };
299    }
300    let lambda = -(physical_rate / threshold).ln();
301    ScalingLaw { name: "Logical error scaling (below threshold)".into(), exponent: lambda, prefactor: 0.1 }
302}
303
304// == Tests ===================================================================
305
306#[cfg(test)]
307mod tests {
308    use super::*;
309
310    fn make_plant(d: u32, p: f64) -> QuantumPlant {
311        QuantumPlant { code_distance: d, physical_error_rate: p, num_data_qubits: d * d, coherence_time_ns: 100_000 }
312    }
313    fn make_controller(lat: u64, tp: f64, acc: f64) -> ClassicalController {
314        ClassicalController { decode_latency_ns: lat, decode_throughput: tp, accuracy: acc }
315    }
316    fn make_loop(d: u32, p: f64, lat: u64) -> QecControlLoop {
317        QecControlLoop { plant: make_plant(d, p), controller: make_controller(lat, 1e6, 0.99), state: ControlState::new() }
318    }
319
320    #[test] fn test_control_state_new() {
321        let s = ControlState::new();
322        assert_eq!(s.logical_error_rate, 0.0); assert_eq!(s.error_backlog, 0.0);
323        assert_eq!(s.rounds_decoded, 0); assert_eq!(s.total_latency_ns, 0);
324    }
325    #[test] fn test_control_state_default() { assert_eq!(ControlState::default().rounds_decoded, 0); }
326
327    #[test] fn test_syndrome_period_scales() {
328        assert!(syndrome_period_ns(3) < syndrome_period_ns(5));
329        assert!(syndrome_period_ns(5) < syndrome_period_ns(7));
330    }
331    #[test] fn test_syndrome_period_d3() { assert_eq!(syndrome_period_ns(3), 360); }
332
333    #[test] fn test_stable_loop() {
334        let c = analyze_stability(&make_loop(5, 0.001, 100));
335        assert!(c.is_stable); assert!(c.margin > 0.0); assert!(c.convergence_rate > 0.0);
336    }
337    #[test] fn test_unstable_loop() {
338        let c = analyze_stability(&make_loop(3, 0.001, 1000));
339        assert!(!c.is_stable); assert!(c.margin < 0.0);
340    }
341    #[test] fn test_stability_critical_latency() {
342        assert_eq!(analyze_stability(&make_loop(5, 0.001, 100)).critical_latency_ns, syndrome_period_ns(5));
343    }
344    #[test] fn test_stability_zero_decode() {
345        let c = analyze_stability(&make_loop(3, 0.001, 0));
346        assert!(c.is_stable); assert!(c.margin.is_infinite());
347    }
348
349    #[test] fn test_max_stable_fast() { assert!(max_stable_distance(&make_controller(100, 1e7, 0.99), 0.001) >= 3); }
350    #[test] fn test_max_stable_slow() { assert!(max_stable_distance(&make_controller(10_000, 1e5, 0.99), 0.001) >= 3); }
351    #[test] fn test_max_stable_above_thresh() { assert_eq!(max_stable_distance(&make_controller(100, 1e7, 0.99), 0.5), 3); }
352
353    #[test] fn test_min_throughput_d3() {
354        let tp = min_throughput(&make_plant(3, 0.001));
355        assert!(tp > 2e6 && tp < 3e6);
356    }
357    #[test] fn test_min_throughput_ordering() {
358        assert!(min_throughput(&make_plant(3, 0.001)) > min_throughput(&make_plant(5, 0.001)));
359    }
360
361    #[test] fn test_optimize_basic() {
362        let b = ResourceBudget { total_physical_qubits: 10_000, classical_cores: 8, classical_clock_ghz: 3.0, total_time_budget_us: 1_000 };
363        let a = optimize_allocation(&b, 0.001, 1);
364        assert!(!a.is_empty());
365        for w in a.windows(2) { assert!(w[0].pareto_score >= w[1].pareto_score); }
366    }
367    #[test] fn test_optimize_min_logical() {
368        let b = ResourceBudget { total_physical_qubits: 100, classical_cores: 4, classical_clock_ghz: 2.0, total_time_budget_us: 1_000 };
369        for a in &optimize_allocation(&b, 0.001, 5) { assert!(a.logical_qubits >= 5); }
370    }
371    #[test] fn test_optimize_insufficient() {
372        let b = ResourceBudget { total_physical_qubits: 5, classical_cores: 1, classical_clock_ghz: 1.0, total_time_budget_us: 100 };
373        assert!(optimize_allocation(&b, 0.001, 1).is_empty());
374    }
375    #[test] fn test_optimize_zero_cores() {
376        let b = ResourceBudget { total_physical_qubits: 10_000, classical_cores: 0, classical_clock_ghz: 0.0, total_time_budget_us: 1_000 };
377        assert!(optimize_allocation(&b, 0.001, 1).is_empty());
378    }
379
380    #[test] fn test_latency_budget_d3() {
381        let lb = plan_latency_budget(3, 100);
382        assert_eq!(lb.syndrome_extraction_ns, 360); assert_eq!(lb.decode_ns, 100);
383        assert_eq!(lb.correction_ns, 20); assert_eq!(lb.total_round_ns, 480); assert_eq!(lb.slack_ns, 240);
384    }
385    #[test] fn test_latency_budget_negative_slack() { assert!(plan_latency_budget(3, 1000).slack_ns < 0); }
386    #[test] fn test_latency_budget_scales() {
387        assert!(plan_latency_budget(7, 100).syndrome_extraction_ns > plan_latency_budget(3, 100).syndrome_extraction_ns);
388    }
389
390    #[test] fn test_sim_stable() {
391        let t = simulate_control_loop(&make_loop(5, 0.001, 100), 100, 42);
392        assert_eq!(t.rounds.len(), 100); assert!(t.converged); assert!(t.max_backlog < 50.0);
393    }
394    #[test] fn test_sim_unstable() {
395        let t = simulate_control_loop(&make_loop(3, 0.3, 1000), 200, 42);
396        assert_eq!(t.rounds.len(), 200); assert!(t.max_backlog > 0.0);
397    }
398    #[test] fn test_sim_zero_rounds() {
399        let t = simulate_control_loop(&make_loop(3, 0.001, 100), 0, 42);
400        assert!(t.rounds.is_empty()); assert_eq!(t.final_logical_error_rate, 0.0); assert!(t.converged);
401    }
402    #[test] fn test_sim_deterministic() {
403        let t1 = simulate_control_loop(&make_loop(5, 0.01, 200), 50, 123);
404        let t2 = simulate_control_loop(&make_loop(5, 0.01, 200), 50, 123);
405        for (a, b) in t1.rounds.iter().zip(t2.rounds.iter()) {
406            assert_eq!(a.errors_this_round, b.errors_this_round);
407            assert_eq!(a.errors_corrected, b.errors_corrected);
408        }
409    }
410    #[test] fn test_sim_zero_error_rate() {
411        let t = simulate_control_loop(&make_loop(5, 0.0, 100), 50, 99);
412        assert!(t.converged); assert_eq!(t.final_logical_error_rate, 0.0);
413        for s in &t.rounds { assert_eq!(s.errors_this_round, 0); }
414    }
415    #[test] fn test_sim_snapshot_fields() {
416        let t = simulate_control_loop(&make_loop(3, 0.01, 100), 10, 7);
417        for (i, s) in t.rounds.iter().enumerate() {
418            assert_eq!(s.round, i as u64); assert!(s.errors_corrected <= s.errors_this_round);
419            assert!(s.decode_latency_ns > 0);
420        }
421    }
422
423    #[test] fn test_scaling_uf() { let l = classical_overhead_scaling("union_find"); assert_eq!(l.exponent, 1.0); assert!(l.name.contains("Union-Find")); }
424    #[test] fn test_scaling_mwpm() { assert_eq!(classical_overhead_scaling("mwpm").exponent, 3.0); }
425    #[test] fn test_scaling_neural() { let l = classical_overhead_scaling("neural"); assert_eq!(l.exponent, 1.0); assert!(l.prefactor > 1.0); }
426    #[test] fn test_scaling_unknown() { let l = classical_overhead_scaling("custom"); assert_eq!(l.exponent, 2.0); assert!(l.name.contains("custom")); }
427
428    #[test] fn test_logical_below() { let l = logical_error_scaling(0.001, 0.01); assert!(l.exponent > 0.0); assert_eq!(l.prefactor, 0.1); }
429    #[test] fn test_logical_above() { let l = logical_error_scaling(0.05, 0.01); assert_eq!(l.exponent, 0.0); assert_eq!(l.prefactor, 1.0); }
430    #[test] fn test_logical_at() { assert_eq!(logical_error_scaling(0.01, 0.01).exponent, 0.0); }
431    #[test] fn test_logical_zero_rate() { assert_eq!(logical_error_scaling(0.0, 0.01).exponent, 0.0); }
432    #[test] fn test_logical_zero_thresh() { assert_eq!(logical_error_scaling(0.001, 0.0).exponent, 0.0); }
433}
ruqu_core/control_theory.rs

ruqu_core/
control_theory.rs
