update: add measure model

results/a.txt

+34542 epochs: error at X
+68448 epochs: error at X
+96469 epochs: error at X
+114284 epochs: error at X
+46796 epochs: error at X
+190354 epochs: error at X
+232556 epochs: error at Z
+119949 epochs: error at X
+317618 epochs: error at X
+322432 epochs: error at X
+365193 epochs: error at X
+118307 epochs: error at Z
+319066 epochs: error at X
+23335 epochs: error at Z
+23788 epochs: error at X
+120609 epochs: error at Z
+360194 epochs: error at X
+753124 epochs: error at X
+546134 epochs: error at X
+671985 epochs: error at X
+2650 epochs: error at X
+257045 epochs: error at Z
+453623 epochs: error at Z
+999008 epochs: error at Z
+1219310 epochs: error at X
+267614 epochs: error at X
+1262107 epochs: error at X
+213224 epochs: error at X
+1096497 epochs: error at X
+198847 epochs: error at X
+15685 epochs: error at X
+1629810 epochs: error at Z
+1208869 epochs: error at X
+153543 epochs: error at X
+588797 epochs: error at X
+416652 epochs: error at Z
+697335 epochs: error at X
+315546 epochs: error at Z
+1598568 epochs: error at X
+33362 epochs: error at X
+264370 epochs: error at X
+1482308 epochs: error at Z
+1178478 epochs: error at X
+402478 epochs: error at X
+50241 epochs: error at X
+455593 epochs: error at X
+1123066 epochs: error at X
+611230 epochs: error at X
+721319 epochs: error at X
+2761424 epochs: error at X
+615442 epochs: error at Z
+941473 epochs: error at X
+1424981 epochs: error at Z
+614046 epochs: error at Z
+313860 epochs: error at X
+198051 epochs: error at X
+37496 epochs: error at X
+1666075 epochs: error at X
+12609 epochs: error at X
+103780 epochs: error at X
+147037 epochs: error at X
+306449 epochs: error at X
+1453734 epochs: error at X
+738236 epochs: error at X
+8596 epochs: error at X
+2099605 epochs: error at Z
+1156852 epochs: error at Z
+1202205 epochs: error at X
+927099 epochs: error at Z
+201906 epochs: error at X
+1651736 epochs: error at Z
+523076 epochs: error at X
+281169 epochs: error at X
+477685 epochs: error at X
+24024 epochs: error at X
+375126 epochs: error at Z
+183999 epochs: error at X
+88460 epochs: error at X
+261486 epochs: error at X
+879009 epochs: error at X
+134891 epochs: error at Z
+1729346 epochs: error at X
+28962 epochs: error at X
+211918 epochs: error at X
+809288 epochs: error at X
+389620 epochs: error at X
+669960 epochs: error at Z
+39611 epochs: error at X
+2152525 epochs: error at X
+753842 epochs: error at X
+782729 epochs: error at X
+3079475 epochs: error at Z
+801006 epochs: error at X
+2069827 epochs: error at X
+1858144 epochs: error at Z
+1981299 epochs: error at X
+1985987 epochs: error at X
+2427930 epochs: error at X
+2587372 epochs: error at X
+3275771 epochs: error at X
+average resist epoch: 732085.87

src/controller.rs

+use crate::model::{SurfaceCodeModel, SurfModelClass};
+use crate::error::{NoiseChannel, ErrorModel, MeasureNoiseModel};
+use crate::utils::epoch;
 
+
+pub fn naive_runner(xn: usize, zm: usize, p: f64) -> usize {
+    let mut model = SurfaceCodeModel::new(xn, zm, SurfModelClass::Surf1);
+    let channel = NoiseChannel::new_depolarizing_channel(p);
+    let error_model = ErrorModel::new(channel, None, 0.0);
+
+    let mut epochs = 0;
+    loop {
+        epochs = epochs + 1;
+        if let Some(stablizer) = epoch(&mut model, &error_model) {
+            println!("{} epochs: error at {}", epochs, stablizer);
+            break;
+        }
+    }
+    epochs
+}
+
+pub fn test_runner1(xn: usize, zm: usize, p: f64) -> usize {
+    let mut model = SurfaceCodeModel::new(xn, zm, SurfModelClass::Surf1);
+    let channel = NoiseChannel::new_depolarizing_channel(p);
+    let pm_total = 0.000;
+    let measure_model = MeasureNoiseModel::new(vec![1.0 - pm_total], vec![pm_total]);
+    let error_model = ErrorModel::new(channel, Some(measure_model), 0.0);
+
+    let mut epochs = 0;
+    loop {
+        epochs = epochs + 1;
+        if let Some(stablizer) = epoch(&mut model, &error_model) {
+            println!("{} epochs: error at {}", epochs, stablizer);
+            break;
+        }
+    }
+    epochs
+}

