new type for the probs

src/clap_app.rs

@@ -3,12 +3,13 @@ use clap::{ArgAction, Args, Parser, Subcommand};
 use hsc::{
     process::{ProcessOptions, Snapshot},
     subclone::{from_mean_std_to_shape_scale, Fitness},
+    ProbsPerYear,
 };
 use num_traits::{Float, NumCast};
 use sosa::{IterTime, NbIndividuals, Options};
 use std::{collections::VecDeque, ops::RangeInclusive, path::PathBuf};
 
-use crate::{AppOptions, Probs, SimulationOptionsExp, SimulationOptionsMoran};
+use crate::{AppOptions, SimulationOptionsExp, SimulationOptionsMoran};
 
 #[derive(Clone, Debug)]
 pub enum Parallel {
@@ -302,12 +303,12 @@ impl Cli {
                     save_sfs_only: cli.save_sfs_only,
                     save_population: cli.save_population,
                     tau: moran.tau,
-                    probs: Probs {
+                    probs_per_year: ProbsPerYear {
                         mu_background: moran.mu_background,
                         mu_division: moran.mu_division,
                         mu: moran.mu,
-                        asymmetric: moran.asymmetric,
                     },
+                    asymmetric: moran.asymmetric,
                 };
                 (options_moran, None)
             }
@@ -326,12 +327,12 @@ impl Cli {
                     save_sfs_only: cli.save_sfs_only,
                     save_population: cli.save_population,
                     tau: exp_moran.moran.tau,
-                    probs: Probs {
+                    probs_per_year: ProbsPerYear {
                         mu_background: exp_moran.moran.mu_background,
                         mu_division: exp_moran.moran.mu_division,
                         mu: exp_moran.moran.mu,
-                        asymmetric: exp_moran.moran.asymmetric,
                     },
+                    asymmetric: exp_moran.moran.asymmetric,
                 };
                 let options_exponential = SimulationOptionsExp {
                     gillespie_options: Options {
@@ -344,12 +345,12 @@ impl Cli {
                         verbosity,
                     },
                     tau: exp_moran.exponential.tau_exp,
-                    probs: Probs {
+                    probs_per_year: ProbsPerYear {
                         mu_background: exp_moran.exponential.mu_background_exp,
                         mu_division: exp_moran.exponential.mu_division_exp,
                         mu: exp_moran.exponential.mu_exp,
-                        asymmetric: exp_moran.exponential.asymmetric_exp,
                     },
+                    asymmetric: exp_moran.exponential.asymmetric_exp,
                 };
                 (options_moran, Some(options_exponential))
             }

src/lib.rs

@@ -40,6 +40,84 @@ pub const MAX_SUBCLONES: usize = 1200;
 /// exponential growing phase.
 pub const TIME_AT_BIRTH: f32 = 9. / 12.;
 
+#[derive(Clone, Debug)]
+pub struct ProbsPerYear {
+    /// Arrival rate of neutral background mutations per year
+    pub mu_background: f32,
+    /// Arrival rate of neutral "divisional" mutations per year
+    pub mu_division: f32,
+    /// Arrival rate of fit mutants per year
+    pub mu: f32,
+}
+
+/// Probabilities used in the simulations.
+#[derive(Clone, Debug)]
+pub enum Probs {
+    Asymmetric {
+        /// Arrival rate of fit mutants per cell per division
+        u: f32,
+        /// Probabilities per year
+        probs_per_year: ProbsPerYear,
+        /// Probability of asymmetric division per cell per division
+        asymmetric: f32,
+    },
+    Symmetric {
+        /// Arrival rate of fit mutants per cell per division
+        u: f32,
+        /// Probabilities per year
+        probs_per_year: ProbsPerYear,
+    },
+}
+
+impl Probs {
+    pub fn new(
+        mu_background: f32,
+        mu_division: f32,
+        mu: f32,
+        asymmetric: f32,
+        cells: u64,
+        verbosity: u8,
+    ) -> Probs {
+        //! ## Panics
+        //! Panics when `mu` is greater than `cells` or when asymmetric is not
+        //! within interval of 0 and 1.
+        assert!(mu <= cells as f32);
+        let probs_per_year = ProbsPerYear {
+            mu_background,
+            mu_division,
+            mu,
+        };
+        let mut u = mu / (cells as f32);
+        let probs = if (asymmetric - 0.).abs() > f32::EPSILON {
+            Probs::Asymmetric {
+                u,
+                probs_per_year,
+                asymmetric,
+            }
+        } else {
+            u *= 0.5;
+            Probs::Symmetric { u, probs_per_year }
+        };
+        if verbosity > 0 {
+            println!("probs {:#?}", probs);
+        }
+        dbg!(u);
+        assert!((0f32..1.).contains(&u), "Invalid u: u>=0 and u<1");
+        assert!(
+            (0f32..=1.).contains(&asymmetric),
+            "Invalid asymmetric: asymmetric>=0 and asymmetric<=1"
+        );
+        probs
+    }
+
+    pub fn is_asymmetric(&self) -> bool {
+        match self {
+            Probs::Symmetric { .. } => false,
+            Probs::Asymmetric { .. } => true,
+        }
+    }
+}
+
 pub fn write2file<T: std::fmt::Display>(
     data: &[T],
     path: &Path,
@@ -78,6 +156,7 @@ pub fn write2file<T: std::fmt::Display>(
 
 #[cfg(test)]
 mod tests {
+    use super::*;
     use quickcheck::{Arbitrary, Gen};
     use std::num::NonZeroU8;
 
@@ -90,4 +169,46 @@ mod tests {
             LambdaFromNonZeroU8(lambda.get() as f32)
         }
     }
+
+    #[test]
+    #[should_panic]
+    fn panic_asymmetric_neg_cells_test() {
+        Probs::new(1.1, 1.1, 0.1, -0.1, 10, 0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn panic_asymmetric_inf_cells_test() {
+        Probs::new(1.1, 1.1, 0.1, f32::INFINITY, 10, 0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn panic_asymmetric_nan_cells_test() {
+        Probs::new(1.1, 1.1, 0.1, f32::NAN, 10, 0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn panic_mu_gr_cells_test() {
+        Probs::new(1.1, 1.1, 12., 0., 10, 0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn panic_mu_neg_cells_test() {
+        Probs::new(1.1, 1.1, -0.1, 0., 10, 0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn panic_mu_inf_cells_test() {
+        Probs::new(1.1, 1.1, f32::INFINITY, 0., 10, 0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn panic_mu_nan_cells_test() {
+        Probs::new(1.1, 1.1, f32::NAN, 0., 10, 0);
+    }
 }