#monte-carlo #simulation #time #stochastic #ssa #individual #numbers

sosa

Stochastic simulation algorithm (SSA) with a Monte-Carlo generating method

8 stable releases

3.0.3 Sep 4, 2023
3.0.2 Sep 1, 2023
3.0.1 Jul 20, 2023
2.0.1 Jul 20, 2023
1.0.1 Apr 7, 2023

#1552 in Algorithms

Apache-2.0 OR MIT

20KB
366 lines

The stochastic simulation algorithm (SSA) with a Monte-Carlo generating method.

Example

sosa allows using the SSA with agents that carry some individual proprieties evolving over time.

Consider for example human cells which reproduce asexually and are thought to acquire new point mutations in their genome upon cell division. We could be interested in tracking the evolution in the number of mutations over time, as cells reproduce. Moreover, cells can reproduce at different rates on average, e.g. cells carrying special mutations can reproduce faster compared to other cells. In this case, we can use sosa to perform SSA and at the same time track those mutations over time taking into account the different proliferation rates.

Note that if we are just interested in tracking the number of individuals over time, without taking into consideration the indiviual proprities of the agents, then rebop should be used instead of sosa. Find the next reaction in the system according to a Monte-Carlo generating method.

Returns

None if the maximal number of iterations have been reached or there aren't any individuals left in the total population, see SimState.

Panics

If no inidividuals left or all computed times are not normal. Compute the Gillepsie-time for all reactions. The Gillespie-time is defined as:

-ln(1 - r) / (population[i] * rates[i]) for i 0..N

where r is a random number and rates is self.0.

Returns

  • error when all computed times are infinity or 0
  • array of computed times otherwise Generates a random waiting time using the exponential waiting time with parameter lambda of Poisson StochasticProcess.

Returns

  • a waiting time of 0 if lambda is infinity,
  • a random exponential waiting time if lambda f32::is_normal,
  • infinity otherwise.
use rand_chacha::{ChaCha8Rng, rand_core::SeedableRng};

let mut rng = ChaCha8Rng::seed_from_u64(1u64);

let lambda_gr_than_zero = 0.1_f32;
assert!(exprand(lambda_gr_than_zero, &mut rng).is_sign_positive());

let lambda_zero = 0_f32;
assert!(exprand(lambda_zero, &mut rng).is_infinite());

let lambda_inf = f32::INFINITY;
assert!((exprand(lambda_inf, &mut rng) - 0.).abs() < f32::EPSILON);

Panics

When lambda is negative.

use rand_chacha::{ChaCha8Rng, rand_core::SeedableRng};

let mut rng = ChaCha8Rng::seed_from_u64(1u64);

let lambda_neg = -0.1_f32;
exprand(lambda_neg, &mut rng);

Write vector of float into new file with a precision of 4 decimals. Write NAN if the slice to write to file is empty.

Dependencies

~1.8–2.8MB
~57K SLoC