src/error.rs

@@ -38,24 +38,136 @@ impl NoiseChannel {
     }
 }
 
-#[derive(Debug, Clone, Copy)]
+#[derive(Debug, Clone)]
+pub struct MeasureNoiseModel {
+    pub success: NoiseChannelN,
+    pub failure: NoiseChannelN
+}
+
+impl MeasureNoiseModel {
+    pub fn new(success: Vec<f64>, failure: Vec<f64>) -> MeasureNoiseModel {
+        MeasureNoiseModel {
+            success: NoiseChannelN::new(4, success),
+            failure: NoiseChannelN::new(4, failure)
+        }
+    }
+
+    pub fn is_success(&self, rng: &mut ThreadRng) -> bool {
+        let a = rng.gen_range(0.0..1.0);
+        let p_succ = self.success.prob_sum;
+        return a <= p_succ
+    }
+}
+
+
+#[derive(Debug, Clone)]
 pub struct ErrorModel {
     /// error for time evolution
     pub decay: NoiseChannel,
     /// error model for noise measurement:
-    /// a perfect measurement tailed with a noisy channel
-    pub post_measure: NoiseChannel,
+    pub measure: Option<MeasureNoiseModel>,
     pub mesure_error_rate: f64,
 }
 
 impl ErrorModel {
-    pub fn new(decay: NoiseChannel, post_measure: NoiseChannel, mesure_error_rate: f64) -> Self {
+    pub fn new(decay: NoiseChannel, post_measure: Option<MeasureNoiseModel>, mesure_error_rate: f64) -> Self {
         Self {
             decay,
-            post_measure,
+            measure: post_measure,
             mesure_error_rate,
         }
     }
 }
 
 // TODO channel for multi qubit
+#[derive(Debug, Clone)]
+pub struct NoiseChannelN {
+    n: usize,
+    /// did not sum over different permutations.
+    /// probs in order: IIII, IIIX, IIIY, IIIZ, IIXI, IIXX, ...
+    probs: Vec<f64>,
+    prob_sum: f64,
+}
+
+impl NoiseChannelN {
+    pub fn new(n: usize, probs: Vec<f64>) -> NoiseChannelN {
+        let prob_sum = probs.iter().sum();
+        NoiseChannelN {
+            n,
+            probs,
+            prob_sum,
+        }
+    }
+    pub fn get_rand_op(&self, rng: &mut ThreadRng) -> Vec<u8> {
+        // *qubit = self.apply_with_rng(*qubit, rng);
+        let n: f64 = rng.gen_range(0.0..self.prob_sum);
+        self.get_op(n)
+    }
+    /// determine which kraus operator it belong
+    pub fn get_op(&self, n: f64) -> Vec<u8> {
+        // *qubit = self.apply_with_rng(*qubit, rng);
+        let mut p_sum_l: f64 = 0.0;
+        let mut p_sum_r: f64 = self.probs[0];
+        let mut counter = 0;
+        while counter < self.probs.len() - 1 {
+            if p_sum_l <= n && n < p_sum_r {
+                break;
+            }
+            p_sum_l = p_sum_l + self.probs[counter];
+            p_sum_r = p_sum_r + self.probs[counter + 1];
+            counter = counter + 1;
+        }
+        let digits = to_quanternary(counter, self.n);
+        let mut qubits = vec![0; self.n];
+        for i in 0..qubits.len() {
+            qubits[i] ^= match digits[i] {
+                0 => 0,
+                1 => BIT_X,
+                2 => (BIT_X + BIT_Z),
+                3 => BIT_Z,
+                _ => 0,
+            }
+        }
+        qubits
+    }
+}
+
+fn to_digits(num: usize, d: usize, max_digit: usize) -> Vec<usize> {
+    let mut result = Vec::new();
+    let mut rest = num;
+    while rest > 0 {
+        let a = rest % d;
+        rest = rest / d;
+        result.push(a);
+    }
+
+    if result.len() < max_digit {
+        result.append(&mut vec![0; max_digit - result.len()]);
+    }
+
+    result.reverse();
+    result
+}
+
+fn to_quanternary(num: usize, max_digit: usize) -> Vec<usize> {
+    to_digits(num, 4, max_digit)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn can_convert_n_to_digits() {
+        let n = 14;
+        let d = 4;
+        let result_expected = vec![0, 0, 0, 3, 2];
+        let result = to_digits(n, d, 5);
+        assert_eq!(result_expected, result);
+        let n = 284;
+        let d = 4;
+        let result_expected = vec![1, 0, 1, 3, 0];
+        let result = to_digits(n, d, 3);
+        assert_eq!(result_expected, result);
+    }
+}

src/main.rs

@@ -7,25 +7,7 @@ mod utils;
 use clap::{App, Arg};
 use rayon::prelude::*;
 
-use error::{ErrorModel, NoiseChannel};
-use model::SurfaceCodeModel;
-use utils::epoch;
-
-fn naive_runner(xn: usize, zm: usize, p: f64) -> usize {
-    let mut model = SurfaceCodeModel::new(xn, zm);
-    let channel = NoiseChannel::new_depolarizing_channel(p);
-    let error_model = ErrorModel::new(channel, channel, 0.0);
-
-    let mut epochs = 0;
-    loop {
-        epochs = epochs + 1;
-        if let Some(stablizer) = epoch(&mut model, &error_model) {
-            println!("{} epochs: error at {}", epochs, stablizer);
-            break;
-        }
-    }
-    epochs
-}
+use controller::*;
 
 fn main() {
     std::env::set_var("LD_LIBRARY_PATH", "bin/");
@@ -86,7 +68,7 @@ fn main() {
 
     let result: Vec<usize> = (0..rounds)
         .into_par_iter()
-        .map(|_| naive_runner(xn, zm, p))
+        .map(|_| test_runner1(xn, zm, p))
         .collect();
 
     let average = result.iter().sum::<usize>() as f64 / (rounds as f64);

src/model/codestate.rs

@@ -36,7 +36,7 @@ pub struct CodeState {
     pub measures_z: Vec<bool>,
 
     /// for the sake of efficiency
-    rng: ThreadRng
+    pub rng: ThreadRng
 }
 
 impl CodeState {

src/model/geometry.rs

@@ -4,6 +4,7 @@ use super::codestate::*;
 
 mod coordinate;
 mod grid_geometry;
+mod toric_geometry;
 
 pub use coordinate::Scale;
 pub use grid_geometry::GridGeometry;

src/model/geometry/coordinate.rs

+use crate::model::model::SurfModelClass;
 use crate::model::Stabilizer;
 
 #[derive(Debug, Copy, Clone)]
@@ -6,15 +7,46 @@ pub struct Scale {
     pub zn: usize,
     pub xm: usize,
     pub xn: usize,
+
+    pub class: SurfModelClass,
 }
 
 impl Scale {
-    pub fn new(xn: usize, zm: usize) -> Scale {
+    pub fn new(xn: usize, zm: usize, class: SurfModelClass) -> Scale {
+        let zn = match class {
+            SurfModelClass::Surf1 => xn + 1,
+            SurfModelClass::Toric => xn,
+        };
+        let xm = match class {
+            SurfModelClass::Surf1 => zm + 1,
+            SurfModelClass::Toric => zm,
+        };
         Scale {
             zm,
-            zn: xn + 1,
-            xm: zm + 1,
+            zn,
+            xm,
             xn,
+            class,
+        }
+    }
+
+    pub fn get_data_num(&self) -> usize {
+        match self.class {
+            SurfModelClass::Surf1 => self.zm * (2 * self.xn + 1) + self.xn + 1,
+            SurfModelClass::Toric => 2 * self.zm * self.zn,
+        }
+    }
+
+    pub fn get_measure_num(&self, stabilizer: Stabilizer) -> usize {
+        match self.class {
+            SurfModelClass::Surf1 => match stabilizer {
+                Stabilizer::X => self.xn * self.xm,
+                Stabilizer::Z => self.zn * self.zm,
+            },
+            SurfModelClass::Toric => {
+                // for toric code, it's all the same
+                self.zn * self.zm
+            }
         }
     }
 }
@@ -103,7 +135,7 @@ mod tests {
 
     #[test]
     fn from_data_index_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = Scale::new(5, 5, SurfModelClass::Surf1);
         let data_index_1 = 15;
         let data_coordinate_expected_1 = Coordinate::new(8, 2);
         assert_eq!(
@@ -133,7 +165,7 @@ mod tests {
 
     #[test]
     fn from_measure_index_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = Scale::new(5, 5, SurfModelClass::Surf1);
         let measure_index_1 = 8;
         let measure_coordinate_expected_1 = Coordinate::new(4, 3);
         assert_eq!(
@@ -150,7 +182,7 @@ mod tests {
 
     #[test]
     fn to_data_index_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = Scale::new(5, 5, SurfModelClass::Surf1);
         let coordinate_1 = Coordinate::new(8, 2);
         let data_index_expected_1 = 15;
         assert_eq!(coordinate_1.to_data_index(&scale), data_index_expected_1);
@@ -166,7 +198,7 @@ mod tests {
     }
     #[test]
     fn to_measure_index_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = Scale::new(5, 5, SurfModelClass::Surf1);
 
         let coordinate_1 = Coordinate::new(4, 3);
         let measure_index_expected_1 = 8;

src/model/geometry/grid_geometry.rs

@@ -220,6 +220,10 @@ impl GridGeometry {
 mod tests {
     use super::*;
 
+    fn gen_scale(xn: usize, zm: usize) -> Scale {
+        Scale::new(xn, zm, crate::model::SurfModelClass::Surf1)
+    }
+
     fn vecs_eq(vec1: &[usize], vec2: &[usize]) -> bool {
         for ele in vec1.into_iter() {
             if !vec2.contains(ele) {
@@ -236,7 +240,7 @@ mod tests {
 
     #[test]
     fn get_measureq_datas_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = gen_scale(5, 5);
         let geometry = GridGeometry::new(Stabilizer::X, scale);
         let measure = 4;
         let datas_expected = vec![4, 9, 10, 15];
@@ -247,7 +251,7 @@ mod tests {
 
     #[test]
     fn get_dataq_measures_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = gen_scale(5, 5);
         let geometry = GridGeometry::new(Stabilizer::X, scale);
         let data = 16;
         let measures_expected = vec![5, 11];
@@ -258,7 +262,7 @@ mod tests {
 
     #[test]
     fn get_path_between_measures_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = gen_scale(5, 5);
         let geometry = GridGeometry::new(Stabilizer::X, scale);
         let c1 = Coordinate::new(8, 3);
         let c2 = Coordinate::new(2, 7);
@@ -272,7 +276,7 @@ mod tests {
 
     #[test]
     fn get_path_to_boundary_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = gen_scale(5, 5);
         let geometry = GridGeometry::new(Stabilizer::X, scale);
         {
             let c = Coordinate::new(8, 3);
@@ -307,7 +311,7 @@ mod tests {
 
     #[test]
     fn get_weight_boundary_is_correct() {
-        let scale = Scale::new(5, 5);
+        let scale = gen_scale(5, 5);
         let geometry = GridGeometry::new(Stabilizer::X, scale);
         {
             let c = Coordinate::new(8, 3);

src/model/geometry/toric_geometry.rs

+use petgraph::{Graph, Undirected};
+
+use crate::model::Stabilizer;
+
+use super::coordinate::Coordinate;
+use super::{Geometry, Scale};
+
+/// x ⋅ x ⋅ x ⋅ x · x · -
+/// ⋅ z ⋅ z ⋅ z ⋅ z · z
+/// x ⋅ x ⋅ x ⋅ x · x ·
+/// ⋅ z ⋅ z ⋅ z ⋅ z · z
+/// x ⋅ x ⋅ x ⋅ x · x ·
+/// ⋅ z ⋅ z ⋅ z ⋅ z · z N
+/// x ⋅ x ⋅ x ⋅ x · x ·
+/// ⋅ z ⋅ z ⋅ z ⋅ z · z
+/// x ⋅ x ⋅ x ⋅ x · x ·
+/// ⋅ z ⋅ z ⋅ z ⋅ z · z -
+/// |         M       |
+#[derive(Debug)]
+pub struct ToricGeometry {
+    pub stabilizer: Stabilizer,
+    pub scale: Scale,
+    m: usize,
+    n: usize,
+}
+
+fn abs_sub(op1: usize, op2: usize) -> usize {
+    if op1 < op2 {
+        op2 - op1
+    } else {
+        op1 - op2
+    }
+}
+
+impl Geometry for ToricGeometry {
+    fn new(stabilizer: Stabilizer, scale: Scale) -> Self {
+        let m = scale.xm;
+        let n = scale.xn;
+        ToricGeometry {
+            stabilizer,
+            scale,
+            m,
+            n,
+        }
+    }
+
+    fn get_measureq_datas(&self, measure: usize) -> Vec<usize> {
+        let cm = Coordinate::from_measure_index(self.stabilizer, measure, &self.scale);
+
+        let x: isize = cm.x as isize;
+        let y: isize = cm.y as isize;
+        let mut adjacant_datas: Vec<usize> = Vec::new();
+        if let Some(c) = Coordinate::new_restricted(x - 1, y, &self.scale) {
+            adjacant_datas.push(c.to_data_index(&self.scale));
+        }
+        if let Some(c) = Coordinate::new_restricted(x + 1, y, &self.scale) {
+            adjacant_datas.push(c.to_data_index(&self.scale));
+        }
+        if let Some(c) = Coordinate::new_restricted(x, y - 1, &self.scale) {
+            adjacant_datas.push(c.to_data_index(&self.scale));
+        }
+        if let Some(c) = Coordinate::new_restricted(x, y + 1, &self.scale) {
+            adjacant_datas.push(c.to_data_index(&self.scale));
+        }
+        adjacant_datas
+    }
+
+    fn get_dataq_measures(&self, data: usize) -> Vec<usize> {
+        let cd = Coordinate::from_data_index(data, &self.scale);
+
+        let x: isize = cd.x as isize;
+        let y: isize = cd.y as isize;
+        let mut adjacant_measures: Vec<usize> = Vec::new();
+        if let Some(c) = Coordinate::new_restricted(x - 1, y, &self.scale) {
+            if let Some(mesure_index) = c.to_measure_index(self.stabilizer, &self.scale) {
+                adjacant_measures.push(mesure_index);
+            }
+        }
+        if let Some(c) = Coordinate::new_restricted(x + 1, y, &self.scale) {
+            if let Some(mesure_index) = c.to_measure_index(self.stabilizer, &self.scale) {
+                adjacant_measures.push(mesure_index);
+            }
+        }
+        if let Some(c) = Coordinate::new_restricted(x, y - 1, &self.scale) {
+            if let Some(mesure_index) = c.to_measure_index(self.stabilizer, &self.scale) {
+                adjacant_measures.push(mesure_index);
+            }
+        }
+        if let Some(c) = Coordinate::new_restricted(x, y + 1, &self.scale) {
+            if let Some(mesure_index) = c.to_measure_index(self.stabilizer, &self.scale) {
+                adjacant_measures.push(mesure_index);
+            }
+        }
+        adjacant_measures
+    }
+
+    fn get_path_between_measures(&self, measure1: usize, measure2: usize) -> Vec<usize> {
+        if measure1 == measure2 {
+            return self.get_path_to_boundary(measure1);
+        }
+        let cm1 = Coordinate::from_measure_index(self.stabilizer, measure1, &self.scale);
+        let cm2 = Coordinate::from_measure_index(self.stabilizer, measure2, &self.scale);
+
+        let mut path: Vec<usize> = Vec::new();
+        let mut i = cm1.x;
+        let xl = cm1.x < cm2.x;
+        while i != cm2.x {
+            let x = if xl { i + 1 } else { i - 1 };
+            let c = Coordinate::new(x, cm1.y);
+            let data_index = c.to_data_index(&self.scale);
+            path.push(data_index);
+
+            i = if xl { i + 2 } else { i - 2 }
+        }
+        i = cm1.y;
+        let yl = cm1.y < cm2.y;
+        while i != cm2.y {
+            let y = if yl { i + 1 } else { i - 1 };
+            let c = Coordinate::new(cm2.x, y);
+            let data_index = c.to_data_index(&self.scale);
+            path.push(data_index);
+
+            i = if yl { i + 2 } else { i - 2 }
+        }
+        path
+    }
+
+    fn get_weight_between_measures(&self, measure1: usize, measure2: usize) -> usize {
+        if measure1 == measure2 {
+            return self.get_weight_to_boundary(measure1);
+        }
+        let cm1 = Coordinate