Draft of example 1: simple infection model

examples/basic_infection/README.md

@@ -0,0 +1,46 @@
+# Example - Ixa - Simple transmission model
+
+Infection process in a homogeneous closed population with infection process independently from contact process and disease progression.
+
+- Specifications
+
+  - Time-varying infectiousness
+  - Disease progression: symptom, including hospitalizations.
+
+- Contact patterns
+  - Homogeneous but contacts should be reduced when hospitalized
+  - Closed population
+
+- Outputs
+  - Disease and infection status: incidence and in periodic reports that represent the number of people on each status.
+
+## Motivation
+This example shows how to implement a simple transmission model where symptoms and infection are treated separately. This model requires some basic characteristics from Ixa:
+
+  - person properties: infection and symptom status are specified in person properties.
+  - random number generator: transmission requires the evaluation and generation of random numbers.
+  - input/output handling
+  - basic reports: periodic report where reports are scheduled to output prevalence values of person properties, incidence report for disease status where changes in disease status are outputed.
+  - contact manager: a way to keep track of contacts for each person.
+  - plans: plans are necessary to schedule events related to infection progression, symptom progression, transmission, and reports.
+
+## Model description
+This model reproduces transmission dynamics of a pathogen in a closed homogeneous population. Individuals are modeled as perons with two properties related to their infection and health status. Infection status specifies the course of infection for an individual exposed to the pathogen. An individual is susceptible to infection. When a susceptible individual is exposed to the pathogen through contact with an infectious person, the susceptible individual is considered infected with the pathogen but unable to infect others. After the latent period, the infected individual becomes infectious and is able to infect others through contact. The number of new infections is defined in this example by the basic reproductive number.
+
+Infected individuals may develop symptoms after exposure based on the incubation period. Those who develop symptoms could require hospitalization, which is defined in the model by a probability of hospitalization, time from symptom onset to hospitalization, and hospitalization duration. During hospitalization, the contacts from an individual should be restricted. In this simple example, and because this is a homogeneous population, we assume that hospitalization effectively reduces all contacts to zero.
+
+
+## How to run the model
+To run the model:
+
+`cargo build -r`
+
+Run a single example scenario:
+
+`target/release/cfa-eosim-facemask -i test/input/config.yaml -o test/output/`
+
+Run multiple example scenarios with 4 threads:
+
+`target/release/cfa-eosim-facemask -i test/input/config_multi.yaml -o test/output/ -t 4`
+
+Inspect the output in R using the script `test/output/plot_output.R`

examples/basic_infection/src/lib.rs

@@ -0,0 +1 @@
+pub mod sir;

examples/basic_infection/src/main.rs

@@ -0,0 +1,239 @@
+use std::{fs::File, path::Path};
+
+use cfa_eosim_facemask::sir::{
+    death_manager::DeathManager,
+    global_properties::{
+        DeathRate, HospitalizationDelay, HospitalizationDuration, IncubationPeriod,
+        InfectiousPeriod, InitialInfections, LatentPeriod, MaxDays, Population,
+        ProbabilityHospitalized, ProbabilitySymptoms, SymptomaticPeriod, R0,
+    },
+    infection_manager::InfectionManager,
+    infection_seeder::InfectionSeeder,
+    periodic_report::{PeriodicReport, PeriodicStatus},
+    person_property_report::{PersonPropertyChange, PersonPropertyReport},
+    population_loader::PopulationLoader,
+    transmission_manager::TransmissionManager,
+};
+use clap::Parser;
+use eosim::{
+    context::Context,
+    global_properties::GlobalPropertyContext,
+    random::RandomContext,
+    reports::{get_channel_report_handler, Report, ReportsContext},
+};
+use serde_derive::{Deserialize, Serialize};
+use threadpool::ThreadPool;
+use tokio::runtime::Handle;
+use tokio::sync::mpsc::{self, Sender};
+
+#[derive(Debug, Parser)]
+struct SirArgs {
+    /// Input config file
+    #[arg(short, long)]
+    input: String,
+    /// Output directory
+    #[arg(short, long)]
+    output: String,
+    /// Number of threads
+    #[arg(short, long, default_value_t = 1)]
+    threads: u8,
+}
+
+#[derive(Debug, Serialize, Deserialize, Clone, Copy)]
+struct Parameters {
+    population: usize,
+    r0: f64,
+    infectious_period: f64,
+    latent_period: f64,
+    incubation_period: f64,
+    max_days: usize,
+    probability_symptoms: f64,
+    symptomatic_period: f64,
+    hospitalization_duration: f64,
+    probability_hospitalized: f64,
+    hospitalization_delay: f64,
+    initial_infections: usize,
+    random_seed: u64,
+    death_rate: f64,
+}
+
+#[derive(Debug, Serialize, Deserialize, Clone, Copy)]
+struct Scenario {
+    scenario: usize,
+}
+
+#[derive(Debug, Serialize, Deserialize)]
+#[serde(untagged)]
+enum Config {
+    Single(Parameters),
+    Multiple(Vec<Parameters>),
+}
+
+fn setup_context(context: &mut Context, parameters: &Parameters) {
+    // Set up parameters in simulation
+    context.set_global_property_value::<Population>(parameters.population);
+    context.set_global_property_value::<R0>(parameters.r0);
+    context.set_global_property_value::<InfectiousPeriod>(parameters.infectious_period);
+    context.set_global_property_value::<LatentPeriod>(parameters.latent_period);
+    context.set_global_property_value::<IncubationPeriod>(parameters.incubation_period);
+    context.set_global_property_value::<SymptomaticPeriod>(parameters.symptomatic_period);
+    context
+        .set_global_property_value::<HospitalizationDuration>(parameters.hospitalization_duration);
+    context.set_global_property_value::<HospitalizationDelay>(parameters.hospitalization_delay);
+    context
+        .set_global_property_value::<ProbabilityHospitalized>(parameters.probability_hospitalized);
+    context.set_global_property_value::<ProbabilitySymptoms>(parameters.probability_symptoms);
+    context.set_global_property_value::<InitialInfections>(parameters.initial_infections);
+    context.set_global_property_value::<DeathRate>(parameters.death_rate);
+    context.set_global_property_value::<MaxDays>(parameters.max_days);
+
+    // Set up RNG
+    context.set_base_random_seed(parameters.random_seed);
+
+    // Add reports
+    context.add_component::<PersonPropertyReport>();
+    context.add_component::<PeriodicReport>();
+
+    // Add model components
+    context.add_component::<PopulationLoader>();
+    context.add_component::<InfectionManager>();
+    context.add_component::<TransmissionManager>();
+    context.add_component::<InfectionSeeder>();
+    context.add_component::<DeathManager>();
+}
+
+pub fn get_bounded_channel_report_handler<T: Report, S>(
+    sender: Sender<(S, T::Item)>,
+    id: S,
+) -> impl FnMut(T::Item) + 'static
+where
+    T::Item: serde::Serialize + Send + 'static,
+    S: serde::Serialize + Send + Copy + 'static,
+{
+    move |item| {
+        let sender = sender.clone();
+        let id = id;
+        futures::executor::block_on(async move {
+            if let Err(e) = sender.send((id, item)).await {
+                panic!("Receiver being closed, failed to send item: {:?}", e);
+            }
+        });
+    }
+}
+
+fn run_single_threaded(parameters_vec: Vec<Parameters>, output_path: &Path) {
+    let output_file = File::create(output_path.join("person_property_report.csv"))
+        .expect("Could not create output file.");
+    let output_periodic_file = File::create(output_path.join("periodic_report.csv"))
+        .expect("Could not create output periodic file.");
+    for (scenario, parameters) in parameters_vec.iter().enumerate() {
+        let mut writer_builder = csv::WriterBuilder::new();
+        // Don't re-write the headers
+        if scenario > 0 {
+            writer_builder.has_headers(false);
+        }
+        let mut writer = writer_builder.from_writer(
+            output_file
+                .try_clone()
+                .expect("Could not write to output file"),
+        );
+        let mut periodic_writer = writer_builder.from_writer(
+            output_periodic_file
+                .try_clone()
+                .expect("could not write to output file for periodic report"),
+        );
+        // Set up and execute context
+        let mut context = Context::new();
+        context.set_report_item_handler::<PersonPropertyReport>(move |item| {
+            if let Err(e) = writer.serialize((Scenario { scenario }, item)) {
+                eprintln!("{}", e);
+            }
+        });
+        context.set_report_item_handler::<PeriodicReport>(move |item| {
+            if let Err(e) = periodic_writer.serialize((Scenario { scenario }, item)) {
+                eprintln!("{}", e);
+            }
+        });
+        setup_context(&mut context, parameters);
+        context.execute();
+        println!("Scenario {} completed", scenario);
+    }
+}
+
+async fn run_multi_threaded(parameters_vec: Vec<Parameters>, output_path: &Path, threads: u8) {
+    let output_file = File::create(output_path.join("person_property_report.csv"))
+        .expect("Could not create output file.");
+    let output_periodic_file = File::create(output_path.join("periodic_report.csv"))
+        .expect("Could not create output periodic file.");
+    let pool = ThreadPool::new(threads.into());
+    let (sender, mut receiver) = mpsc::channel::<(Scenario, PersonPropertyChange)>(100000);
+    let (periodic_sender, mut periodic_receiver) =
+        mpsc::channel::<(Scenario, PeriodicStatus)>(100000);
+
+    let handle = Handle::current();
+
+    for (scenario, parameters) in parameters_vec.iter().enumerate() {
+        let sender = sender.clone();
+        let periodic_sender = periodic_sender.clone();
+        let parameters = *parameters;
+        let handle = handle.clone();
+        pool.execute(move || {
+            let _guard = handle.enter();
+            // Set up and execute context
+            let mut context = Context::new();
+            context.set_report_item_handler::<PersonPropertyReport>(
+                get_bounded_channel_report_handler::<PersonPropertyReport, Scenario>(
+                    sender,
+                    Scenario { scenario },
+                ),
+            );
+            context.set_report_item_handler::<PeriodicReport>(
+                get_bounded_channel_report_handler::<PeriodicReport, Scenario>(
+                    periodic_sender,
+                    Scenario { scenario },
+                ),
+            );
+            setup_context(&mut context, &parameters);
+            context.execute();
+            println!("Scenario {} completed", scenario);
+        });
+    }
+    drop(sender);
+    drop(periodic_sender);
+
+    // Write output from main thread
+    let mut person_property_writer = csv::Writer::from_writer(output_file);
+    let mut periodic_writer = csv::Writer::from_writer(output_periodic_file);
+    loop {
+        tokio::select! {
+            Some(item) = receiver.recv() => {
+                person_property_writer.serialize(item).unwrap();
+            },
+            Some(item) = periodic_receiver.recv() => {
+                periodic_writer.serialize(item).unwrap();
+            },
+            else => break,
+        }
+    }
+}
+
+#[tokio::main]
+async fn main() {
+    // Parse args and load parameters
+    let args = SirArgs::parse();
+    let config_file = File::open(&args.input)
+        .unwrap_or_else(|_| panic!("Could not open config file: {}", args.input));
+    let config: Config = serde_yaml::from_reader(config_file).expect("Could not parse config file");
+    let output_path = Path::new(&args.output);
+
+    match config {
+        Config::Single(parameters) => run_single_threaded(vec![parameters], output_path),
+        Config::Multiple(parameters_vec) => {
+            if args.threads <= 1 {
+                run_single_threaded(parameters_vec, output_path)
+            } else {
+                run_multi_threaded(parameters_vec, output_path, args.threads).await;
+            }
+        }
+    }
+}

