Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Births and deaths closes issue #51 #65

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
Open

Conversation

confunguido
Copy link
Collaborator

@confunguido confunguido commented Oct 24, 2024

This PR provides an example of births, deaths, and aging. This example highlights the need for a tool to look up people by person properties (1+) as well as how to sample from those.

Being alive
To handle births and deaths, I created a person property named Alive. This is set to true when a person is added to the simulation and to false when they die.

Aging
I defined a person property as Birth, instead of age to determine someone's current risk category, which is based on age. For instance, in 1+ years time, someone may switch from one group to the other, and that should be considered. Age and Birth here are defined in years but as f64, so that would account for 1.5 years).

fn get_person_age(&mut context, person_id: PersonId) -> f64 {
  let birth_time = self.get_person_property(person_id, Birth);
 (context.get_current_time() - birth_time) / 365.0
}

How to get current alive population
For every operation that requires sampling an individual, estimating current population, etc., we need to filter by people who are alive. I created a function called get_population_by_property which finds the total number of people by a specified property.

fn get_population_by_property<T: PersonProperty + 'static> (&mut self, property: T, value: T::Value) -> usize {
  let mut population_counter = 0;
  for i in 0..self.get_current_population() {
     let person = self.get_person_id(i);
     if self.get_person_property(person_id, property) == value { 
        population_counter += 1;
     } 
  }
  population_counter
}

Should alive be treated differently when filtering people

I wonder whether we should make being alive a special case for some of these filters. For instance, if I want to sample someone to infect, it is obvious I want someone alive, and if the FOI depends on the current population, it is also obvious that it requires the current alive population. Another example of this is periodic reports, where we want to get the number of people per person property, but we want to get those alive. Of course, not all models include aging, so maybe this isn't an issue.

Canceling plans after a dead event
I implemented a rough approach to store infection events and cancel them after a person dies in the simulation. This is done in the infection manager given that it's the only place where people may need to cancel scheduled plans after death.

@confunguido confunguido force-pushed the gce_births_deaths branch 2 times, most recently from cf367a7 to b24cdd1 Compare October 26, 2024 00:15
@confunguido confunguido marked this pull request as ready for review October 26, 2024 00:21
Copy link
Collaborator

@k88hudson-cfa k88hudson-cfa left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few questions, mostly this looks pretty reasonable to me – the age / age group varying with time will make it tricky for indexing to be useful for a lot of stuff here

For a more general API plan cancelling – we could consider something like cancel_plans_for_person (which would require you to also register it, I think), but what you're doing here doesn't seem too bad. What are your thoughts after implementing?

define_person_property_with_default!(InfectionStatusType, InfectionStatus, InfectionStatus::S);

define_person_property!(Birth, f64);
define_person_property!(Alive, bool);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems like an appropriate place to use define_person_property_with_default!, no?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, good point. I guess people by default are alive XD

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@k88hudson-cfa I have changed the way age is implemented. Now is back to being u8 and age scheduled to update every year. For canceling plans, I think creating something like cancel_plans_for_person would be fine. I think it wasn't too hard to keep track of the plans but made me think of how to keep track of which plans specifically need to be canceled when a person is removed from the simulation.


pub trait ContextPopulationExt {
fn create_new_person(&mut self, birth_time: f64) -> PersonId;
fn attempt_death(&mut self, person_id: PersonId);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why attempt if this is just setting death?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Honestly, this is nothing else that me being unable to write something like kill_person. I will change that.

schedule_death(context);
});
}
// Cancel all plans
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you mean to leave this as TODO?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not really. I think I realized later that I don't actually know who should take care of canceling all plans. I ended up implementing it in infection_manager and just registering for death events. I think that approach isn't too bad, so I just forgot to remove from population manager.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok thanks. Yeah, just remove the comment then

});
}

