Improvements, fixes, and tests for kTAM model

pull request #4

.github/workflows/ci.yml

@@ -93,11 +93,10 @@ jobs:
       - name: Install maturin and pytest
         run: |
           source venv/bin/activate
-          pip install maturin pytest pytest-cov
+          pip install maturin pytest pytest-cov hypothesis
       # - name: Install X components
       #   run: sudo apt-get install -y gdb libxinerama-dev libxcursor-dev libxfixes-dev libxft-dev libpango1.0-dev libpangoxft-1.0-0 libpangocairo-1.0-0
       - name: Run coverage #           cargo clippy --features python,use_rayon
-        continue-on-error: true
         run: |
           source venv/bin/activate
           cd py-rgrow

.gitignore

@@ -27,4 +27,5 @@ rgrow/target
 Cargo.lock
 .coverage
 coverage.xml
-coverage.lcov
+coverage.lcov
+**/.hypothesis/*

CHANGELOG.md

@@ -1,6 +1,7 @@
 # 0.14.1
 
 - Fix dimer formation rates in kTAM model.
+- Fix k_f specification in kTAM model.
 
 # 0.14.0
 

Cargo.toml

@@ -15,7 +15,7 @@ pyo3-polars = "^0.14"
 polars = {version = "^0.40", features = ["lazy"]}
 
 [workspace.package]
-version = "0.14.0"
+version = "0.14.1"
 authors = ["Constantine Evans <[email protected]>"]
 edition = "2021"
 repository = "https://github.com/cgevans/rgrow"

py-rgrow/pyproject.toml

@@ -3,6 +3,13 @@ name = "rgrow"
 dependencies = ["numpy ~= 1.26", "attrs ~= 23.2", "matplotlib ~= 3.8.2"]
 requires-python = ">=3.9"
 
+[project.optional-dependencies]
+
+testing = [
+    "pytest",
+    "hypothesis"
+]
+
 [build-system]
 requires = ["maturin~=1.4"]
 build-backend = "maturin"

py-rgrow/rgrow/rgrow.pyi

@@ -13,6 +13,8 @@ class TileShape(Enum):
 
 class EvolveOutcome(object): ...
 
+class NeededUpdate(object): ...
+
 class State(object):
     @property
     def canvas_view(self) -> np.ndarray: ...
@@ -67,7 +69,7 @@ class System(object):
     def tile_names(self) -> list[str]: ...
     @property
     def tile_colors(self) -> np.ndarray: ...
-    def update_all(self, state: State, needed) -> None: ...
+    def update_all(self, state: State, needed: NeededUpdate | None = None) -> None: ...
     def plot_canvas(
         self,
         state: State,
@@ -97,6 +99,7 @@ class FFSStateRef(object):
     def time(self) -> float: ...
     def tracking_copy(self) -> np.ndarray: ...
     def clone_state(self) -> State: ...
+    def rate_at_point(self, point: tuple[int, int]) -> float: ...
 
 
 class FFSLevelRef(object):

py-rgrow/tests/test_sanity_ktam.py

@@ -0,0 +1,119 @@
+from rgrow import Tile, TileSet, Bond  # noqa: F841
+import pytest # noqa: F401
+from pytest import approx
+import hypothesis.strategies as st
+from hypothesis import given
+import numpy as np
+
+
+@given(
+    gse=st.floats(4, 12),
+    concs_nM=st.floats(1, 10000),
+    alpha=st.floats(-10, 10),
+    kf=st.floats(1e3, 1e9),
+    stoic=st.floats(0.1, 5),
+    bond_strength=st.floats(0.5, 3)
+)
+def test_basic_rates_ktam(gse, concs_nM, alpha, kf, stoic, bond_strength):
+    gmc = alpha-np.log(concs_nM / 1e9)
+
+    ts = TileSet([
+        Tile(["a", "b", "c", "d"], stoic = stoic),
+        Tile(["e", "d", "g", "b"]),
+        ],
+        [Bond("b", bond_strength)],
+        [("e", "g", 1)],
+        kf=kf,
+        alpha=alpha,
+        gse=gse,
+        gmc=gmc,
+        )
+
+    sys, state = ts.create_system_and_state()
+
+    cv = state.canvas_view
+
+    cv[5,5] = 1
+    cv[5,6] = 2
+    cv[6,6] = 2
+
+    sys.update_all(state)
+
+    dimers = sys.calc_dimers()
+
+    for d in dimers:
+        if (d.t1 == 1) and (d.t2 == 1):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2 * stoic**2)
+        elif (d.t1 == 1) and (d.t2 == 2):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2 * stoic)
+        elif (d.t1 == 2) and (d.t2 == 1):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2 * stoic)
+        elif (d.t1 == 2) and (d.t2 == 2):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2)
+        else:
+            raise ValueError("Unexpected Dimer: ", d)
+
+    assert state.rate_at_point((5,5)) == approx(kf * np.exp(-bond_strength * gse + alpha))
+
+    assert state.rate_at_point((5,6)) == approx(kf * np.exp(-(1+bond_strength) * gse + alpha))
+
+    assert state.rate_at_point((6,6)) == approx(kf * np.exp( -gse + alpha))
+
+    assert state.rate_at_point((6,5)) == approx(kf * stoic * concs_nM / 1e9)
+
+
+@given(
+    gse=st.floats(4, 12),
+    concs_nM=st.floats(1, 10000),
+    alpha=st.floats(-10, 10),
+    kf=st.floats(1e3, 1e9),
+    stoic=st.floats(0.1, 5),
+    bond_strength=st.floats(0.5, 3)
+)
+def test_basic_rates_oldktam(gse, concs_nM, alpha, kf, stoic, bond_strength):
+    gmc = alpha-np.log(concs_nM / 1e9)
+
+    ts = TileSet([
+        Tile(["a", "b", "c", "d"], stoic = stoic),
+        Tile(["e", "d", "g", "b"]),
+        ],
+        [Bond("b", bond_strength)],
+        [("e", "g", 1)],
+        kf=kf,
+        alpha=alpha,
+        gse=gse,
+        gmc=gmc,
+        model='oldktam'
+        )
+
+    sys, state = ts.create_system_and_state()
+
+    cv = state.canvas_view
+
+    cv[5,5] = 1
+    cv[5,6] = 2
+    cv[6,6] = 2
+
+    sys.update_all(state)
+
+    dimers = sys.calc_dimers()
+
+    for d in dimers:
+        if (d.t1 == 1) and (d.t2 == 1):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2 * stoic**2)
+        elif (d.t1 == 1) and (d.t2 == 2):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2 * stoic)
+        elif (d.t1 == 2) and (d.t2 == 1):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2 * stoic)
+        elif (d.t1 == 2) and (d.t2 == 2):
+            assert d.formation_rate == approx(kf * (concs_nM / 1e9)**2)
+        else:
+            raise ValueError("Unexpected Dimer: ", d)
+
+    assert state.rate_at_point((5,5)) == approx(kf * np.exp(-bond_strength * gse + alpha))
+
+    assert state.rate_at_point((5,6)) == approx(kf * np.exp(-(1+bond_strength) * gse + alpha))
+
+    assert state.rate_at_point((6,6)) == approx(kf * np.exp( -gse + alpha))
+
+    assert state.rate_at_point((6,5)) == approx(kf * stoic * concs_nM / 1e9)

rgrow/src/models/ktam.rs

@@ -3,10 +3,11 @@
 ///
 /// Implementation notes:
 ///
-/// - Concentration is stored *in nM*, as this is most likely to result in reasonable values.
+/// - Concentration is stored in M.
 /// - Concentration, rather than stoichiometry and $G_{mc}$ values, are stored internally.
 /// - As with Xgrow, duples are treated as two tiles, where the right/bottom tile is "fake".
 /// - Unlike Xgrow, there is no *requirement* that there be a seed.
+/// - Dimer detachment is not currently implemented.
 use super::ktam_fission::*;
 use crate::{
     base::{GrowError, RgrowError},
@@ -27,10 +28,16 @@ use std::collections::HashMap;
 
 use crate::base::{Glue, Tile};
 
+/// Concentration (M)
 type Conc = f64;
 type Strength = f64;
+
+/// Rate per concentration (M/s)
 type RatePerConc = f64;
+
 type Energy = f64;
+
+/// Rate (1/s)
 type Rate = f64;
 
 trait NonZero {
@@ -45,8 +52,8 @@ impl NonZero for Tile {
 
 const FAKE_EVENT_RATE: f64 = 1e-20;
 
-fn energy_exp_times_u0(x: f64) -> Conc {
-    1.0e9 * x.exp()
+fn energy_exp_times_u0(x: Energy) -> Conc {
+    x.exp()
 }
 
 #[derive(Clone, Debug, Serialize, Deserialize)]
@@ -75,7 +82,7 @@ impl Default for TileShape {
 pub struct KTAM {
     /// Tile names, as strings.  Only used for reference.
     pub tile_names: Vec<String>,
-    /// Tile concentrations, actual (not modified by alpha/Gse/etc) in nM.
+    /// Tile concentrations, actual (not modified by alpha/Gse/etc) in M.
     pub tile_concs: Array1<Conc>,
     /// Glues (by number) on tile edges.
     pub tile_edges: Array2<Glue>,
@@ -85,8 +92,11 @@ pub struct KTAM {
     /// Strengths of links between different glues (eg, glue 1 binding to
     /// glue 2).  Should be symmetric.  Will be added with glue_strengths.
     pub glue_links: Array2<Strength>,
+    /// kTAM $G_{se}$ value (unitless, positive is favorable)
     pub g_se: Energy,
+    /// kTAM $\alpha$ value (unitless, positive is favorable)
     pub alpha: Energy,
+    /// Rate constant for monomer attachment events, in M/s.
     pub kf: RatePerConc,
     pub double_to_right: Array1<Tile>,
     pub double_to_bottom: Array1<Tile>,
@@ -690,25 +700,25 @@ impl SystemWithDimers for KTAM {
 
         for ((t1, t2), e) in self.energy_ns.indexed_iter() {
             if *e > 0. {
-                let biconc = self.tile_concs[t1] * self.tile_concs[t2] / 1.0e9;
+                let biconc = self.tile_concs[t1] * self.tile_concs[t2];
                 dvec.push(DimerInfo {
                     t1: t1 as Tile,
                     t2: t2 as Tile,
                     orientation: Orientation::NS,
-                    formation_rate: self.kf * biconc / 1e9, // FIXME: 1e9 because we're using nM for concs
+                    formation_rate: self.kf * biconc,
                     equilibrium_conc: biconc * f64::exp(*e - self.alpha),
                 });
             }
         }
 
         for ((t1, t2), e) in self.energy_we.indexed_iter() {
             if *e > 0. {
-                let biconc = self.tile_concs[t1] * self.tile_concs[t2] / 1.0e9;
+                let biconc = self.tile_concs[t1] * self.tile_concs[t2];
                 dvec.push(DimerInfo {
                     t1: t1 as Tile,
                     t2: t2 as Tile,
                     orientation: Orientation::WE,
-                    formation_rate: self.kf * biconc / 1e9, // FIXME: 1e9 because we're using nM for concs
+                    formation_rate: self.kf * biconc,
                     equilibrium_conc: biconc * f64::exp(*e - self.alpha),
                 });
             }
@@ -766,7 +776,7 @@ impl KTAM {
             glue_links: Array2::zeros((nglues + 1, nglues + 1)),
             g_se: (9.),
             alpha: (0.),
-            kf: (1e-3),
+            kf: (1e6),
             double_to_right: Array1::zeros(ntiles + 1),
             double_to_bottom: Array1::zeros(ntiles + 1),
             seed: Seed::None(),
@@ -815,29 +825,32 @@ impl KTAM {
         g_se: f64,
         g_mc: f64,
         alpha: Option<f64>,
-        _k_f: Option<f64>,
+        k_f: Option<f64>,
         seed: Option<Seed>,
         fission_handling: Option<FissionHandling>,
-        _chunk_handling: Option<ChunkHandling>,
-        _chunk_size: Option<ChunkSize>,
+        _chunk_handling: Option<ChunkHandling>, // Fixme
+        _chunk_size: Option<ChunkSize>, // Fixme
         tile_names: Option<Vec<String>>,
         tile_colors: Option<Vec<[u8; 4]>>,
     ) -> Self {
         let ntiles = tile_stoics.len() as Tile;
 
-        tile_stoics.map_inplace(|x| *x *= 1.0e9 * (-g_mc + alpha.unwrap_or(0.)).exp());
 
-        let mut ktam = Self::new_sized(tile_stoics.len() as Tile - 1, glue_strengths.len() - 1);
+        let mut ktam = Self::new_sized(ntiles - 1, glue_strengths.len() - 1);
 
-        ktam.tile_concs = tile_stoics;
         ktam.tile_edges = tile_edges;
         ktam.glue_strengths = glue_strengths;
         ktam.g_se = g_se;
         ktam.alpha = alpha.unwrap_or(ktam.alpha);
+        tile_stoics.map_inplace(|x| *x *= (-g_mc + ktam.alpha).exp());
+        ktam.tile_concs = tile_stoics;
         ktam.seed = seed.unwrap_or(ktam.seed);
-        //ktam.tile_colors = tile_colors.unwrap_or(ktam.tile_colors);
         ktam.tile_names = tile_names.unwrap_or(ktam.tile_names);
 
+        ktam.kf = k_f.unwrap_or(ktam.kf);
+
+
+
         ktam.tile_colors = match tile_colors {
             Some(tc) => tc,
             None => {

rgrow/src/python.rs

@@ -2,12 +2,12 @@ use std::ops::DerefMut;
 use std::time::Duration;
 
 use crate::base::{NumEvents, NumTiles, RustAny, Tile};
-use crate::canvas::Canvas;
+use crate::canvas::{Canvas, PointSafeHere};
 use crate::ffs::{BoxedFFSResult, FFSRunConfig, FFSStateRef};
 use crate::ratestore::RateStore;
 use crate::state::{ClonableState, StateEnum, StateStatus, TrackerData};
 use crate::system::{
-    DynSystem, EvolveBounds, EvolveOutcome, NeededUpdate, SystemEnum, TileBondInfo,
+    DimerInfo, DynSystem, EvolveBounds, EvolveOutcome, NeededUpdate, SystemEnum, SystemWithDimers, TileBondInfo
 };
 use crate::tileset::CanvasType;
 use ndarray::Array2;
@@ -48,6 +48,10 @@ impl PyState {
         Ok(PyArray2::from_array_bound(py, &ra))
     }
 
+    pub fn rate_at_point(&self, point: (usize, usize)) -> f64 {
+        self.0.rate_at_point(PointSafeHere(point)) // FIXME: check on bounds
+    }
+
     pub fn tracking_copy(
         this: &Bound<Self>,
     ) -> PyResult<RustAny> {
@@ -235,6 +239,10 @@ impl PySystem {
         }
     }
 
+    fn calc_dimers(&self) -> Vec<DimerInfo> {
+        self.0.calc_dimers()
+    }
+
     fn calc_mismatch_locations<'py>(
         this: &Bound<'py, Self>,
         state: PyStateOrRef,
@@ -292,6 +300,7 @@ impl PySystem {
         Ok(RustAny(self.0.get_param(param_name)?))
     }
 
+    #[pyo3(signature = (state, needed = &NeededUpdate::All))]
     fn update_all(&self, state: &mut PyState, needed: &NeededUpdate) {
         self.0.update_all(&mut state.0, needed)
     }