examples/basic_infection/src/sir/death_manager.rs

@@ -0,0 +1,44 @@
+use eosim::{
+    context::{Component, Context},
+    global_properties::GlobalPropertyContext,
+    people::PersonId,
+    person_properties::PersonPropertyContext,
+    random::RandomContext,
+};
+use rand::Rng;
+
+use super::{global_properties::DeathRate, person_properties::DiseaseStatus};
+
+pub struct DeathManager {}
+
+impl Component for DeathManager {
+    fn init(context: &mut Context) {
+        context
+            .observe_person_property_changes::<DiseaseStatus>(handle_person_disease_status_change);
+    }
+}
+
+eosim::define_random_id!(DeathRandomId);
+
+pub fn handle_person_disease_status_change(
+    context: &mut Context,
+    person_id: PersonId,
+    _: DiseaseStatus,
+) {
+    let disease_status = context.get_person_property_value::<DiseaseStatus>(person_id);
+    if matches!(disease_status, DiseaseStatus::I) {
+        schedule_death_check(context, person_id);
+    }
+}
+
+pub fn schedule_death_check(context: &mut Context, person_id: PersonId) {
+    let death_rate = *context
+        .get_global_property_value::<DeathRate>()
+        .expect("Death Rate not specified");
+    let mut rng = context.get_rng::<DeathRandomId>();
+    let should_die = rng.gen::<f64>() < death_rate;
+    drop(rng);
+    if should_die {
+        context.set_person_property_value::<DiseaseStatus>(person_id, DiseaseStatus::D);
+    }
+}

examples/basic_infection/src/sir/global_properties.rs

@@ -0,0 +1,25 @@
+eosim::define_global_property!(R0, f64);
+
+eosim::define_global_property!(InfectiousPeriod, f64);
+
+eosim::define_global_property!(LatentPeriod, f64);
+
+eosim::define_global_property!(SymptomaticPeriod, f64);
+
+eosim::define_global_property!(IncubationPeriod, f64);
+
+eosim::define_global_property!(HospitalizationDuration, f64);
+
+eosim::define_global_property!(ProbabilityHospitalized, f64);
+
+eosim::define_global_property!(HospitalizationDelay, f64);
+
+eosim::define_global_property!(ProbabilitySymptoms, f64);
+
+eosim::define_global_property!(Population, usize);
+
+eosim::define_global_property!(InitialInfections, usize);
+
+eosim::define_global_property!(DeathRate, f64);
+
+eosim::define_global_property!(MaxDays, usize);