context.add_plan(1.1, move |context| {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why 1.1?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wanted to kill right after infection so that recovery would have to be canceled. It could also be at 1.0, I guess.

<T as PersonProperty>::Value: PartialEq,
{
let mut people_vec = Vec::<PersonId>::new();
for i in 0..self.get_current_population() {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would probably use filter_map here but this is perfectly fine

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Noted. I'm not sure I should work too much on this sample since it will be replaced once we land the PR to query people, right?

Copy link
Contributor

@ChiragKumar9 ChiragKumar9 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The only change I am tied to is the comment about the mathematical inexactness of this simulation compared to what is described in the readme.

#[allow(clippy::cast_precision_loss)]
let next_attempt_time = context.get_current_time()
+ context.sample_distr(TransmissionRng, Exp::new(foi).unwrap())
/ population_size as f64;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This needs a note that what is modeled here is not exact to what is described in the readme. In particular, the population changes over time, so the infection attempt times must also change with the population size. Although this expression uses the current population size, it is not updating to changes in the population size. A rigorous treatment of this would require rejection sampling.

@@ -0,0 +1,143 @@
# Birth and death processes

This example describes the process of loading a population, adding new births and deaths. It has a simple sir model that represents infections from a constant force of infection that differs by age category (0,2m, 6m, 65y).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
This example describes the process of loading a population, adding new births and deaths. It has a simple sir model that represents infections from a constant force of infection that differs by age category (0,2m, 6m, 65y).
This example describes the process of loading a population, adding new births and deaths. It has a simple disease model that represents infections from a constant force of infection that differs by age category (0,2m, 6m, 65y).

The exposure process is initialized by scheduling an initial infection attempt for each age group. This depends on the force of infection defined for each age group, and the current size of the group, which only includes people who are alive. Individuals who are infected by the pathogen may recover based on a recovery period (`t + infection_period`) specified at initialization. These recovery times are scheduled at the time of infection as a `plan`. The `id` for this plan is saved in a data structure that holds the current plans for recovery scheduled for each individual. Once the person recovers, these plans are removed from the data structure. However, if the person dies before recovering, these plans are canceled at the time of death. The infection status of recovered individuals remains as recovered for the rest of the simulation, which will stop when the plan queue is empty.

# Population manager
At initialization, the population manager adds people to the simulation as defined in the input parameters file, and initializes each person with two person properties: `Birth` time and `Alive` status. Birth time is estimated as a random number between 0 and -100 to represent ages from 0 to 100 years. To trait `get_person_age` can return a person's age based on their time of birth.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
At initialization, the population manager adds people to the simulation as defined in the input parameters file, and initializes each person with two person properties: `Birth` time and `Alive` status. Birth time is estimated as a random number between 0 and -100 to represent ages from 0 to 100 years. To trait `get_person_age` can return a person's age based on their time of birth.
At initialization, the population manager adds people to the simulation as defined in the input parameters file, and initializes each person with two person properties: `Birth` time and `Alive` status. Birth time is estimated as a random number between 0 and -100 to represent ages from 0 to 100 years. The method `get_person_age` (implemented as a trait extension on `context`) can return a person's age based on their time of birth.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also think an explanation of why a trait extension is used here would be valuable as it is likely the user's first time of seeing this pattern in their ixa learning.

```

## Age Groups
Age groups are defined in the population manager as an `enum`. These groups are determined for each person using the trait `get_person_age_group`. This function estimates the current age group based on the time of the simulation and the time of birth. In this example, the force of infection varies for each of the three age groups defined below. Hence, a hash map contains the force of infection for each of these groups and is saved as a global property.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Age groups are defined in the population manager as an `enum`. These groups are determined for each person using the trait `get_person_age_group`. This function estimates the current age group based on the time of the simulation and the time of birth. In this example, the force of infection varies for each of the three age groups defined below. Hence, a hash map contains the force of infection for each of these groups and is saved as a global property.
Age groups are defined in the population manager as an `enum`. These groups are determined for each person using the method `get_person_age_group`. This function estimates the current age group based on the time of the simulation and the time of birth. In this example, the force of infection varies for each of the three age groups defined below. Hence, a hash map contains the force of infection for each of these groups and is saved as a global property.

use ixa::define_rng;
use ixa::global_properties::ContextGlobalPropertiesExt;
use ixa::people::{ContextPeopleExt, PersonId, PersonPropertyChangeEvent};
use ixa::plan;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like we are adopting a convention of importing the exact struct needed, so should this be use ixa::plan::Id?

output_df <- read_csv(file.path(dir, "incidence.csv")) |>
dplyr::filter(infection_status == "I") |>
group_by(time, age_group) |>
mutate(inf = n()) |>
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should be a summarize. Imagine you have two infections that occur at the same time -- you want to count those each once, not twice.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah that makes sense. Though, unlikely with this approach.

pub fn init(context: &mut Context) {
let parameters = context.get_global_property_value(Parameters).clone();

let foi_map = parameters
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not do this in the parameters_loader?

Copy link
Collaborator

@k88hudson-cfa k88hudson-cfa left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This mostly looks good and I like the age incrementing approach – just have a question / alternative idea for scheduling by birthday (rather than by individual) I think we should discuss.

Looks like @ChiragKumar9 has some comments too

context.set_person_property(person_id, Age, prev_age + 1);
let next_age_event = context.get_current_time() + 365.0;
context.add_plan(next_age_event, move |context| {
schedule_aging(context, person_id);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't it kind of weird that everyone in the simulation gets a birthday of t = 0 (well, -365*Age) except people who are born after? I defer to you on what is more natural on the simulation side of things but I'd probably do this a bit differently....

Rather than scheduling aging events for each person, I'd probably index by birthday and then check who's age needs to be incremented. Alternatively, you could assume everyone's birthday is t=0 and just increment everyone's age every 365 days. What do you think?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I modified this by adding a random birthday from 0 to 365 when creating a new person at the beginning of the simulation. Let me know if this looks good for you. I addressed @ChiragKumar9 comments as well but let me know if I missed anything.

@k88hudson-cfa
Copy link
Collaborator

LGTM! @ChiragKumar9 needs to approve to merge

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants