A benchmark and slight optimization for SDC

.gitignore

@@ -28,4 +28,5 @@ Cargo.lock
 .coverage
 coverage.xml
 coverage.lcov
-**/.hypothesis/*
+**/.hypothesis/*
+rgrow/profile.json

Cargo.toml

@@ -25,4 +25,8 @@ categories = ["science", "simulation"]
 
 [profile.release]
 # debug = true
-# lto = true
+# lto = true
+
+[profile.profiling]
+inherits = "release"
+debug = true

rgrow/Cargo.toml

@@ -18,6 +18,10 @@ criterion = { version = "0.5", features = ["html_reports"] }
 name = "sierpinski"
 harness = false
 
+[[bench]]
+name = "sdc"
+harness = false
+
 # [[bench]]
 # name = "ratestore"
 # harness = false

rgrow/benches/sdc.rs

@@ -0,0 +1,97 @@
+use std::collections::HashMap;
+
+use rgrow::{canvas::CanvasSquare, models::sdc1d::{RefOrPair, SDCParams, SDC}, state::{NullStateTracker, QuadTreeState, StateStatus, StateWithCreate}, system::{DynSystem, EvolveBounds, NeededUpdate, System}};
+use criterion::{black_box, criterion_group, criterion_main, Criterion};
+
+fn make_params(c: &mut Criterion) -> SDCParams {
+    let mut scaffold = vec![None; 2];
+    scaffold.extend((1..20).map(|x| Some(format!("s{}*", x))));
+    scaffold.extend((0..2).map(|x| None));
+
+    let mut glue_dg_s = HashMap::new();
+    for s in 1..20 {
+        glue_dg_s.insert(RefOrPair::Ref(format!("s{}", s)), (-27.3, -0.486));
+    }
+
+    for s in "abcdefghijklmnopqrst".chars() {
+        glue_dg_s.insert(RefOrPair::Ref(format!("{}", s)), (-15.36, -0.2366));
+    }
+
+    for s in "jklmnopqrst".chars() {
+        glue_dg_s.insert(RefOrPair::Ref(format!("e{}", s)), (-15.36, -0.2366));
+    }
+
+    let params = SDCParams {
+        tile_glues: [
+            ["a", "s1", "b"],
+            ["b*", "s2", "c"],
+            ["c*", "s3", "d"],
+            ["d*", "s4", "e"],
+            ["e*", "s5", "f"],
+            ["f*", "s6", "g"],
+            ["g*", "s7", "h"],
+            ["h*", "s8", "i"],
+            ["i*", "s9", "j"],
+            ["j*", "s10", "k"],
+            ["k*", "s11", "l"],
+            ["l*", "s12", "m"],
+            ["m*", "s13", "n"],
+            ["n*", "s14", "o"],
+            ["o*", "s15", "p"],
+            ["p*", "s16", "q"],
+            ["q*", "s17", "r"],
+            ["r*", "s18", "s"],
+            ["s*", "s19", "t"],
+            ["ej*", "s10", "ek"],
+            ["ek*", "s11", "el"],
+            ["el*", "s12", "em"],
+            ["em*", "s13", "en"],
+            ["en*", "s14", "eo"],
+            ["eo*", "s15", "ep"],
+            ["ep*", "s16", "eq"],
+            ["eq*", "s17", "er"],
+            ["er*", "s18", "es"],
+            ["es*", "s19", "et"],
+        ].map(|x| x.map(|y| Some(y.to_string())).to_vec()).to_vec(),
+        tile_concentration: vec![1000.0e-9; 30],
+        tile_names: (1..30).map(|x| Some(format!("tile{}", x))).collect(),
+        tile_colors: (1..30).map(|x| None).collect(),
+        scaffold: rgrow::models::sdc1d::SingleOrMultiScaffold::Single(scaffold),
+        glue_dg_s,
+        k_f: 1e6,
+        k_n: 1e5,
+        k_c: 1e4,
+        temperature: 70.0,
+    };
+
+    params
+}
+
+fn bench_hold(c: &mut Criterion) {
+    let params = make_params(c);
+    let mut sdc = SDC::from_params(params);
+
+    let mut state = QuadTreeState::<CanvasSquare, NullStateTracker>::empty((64, 24)).unwrap();
+    let bounds = EvolveBounds::default().for_events(100);
+
+    System::update_all(&mut sdc, &mut state, &NeededUpdate::All);
+
+    // c.bench_function("evolve100_sdc_at_temp", |b| b.iter(|| {
+    //     System::evolve(&sdc, &mut state, bounds).unwrap();
+    //     assert!(state.n_tiles() > 0);
+    // }));
+
+    c.bench_function("unistep_sdc_at_temp", |b| b.iter(|| {
+        System::take_single_step(&sdc, &mut state, 1e6);
+    }));
+
+
+    c.bench_function("update_sdc_temps", |b| b.iter(|| {
+        sdc.change_temperature_to(50.0);
+        System::update_all(&mut sdc, &mut state, &NeededUpdate::NonZero);
+    }));
+
+}
+
+criterion_group!(benches, bench_hold);
+criterion_main!(benches);

rgrow/src/models/sdc1d.rs

@@ -196,21 +196,36 @@ impl SDC {
     }
 
     fn update_monomer_point<S: State>(&self, state: &mut S, scaffold_point: &PointSafe2) {
-        let points = [
-            state.move_sa_w(*scaffold_point),
-            state.move_sa_e(*scaffold_point),
-            PointSafeHere(scaffold_point.0),
-        ]
-        .map(|point| (point, self.event_rate_at_point(state, point)));
-
-        state.update_multiple(&points);
+        let here = PointSafeHere(scaffold_point.0);
+        let w = state.move_sa_w(*scaffold_point);
+        let e = state.move_sa_e(*scaffold_point);
+
+        match (state.v_sh(w), state.v_sh(e)) {
+            (0, 0) => state.update_point(here, self.event_rate_at_point(state, here)),
+            (0, _) => {
+                state.update_multiple(
+                    &[here, e].map(|point| (point, self.event_rate_at_point(state, point))),
+                );
+            }
+            (_, 0) => {
+                state.update_multiple(
+                    &[here, w].map(|point| (point, self.event_rate_at_point(state, point))),
+                );
+            }
+            (_, _) => {
+                state.update_multiple(
+                    &[here, w, e].map(|point| (point, self.event_rate_at_point(state, point))),
+                );
+            }
+        }
     }
 
     /// Fill the energy_bonds array
     fn fill_energy_array(&mut self) {
         let num_of_strands = self.strand_names.len();
         // For each *possible* pair of strands, calculate the energy bond
-        for strand_f in 1..(num_of_strands as usize) { // 1: no point in calculating for 0
+        for strand_f in 1..(num_of_strands as usize) {
+            // 1: no point in calculating for 0
             let (f_west_glue, f_btm_glue, f_east_glue) = {
                 let glues = self.glues.row(strand_f);
                 (
@@ -244,8 +259,11 @@ impl SDC {
                 continue;
             }
 
-            let b_inverse = if f_btm_glue % 2 == 1 { f_btm_glue + 1 } else { f_btm_glue - 1 };
-
+            let b_inverse = if f_btm_glue % 2 == 1 {
+                f_btm_glue + 1
+            } else {
+                f_btm_glue - 1
+            };
 
             // Calculate the binding strength of the strand with the scaffold
             self.scaffold_energy_bonds[strand_f] = self.glue_links[(f_btm_glue, b_inverse)];