//! This module contains functions to simulate population demographics (not including genetics) using forward-direction population-level simulations. Populations are represented by matrices and vectors containing demographic and behavioral information.
/// This struct represents a population by a (typically integer) vector. Each value of the vector represents the number of individuals in a lifestage present in the population. For example, a population with 15 hatchlings, 8 juveniles, and 30 adults could be represented by this vector: `[40, 20, 100]`. This struct is meant to contain this type of information. The data is stored as f64 (floating point) values to accommodate conditions when decimal populations are desirable and facilitate calculations that may not return integer values.
use rand::distributions::Uniform;
use rand::{thread_rng, Rng};
use statrs::statistics::{self, Statistics};
#[derive(Clone, Debug)]
pub struct PopulationVector {
vector: Vec<f64>,
lifestage_count: u8,
}
impl PopulationVector {
/// Create a new Population Vector instance by inputting a vector containing f64 values.
pub fn new(vector: Vec<f64>) -> PopulationVector {
PopulationVector {
lifestage_count: vector.len() as u8,
vector: vector,
}
}
/// Return the value stored at a specifc index in the Population Vector based on inputed
/// integer (u32). The first value is 0.
pub fn get_value_at_index(&self, index: u32) -> Option<&f64> {
self.vector.get(index as usize)
}
// Return full vector stored in the Population Vector as a `Vec<f64>`.
pub fn get_vector(&self) -> &Vec<f64> {
return &self.vector;
}
// Return the number of items stored in the Population Vector instance. This is used to prevent
// errors in calculations that require matching vector/matrix lengths.
pub fn get_lifestage_count(&self) -> u8 {
return self.lifestage_count;
}
pub fn total_individuals(&self) -> f64 {
return self.vector.iter().sum();
}
}
/// This struct represents the likelihood of different lifestages of an organism to survive, grow,
/// and reproduce to the next lifestage (or stay in the same lifestage) by storing decimal (f64) values in a matrix (a vector of vectors, denoted
/// `Vec<Vec<f64>>`). Each sub-vector (row) represents a lifestage, and each item (column)
/// represents the yearly proportion of individuals recruited to that lifestage. The first row typically represents newborns, seedlings, etc.
///
/// ## Examples
/// Index [1][1] is the number of lifestage 1 individuals that were in lifestage 1 the
/// previous year.
///
/// Index [1][2] is the number of lifestage 1 indidviduals that were in lifestage 2 the year
/// before.
///
/// Indeg [2][1] is the number of lifestage 2 individuals who were lifetage 1 individuals the
/// previous year.
///
/// ...etc...
///
/// An example matrix could look like this:
/// [0][0][0.1][0.2]
/// [0.6][0][0][0]
/// [0][0.8][0][0]
/// [0][0][0.8][0.94]
pub struct PopulationMatrix {
matrix: Vec<Vec<f64>>,
lifestage_count: u8,
}
impl PopulationMatrix {
/// This function builds a Population Matrix from a square vector of vectors (Vec<Vec<f64>>), ensuring that it contains a consistent
/// number of lifestages across all inputted Lifestage Survival Vectors and in the number of
/// inputted Lifestage Survival Vectors. If these conditions are not met, it will return an
/// error message.
pub fn build(input: Vec<Vec<f64>>) -> Result<PopulationMatrix, &'static str> {
if input.len() == input[0].len() {
for count in 1..input.len() {
if input[count].len() != input[count - 1].len() {
return Err("All sub-vectors must be of matching lengths to construct a population matrix.");
}
}
return Ok(PopulationMatrix {
lifestage_count: input.len() as u8,
matrix: input,
});
} else {
return Err("Number of items in lifestages must match number of inputted sub-vectors.");
}
}
/// Returns the number of listages represented in the Population Matrix, useful for calculations requiring
/// matching numbers of lifestages.
pub fn get_lifestage_count(&self) -> u8 {
return self.lifestage_count;
}
/// Returns the full matrix of matrices stored in the Population Matrix.
pub fn get_matrix(&self) -> &Vec<Vec<f64>> {
return &self.matrix;
}
/// Given an input of a PopulationMatrix and a PopulationVector with the same number of items
/// in their `matrix` and `vector` values respectively, the function will return a
/// PopulationMatrix.
/// ##Use
/// This function takes a population matrix (formatted as a PopulationMatrix struct) and a population vector (formatted as a PopulationVector struct) as an imput, calculating the sum of the components of each column of the vector within the population matrix and muliplying it by the value of the population vector at the same index. This is essentially a column-wise matrix-vector product. It then returns a new PopulationVector.
///
/// This can be visually thought of as such:
///
///[ 40] [0 , 0 , 0.1 ] [ 10]
///[ 20] x [0.6, 0.8, 0 ] = [ 40]
///[100] [0 , 0.8, 0.95] [111]
///
///This calculation is common in Population Variability Analysis (PVA) wherein each row (LifeStageSurvivalVector) represents the probability of recruitment into that life stage over the course of a year. By multiplying a matrix of these probabilities by a vector containing the current population, a researcher can estimate the following year's population.
///
/// ## Errors
/// This function will return an Err('static str') if the number of rows or items within rows in the matrix is not equal to the number of items in the population vector.
///
/// Although this is theoretically impossible, the program could also panic if it recieves an out-of-bounds index request for the population vector. However, the function checks for this earlier in order to return a useful error code and prevent other mistakes, so should never occur.
/// # Examples
/// ```
/// use ecolysis_core::populations::population_level_simulation::{PopulationMatrix, PopulationVector}; // import relevant structs
/// let popvector = PopulationVector::new(vec![150.0, 200.0, 33.0]); // create a population vector type
/// let popmatrix = PopulationMatrix::build(vec![
/// vec![0.25, 0.001, 0.75],
/// vec![0.3, 0.4346, 0.002],
/// vec![0.98, 0.66, 0.161]
/// ]).unwrap(); // create a Population Matrix type
/// let new_popvector = popmatrix.project_vector(&popvector); // project the vector by the matrix
/// println!("{:?}", new_popvector.unwrap().get_vector()); // print results
/// ```
pub fn project_vector(
&self,
vector: &PopulationVector,
) -> Result<PopulationVector, &'static str> {
if self.lifestage_count != vector.get_lifestage_count() {
return Err(
"the length of inputted population matrix and population vector do not match.",
);
}
let mut new_population_vector: Vec<f64> = vec![];
let mut storage: f64 = 0.0;
for lifestage in &self.matrix {
for (count2, item) in lifestage.iter().enumerate() {
storage += item
* vector
.get_value_at_index(count2 as u32)
.expect("Index out of bounds.");
}
new_population_vector.push(storage);
storage = 0.0;
}
Ok(PopulationVector::new(new_population_vector))
}
}
/// The PvaDeterministicPopulation struct stores population data for deterministic PVA models, allowing PVA operations to be performed by simply calling
/// functions on an instance. It contains the following:
/// - A Population Vector representing the initial population size.
/// - A Population Matrix contains data on the survival rates
/// and recruitment rates of verious lifestages.
pub struct PvaDeterministicPopulation {
initial_population: PopulationVector,
projection_matrix: PopulationMatrix,
}
impl PvaDeterministicPopulation {
/// Return a Result enum containing a new PvaDeterministicPopulation instance with the input of a Population Vector and a Population Matrix. /// one Population Matrix in the latter vector.
/// # Errors
/// Will return `Err<'static str>` if the lengths of the Population Vector the Matrix do not match.
pub fn build(
initial_population: PopulationVector,
matrix: PopulationMatrix,
) -> Result<PvaDeterministicPopulation, &'static str> {
let expected_lifestage_length = initial_population.get_lifestage_count();
if matrix.get_lifestage_count() != expected_lifestage_length {
return Err("Population vector size does not match matrices.");
}
Ok(PvaDeterministicPopulation {
initial_population: initial_population,
projection_matrix: matrix,
})
}
/// Return a Result enum containing a new PvaDeterministicPopulation instance with the input of a vector containing f64 values (a population vector) and a square set of vector of vectors containing f64 values (a population matrix).
/// # Errors
/// Will return `Err<'static str>` if the lengths of the Population Vector the Matrix do not match.
/// ```
///use ecolysis_core::populations::population_level_simulation::PvaDeterministicPopulation;
///let new_population = PvaDeterministicPopulation::build_from_vectors(vec![12.0, 55.0, 172.0],
///vec![
///vec![0.0, 0.0, 0.9],
///vec![0.6, 0.9, 0.0],
///vec![0.0, 0.95, 0.99]
///]);
/// ```
pub fn build_from_vectors(
initial_population: Vec<f64>,
matrix: Vec<Vec<f64>>,
) -> Result<PvaDeterministicPopulation, &'static str> {
PvaDeterministicPopulation::build(
PopulationVector::new(initial_population),
PopulationMatrix::build(matrix)?,
)
}
/// Return a Result enum containing a PvaDeterministicOutput type that holds the output of a determinisitc simulation
/// given the number of simulation steps to perform (as a u32). This function performs the
/// actual simulation.
///
/// The function may panic if the Populatiom Matrix is not square or the
/// lengths of Population Vector and Population Matrix do not match, although this situation should
/// be prevented by checks when building a PVA Population instance.
/// ```
///use ecolysis_core::populations::population_level_simulation::PvaDeterministicPopulation;
///let new_population = PvaDeterministicPopulation::build_from_vectors(vec![12.0, 55.0, 172.0],
///vec![
///vec![0.0, 0.0, 0.9],
///vec![0.6, 0.9, 0.0],
///vec![0.0, 0.95, 0.99]
///]).unwrap();
///let simulation_output = new_population.deterministic_projection(100);
///simulation_output.print_output();
/// ```
pub fn deterministic_projection(&self, steps: u32) -> PvaDeterministicOutput {
let mut active_vector = self.initial_population.clone();
let mut result: Vec<PopulationVector> = Vec::new();
for _ in 1..=steps {
active_vector = self.projection_matrix.project_vector(&active_vector).expect("This error should not be possible. Mismatched Vector and Matrix lengths, or non-square Matrix. Please file a bug report.");
result.push(active_vector.clone());
}
return PvaDeterministicOutput::new(result);
}
}
/// This enum stores the output of Population Viability Analysis operations performed by the PVA
/// Deterministic Population struct.
pub struct PvaDeterministicOutput {
result: Vec<PopulationVector>,
}
impl PvaDeterministicOutput {
// Create a new PvaDeterministicOutput struct from a vector of PopulationVectors
// (Vec<PopulationVector>).
pub fn new(simulation_output: Vec<PopulationVector>) -> PvaDeterministicOutput {
return PvaDeterministicOutput {
result: simulation_output,
};
}
/// Print a CSV containing the output of each simulation step to the console.
pub fn print_output(&self) {
let mut string = String::new();
for (counti, i) in self.result.iter().enumerate() {
for (countj, j) in i.get_vector().iter().enumerate() {
string.push_str(&j.to_string());
if countj + 1 < i.get_vector().len() {
string.push_str(", ");
}
}
if counti + 1 < self.result.len() {
string.push_str("\n");
}
}
println!("{}", string);
}
/// Return a vector of Population Vectors representing all the data from
/// each step of the simulation for a determinisitc model. Each item is the output of an
/// iteration of the simulation. The first item is the first iteration, the last itemn is
/// the last iteration of the simulation.
pub fn return_typed_output(&self) -> &Vec<PopulationVector> {
&self.result
}
/// Return a vecotr of vectors, containing floating point values (Vec<Vec<f64>>), representing
/// all the data from each step of the simulation for a determinisic model. Each item in the outer vector is the output of an
/// iteration of the simulation. The first item is the first iteration, the last itemn is the
/// last iteration of the simulation. Each of the sub-vectors is a de-typed population vector,
/// representing the demographics of a population.
pub fn return_numerical_output(&self) -> Vec<Vec<f64>> {
let mut num_vec: Vec<Vec<f64>> = Vec::new();
for i in &self.result {
num_vec.push(i.get_vector().clone());
}
num_vec
}
}
pub struct PvaStochasticPopulation {
initial_population: PopulationVector,
projection_matrices: Vec<PopulationMatrix>,
}
impl PvaStochasticPopulation {
/// Return a Result enum containing a new PvaStochasticPopulation instance with the input of a Population Vector and a vector of Population Matrices.
/// # Errors
/// Will return `Err<'static str>` if the lengths of the Population Vector and the Matrices do not match.
pub fn build(
initial_population: PopulationVector,
matrices: Vec<PopulationMatrix>,
) -> Result<PvaStochasticPopulation, &'static str> {
let expected_lifestage_length = initial_population.get_lifestage_count();
for matrix in &matrices {
if matrix.get_lifestage_count() != expected_lifestage_length {
return Err("Population vector size does not match matrices.");
}
}
Ok(PvaStochasticPopulation {
initial_population: initial_population,
projection_matrices: matrices,
})
}
/// Return a Result enum containing a new PvaStochasticPopulation instance with the input of a vector containing f64 values (a population vector) and list of square sets of vector of vectors containing f64 values (a vector of population matrices).
/// # Errors
/// Will return `Err<'static str>` if the lengths of the Population Vector the Matrices do not match.
/// ```
///use ecolysis_core::populations::population_level_simulation::PvaStochasticPopulation;
///let new_population = PvaStochasticPopulation::build_from_vectors(vec![12.0, 55.0, 172.0],
///vec![
///vec![
///vec![0.0, 0.0, 0.9],
///vec![0.6, 0.9, 0.0],
///vec![0.0, 0.95, 0.99]
///],
///vec![
///vec![0.0, 0.0, 0.7],
///vec![0.7, 0.4, 0.0],
///vec![0.0, 0.93, 0.95]
///],
///vec![
///vec![0.0, 0.0, 0.8],
///vec![0.7, 0.91, 0.0],
///vec![0.0, 0.93, 0.98]
///],
///]);
/// ```
pub fn build_from_vectors(
initial_population: Vec<f64>,
matrices: Vec<Vec<Vec<f64>>>,
) -> Result<PvaStochasticPopulation, &'static str> {
let mut formatted_matrices: Vec<PopulationMatrix> = vec![];
for matrix in matrices {
let temp_matrix = match PopulationMatrix::build(matrix) {
Ok(valid_matrix) => valid_matrix,
Err(e) => return Err(e),
};
formatted_matrices.push(temp_matrix)
}
PvaStochasticPopulation::build(
PopulationVector::new(initial_population),
formatted_matrices,
)
}
pub fn stochastic_projection(&self, steps: u32, iterations: u32) -> PvaStochasticOutput {
let mut active_vector = self.initial_population.clone();
let mut step_result: Vec<PopulationVector> = Vec::new();
let mut iteration_result: Vec<Vec<PopulationVector>> = Vec::new();
let mut rng = thread_rng();
let rand_generator = Uniform::new_inclusive(1, self.projection_matrices.len() as u32 - 1);
for _ in 1..=iterations {
for _ in 1..=steps {
active_vector = self.projection_matrices[rng.sample(rand_generator) as usize].project_vector(&active_vector).expect("This error should not be possible. Mismatched Vector and Matrix lengths, or non-square Matrix. Please file a bug report.");
step_result.push(active_vector.clone());
}
iteration_result.push(step_result.clone());
step_result = vec![];
}
return PvaStochasticOutput::new(iteration_result);
}
}
pub struct PvaStochasticOutput {
result: Vec<Vec<PopulationVector>>,
steps: u32,
iterations: u32,
}
impl PvaStochasticOutput {
pub fn new(simulation_output: Vec<Vec<PopulationVector>>) -> PvaStochasticOutput {
return PvaStochasticOutput {
steps: simulation_output[1].len() as u32,
iterations: simulation_output.len() as u32,
result: simulation_output,
};
}
/// Refactors result data to be organized by time step, in the form of a 2-dimensional vector
/// containing population vectors. The output is structured such that the out outer layer is
/// the time-step of the simulation and the inner vector is the simulation iteration.
fn organize_by_time_step(&self) -> Vec<Vec<PopulationVector>> {
let mut new_result: Vec<Vec<PopulationVector>> = Vec::new();
let mut temp_vec: Vec<PopulationVector> = Vec::new();
for step in 0..self.steps {
for item in &self.result {
temp_vec.push(item[step as usize].clone());
}
new_result.push(temp_vec);
temp_vec = vec![];
}
new_result
}
/// Calculates the average population across simulation iterations, expressed as an f64 value, at each time step in the siulation.
pub fn result_mean_total_population(&self) -> Vec<f64> {
let working_result = self.organize_by_time_step();
let mut mean_vec: Vec<f64> = Vec::new();
for step in working_result {
mean_vec.push(step.iter().map(|item| item.total_individuals()).mean());
}
mean_vec
}
/// Returns a specific iteration of the simulation using its index.
pub fn iteration_output_by_index(&self, index: usize) -> &Vec<PopulationVector> {
return &self.result[index];
}
/// This function exports simulation output while retaining the type structure: a 2 dimensional
/// matrix of vectors containing Population Vector types. The outer vector is the simulation
/// iteration, while the inner vector is the time-step.
pub fn result_typed_output(&self) -> &Vec<Vec<PopulationVector>> {
return &self.result;
}
/// This function converts all the output data into a 3-dimensional nested vector of f64
/// values.
pub fn result_vec_output(&self) -> Vec<Vec<&Vec<f64>>> {
return self
.result
.iter()
.map(|item1| {
item1
.iter()
.map(|iter2| iter2.get_vector())
.collect::<Vec<_>>()
})
.collect::<Vec<_>>();
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn matrix_multiplication() {
let popvector = PopulationVector::new(vec![40.0, 20.0, 100.0]);
let mut lifestage_recruit: Vec<Vec<f64>> = vec![vec![0.0, 0.0, 0.1]];
lifestage_recruit.push(vec![0.6, 0.8, 0.0]);
lifestage_recruit.push(vec![0.0, 0.8, 0.95]);
let popmatrix =
PopulationMatrix::build(lifestage_recruit).expect("Invalid population matrix."); // Calculated vectors for tests: vec![vec![37.5, 0.15, 1500], vec![20000, 86.92, 400], vec![2178, 3861, 132]];
let result: Vec<f64> = vec![10.0, 40.0, 111.0];
assert_eq!(
&result,
&popmatrix
.project_vector(&popvector)
.unwrap()
.get_vector()
.into_iter()
.map(|x| { (x * 10.0).round() / 10.0 }) // Rounding is necessary to get rid of floating point errors.
.collect::<Vec<_>>(),
);
}
#[test]
fn matrix_invalid_matrix_length() {
assert!(PopulationMatrix::build(vec![vec![0.5, 0.7, 0.3], vec![0.1, 0.11, 0.6]]).is_err());
}
#[test]
fn matrix_unmatched_lifestage_lengths() {
assert!(PopulationMatrix::build(vec![
vec![0.5, 0.7, 0.3],
vec![0.1, 0.11, 0.6],
vec![0.2, 0.91]
])
.is_err());
}
#[test]
fn pva_simple_matrix_projection_test() {
let population_vec = PopulationVector::new(vec![40.0, 20.0, 100.0]);
let matrix = PopulationMatrix::build(vec![
vec![0.0, 0.0, 0.1],
vec![0.6, 0.8, 0.0],
vec![0.0, 0.8, 0.95],
])
.unwrap();
let population = PvaDeterministicPopulation::build(population_vec, matrix).unwrap();
let result = population.deterministic_projection(8);
result.print_output();
let correct_result = vec![24.9, 50.8, 273.5];
let mut temp_vec: Vec<f64> = Vec::new();
let mut clean_output: Vec<Vec<f64>> = Vec::new();
for i in result.return_typed_output() {
for j in i.get_vector() {
temp_vec.push(((j * 10.0 as f64).round()) / 10.0);
}
clean_output.push(temp_vec);
temp_vec = vec![];
}
assert_eq!(correct_result, clean_output[clean_output.len() - 1])
}
#[test]
fn pva_stochastic_output_summary_mean_test() {
let test_sim_output = PvaStochasticOutput::new(vec![
vec![
PopulationVector::new(vec![1.0, 2.0, 1.0]),
PopulationVector::new(vec![5.0, 25.0, 2.0]),
],
vec![
PopulationVector::new(vec![1.0, 1.0, 1.0]),
PopulationVector::new(vec![3.0, 9.0, 10.0]),
],
vec![
PopulationVector::new(vec![3.0, 1.0, 1.0]),
PopulationVector::new(vec![8.0, 8.0, 8.0]),
],
]);
assert_eq!(
test_sim_output.result_mean_total_population(),
vec![4.0, 26.0]
);
}
#[test]
fn pva_stochastic_simulation_test() {
let population_vec = PopulationVector::new(vec![40.0, 20.0, 100.0]);
let matrices: Vec<PopulationMatrix> = vec![
PopulationMatrix::build(vec![
vec![0.0, 0.0, 0.1],
vec![0.6, 0.8, 0.0],
vec![0.0, 0.8, 0.95],
])
.unwrap(),
PopulationMatrix::build(vec![
vec![0.0, 0.0, 0.2],
vec![0.5, 0.9, 0.0],
vec![0.0, 0.7, 0.99],
])
.unwrap(),
PopulationMatrix::build(vec![
vec![0.0, 0.0, 0.05],
vec![0.4, 0.5, 0.0],
vec![0.0, 0.3, 0.7],
])
.unwrap(),
PopulationMatrix::build(vec![
vec![0.0, 0.0, 0.2],
vec![0.5, 0.75, 0.0],
vec![0.0, 0.6, 0.9],
])
.unwrap(),
PopulationMatrix::build(vec![
vec![0.0, 0.0, 0.5],
vec![0.9, 0.9, 0.0],
vec![0.0, 0.95, 0.99],
])
.unwrap(),
];
let population = PvaStochasticPopulation::build(population_vec, matrices).unwrap();
let binding = population.stochastic_projection(100, 10);
let simulation_output = binding.result_vec_output();
for (count_i, i) in simulation_output.iter().enumerate() {
for (count_j, j) in simulation_output.iter().enumerate() {
if count_i != count_j {
//println!("{:?}", i);
//println!("{:?}", j);
assert_ne!(i, j);
}
}
}
}
}