Skip to content

Commit

Permalink
v0.22.0
Browse files Browse the repository at this point in the history
  • Loading branch information
@fraterenz
fraterenz committed Dec 15, 2023
1 parent 90e1e5a commit ab15ac1
Show file tree
Hide file tree
Showing 8 changed files with 168 additions and 246 deletions.
3 changes: 3 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
# Changelog

## 0.22.0
change the way we save the data at the end of the simulations to simplify subsampling.

## 0.21.0
Assign idx such that it corresponds to the seed used.

Expand Down
6 changes: 3 additions & 3 deletions Cargo.lock
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

4 changes: 2 additions & 2 deletions Cargo.toml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "ecdna-evo"
version = "0.21.0"
version = "0.22.0"
edition = "2021"
repository = "https://github.com/fraterenz/ecdna-evo"
description = "Evolutionary models of extra-chromosomal DNA (ecDNA)"
Expand All @@ -15,11 +15,11 @@ anyhow = "1.0"
chrono = "0.4"
clap = { version = "4", features = ["derive"] }
csv = "1.1"
ecdna-lib = "3.0.2"
indicatif = { version = "*", features = ["rayon"] }
serde = { version = "1.0.163", features = ["derive"] }
serde_json = "1.0"
sosa = "2"
ecdna-lib = "2.0.6"
rand = { version = "0.8.0", features = ["small_rng"] }
rand_chacha = "0.3.1"
rand_distr = "0.4.0"
Expand Down
196 changes: 91 additions & 105 deletions src/app.rs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
use anyhow::Context;
use ecdna_evo::{distribution::SamplingStrategy, RandomSampling, ToFile};
use ecdna_lib::distribution::EcDNADistribution;
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use sosa::{
Expand All @@ -10,13 +11,39 @@ use std::{
path::{Path, PathBuf},
};

use crate::NB_RESTARTS;
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 {
Expand All @@ -42,128 +69,87 @@ impl Dynamics {
+ RandomSampling,
REACTION: std::fmt::Debug,
{
let run_helper =
|sampling: Option<(NbIndividuals, SamplingStrategy)>,
path2dir: &Path,
state: &mut CurrentState<NB_REACTIONS>,
process: &mut P| {
let stream = idx as u64 * NB_RESTARTS;
let mut rng = ChaCha8Rng::seed_from_u64(self.seed);
rng.set_stream(stream);

// clone the initial state in case we restart
let stop_reason = simulate(
state,
rates,
possible_reactions,
process,
&self.options,
&mut rng,
);

// restart, this can happen when high death rates
// TODO: bug, see issue #104 and also biases the results since
// we are kind of conditionning upon survival
//
// while (stop == StopReason::NoIndividualsLeft
// || stop == StopReason::AbsorbingStateReached)
// && j <= NB_RESTARTS
// {
// // restart process as well
// process_copy = process.clone();
// let stream = idx as u64 * NB_RESTARTS + self.seed + j;
// if self.options.verbosity > 1 {
// println!(
// "Restarting with stream {} because {:#?}",
// stream, stop
// );
// }
// rng.set_stream(stream);

// // clone the initial state in case we restart
// stop = simulate(
// &mut initial_state.clone(),
// rates,
// possible_reactions,
// &mut process_copy,
// &self.options,
// &mut rng,
// );
// j += 1;
// }

if self.options.verbosity > 1 {
println!("stop reason: {:#?}", stop_reason);
}
if self.options.verbosity > 1 {
println!("{:#?}", process);
}

if let Some(sampling) = sampling {
if self.options.verbosity > 1 {
println!("subsample with {} cells", sampling.0);
}
process
.save(&path2dir.join("before_subsampling"), rates, idx)
.with_context(|| {
format!(
"Cannot save run {} before subsampling",
idx
)
})
.unwrap();
let stream = idx as u64;
let mut rng = ChaCha8Rng::seed_from_u64(self.seed);
rng.set_stream(stream);

if self.options.verbosity > 0 {
println!("Subsampling");
}
process.random_sample(&sampling.1, sampling.0, rng);
// after subsampling we need to update the state
process.update_state(state);
}
};
let stop_reason = simulate(
initial_state,
rates,
possible_reactions,
&mut process,
&self.options,
&mut rng,
);

if self.options.verbosity > 1 {
println!("{:#?}", process);
if self.options.verbosity > 0 {
println!("stop reason: {:#?}", stop_reason);
}

if let Some(sampling) = sampling {
let path2dir =
&self.path2dir.join((sampling.at.len() - 1).to_string());
for (i, sample_at) in sampling.at.iter().enumerate() {
// save only at last sample
let last_sample = i == sampling.at.len() - 1;
run_helper(
Some((*sample_at, sampling.strategy)),
path2dir,
initial_state,
&mut process,
);
// happens with birth-death processes
let no_individuals_left =
initial_state.population.iter().sum::<u64>() == 0u64;
if last_sample || no_individuals_left {
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(&path2dir.join("after_subsampling"), rates, idx)
.save(
&self.path2dir.join(format!(
"{}cells/before_subsampling",
sample_at
)),
rates,
idx,
)
.with_context(|| {
format!(
"Cannot save run {} for timepoint {}",
idx, i
"Cannot save run {} before subsampling",
idx
)
})
.unwrap();
}
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!("Subsampling");
}
}
} else {
if self.options.verbosity > 1 {
println!("Nosubsampling");
}
run_helper(None, &self.path2dir, initial_state, &mut process);
process
.save(&self.path2dir, rates, idx)
.with_context(|| format!("Cannot save run {}", idx))
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(())
}
Expand Down
Loading

0 comments on commit ab15ac1

Please sign in to comment.