igraph/misc/epidemics.rs
1//! SIR epidemics (`igraph_epidemics.h`).
2//!
3//! `igraph_sir_init` / `igraph_sir_destroy` are not wrapped: they are memory
4//! plumbing. [`SirRun`] owns the four vectors of a run and frees them on drop.
5
6use crate::{error::Result, ffi::*, igraph_call};
7use std::mem::MaybeUninit;
8
9/// The outcome of one run of the SIR epidemic model (the owned counterpart of
10/// `igraph_sir_t`), as returned by [`Graph::sir`](crate::Graph::sir).
11///
12/// The four vectors have the same length: entry `k` describes the state of
13/// the population right after the `k`-th event (an infection or a recovery).
14/// Entry `0` is the initial state, at time `0`, with a single infected
15/// individual; the last entry is the moment the last infected individual
16/// recovers. At every step `susceptible + infected + recovered` equals the
17/// number of vertices of the graph.
18#[derive(Debug, Clone, PartialEq)]
19pub struct SirRun {
20 /// The (continuous) times of the events, starting at `0.0`, non-decreasing.
21 pub times: Vec<f64>,
22 /// Number of susceptible individuals after each event.
23 pub susceptible: Vec<i64>,
24 /// Number of infected individuals after each event.
25 pub infected: Vec<i64>,
26 /// Number of recovered individuals after each event.
27 pub recovered: Vec<i64>,
28}
29
30impl SirRun {
31 /// Number of recorded states (events plus the initial state).
32 pub fn len(&self) -> usize {
33 self.times.len()
34 }
35
36 /// Whether no state was recorded (never the case for runs produced by igraph).
37 pub fn is_empty(&self) -> bool {
38 self.times.is_empty()
39 }
40
41 /// The time at which the epidemic died out (the time of the last event).
42 pub fn duration(&self) -> f64 {
43 self.times.last().copied().unwrap_or(0.0)
44 }
45
46 /// The *final size* of the epidemic: how many individuals were ever
47 /// infected (they have all recovered at the end).
48 pub fn final_size(&self) -> i64 {
49 self.recovered.last().copied().unwrap_or(0)
50 }
51
52 /// The largest number of simultaneously infected individuals.
53 pub fn peak_infected(&self) -> i64 {
54 self.infected.iter().copied().max().unwrap_or(0)
55 }
56
57 /// Iterates over `(time, susceptible, infected, recovered)` tuples.
58 pub fn states(&self) -> impl Iterator<Item = (f64, i64, i64, i64)> + '_ {
59 (0..self.len()).map(move |k| {
60 (
61 self.times[k],
62 self.susceptible[k],
63 self.infected[k],
64 self.recovered[k],
65 )
66 })
67 }
68}
69
70impl From<igraph_sir_t> for SirRun {
71 fn from(sir: igraph_sir_t) -> Self {
72 // The fields are owned igraph vectors: they are freed when `sir` drops.
73 Self {
74 times: sir.times.to_vec(),
75 susceptible: sir.no_s.to_vec(),
76 infected: sir.no_i.to_vec(),
77 recovered: sir.no_r.to_vec(),
78 }
79 }
80}
81
82/// Owns an `igraph_vector_ptr_t` filled by `igraph_sir` with heap allocated
83/// `igraph_sir_t` items, and frees everything (items included) on drop.
84struct SirList(igraph_vector_ptr_t);
85
86impl SirList {
87 fn new() -> Self {
88 crate::error::ensure_init();
89 let mut raw = MaybeUninit::<igraph_vector_ptr_t>::uninit();
90 crate::error::check(unsafe { igraph_vector_ptr_init(raw.as_mut_ptr(), 0) })
91 .expect("igraph failed to allocate a pointer vector");
92 Self(unsafe { raw.assume_init() })
93 }
94
95 /// Moves the runs out of the list, leaving null pointers behind.
96 fn take_runs(&mut self) -> Vec<SirRun> {
97 let n = unsafe { igraph_vector_ptr_size(&self.0) } as usize;
98 let mut runs = Vec::with_capacity(n);
99 for k in 0..n {
100 let slot = unsafe { self.0.stor_begin.add(k) };
101 let item = unsafe { *slot } as *mut igraph_sir_t;
102 if item.is_null() {
103 continue;
104 }
105 // Take ownership of the four vectors, then free the C allocation
106 // holding the struct itself (allocated with IGRAPH_CALLOC).
107 let sir = unsafe { std::ptr::read(item) };
108 unsafe {
109 igraph_free(item.cast());
110 *slot = std::ptr::null_mut();
111 }
112 runs.push(SirRun::from(sir));
113 }
114 runs
115 }
116}
117
118impl Drop for SirList {
119 fn drop(&mut self) {
120 // Frees any run not moved out (e.g. never reached on error paths,
121 // where igraph already nulled them), then the pointer vector itself.
122 drop(self.take_runs());
123 unsafe { igraph_vector_ptr_destroy(&mut self.0) };
124 }
125}
126
127impl igraph_t {
128 /// Runs `num_simulations` stochastic SIR (susceptible–infected–recovered)
129 /// epidemics on the graph.
130 ///
131 /// Each individual (vertex) is susceptible, infected or recovered;
132 /// recovered individuals are immune. A susceptible vertex with `n`
133 /// infected neighbors becomes infected at rate `n * beta`, an infected one
134 /// recovers at rate `gamma` (both are rates of exponential
135 /// distributions, so the model runs in continuous time, as a Gillespie
136 /// simulation). Every simulation starts with a single, uniformly chosen,
137 /// infected vertex and stops when no infected vertex is left. It uses the
138 /// thread's default random number generator: seed it with
139 /// [`rng::seed`](crate::rng::seed) for reproducible runs.
140 ///
141 /// Edge directions are ignored (with a warning) for directed graphs, so a
142 /// directed graph with a pair of opposite edges `u → v`, `v → u` counts
143 /// as a multigraph and is rejected.
144 /// Seeded runs are reproducible, also when several threads simulate at
145 /// the same time: each thread has its own default generator.
146 ///
147 /// See also [`Graph::famous`](crate::Graph::famous) and the random graph models of
148 /// [`crate::games`] for contact networks, and
149 /// [`set_interruption_handler`](crate::misc::set_interruption_handler) to
150 /// cancel long simulations.
151 ///
152 /// Binds [`igraph_sir`](https://igraph.org/c/html/latest/igraph-Processes.html#igraph_sir).
153 /// Time complexity: `O(num_simulations * (|V| + |E| log |V|))`.
154 ///
155 /// # Errors
156 /// [`ErrorKind::InvalidValue`](crate::ErrorKind::InvalidValue) when the
157 /// graph has no vertices or is not simple, when `beta < 0`, `gamma <= 0`
158 /// or `num_simulations == 0`. The computation can be stopped by an
159 /// [interruption handler](crate::misc::set_interruption_handler)
160 /// ([`ErrorKind::Interrupted`](crate::ErrorKind::Interrupted)).
161 ///
162 /// # Examples
163 ///
164 /// ```
165 /// use igraph::prelude::*;
166 ///
167 /// rng::seed(42)?;
168 /// let ring = Graph::ring(10, false, false, true)?;
169 /// let runs = ring.sir(2.0, 1.0, 5)?;
170 /// assert_eq!(runs.len(), 5);
171 /// for run in &runs {
172 /// for (_, s, i, r) in run.states() {
173 /// assert_eq!(s + i + r, 10); // the population is conserved
174 /// }
175 /// assert_eq!(*run.infected.last().unwrap(), 0); // the epidemic dies out
176 /// }
177 /// # Ok::<(), igraph::Error>(())
178 /// ```
179 ///
180 /// Epidemics spread further on the karate club network when the
181 /// infection rate grows:
182 ///
183 /// ```
184 /// use igraph::prelude::*;
185 ///
186 /// rng::seed(1)?;
187 /// let club = Graph::famous("Zachary")?;
188 /// let mean_final_size = |beta: f64| -> igraph::Result<f64> {
189 /// let runs = club.sir(beta, 1.0, 300)?;
190 /// Ok(runs.iter().map(|r| r.final_size() as f64).sum::<f64>() / 300.0)
191 /// };
192 /// assert!(mean_final_size(0.05)? < mean_final_size(1.0)?);
193 /// # Ok::<(), igraph::Error>(())
194 /// ```
195 pub fn sir(&self, beta: f64, gamma: f64, num_simulations: usize) -> Result<Vec<SirRun>> {
196 let mut list = SirList::new();
197 igraph_call!(igraph_sir(
198 self,
199 beta,
200 gamma,
201 num_simulations as igraph_int_t,
202 &mut list.0
203 ))?;
204 Ok(list.take_runs())
205 }
206}