diff --git a/CHANGELOG.md b/CHANGELOG.md index 14adae7..b0ad191 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # Changelog +## 0.23.0 +Huge update to remove useless stuff and keep only the content relevant for the PhD thesis. +Remove `abc` in rust, do it in python. + ## 0.22.0 change the way we save the data at the end of the simulations to simplify subsampling. diff --git a/Cargo.lock b/Cargo.lock index ad863f8..9954b31 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -400,7 +400,7 @@ dependencies = [ [[package]] name = "ecdna-evo" -version = "0.22.0" +version = "0.23.0" dependencies = [ "anyhow", "chrono", @@ -942,9 +942,9 @@ dependencies = [ [[package]] name = "sosa" -version = "2.0.0" +version = "3.0.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bf505cb5c17fadbdf3303e046049c6c03a89f7fb956696fd852f554b6a4fe6df" +checksum = "9703c29f4ff2a2a5d53ee1d6360b17e569a8315812d96b8b9271d7ac064a3ca1" dependencies = [ "anyhow", "rand", diff --git a/Cargo.toml b/Cargo.toml index 6d19ef5..ef9eac4 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "ecdna-evo" -version = "0.22.0" +version = "0.23.0" edition = "2021" repository = "https://github.com/fraterenz/ecdna-evo" description = "Evolutionary models of extra-chromosomal DNA (ecDNA)" @@ -19,7 +19,7 @@ ecdna-lib = "3.0.2" indicatif = { version = "*", features = ["rayon"] } serde = { version = "1.0.163", features = ["derive"] } serde_json = "1.0" -sosa = "2" +sosa = "3.0.3" rand = { version = "0.8.0", features = ["small_rng"] } rand_chacha = "0.3.1" rand_distr = "0.4.0" diff --git a/src/app.rs b/src/app.rs deleted file mode 100644 index 4fdfddd..0000000 --- a/src/app.rs +++ /dev/null @@ -1,156 +0,0 @@ -use anyhow::Context; -use ecdna_evo::{distribution::SamplingStrategy, RandomSampling, ToFile}; -use ecdna_lib::distribution::EcDNADistribution; -use rand::SeedableRng; -use rand_chacha::ChaCha8Rng; -use sosa::{ - simulate, AdvanceStep, CurrentState, NbIndividuals, Options, ReactionRates, -}; -use std::{ - fmt::Debug, - path::{Path, PathBuf}, -}; - -pub fn save_distribution( - path2dir: &Path, - distribution: &EcDNADistribution, - rates: &[f32], - id: usize, - verbosity: u8, -) -> anyhow::Result<()> { - let mut filename = match rates { - [b0, b1] => format!("{}b0_{}b1_{}d0_{}d1_{}idx", b0, b1, 0, 0, id), - [b0, b1, d0, d1] => { - format!("{}b0_{}b1_{}d0_{}d1_{}idx", b0, b1, d0, d1, id) - } - _ => unreachable!(), - }; - filename = filename.replace('.', "dot"); - distribution.save( - &path2dir - .join(format!( - "{}cells/ecdna/{}", - distribution.compute_nplus() + distribution.get_nminus(), - filename - )) - .with_extension("json"), - verbosity, - ) -} - -pub struct Dynamics { - pub seed: u64, - pub path2dir: PathBuf, - pub max_cells: NbIndividuals, - pub options: Options, - pub save_before_subsampling: bool, -} - -pub struct Sampling { - pub at: Vec, - pub strategy: SamplingStrategy, -} - -impl Dynamics { - pub fn run( - &self, - idx: usize, - mut process: P, - initial_state: &mut CurrentState, - rates: &ReactionRates, - possible_reactions: &[REACTION; NB_REACTIONS], - sampling: &Option, - ) -> anyhow::Result<()> - where - P: AdvanceStep - + Clone - + Debug - + ToFile - + RandomSampling, - REACTION: std::fmt::Debug, - { - if self.options.verbosity > 1 { - println!("{:#?}", process); - } - - let stream = idx as u64; - let mut rng = ChaCha8Rng::seed_from_u64(self.seed); - rng.set_stream(stream); - - let stop_reason = simulate( - initial_state, - rates, - possible_reactions, - &mut process, - &self.options, - &mut rng, - ); - - if self.options.verbosity > 0 { - println!("stop reason: {:#?}", stop_reason); - } - - if let Some(sampling) = sampling { - for (i, sample_at) in sampling.at.iter().enumerate() { - // happens with birth-death processes - let no_individuals_left = - initial_state.population.iter().sum::() == 0u64; - if no_individuals_left { - if self.options.verbosity > 0 { - println!("do not subsample, early stopping, due to no cells left"); - } - break; - } - if self.options.verbosity > 0 { - println!( - "subsample {}th timepoint with {} cells from the ecDNA distribution with {:#?} cells", - i, sample_at, process - ); - } - if self.save_before_subsampling { - process - .save( - &self.path2dir.join(format!( - "{}cells/before_subsampling", - sample_at - )), - rates, - idx, - ) - .with_context(|| { - format!( - "Cannot save run {} before subsampling", - idx - ) - }) - .unwrap(); - } - - if self.options.verbosity > 0 { - println!("Subsampling"); - } - let distribution = process - .random_sample(*sample_at, &mut rng) - .with_context(|| { - format!( - "wrong number of cells to subsample? {}", - sample_at - ) - }) - .unwrap(); - save_distribution( - &self.path2dir, - &distribution, - rates.0.as_ref(), - idx, - self.options.verbosity, - ) - .with_context(|| "cannot save the ecDNA distribution") - .unwrap(); - } - } else if self.options.verbosity > 1 { - println!("Nosubsampling"); - }; - Ok(()) - } -} diff --git a/src/clap_app.rs b/src/clap_app.rs index fcb588d..dbb915a 100644 --- a/src/clap_app.rs +++ b/src/clap_app.rs @@ -1,32 +1,34 @@ -use anyhow::Context; +use anyhow::{ensure, Context}; use clap::{ArgAction, Parser, ValueEnum}; use ecdna_evo::{ - distribution::{EcDNADistribution, SamplingStrategy, Sigma}, + distribution::EcDNADistribution, segregation::{ Binomial, BinomialNoNminus, BinomialNoUneven, Deterministic, Segregate, }, + Snapshot, }; use rand::Rng; -use sosa::{NbIndividuals, Options}; -use std::{collections::HashMap, path::PathBuf}; - -use crate::{ - app::{Dynamics, Sampling}, - SimulationOptions, MAX_ITER, +use sosa::{IterTime, NbIndividuals, Options}; +use std::{ + collections::{HashMap, VecDeque}, + path::PathBuf, }; +use crate::{SimulationOptions, MAX_ITER}; + +#[derive(Debug)] pub enum Parallel { False, True, Debug, } -#[derive(Debug, Parser)] // requires `derive` feature +#[derive(Debug, Parser)] #[command(name = "Dynamics")] #[command( version, - about = "Mathematical modelling of the ecDNA dynamics", // TODO - long_about = "Study the effect of the random segregation on the ecDNA dynamics using a stochastic simulation algorithm (SSA) aka Gillespie algorithm" + about = "Agent-based modelling of the ecDNA dynamics", + long_about = "Study the effect of the random segregation and positive selection on the ecDNA dynamics using a stochastic simulation algorithm (SSA) aka Gillespie algorithm" )] pub struct Cli { /// The ecDNA segregation type @@ -51,9 +53,9 @@ pub struct Cli { /// process. #[arg(long, value_name = "RATE")] d1: Option, - /// Number of cells to simulate - #[arg(long, short, default_value_t = 100000, conflicts_with = "debug")] - cells: NbIndividuals, + /// Number of years to simulate + #[arg(long, short, default_value_t = 5, conflicts_with = "debug")] + years: NbIndividuals, /// Seed for reproducibility #[arg(long, default_value_t = 26)] seed: u64, @@ -63,40 +65,6 @@ pub struct Cli { /// Run sequentially each run instead of using rayon for parallelisation #[arg(short, long, action = ArgAction::SetTrue, default_value_t = false, conflicts_with = "debug")] sequential: bool, - /// Whether to track over simulations the evolution of the cells with and - /// without any ecDNAs. - #[arg(short, long, action = ArgAction::SetTrue, default_value_t = false)] - nplus_nminus: bool, - /// Whether to track over simulations the evolution of the mean ecDNA - /// copies in the tumour population - #[arg(short, long, action = ArgAction::SetTrue, default_value_t = false)] - mean: bool, - /// Whether to track over simulations the evolution of the variance ecDNA - /// copies in the tumour population - #[arg(long, action = ArgAction::SetTrue, default_value_t = false)] - variance: bool, - /// Whether to track over simulations the evolution of the entropy ecDNA - /// copies in the tumour population - #[arg(long, action = ArgAction::SetTrue, default_value_t = false)] - entropy: bool, - /// Save the data also just before subsampling - #[arg(long, action = ArgAction::SetTrue, default_value_t = false, requires = "sample")] - sample_save_before: bool, - /// The number of cells kept after subsampling. - /// - /// To take the full distribution pass 0 as argument, e.g. when we want to - /// sample with a custom strategy the full ecDNA distribution (see - /// `--sampling-strategy`) - #[arg(long, num_args = 0.., value_name = "CELLS")] - sample: Option>, - #[arg( - long, - value_enum, - requires = "sample", - required = false, - default_value_t - )] - sample_strategy: SamplingStrategyArg, /// Path to store the results of the simulations #[arg(value_name = "DIR")] path: PathBuf, @@ -118,14 +86,123 @@ pub struct Cli { #[arg(short, long, default_value_t = 12, conflicts_with = "debug")] /// Number of independent realisations to simulate of the same stochastic process runs: usize, + #[arg(long, requires = "subsamples", value_delimiter = ',', require_equals = true, num_args = 0..)] + /// Snapshots to take to save the simulation, requires `subsamples`. + /// + /// The combination of `snapshots` with `subsamples` gives four different + /// behaviours: + /// + /// 1. when `snapshots.len() = 1` and `subsamples.len() = 1`: subsample once with the number of cells corresponding to `snapshots[0]` + /// + /// 2. when `snapshots.len() > 1` and `subsamples.len() = 1`: for every `s` in `snapshots`, subsample with the number of cells corresponding to `snapshots[0]` + /// + /// 3. when `snapshots.len() = 1` and `subsamples.len() > 1`: for every `c` in `subsamples`, subsample once with the number of cells corresponding to `c` + /// + /// 4. when `snapshots.len() > 1` and `subsamples.len() > 1`: for every pair `(s, c)` in `snapshots.zip(subsamples)`, subsample at time `s` with `c` cells + snapshots: Option>, + /// Number of cells to subsample before saving the measurements, leave + /// empty when no subsample is needed. + /// If subsampling is performed, the measurements of the whole population + /// will also be saved. + /// + /// See help for `snapshots` for more details. + #[arg(long, requires = "snapshots", num_args = 0.., value_delimiter = ',', require_equals = true)] + subsamples: Option>, #[arg(short, long, action = clap::ArgAction::Count, conflicts_with = "debug", default_value_t = 0)] - verbose: u8, + verbosity: u8, +} + +fn build_snapshots_from_time(n_snapshots: usize, time: f32) -> Vec { + let dx = time / ((n_snapshots - 1) as f32); + let mut x = vec![0.; n_snapshots]; + for i in 1..n_snapshots - 1 { + x[i] = x[i - 1] + dx; + } + + x.shrink_to_fit(); + x[n_snapshots - 1] = time; + x } impl Cli { - pub fn build() -> SimulationOptions { + pub fn build() -> anyhow::Result { let cli = Cli::parse(); + let (years, verbosity, parallel, runs) = if cli.debug { + (2, u8::MAX, Parallel::Debug, 1) + } else if cli.sequential { + (cli.years, cli.verbosity, Parallel::False, cli.runs) + } else { + (cli.years, cli.verbosity, Parallel::True, cli.runs) + }; + + let mut snapshots = match (cli.subsamples, cli.snapshots) { + (Some(sub), Some(snap)) => { + match (&sub[..], &snap[..]) { + // subsample `unique_sub` once at `unique_snap` time + ([unique_sub], [unique_snap]) => VecDeque::from_iter( + [(unique_sub, unique_snap)].into_iter().map( + |(&cells2sample, &time)| Snapshot { + cells2sample, + time, + }, + ), + ), + // subsample `unique_sub` at several `snaps` times + ([unique_sub], snaps) => VecDeque::from_iter( + snaps.iter().zip(std::iter::repeat(unique_sub)).map( + |(&time, &cells2sample)| Snapshot { + cells2sample, + time, + }, + ), + ), + // subsample with different cells `unique_sub` once at `unique_snap` time + (subs, [unique_snap]) => VecDeque::from_iter( + subs.iter().zip(std::iter::repeat(unique_snap)).map( + |(&cells2sample, &time)| Snapshot { + cells2sample, + time, + }, + ), + ), + // subsample with many cells `subs` at specific `snaps` + (subs, snaps) => { + ensure!( + subs.len() == snaps.len(), + "the lenght of snapshots do not match the lenght of subsamples" + ); + VecDeque::from_iter(subs.iter().zip(snaps).map( + |(&cells2sample, &time)| Snapshot { + cells2sample, + time, + }, + )) + } + } + } + (None, None) => VecDeque::from_iter( + build_snapshots_from_time( + 10usize, + if years > 1 { years as f32 - 1. } else { 0.9 }, + ) + .into_iter() + .map(|t| Snapshot { cells2sample: 10, time: t }), + ), + _ => unreachable!(), + }; + + snapshots.make_contiguous(); + snapshots + .as_mut_slices() + .0 + .sort_by(|s1, s2| s1.time.partial_cmp(&s2.time).unwrap()); + + ensure!( + snapshots.iter().all(|s| s.time < cli.years as f32), + "times to take `snapshots` must be smaller than total `years`" + ); + let segregation = cli.segregation.into(); // if both d0, d1 are either unset or equal to 0, pure birth, @@ -140,120 +217,42 @@ impl Cli { None => (false, 0f32), }; let is_birth_death = is_birth_death_d0 | is_birth_death_d1; - let cells = if cli.debug { 10 } else { cli.cells }; - let iterations = if is_birth_death { - MAX_ITER * cells as usize - } else { - cells as usize - }; // if no initial distribution, start with 1 cell with 1 ecDNA copy let distribution = match cli.initial { - Some(path) => EcDNADistribution::load(&path, iterations) - .with_context(|| { - format!( - "Cannot load the ecDNA distribution from {:#?}", - cli.path - ) - }) - .unwrap(), + Some(path) => EcDNADistribution::load( + &path, MAX_ITER, // TODO double check this + ) + .with_context(|| { + format!( + "Cannot load the ecDNA distribution from {:#?}", + cli.path + ) + }) + .unwrap(), None => EcDNADistribution::new( HashMap::::from([(1, 1)]), - iterations, + MAX_ITER, // TODO double check this ), }; - let verbose = if cli.debug { u8::MAX } else { cli.verbose }; - - let (parallel, runs) = if cli.debug { - (Parallel::Debug, 1) - } else if cli.sequential { - (Parallel::False, cli.runs) - } else { - (Parallel::True, cli.runs) - }; let process_type = { - match is_birth_death { - true => match cli.mean { - true => { - if cli.variance { - if cli.entropy { - ProcessType::BirthDeath( - BirthDeathType::MeanVarianceEntropy, - ) - } else { - ProcessType::BirthDeath( - BirthDeathType::MeanVariance, - ) - } - } else { - ProcessType::BirthDeath(BirthDeathType::Mean) - } - } - false => { - if cli.nplus_nminus { - ProcessType::BirthDeath( - BirthDeathType::NMinusNPlus, - ) - } else { - ProcessType::BirthDeath(BirthDeathType::BirthDeath) - } - } - }, - false => match cli.mean { - true => ProcessType::PureBirth(PureBirthType::Mean), - false => { - if cli.nplus_nminus { - ProcessType::PureBirth(PureBirthType::NMinusNPlus) - } else { - ProcessType::PureBirth(PureBirthType::PureBirth) - } - } - }, + if is_birth_death { + ProcessType::BirthDeath + } else { + ProcessType::PureBirth } }; - let sampling = if let Some(sampling_at) = cli.sample { - let sampling_strategy = match cli.sample_strategy { - SamplingStrategyArg::Uniform => SamplingStrategy::Uniform, - SamplingStrategyArg::Gaussian => { - SamplingStrategy::Gaussian(Sigma::new(1.)) - } - }; - let (at, strategy) = { - let strategy = sampling_strategy; - // 0 means take the full ecDNA distribution which does not make - // sense with a uniform distribution - if sampling_at.len() == 1 && sampling_at[0] == 0 { - assert_ne!( - strategy, - SamplingStrategy::Uniform, - "It doesn't make sense to take the full ecDNA distribution as sample with an uniform sampling strategy" - ); - (vec![cells], strategy) - } else { - (sampling_at, strategy) - } - }; - Some(Sampling { at, strategy }) - } else { - assert!(cli.sample.is_none()); - None - }; - - SimulationOptions { - simulation: Dynamics { - seed: cli.seed, - max_cells: cells, - options: Options { - max_iter: iterations, - init_iter: 0, - max_cells: cells, - verbosity: verbose, - }, - path2dir: cli.path, - save_before_subsampling: cli.sample_save_before, + let options = SimulationOptions { + seed: cli.seed, + options: Options { + max_iter_time: IterTime { iter: MAX_ITER, time: years as f32 }, + init_iter: 0, + max_cells: 1_000_000, + verbosity, }, + path2dir: cli.path, process_type, b0: cli.b0, b1: cli.b1, @@ -262,10 +261,15 @@ impl Cli { distribution, parallel, segregation, - sampling, + snapshots, growth: cli.growth, runs, + }; + if verbosity > 0 { + println!("running sims with {:#?}", options); } + + Ok(options) } } @@ -348,26 +352,10 @@ impl std::fmt::Display for GrowthOptions { } } -#[derive(Clone, Copy)] +#[derive(Clone, Copy, Debug)] pub enum ProcessType { - PureBirth(PureBirthType), - BirthDeath(BirthDeathType), -} - -#[derive(Clone, Copy)] -pub enum PureBirthType { PureBirth, - NMinusNPlus, - Mean, -} - -#[derive(Clone, Copy)] -pub enum BirthDeathType { BirthDeath, - NMinusNPlus, - Mean, - MeanVariance, - MeanVarianceEntropy, } #[derive( diff --git a/src/dynamics.rs b/src/dynamics.rs deleted file mode 100644 index d14cae0..0000000 --- a/src/dynamics.rs +++ /dev/null @@ -1,224 +0,0 @@ -//! The ecDNA data model. -use std::path::Path; - -use crate::distribution::EcDNADistribution; -use anyhow::Context; -use sosa::{write2file, NbIndividuals}; - -/// Some summary statistics for the [`EcDNADistribution`]. -#[derive(Debug, Clone)] -pub struct EcDNASummary { - pub mean: f32, - pub frequency: f32, - pub entropy: f32, -} - -/// The quantities of interest that evolve over time for the ecDNA problem. -/// -/// Those quantities are not save for every new reaction (i.e. every new -/// birth/death event), instead for snapshots taken at `timepoints`. -#[derive(Debug, Clone)] -pub struct EcDNADynamics { - /// The number of cells with ecDNAs for each iteration. - nplus: Vec, - /// The number of cells without ecDNAs for each iteration. - nminus: Vec, - /// The times at which the dynamics will be registered, used as a stack. - timepoints: Vec, - /// The timepoints that will be saved, not used as a stack but const from - /// the start to the end of the simulation - timepoints2save: Vec, - /// The time for the **current** iteration. - pub time: f32, - /// The ecDNA distribution for the **current** iteration. - pub distribution: EcDNADistribution, -} - -/// When saving the dynamics, we look at some fixed timepoints such that all -/// simulations will share the same time. -/// To do this, we bin the continous time variable into fixed-sized bins and -/// look whether the new reaction's time is smaller or equal to the timepoint. -/// It's equal, then we save the new dynamic with `SaveNewOne`. -/// If it's greater than the current timepoints (timepoints get removed once -/// used), then we save the dynamic using the one of the previous iteration with -/// `SavePreviousOne`. -/// Else, we dont save `DoNotSave`. -pub enum SaveDynamic { - /// Save the previous value for the dynamic `times` times. - SavePreviousOne { - times: usize, - }, - SaveNewOne, - DoNotSave, -} - -impl EcDNADynamics { - pub fn new( - distribution: EcDNADistribution, - iterations: usize, - timepoints: &[f32], - initial_time: f32, - ) -> Self { - // an ecDNA distribution as always an entry for the nminus cells - let mut nminus = Vec::with_capacity(iterations); - nminus.push(*distribution.get_nminus()); - let mut nplus = Vec::with_capacity(iterations); - nplus.push(distribution.compute_nplus()); - - let mut timepoints2save = timepoints.to_owned(); - // sort them to use vec as stack - timepoints2save.sort_by(|a, b| a.partial_cmp(b).unwrap()); - let timepoints = timepoints2save.clone().into_iter().rev().collect(); - // sort them decreasing order such that FILO (first the smallest) - - Self { - nplus, - nminus, - distribution, - timepoints, - time: initial_time, - timepoints2save: timepoints2save.to_vec(), - } - } - - pub fn update_nplus_nminus( - &mut self, - nplus: NbIndividuals, - nminus: NbIndividuals, - ) { - self.nplus.push(nplus); - self.nminus.push(nminus); - } - - pub fn is_time_to_save( - &mut self, - time: f32, - verbosity: u8, - ) -> SaveDynamic { - //! Is it time to save based on `time` and on the timepoints? - //! - //! This is true when the time is smaller than the timepoint and close - //! enough to it (due to float precision). - if let Some(current_t) = self.get_current_timepoint() { - let time2save = (current_t - time).abs() < f32::EPSILON; - if time2save { - if verbosity > 1 { - println!( - "saving the new dynamic at timepoint {} for time {}", - current_t, time - ); - } - self.register_timepoint(); - return SaveDynamic::SaveNewOne; - } - if &time >= current_t { - if verbosity > 1 { - println!( - "saving the previous dynamic at timepoint {} for time {}", - current_t, time - ); - } - // nb of times to repeat the dynamic - let tot_nb_bins = self.timepoints.len(); - self.timepoints.retain_mut(|&mut t| t >= time); - let nb_skipped_bins = tot_nb_bins - self.timepoints.len(); - if verbosity > 1 { - println!( - "saving the previous dynamic {} times", - nb_skipped_bins - ); - } - return SaveDynamic::SavePreviousOne { - times: nb_skipped_bins, - }; - } - - if verbosity > 1 { - println!( - "do not save at timepoint {} at time {}", - current_t, time - ); - } - }; - SaveDynamic::DoNotSave - } - - fn get_current_timepoint(&self) -> Option<&f32> { - self.timepoints.last() - } - - fn register_timepoint(&mut self) { - self.timepoints.pop().expect("timepoint already registered"); - } - - pub fn save( - &self, - path2dir: &Path, - id: &str, - verbosity: u8, - ) -> anyhow::Result<()> { - let mut nplus = path2dir.join("nplus").join(id); - nplus.set_extension("csv"); - if verbosity > 1 { - println!("Saving nplus to {:#?}", nplus); - } - write2file(&self.nplus, &nplus, None, false) - .with_context(|| "Cannot save nplus".to_string())?; - - let mut nminus = path2dir.join("nminus").join(id); - nminus.set_extension("csv"); - if verbosity > 1 { - println!("Saving nminus to {:#?}", nminus); - } - write2file(&self.nminus, &nminus, None, false) - .with_context(|| "Cannot save nminus".to_string())?; - - let mut ecdna = path2dir.join("ecdna").join(id); - ecdna.set_extension("json"); - self.distribution.save(&ecdna, verbosity)?; - - let mut times = path2dir.join("times").join(id); - times.set_extension("csv"); - if verbosity > 1 { - println!("Saving times to {:#?}", times); - } - write2file(&self.timepoints2save, ×, None, false) - .with_context(|| "Cannot save times".to_string())?; - - // the absolute gillespie time at the last iteration - let mut gillespie_time = path2dir.join("gillespie_time").join(id); - gillespie_time.set_extension("csv"); - if verbosity > 1 { - println!("Saving gillespie time to {:#?}", gillespie_time); - } - write2file(&[self.time], &gillespie_time, None, false) - .with_context(|| "Cannot save times".to_string())?; - Ok(()) - } -} - -#[cfg(test)] -mod tests { - - use quickcheck_macros::quickcheck; - - use crate::test_util::NonEmptyDistribtionWithNPlusCells; - - use super::*; - - #[quickcheck] - fn ecdna_new_test( - time: u8, - distribution: NonEmptyDistribtionWithNPlusCells, - ) -> bool { - let ecdna = EcDNADynamics::new( - distribution.clone().0, - 1, - &[time as f32, time as f32 + 0.4], - 0., - ); - ecdna.nplus.last().unwrap() == &distribution.0.compute_nplus() - && ecdna.nminus.last().unwrap() == distribution.0.get_nminus() - && ecdna.timepoints == vec![time as f32 + 0.4, time as f32] - } -} diff --git a/src/lib.rs b/src/lib.rs index ece78f8..122cbb8 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,38 +1,48 @@ -//! The two-type/k-type ecDNA simulation problem modelling the effect of the -//! random segregation on the ecDNA dynamics. -/// The ecDNA dynamics such as the mean in function of time. -pub mod dynamics; -/// The available ecDNA processes. +//! A k-type population model to simulate the effect of the random segregation +//! and positive selection on the ecDNA dynamics +/// The pure-birth and birth-death processes simulating the ecDNA dynamics. pub mod process; -/// EcDNA growth models +/// EcDNA growth models. pub mod proliferation; -/// EcDNA segregation models. pub mod segregation; -use std::path::Path; - -use ecdna_lib::distribution::EcDNADistribution; -pub use ecdna_lib::{abc, distribution, DNACopy}; -use rand::Rng; -use sosa::{NbIndividuals, ReactionRates}; - -/// Save the results of the simulations. -pub trait ToFile { - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates, - id: usize, - ) -> anyhow::Result<()>; +use std::{collections::VecDeque, path::PathBuf}; + +pub use ecdna_lib::{distribution, DNACopy}; + +#[derive(Debug, Clone)] +pub struct Snapshot { + /// The number of cells to subsample + pub cells2sample: usize, + /// The time at which we subsample + pub time: f32, +} + +#[derive(Debug, Clone)] +pub struct SavingOptions { + pub snapshots: VecDeque, + pub path2dir: PathBuf, + pub filename: String, +} + +pub fn create_filename_birth_death(rates: &[f32; 4], id: usize) -> String { + format!( + "{}b0_{}b1_{}d0_{}d1_{}idx", + rates[0].to_string().replace('.', "dot"), + rates[1].to_string().replace('.', "dot"), + rates[2].to_string().replace('.', "dot"), + rates[3].to_string().replace('.', "dot"), + id, + ) } -/// Undersample a process by randomly reducing the number of individuals. -pub trait RandomSampling { - fn random_sample( - &self, - nb_individuals: NbIndividuals, - rng: &mut impl Rng, - ) -> anyhow::Result; +pub fn create_filename_pure_birth(rates: &[f32; 2], id: usize) -> String { + format!( + "{}b0_{}b1_0d0_0d1_{}idx", + rates[0].to_string().replace('.', "dot"), + rates[1].to_string().replace('.', "dot"), + id, + ) } #[cfg(test)] diff --git a/src/main.rs b/src/main.rs index dac569b..513f242 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,35 +1,31 @@ -use app::Dynamics; +use std::{collections::VecDeque, path::PathBuf}; + use chrono::Utc; use clap_app::{GrowthOptions, Segregation}; use ecdna_evo::{ + create_filename_birth_death, create_filename_pure_birth, distribution::EcDNADistribution, - process::{ - BirthDeath, BirthDeathMean, BirthDeathMeanVariance, - BirthDeathMeanVarianceEntropy, BirthDeathNMinusNPlus, EcDNAEvent, - PureBirth, PureBirthMean, PureBirthNMinusNPlus, - }, - proliferation::{EcDNADeath, Exponential}, + process::{BirthDeath, EcDNAEvent, PureBirth}, + proliferation::{CellDeath, Exponential}, + SavingOptions, Snapshot, }; use indicatif::ParallelProgressIterator; +use rand::SeedableRng; +use rand_chacha::ChaCha8Rng; use rayon::prelude::{IntoParallelIterator, ParallelIterator}; -use sosa::{CurrentState, ReactionRates}; +use sosa::{simulate, CurrentState, Options, ReactionRates}; -use crate::{ - app::Sampling, - clap_app::{BirthDeathType, Cli, Parallel, ProcessType, PureBirthType}, -}; +use crate::clap_app::{Cli, Parallel, ProcessType}; -mod app; mod clap_app; /// The number of max iterations that we max simulate compared to the cells. /// The max number of iterations will be MAX_ITER * cells. pub const MAX_ITER: usize = 10; +#[derive(Debug)] pub struct SimulationOptions { - simulation: Dynamics, parallel: Parallel, - sampling: Option, process_type: ProcessType, b0: f32, b1: f32, @@ -39,219 +35,135 @@ pub struct SimulationOptions { segregation: Segregation, growth: GrowthOptions, runs: usize, + seed: u64, + path2dir: PathBuf, + options: Options, + snapshots: VecDeque, } fn main() { - let app = Cli::build(); + let app = Cli::build().expect("cannot construct the app"); let proliferation = match app.growth { GrowthOptions::Constant => todo!(), GrowthOptions::Exponential => Exponential {}, }; println!("{} Starting the simulation", Utc::now()); - let timepoints: [f32; 300] = std::array::from_fn(|i| i as f32 * 0.1 + 0.1); - let my_closure = |idx| match app.process_type { - ProcessType::PureBirth(p_type) => { - let mut initial_state = CurrentState { - population: [ - *app.distribution.get_nminus(), - app.distribution.compute_nplus(), - ], - }; - let rates = ReactionRates([app.b0, app.b1]); - let reactions = - [EcDNAEvent::ProliferateNMinus, EcDNAEvent::ProliferateNPlus]; - match p_type { - PureBirthType::NMinusNPlus => app - .simulation - .run( - idx, - PureBirthNMinusNPlus::with_distribution( - proliferation, - app.segregation, - app.distribution.clone(), - app.simulation.options.max_iter, - &timepoints, - 0., - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), - PureBirthType::PureBirth => app - .simulation - .run( - idx, - PureBirth::new( - app.distribution.clone(), - proliferation, - app.segregation, - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), - PureBirthType::Mean => app - .simulation - .run( - idx, - PureBirthMean::new( - proliferation, - app.segregation, - app.distribution.clone(), - app.simulation.options.max_iter, - &timepoints, - 0., - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), - }; - } - ProcessType::BirthDeath(p_type) => { - let mut initial_state = CurrentState { - population: [ - *app.distribution.get_nminus(), - app.distribution.compute_nplus(), - *app.distribution.get_nminus(), - app.distribution.compute_nplus(), - ], - }; - let rates = ReactionRates([app.b0, app.b1, app.d0, app.d1]); - let reactions = [ - EcDNAEvent::ProliferateNMinus, - EcDNAEvent::ProliferateNPlus, - EcDNAEvent::DeathNMinus, - EcDNAEvent::DeathNPlus, - ]; - match p_type { - BirthDeathType::BirthDeath => app - .simulation - .run( - idx, - BirthDeath::new( - app.distribution.clone(), - proliferation, - app.segregation, - EcDNADeath, - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), - BirthDeathType::NMinusNPlus => app - .simulation - .run( - idx, - BirthDeathNMinusNPlus::with_distribution( - proliferation, - app.segregation, - app.distribution.clone(), - app.simulation.options.max_iter, - 0., - &timepoints, - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), - BirthDeathType::Mean => app - .simulation - .run( - idx, - BirthDeathMean::new( - proliferation, - app.segregation, - app.distribution.clone(), - app.simulation.options.max_iter, - &timepoints, - 0., - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), - BirthDeathType::MeanVariance => app - .simulation - .run( - idx, - BirthDeathMeanVariance::new( - 0., - proliferation, - app.segregation, - app.distribution.clone(), - app.simulation.options.max_iter, - &timepoints, - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), - BirthDeathType::MeanVarianceEntropy => app - .simulation - .run( - idx, - BirthDeathMeanVarianceEntropy::new( - 0., - proliferation, - app.segregation, - app.distribution.clone(), - app.simulation.options.max_iter, - &timepoints, - app.simulation.options.verbosity, - ) - .unwrap(), - &mut initial_state, - &rates, - &reactions, - &app.sampling, - ) - .unwrap(), + let run_simulations = |idx| { + let stream = idx as u64; + let mut rng = ChaCha8Rng::seed_from_u64(app.seed); + rng.set_stream(stream); + match app.process_type { + ProcessType::PureBirth => { + let mut initial_state = CurrentState { + population: [ + *app.distribution.get_nminus(), + app.distribution.compute_nplus(), + ], + }; + let rates = ReactionRates([app.b0, app.b1]); + let reactions = [ + EcDNAEvent::ProliferateNMinus, + EcDNAEvent::ProliferateNPlus, + ]; + let filename = create_filename_pure_birth(&rates.0, idx); + + let mut process = PureBirth::new( + app.distribution.clone(), + proliferation, + app.segregation, + 0., + SavingOptions { + snapshots: app.snapshots.clone(), + path2dir: app.path2dir.clone(), + filename, + }, + app.options.verbosity, + ) + .unwrap(); + + if app.options.verbosity > 1 { + println!("{:#?}", process); + } + + let stop_reason = simulate( + &mut initial_state, + &rates, + &reactions, + &mut process, + &app.options, + &mut rng, + ); + + if app.options.verbosity > 0 { + println!("stop reason: {:#?}", stop_reason); + } } - } + ProcessType::BirthDeath => { + let mut initial_state = CurrentState { + population: [ + *app.distribution.get_nminus(), + app.distribution.compute_nplus(), + *app.distribution.get_nminus(), + app.distribution.compute_nplus(), + ], + }; + let rates = ReactionRates([app.b0, app.b1, app.d0, app.d1]); + let reactions = [ + EcDNAEvent::ProliferateNMinus, + EcDNAEvent::ProliferateNPlus, + EcDNAEvent::DeathNMinus, + EcDNAEvent::DeathNPlus, + ]; + let filename = create_filename_birth_death(&rates.0, idx); + + let mut process = BirthDeath::new( + app.distribution.clone(), + proliferation, + app.segregation, + CellDeath, + 0., + SavingOptions { + snapshots: app.snapshots.clone(), + path2dir: app.path2dir.clone(), + filename, + }, + app.options.verbosity, + ) + .unwrap(); + if app.options.verbosity > 1 { + println!("{:#?}", process); + } + + let stop_reason = simulate( + &mut initial_state, + &rates, + &reactions, + &mut process, + &app.options, + &mut rng, + ); + + if app.options.verbosity > 0 { + println!("stop reason: {:#?}", stop_reason); + } + } + }; }; std::process::exit({ // start from seed for the array job - let start = (app.simulation.seed * 10) as usize; + let start = (app.seed * 10) as usize; let start_end = start..app.runs + start; + match app.parallel { Parallel::Debug | Parallel::False => { - (start_end).for_each(my_closure) + (start_end).for_each(run_simulations) } Parallel::True => (start_end) .into_par_iter() .progress_count(app.runs as u64) - .for_each(my_closure), + .for_each(run_simulations), } println!("{} End simulation", Utc::now()); 0 diff --git a/src/process.rs b/src/process.rs index 813de34..227e5ac 100644 --- a/src/process.rs +++ b/src/process.rs @@ -1,43 +1,22 @@ -use std::{fs, path::Path}; - use anyhow::ensure; use rand::Rng; -use sosa::{ - write2file, AdvanceStep, CurrentState, NbIndividuals, NextReaction, - ReactionRates, +use sosa::{AdvanceStep, CurrentState, NextReaction}; +use std::{ + collections::VecDeque, + fs, + path::{Path, PathBuf}, }; use crate::{ - distribution::EcDNADistribution, dynamics::SaveDynamic, - segregation::Segregate, RandomSampling, ToFile, + distribution::EcDNADistribution, segregation::Segregate, SavingOptions, + Snapshot, }; use super::{ - dynamics::EcDNADynamics, - proliferation::{EcDNADeath, EcDNAProliferation}, + proliferation::{CellDeath, EcDNAProliferation}, segregation::IsUneven, }; -fn create_filename_birth_death(rates: &[f32; 4], id: usize) -> String { - format!( - "{}b0_{}b1_{}d0_{}d1_{}idx", - rates[0].to_string().replace('.', "dot"), - rates[1].to_string().replace('.', "dot"), - rates[2].to_string().replace('.', "dot"), - rates[3].to_string().replace('.', "dot"), - id, - ) -} - -fn create_filename_pure_birth(rates: &[f32; 2], id: usize) -> String { - format!( - "{}b0_{}b1_0d0_0d1_{}idx", - rates[0].to_string().replace('.', "dot"), - rates[1].to_string().replace('.', "dot"), - id, - ) -} - #[derive(Debug, Clone, Copy)] pub enum EcDNAEvent { ProliferateNPlus, @@ -49,11 +28,34 @@ pub enum EcDNAEvent { SymmetricDifferentiation, } +fn save( + path2dir: &Path, + filename: &str, + time: f32, + distribution: &EcDNADistribution, + verbosity: u8, +) -> anyhow::Result<()> { + //! Helper fn to save the ecDNA distribution + let cells = distribution.compute_nplus() + *distribution.get_nminus(); + let path2file = &path2dir.join(format!("{}cells/ecdna/", cells)); + let mut timepoint = format!("{:.1}", time).replace('.', "dot"); + timepoint.push_str("years"); + let mut path2file = path2file.join(timepoint).join(filename); + path2file = path2file.with_extension("json"); + fs::create_dir_all(path2file.parent().unwrap()) + .expect("Cannot create dir"); + if verbosity > 0 { + println!( + "saving state at time {} with {} cells in {:#?}", + time, cells, path2file + ); + } + distribution.save(&path2file, verbosity)?; + Ok(()) +} + /// The simulation of the [`EcDNADistribution`] according to a pure-birth /// stochastic process. -/// -/// The distribution can be saved during or at the end of the simulation using -/// [`Self::save`]. #[derive(Debug, Clone)] pub struct PureBirth where @@ -63,6 +65,10 @@ where distribution: EcDNADistribution, proliferation: P, segregation: S, + time: f32, + snapshots: VecDeque, + path2dir: PathBuf, + filename: String, verbosity: u8, } @@ -75,10 +81,21 @@ where distribution: EcDNADistribution, proliferation: P, segregation: S, + time: f32, + saving_options: SavingOptions, verbosity: u8, ) -> anyhow::Result { ensure!(!distribution.is_empty()); - Ok(Self { distribution, proliferation, segregation, verbosity }) + Ok(Self { + distribution, + proliferation, + segregation, + time, + path2dir: saving_options.path2dir, + filename: saving_options.filename, + snapshots: saving_options.snapshots, + verbosity, + }) } pub fn get_mut_ecdna_distribution(&mut self) -> &mut EcDNADistribution { @@ -90,46 +107,6 @@ where } } -impl RandomSampling for PureBirth -where - P: EcDNAProliferation, - S: Segregate, -{ - fn random_sample( - &self, - nb_individuals: NbIndividuals, - rng: &mut impl Rng, - ) -> anyhow::Result { - ensure!( - nb_individuals - < self.distribution.get_nminus() - + self.distribution.compute_nplus() - ); - Ok(self.distribution.into_subsampled(nb_individuals, rng)) - } -} - -impl ToFile<2> for PureBirth -where - P: EcDNAProliferation, - S: Segregate, -{ - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<2>, - id: usize, - ) -> anyhow::Result<()> { - //! Save the ecDNA distribution - let path2ecdna = path2dir.join("ecdna"); - fs::create_dir_all(&path2ecdna).expect("Cannot create dir"); - let filename = create_filename_pure_birth(&rates.0, id); - let path2file = path2ecdna.join(filename).with_extension("json"); - self.distribution.save(&path2file, self.verbosity)?; - Ok(()) - } -} - impl AdvanceStep<2> for PureBirth { type Reaction = EcDNAEvent; @@ -138,6 +115,44 @@ impl AdvanceStep<2> for PureBirth { reaction: NextReaction, rng: &mut impl Rng, ) { + while !self.snapshots.is_empty() + && self.snapshots.iter().any(|s| self.time >= s.time) + { + let snapshot = self.snapshots.pop_front().unwrap(); + if self.verbosity > 0 { + println!( + "saving state for timepoint at time {:#?} at simulation's time {} with {} cells", + snapshot.time, self.time, snapshot.cells2sample + ); + } + let cells = (*self.distribution.get_nminus() + + self.distribution.compute_nplus()) + as usize; + + if snapshot.cells2sample == cells || cells < snapshot.cells2sample + { + save( + &self.path2dir, + &self.filename, + snapshot.time, + &self.distribution, + self.verbosity, + ) + .expect("cannot save snapshot"); + } else { + save( + &self.path2dir, + &self.filename, + snapshot.time, + &self + .distribution + .into_subsampled(snapshot.cells2sample as u64, rng), + self.verbosity, + ) + .expect("cannot save snapshot"); + } + } + match reaction.event { EcDNAEvent::ProliferateNPlus => { if let Ok(is_uneven) = self.proliferation.increase_nplus( @@ -175,6 +190,7 @@ impl AdvanceStep<2> for PureBirth { if self.verbosity > 1 { println!("Distribution {:#?}", self.distribution); } + self.time += reaction.time; } fn update_state(&self, state: &mut CurrentState<2>) { @@ -189,340 +205,8 @@ impl AdvanceStep<2> for PureBirth { } } -/// A pure-birth stochastic process tracking the evolution of the -/// [`EcDNADynamics`] over time. -#[derive(Debug, Clone)] -pub struct PureBirthNMinusNPlus { - data: EcDNADynamics, - proliferation: P, - segregation: S, - verbosity: u8, -} - -impl PureBirthNMinusNPlus { - pub fn with_distribution( - proliferation: P, - segregation: S, - distribution: EcDNADistribution, - max_iter: usize, - timepoints: &[f32], - initial_time: f32, - verbosity: u8, - ) -> anyhow::Result { - //! Create a pure-birth process from a non-empty `distribution` - //! tracking the nminus and nplus cells over time. - //! - //! ## Error - //! Returns error when the distribution is empty. - ensure!(!distribution.is_empty()); - Ok(Self { - data: EcDNADynamics::new( - distribution, - max_iter, - timepoints, - initial_time, - ), - proliferation, - segregation, - verbosity, - }) - } -} - -impl RandomSampling for PureBirthNMinusNPlus -where - P: EcDNAProliferation, - S: Segregate, -{ - fn random_sample( - &self, - _nb_individuals: NbIndividuals, - _rng: &mut impl Rng, - ) -> anyhow::Result { - todo!(); - // let (nplus, nminus) = ( - // self.data.distribution.compute_nplus(), - // *self.data.distribution.get_nminus(), - // ); - // self.data.nplus.push(nplus); - // self.data.nminus.push(nminus); - // if self.verbosity > 1 { - // println!( - // "NPlus/NMinus after subsampling {}, {}", - // self.data.nplus.last().unwrap(), - // self.data.nminus.last().unwrap() - // ); - // } - } -} - -impl ToFile<2> for PureBirthNMinusNPlus -where - P: EcDNAProliferation, - S: Segregate, -{ - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<2>, - id: usize, - ) -> anyhow::Result<()> { - fs::create_dir_all(path2dir).expect("Cannot create dir"); - let filename = create_filename_pure_birth(&rates.0, id); - self.data.save(path2dir, &filename, self.verbosity)?; - Ok(()) - } -} - -impl AdvanceStep<2> - for PureBirthNMinusNPlus -{ - type Reaction = EcDNAEvent; - - fn advance_step( - &mut self, - reaction: NextReaction, - rng: &mut impl Rng, - ) { - match reaction.event { - EcDNAEvent::ProliferateNPlus => { - if let Ok(is_uneven) = self.proliferation.increase_nplus( - &mut self.data.distribution, - &self.segregation, - rng, - self.verbosity, - ) { - match is_uneven { - IsUneven::False => { - if self.verbosity > 1 { - println!("IsUneven is false"); - } - } - IsUneven::True => { - if self.verbosity > 1 { - println!("IsUneven is true"); - } - } - IsUneven::TrueWithoutNMinusIncrease => { - if self.verbosity > 1 { - println!( - "IsUneven is true but w/o nminus increase" - ); - } - } - } - } - } - EcDNAEvent::ProliferateNMinus => { - self.proliferation - .increase_nminus(&mut self.data.distribution); - } - _ => unreachable!(), - }; - if self.verbosity > 1 { - println!("Distribution {:#?}", self.data.distribution); - } - - // store absolute time - self.data.time += reaction.time; - - match self.data.is_time_to_save(self.data.time, self.verbosity) { - SaveDynamic::DoNotSave => {} - SaveDynamic::SaveNewOne => { - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - SaveDynamic::SavePreviousOne { times } => { - for _ in 0..times { - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - } - }; - } - - fn update_state(&self, state: &mut CurrentState<2>) { - state.population[0] = *self.data.distribution.get_nminus(); - state.population[1] = self.data.distribution.compute_nplus(); - if self.verbosity > 1 { - println!( - "Update iteration after process update, population: {:#?}", - state.population - ); - } - } -} - -#[derive(Debug, Clone)] -pub struct PureBirthMean { - data: EcDNADynamics, - mean: Vec, - proliferation: P, - segregation: S, - verbosity: u8, -} - -impl PureBirthMean { - pub fn new( - proliferation: P, - segregation: S, - distribution: EcDNADistribution, - max_iter: usize, - timepoints: &[f32], - initial_time: f32, - verbosity: u8, - ) -> anyhow::Result { - ensure!(!distribution.is_empty()); - let mut mean = Vec::with_capacity(max_iter); - mean.push(distribution.compute_mean()); - Ok(Self { - data: EcDNADynamics::new( - distribution, - max_iter, - timepoints, - initial_time, - ), - mean, - proliferation, - segregation, - verbosity, - }) - } -} - -impl RandomSampling for PureBirthMean -where - P: EcDNAProliferation, - S: Segregate, -{ - fn random_sample( - &self, - _nb_individuals: NbIndividuals, - _rng: &mut impl Rng, - ) -> anyhow::Result { - todo!() - } -} - -impl ToFile<2> for PureBirthMean -where - P: EcDNAProliferation, - S: Segregate, -{ - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<2>, - id: usize, - ) -> anyhow::Result<()> { - fs::create_dir_all(path2dir).expect("Cannot create dir"); - - let filename = create_filename_pure_birth(&rates.0, id); - if self.verbosity > 1 { - println!("Saving data in {:#?}", path2dir) - } - - self.data.save(path2dir, &filename, self.verbosity)?; - - let mut mean = path2dir.join("mean").join(filename); - mean.set_extension("csv"); - if self.verbosity > 1 { - println!("Mean data in {:#?}", mean) - } - write2file(&self.mean, &mean, None, false)?; - Ok(()) - } -} - -impl AdvanceStep<2> - for PureBirthMean -{ - type Reaction = EcDNAEvent; - - fn advance_step( - &mut self, - reaction: NextReaction, - rng: &mut impl Rng, - ) { - match reaction.event { - EcDNAEvent::ProliferateNPlus => { - if let Ok(is_uneven) = self.proliferation.increase_nplus( - &mut self.data.distribution, - &self.segregation, - rng, - self.verbosity, - ) { - match is_uneven { - IsUneven::False => { - if self.verbosity > 1 { - println!("IsUneven is false"); - } - } - IsUneven::True => { - if self.verbosity > 1 { - println!("IsUneven is true"); - } - } - IsUneven::TrueWithoutNMinusIncrease => { - if self.verbosity > 1 { - println!( - "IsUneven is true but w/o nminus increase" - ); - } - } - } - } - } - EcDNAEvent::ProliferateNMinus => { - self.proliferation - .increase_nminus(&mut self.data.distribution); - } - _ => unreachable!(), - }; - if self.verbosity > 1 { - println!("Distribution {:#?}", self.data.distribution); - } - - // store absolute time - self.data.time += reaction.time; - - match self.data.is_time_to_save(self.data.time, self.verbosity) { - SaveDynamic::DoNotSave => {} - SaveDynamic::SaveNewOne => { - self.mean.push(self.data.distribution.compute_mean()); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - SaveDynamic::SavePreviousOne { times } => { - let mean = self.data.distribution.compute_mean(); - for _ in 0..times { - self.mean.push(mean); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - } - }; - } - - fn update_state(&self, state: &mut CurrentState<2>) { - state.population[0] = *self.data.distribution.get_nminus(); - state.population[1] = self.data.distribution.compute_nplus(); - } -} - /// The simulation of the [`EcDNADistribution`] according to a birth-death /// stochastic process. -/// -/// The distribution can be saved during or at the end of the simulation using -/// [`Self::save`]. #[derive(Debug, Clone)] pub struct BirthDeath where @@ -532,7 +216,11 @@ where distribution: EcDNADistribution, proliferation: P, segregation: S, - death: EcDNADeath, + death: CellDeath, + time: f32, + snapshots: VecDeque, + path2dir: PathBuf, + filename: String, verbosity: u8, } @@ -545,11 +233,23 @@ where distribution: EcDNADistribution, proliferation: P, segregation: S, - death: EcDNADeath, + death: CellDeath, + time: f32, + saving_options: SavingOptions, verbosity: u8, ) -> anyhow::Result { ensure!(!distribution.is_empty()); - Ok(Self { distribution, proliferation, segregation, death, verbosity }) + Ok(Self { + distribution, + proliferation, + segregation, + death, + time, + path2dir: saving_options.path2dir, + filename: saving_options.filename, + snapshots: saving_options.snapshots, + verbosity, + }) } pub fn get_mut_ecdna_distribution(&mut self) -> &mut EcDNADistribution { @@ -569,6 +269,42 @@ impl AdvanceStep<4> for BirthDeath { reaction: NextReaction, rng: &mut impl Rng, ) { + while !self.snapshots.is_empty() + && self.snapshots.iter().any(|s| self.time >= s.time) + { + let snapshot = self.snapshots.pop_front().unwrap(); + if self.verbosity > 0 { + println!( + "saving state for timepoint at time {:#?} at simulation's time {} with {} cells", + snapshot.time, self.time, snapshot.cells2sample + ); + } + let cells = (*self.distribution.get_nminus() + + self.distribution.compute_nplus()) + as usize; + + if snapshot.cells2sample == cells { + save( + &self.path2dir, + &self.filename, + snapshot.time, + &self.distribution, + self.verbosity, + ) + .expect("cannot save snapshot"); + } else { + save( + &self.path2dir, + &self.filename, + snapshot.time, + &self + .distribution + .into_subsampled(snapshot.cells2sample as u64, rng), + self.verbosity, + ) + .expect("cannot save snapshot"); + } + } match reaction.event { EcDNAEvent::ProliferateNPlus => { if let Ok(is_uneven) = self.proliferation.increase_nplus( @@ -614,6 +350,7 @@ impl AdvanceStep<4> for BirthDeath { if self.verbosity > 1 { println!("Distribution {:#?}", self.distribution); } + self.time += reaction.time; } fn update_state(&self, state: &mut CurrentState<4>) { @@ -624,792 +361,6 @@ impl AdvanceStep<4> for BirthDeath { } } -impl ToFile<4> for BirthDeath { - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<4>, - id: usize, - ) -> anyhow::Result<()> { - //! Save the ecDNA distribution - let path2ecdna = path2dir.join("ecdna"); - let filename = create_filename_birth_death(&rates.0, id); - fs::create_dir_all(&path2ecdna).expect("Cannot create dir"); - let path2file = path2ecdna.join(filename).with_extension("json"); - self.distribution.save(&path2file, self.verbosity)?; - Ok(()) - } -} - -impl RandomSampling for BirthDeath { - fn random_sample( - &self, - nb_individuals: NbIndividuals, - rng: &mut impl Rng, - ) -> anyhow::Result { - ensure!( - nb_individuals - < self.distribution.get_nminus() - + self.distribution.compute_nplus() - ); - Ok(self.distribution.into_subsampled(nb_individuals, rng)) - } -} - -/// A birth-death stochastic process tracking the evolution of the -/// [`EcDNADynamics`] per iterations. -#[derive(Debug, Clone)] -pub struct BirthDeathNMinusNPlus { - data: EcDNADynamics, - proliferation: P, - segregation: S, - death: EcDNADeath, - verbosity: u8, -} - -impl BirthDeathNMinusNPlus { - pub fn with_distribution( - proliferation: P, - segregation: S, - distribution: EcDNADistribution, - max_iter: usize, - initial_time: f32, - timepoints: &[f32], - verbosity: u8, - ) -> anyhow::Result { - ensure!(!distribution.is_empty()); - Ok(Self { - data: EcDNADynamics::new( - distribution, - max_iter, - timepoints, - initial_time, - ), - proliferation, - death: EcDNADeath, - segregation, - verbosity, - }) - } -} -impl RandomSampling for BirthDeathNMinusNPlus -where - P: EcDNAProliferation, - S: Segregate, -{ - fn random_sample( - &self, - _nb_individuals: NbIndividuals, - _rng: &mut impl Rng, - ) -> anyhow::Result { - todo!() - } -} - -impl ToFile<4> for BirthDeathNMinusNPlus -where - P: EcDNAProliferation, - S: Segregate, -{ - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<4>, - id: usize, - ) -> anyhow::Result<()> { - fs::create_dir_all(path2dir).expect("Cannot create dir"); - - let filename = create_filename_birth_death(&rates.0, id); - if self.verbosity > 1 { - println!("Saving data in {:#?}", path2dir) - } - - self.data.save(path2dir, &filename, self.verbosity)?; - - Ok(()) - } -} - -impl AdvanceStep<4> - for BirthDeathNMinusNPlus -{ - type Reaction = EcDNAEvent; - fn advance_step( - &mut self, - reaction: NextReaction, - rng: &mut impl Rng, - ) { - match reaction.event { - EcDNAEvent::ProliferateNPlus => { - if let Ok(is_uneven) = self.proliferation.increase_nplus( - &mut self.data.distribution, - &self.segregation, - rng, - self.verbosity, - ) { - match is_uneven { - IsUneven::False => { - if self.verbosity > 1 { - println!("IsUneven is false"); - } - } - IsUneven::True => { - if self.verbosity > 1 { - println!("IsUneven is true"); - } - } - IsUneven::TrueWithoutNMinusIncrease => { - if self.verbosity > 1 { - println!( - "IsUneven is true but w/o nminus increase" - ); - } - } - } - } - } - EcDNAEvent::ProliferateNMinus => { - self.proliferation - .increase_nminus(&mut self.data.distribution); - } - EcDNAEvent::DeathNMinus => { - self.death.decrease_nminus(&mut self.data.distribution) - } - - EcDNAEvent::DeathNPlus => self.death.decrease_nplus( - &mut self.data.distribution, - rng, - self.verbosity, - ), - - _ => unreachable!(), - }; - if self.verbosity > 1 { - println!("Distribution {:#?}", self.data.distribution); - } - - // store absolute time - self.data.time += reaction.time; - - match self.data.is_time_to_save(self.data.time, self.verbosity) { - SaveDynamic::DoNotSave => {} - SaveDynamic::SaveNewOne => { - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - SaveDynamic::SavePreviousOne { times } => { - for _ in 0..times { - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - } - }; - } - - fn update_state(&self, state: &mut CurrentState<4>) { - state.population[0] = *self.data.distribution.get_nminus(); - state.population[1] = self.data.distribution.compute_nplus(); - state.population[2] = *self.data.distribution.get_nminus(); - state.population[3] = self.data.distribution.compute_nplus(); - if self.verbosity > 1 { - println!( - "Update iteration after process update, population: {:#?}", - state.population - ); - } - } -} - -#[derive(Debug, Clone)] -pub struct BirthDeathMean { - data: EcDNADynamics, - mean: Vec, - proliferation: P, - segregation: S, - death: EcDNADeath, - verbosity: u8, -} - -impl BirthDeathMean { - pub fn new( - proliferation: P, - segregation: S, - distribution: EcDNADistribution, - max_iter: usize, - timepoints: &[f32], - initial_time: f32, - verbosity: u8, - ) -> anyhow::Result { - ensure!(!distribution.is_empty()); - let mut mean = Vec::with_capacity(max_iter); - mean.push(distribution.compute_mean()); - Ok(Self { - data: EcDNADynamics::new( - distribution, - max_iter, - timepoints, - initial_time, - ), - proliferation, - death: EcDNADeath, - mean, - segregation, - verbosity, - }) - } -} - -impl RandomSampling for BirthDeathMean -where - P: EcDNAProliferation, - S: Segregate, -{ - fn random_sample( - &self, - _nb_individuals: NbIndividuals, - _rng: &mut impl Rng, - ) -> anyhow::Result { - todo!() - } -} - -impl ToFile<4> for BirthDeathMean -where - P: EcDNAProliferation, - S: Segregate, -{ - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<4>, - id: usize, - ) -> anyhow::Result<()> { - fs::create_dir_all(path2dir).expect("Cannot create dir"); - - let filename = create_filename_birth_death(&rates.0, id); - if self.verbosity > 1 { - println!("Saving data in {:#?}", path2dir) - } - - self.data.save(path2dir, &filename, self.verbosity)?; - - let mut mean = path2dir.join("mean").join(filename); - mean.set_extension("csv"); - if self.verbosity > 1 { - println!("Mean data in {:#?}", mean) - } - write2file(&self.mean, &mean, None, false)?; - Ok(()) - } -} - -impl AdvanceStep<4> - for BirthDeathMean -{ - type Reaction = EcDNAEvent; - - fn advance_step( - &mut self, - reaction: NextReaction, - rng: &mut impl Rng, - ) { - match reaction.event { - EcDNAEvent::ProliferateNPlus => { - if let Ok(is_uneven) = self.proliferation.increase_nplus( - &mut self.data.distribution, - &self.segregation, - rng, - self.verbosity, - ) { - match is_uneven { - IsUneven::False => { - if self.verbosity > 1 { - println!("IsUneven is false"); - } - } - IsUneven::True => { - if self.verbosity > 1 { - println!("IsUneven is true"); - } - } - IsUneven::TrueWithoutNMinusIncrease => { - if self.verbosity > 1 { - println!( - "IsUneven is true but w/o nminus increase" - ); - } - } - } - } - } - EcDNAEvent::ProliferateNMinus => { - self.proliferation - .increase_nminus(&mut self.data.distribution); - } - EcDNAEvent::DeathNMinus => { - self.death.decrease_nminus(&mut self.data.distribution) - } - EcDNAEvent::DeathNPlus => self.death.decrease_nplus( - &mut self.data.distribution, - rng, - self.verbosity, - ), - _ => unreachable!(), - }; - if self.verbosity > 1 { - println!("Distribution {:#?}", self.data.distribution); - } - - // store absolute time - self.data.time += reaction.time; - - match self.data.is_time_to_save(self.data.time, self.verbosity) { - SaveDynamic::DoNotSave => {} - SaveDynamic::SaveNewOne => { - self.mean.push(self.data.distribution.compute_mean()); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - SaveDynamic::SavePreviousOne { times } => { - let mean = self.data.distribution.compute_mean(); - for _ in 0..times { - self.mean.push(mean); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - } - }; - } - - fn update_state(&self, state: &mut CurrentState<4>) { - state.population[0] = *self.data.distribution.get_nminus(); - state.population[1] = self.data.distribution.compute_nplus(); - state.population[2] = *self.data.distribution.get_nminus(); - state.population[3] = self.data.distribution.compute_nplus(); - } -} - -#[derive(Debug, Clone)] -pub struct BirthDeathMeanVariance { - data: EcDNADynamics, - mean: Vec, - variance: Vec, - proliferation: P, - segregation: S, - death: EcDNADeath, - verbosity: u8, -} - -impl BirthDeathMeanVariance { - pub fn new( - time: f32, - proliferation: P, - segregation: S, - distribution: EcDNADistribution, - max_iter: usize, - timepoints: &[f32], - verbosity: u8, - ) -> anyhow::Result { - ensure!(!distribution.is_empty()); - let mut mean = Vec::with_capacity(max_iter); - mean.push(distribution.compute_mean()); - let mut variance = Vec::with_capacity(max_iter); - variance.push(distribution.compute_variance()); - Ok(Self { - data: EcDNADynamics::new(distribution, max_iter, timepoints, time), - mean, - variance, - proliferation, - death: EcDNADeath, - segregation, - verbosity, - }) - } -} - -impl RandomSampling for BirthDeathMeanVariance -where - P: EcDNAProliferation, - S: Segregate, -{ - fn random_sample( - &self, - _nb_individuals: NbIndividuals, - _rng: &mut impl Rng, - ) -> anyhow::Result { - todo!(); - // let (nplus, nminus) = ( - // self.data.distribution.compute_nplus(), - // *self.data.distribution.get_nminus(), - // ); - // self.data.nplus.push(nplus); - // self.data.nminus.push(nminus); - // self.data.time.push(0f32); - // self.mean.push(self.data.distribution.compute_mean()); - // self.variance.push(self.data.distribution.compute_variance()); - - // if self.verbosity > 1 { - // println!( - // "NPlus/NMinus after subsampling {}, {}", - // self.data.nplus.last().unwrap(), - // self.data.nminus.last().unwrap() - // ); - // } - } -} - -impl ToFile<4> for BirthDeathMeanVariance -where - P: EcDNAProliferation, - S: Segregate, -{ - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<4>, - id: usize, - ) -> anyhow::Result<()> { - fs::create_dir_all(path2dir).expect("Cannot create dir"); - - let filename = create_filename_birth_death(&rates.0, id); - if self.verbosity > 1 { - println!("Saving data in {:#?}", path2dir) - } - - self.data.save(path2dir, &filename, self.verbosity)?; - - let mut mean = path2dir.join("mean").join(filename.clone()); - mean.set_extension("csv"); - if self.verbosity > 1 { - println!("Mean data in {:#?}", mean) - } - write2file(&self.mean, &mean, None, false)?; - - let mut variance = path2dir.join("variance").join(filename); - variance.set_extension("csv"); - if self.verbosity > 1 { - println!("Variance data in {:#?}", path2dir) - } - write2file(&self.variance, &variance, None, false)?; - Ok(()) - } -} - -impl AdvanceStep<4> - for BirthDeathMeanVariance -{ - type Reaction = EcDNAEvent; - - fn advance_step( - &mut self, - reaction: NextReaction, - rng: &mut impl Rng, - ) { - match reaction.event { - EcDNAEvent::ProliferateNPlus => { - if let Ok(is_uneven) = self.proliferation.increase_nplus( - &mut self.data.distribution, - &self.segregation, - rng, - self.verbosity, - ) { - match is_uneven { - IsUneven::False => { - if self.verbosity > 1 { - println!("IsUneven is false"); - } - } - IsUneven::True => { - if self.verbosity > 1 { - println!("IsUneven is true"); - } - } - IsUneven::TrueWithoutNMinusIncrease => { - if self.verbosity > 1 { - println!( - "IsUneven is true but w/o nminus increase" - ); - } - } - } - } - } - EcDNAEvent::ProliferateNMinus => { - self.proliferation - .increase_nminus(&mut self.data.distribution); - } - EcDNAEvent::DeathNMinus => { - self.death.decrease_nminus(&mut self.data.distribution) - } - EcDNAEvent::DeathNPlus => self.death.decrease_nplus( - &mut self.data.distribution, - rng, - self.verbosity, - ), - _ => unreachable!(), - }; - if self.verbosity > 1 { - println!("Distribution {:#?}", self.data.distribution); - } - - // store absolute time - self.data.time += reaction.time; - - match self.data.is_time_to_save(self.data.time, self.verbosity) { - SaveDynamic::DoNotSave => {} - SaveDynamic::SaveNewOne => { - self.mean.push(self.data.distribution.compute_mean()); - self.variance.push(self.data.distribution.compute_variance()); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - SaveDynamic::SavePreviousOne { times } => { - let mean = self.data.distribution.compute_mean(); - let variance = self.data.distribution.compute_variance(); - for _ in 0..times { - self.mean.push(mean); - self.variance.push(variance); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - } - }; - } - - fn update_state(&self, state: &mut CurrentState<4>) { - state.population[0] = *self.data.distribution.get_nminus(); - state.population[1] = self.data.distribution.compute_nplus(); - state.population[2] = *self.data.distribution.get_nminus(); - state.population[3] = self.data.distribution.compute_nplus(); - } -} - -#[derive(Debug, Clone)] -pub struct BirthDeathMeanVarianceEntropy { - data: EcDNADynamics, - mean: Vec, - variance: Vec, - entropy: Vec, - proliferation: P, - segregation: S, - death: EcDNADeath, - verbosity: u8, -} - -impl BirthDeathMeanVarianceEntropy { - pub fn new( - time: f32, - proliferation: P, - segregation: S, - distribution: EcDNADistribution, - max_iter: usize, - timepoints: &[f32], - verbosity: u8, - ) -> anyhow::Result { - ensure!(!distribution.is_empty()); - let mut mean = Vec::with_capacity(max_iter); - mean.push(distribution.compute_mean()); - let mut variance = Vec::with_capacity(max_iter); - variance.push(distribution.compute_variance()); - let mut entropy = Vec::with_capacity(max_iter); - entropy.push(distribution.compute_entropy()); - Ok(Self { - data: EcDNADynamics::new(distribution, max_iter, timepoints, time), - mean, - variance, - entropy, - proliferation, - death: EcDNADeath, - segregation, - verbosity, - }) - } -} - -impl RandomSampling for BirthDeathMeanVarianceEntropy -where - P: EcDNAProliferation, - S: Segregate, -{ - fn random_sample( - &self, - _nb_individuals: NbIndividuals, - _rng: &mut impl Rng, - ) -> anyhow::Result { - todo!(); - // let (nplus, nminus) = ( - // self.data.distribution.compute_nplus(), - // *self.data.distribution.get_nminus(), - // ); - // self.data.nplus.push(nplus); - // self.data.nminus.push(nminus); - // self.data.time.push(0f32); - // self.mean.push(self.data.distribution.compute_mean()); - // self.variance.push(self.data.distribution.compute_variance()); - // self.entropy.push(self.data.distribution.compute_entropy()); - - // if self.verbosity > 1 { - // println!( - // "NPlus/NMinus after subsampling {}, {}", - // self.data.nplus.last().unwrap(), - // self.data.nminus.last().unwrap() - // ); - // } - } -} - -impl ToFile<4> for BirthDeathMeanVarianceEntropy -where - P: EcDNAProliferation, - S: Segregate, -{ - fn save( - &self, - path2dir: &Path, - rates: &ReactionRates<4>, - id: usize, - ) -> anyhow::Result<()> { - fs::create_dir_all(path2dir).expect("Cannot create dir"); - - let filename = create_filename_birth_death(&rates.0, id); - if self.verbosity > 1 { - println!("Saving data in {:#?}", path2dir) - } - - self.data.save(path2dir, &filename, self.verbosity)?; - - let mut mean = path2dir.join("mean").join(filename.clone()); - mean.set_extension("csv"); - if self.verbosity > 1 { - println!("Mean data in {:#?}", mean) - } - write2file(&self.mean, &mean, None, false)?; - - let mut variance = path2dir.join("variance").join(filename.clone()); - variance.set_extension("csv"); - if self.verbosity > 1 { - println!("Variance data in {:#?}", path2dir) - } - write2file(&self.variance, &variance, None, false)?; - - let mut entropy = path2dir.join("entropy").join(filename); - entropy.set_extension("csv"); - if self.verbosity > 1 { - println!("entropy data in {:#?}", path2dir) - } - write2file(&self.entropy, &entropy, None, false)?; - Ok(()) - } -} - -impl AdvanceStep<4> - for BirthDeathMeanVarianceEntropy -{ - type Reaction = EcDNAEvent; - - fn advance_step( - &mut self, - reaction: NextReaction, - rng: &mut impl Rng, - ) { - match reaction.event { - EcDNAEvent::ProliferateNPlus => { - if let Ok(is_uneven) = self.proliferation.increase_nplus( - &mut self.data.distribution, - &self.segregation, - rng, - self.verbosity, - ) { - match is_uneven { - IsUneven::False => { - if self.verbosity > 1 { - println!("IsUneven is false"); - } - } - IsUneven::True => { - if self.verbosity > 1 { - println!("IsUneven is true"); - } - } - IsUneven::TrueWithoutNMinusIncrease => { - if self.verbosity > 1 { - println!( - "IsUneven is true but w/o nminus increase" - ); - } - } - } - } - } - EcDNAEvent::ProliferateNMinus => { - self.proliferation - .increase_nminus(&mut self.data.distribution); - } - EcDNAEvent::DeathNMinus => { - self.death.decrease_nminus(&mut self.data.distribution) - } - EcDNAEvent::DeathNPlus => self.death.decrease_nplus( - &mut self.data.distribution, - rng, - self.verbosity, - ), - _ => unreachable!(), - }; - if self.verbosity > 1 { - println!("Distribution {:#?}", self.data.distribution); - } - - // store absolute time - self.data.time += reaction.time; - - match self.data.is_time_to_save(self.data.time, self.verbosity) { - SaveDynamic::DoNotSave => {} - SaveDynamic::SaveNewOne => { - self.mean.push(self.data.distribution.compute_mean()); - self.variance.push(self.data.distribution.compute_variance()); - self.entropy.push(self.data.distribution.compute_entropy()); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - SaveDynamic::SavePreviousOne { times } => { - let mean = self.data.distribution.compute_mean(); - let variance = self.data.distribution.compute_variance(); - let entropy = self.data.distribution.compute_entropy(); - for _ in 0..times { - self.mean.push(mean); - self.variance.push(variance); - self.entropy.push(entropy); - self.data.update_nplus_nminus( - self.data.distribution.compute_nplus(), - *self.data.distribution.get_nminus(), - ); - } - } - }; - } - - fn update_state(&self, state: &mut CurrentState<4>) { - state.population[0] = *self.data.distribution.get_nminus(); - state.population[1] = self.data.distribution.compute_nplus(); - state.population[2] = *self.data.distribution.get_nminus(); - state.population[3] = self.data.distribution.compute_nplus(); - } -} - #[cfg(test)] mod tests { use super::*; @@ -1420,9 +371,8 @@ mod tests { use quickcheck_macros::quickcheck; #[quickcheck] - fn create_process_nplus_nminus_mean_time_test( + fn create_birth_death_process_test( distr: NonEmptyDistribtionWithNPlusCells, - max_iter: u8, verbosity: u8, time: u8, ) -> bool { @@ -1430,19 +380,26 @@ mod tests { let nplus = distr.0.compute_nplus(); let nminus = *distr.0.get_nminus(); let mean = distr.0.compute_mean(); - let p = BirthDeathMean::new( + let p = BirthDeath::new( + distr.0, Exponential {}, Binomial, - distr.0, - max_iter as usize, - &[1., 2.], + CellDeath, time, + SavingOptions { + snapshots: VecDeque::from([Snapshot { + cells2sample: 10, + time: 2., + }]), + path2dir: PathBuf::default(), + filename: String::from("filename"), + }, verbosity, ) .unwrap(); - p.data.distribution.compute_nplus() == nplus - && p.data.distribution.get_nminus() == &nminus - && p.data.time == time - && p.data.distribution.compute_mean() == mean + p.distribution.compute_nplus() == nplus + && p.distribution.get_nminus() == &nminus + && p.time == time + && p.distribution.compute_mean() == mean } } diff --git a/src/proliferation.rs b/src/proliferation.rs index c542817..c0de0fc 100644 --- a/src/proliferation.rs +++ b/src/proliferation.rs @@ -4,7 +4,7 @@ use std::num::NonZeroU16; use crate::distribution::EcDNADistribution; use crate::segregation::{DNACopySegregating, IsUneven, Segregate}; -/// Update the [`EcDNADistribution`]. +/// Update the [`EcDNADistribution`] according to the random segregation. pub trait EcDNAProliferation { fn increase_nplus( &self, @@ -16,6 +16,8 @@ pub trait EcDNAProliferation { fn increase_nminus(&self, distribution: &mut EcDNADistribution); } +/// The ecDNA dynamics according to an exponential growth model with random +/// segregation. #[derive(Debug, Clone, Copy)] pub struct Exponential {} @@ -115,10 +117,12 @@ impl EcDNAProliferation for Exponential { } } +/// Add cell death as an event in the model to simulate a birth-death process. #[derive(Debug, Clone)] -pub struct EcDNADeath; +pub struct CellDeath; -impl EcDNADeath { +/// Update the ecDNA distribution upon cell death. +impl CellDeath { pub fn decrease_nplus( &self, distribution: &mut EcDNADistribution, @@ -256,7 +260,7 @@ mod tests { seed: u64, distribution: NonEmptyDistribtionWithNPlusCells, ) -> bool { - let ecdna = EcDNADeath; + let ecdna = CellDeath; let mut distribution = distribution.0; let mut rng = ChaCha8Rng::seed_from_u64(seed); let expected_nminus = *distribution.get_nminus(); @@ -271,7 +275,7 @@ mod tests { fn decrease_nminus_test( distribution: NonEmptyDistribtionWithNPlusCells, ) -> bool { - let ecdna = EcDNADeath; + let ecdna = CellDeath; let mut distribution = distribution.0; let expected_nminus = *distribution.get_nminus() - 1; let expected_nplus = distribution.compute_nplus(); diff --git a/src/segregation.rs b/src/segregation.rs index a3866d0..cf8e768 100644 --- a/src/segregation.rs +++ b/src/segregation.rs @@ -241,9 +241,9 @@ mod tests { #[quickcheck] fn try_from_dna_copy_test(copy: DNACopySegregatingGreatherThanOne) { - let copy_segregating: DNACopy = (copy.0).try_into().unwrap(); + let copy_segregating: DNACopy = (copy.0).into(); - assert_eq!(u16::try_from(copy.0).unwrap(), copy_segregating.get()); + assert_eq!(u16::from(copy.0), copy_segregating.get()); } #[derive(Clone, Debug)] @@ -269,7 +269,7 @@ mod tests { Deterministic {}.ecdna_segregation(copies.0, &mut rng, 0); // uneven numbers assert_eq!(k1, k2); - assert_eq!(u16::try_from(copies.0).unwrap() as u64, 2 * k1); + assert_eq!(u16::from(copies.0) as u64, 2 * k1); assert_eq!(is_uneven, IsUneven::False); }