Rc0.0.2 #17
Rc0.0.2 #17
Conversation
stisim/parameters.py
Outdated
| import sciris as sc | ||
| from .settings import options as sso # For setting global options | ||
| from . import misc as ssm | ||
| from .data import loaders as ssdata |
There was a problem hiding this comment.
I know data loading is not set up yet -- perhaps comment this import until this functionality is enabled to avoid import errors.
stisim/people.py
Outdated
| from . import utils as ssu | ||
| from . import settings as sss | ||
| from .data import loaders as ssdata |
There was a problem hiding this comment.
Perhaps comment this import until this functionality is enabled to avoid import errors -- original comment was about this import line in people.py, thought it probably applies to import in parameters.py too.
|
|
||
| class Analyzer: |
There was a problem hiding this comment.
I think eventually we want this to be a class instance like it is in hpvsim and covasim?
|
Looks great to me! |
| # %% Helper functions to create popdicts | ||
|
|
||
| def set_static(new_n, existing_n=0, pars=None, f_ratio=0.5): |
There was a problem hiding this comment.
Suggested docstring:
"""
Set static population characteristics that do not change over time.
This function can be used when adding new births, in which case the existing popsize can be given as `existing_n`.
Arguments:
new_n (int): Number of new individuals to add to the population.
existing_n (int, optional): Number of existing individuals in the population. Default is 0.
pars (dict, optional): Dictionary of parameters. Default is None.
f_ratio (float, optional): Female ratio in the population. Default is 0.5.
Returns:
uid (ndarray, int): unique identifiers for the individuals.
female (ndarray, bool): array indicating whether an individual is female or not.
debut (ndarray, bool): array indicating the debut value for each individual.
"""
| uid = np.arange(existing_n, existing_n+new_n, dtype=sss.default_int) | ||
| female = np.random.choice([True, False], size=new_n, p=f_ratio) | ||
| n_females = len(ssu.true(female)) |
There was a problem hiding this comment.
I think np.count_zeros may be slightly faster than len(...). This is tiny detail, may not make a difference if this function is not called multiple times.
There was a problem hiding this comment.
Agree, seems like np.random.choice([True, False], size=new_n, p=f_ratio), is equivalent to but probably slower than np.random.rand(new_n)>f_ratio
| age_data_prob /= age_data_prob.sum() # Ensure it sums to 1 | ||
| age_bins = ssu.n_multinomial(age_data_prob, n) # Choose age bins | ||
| if dt_round_age: |
There was a problem hiding this comment.
Wondering if this age rounding/correction may be moved to a different place (a standalone function perhaps?), one where dt has already passed validation checks and fulfills the condition that age_data_range[age_bins] / dt are integer numbers, which is the implicit assumption I think.
| pars = sc.objdict() | ||
|
|
||
| # Population parameters | ||
| pars['n_agents'] = 10e3 # Number of agents |
There was a problem hiding this comment.
Following up on the discussion about parameters. Cliff is right that a nested dictionary is a good way to represent parameters and provides a lot of flexibility. The nested dict can be exported as json (sim.export_pars), or parsed and exported in a different way (eg, to excel) via a user-defined function. However what I think we need is a richer description of a parameter, rather than only the value, we could also consider a short, human readable description, units if applicable, valid range (may come handy for calibration later). All these additional parameter attributes can be accommodated in a nested dict, but what would need to change how the rest of the classes access parameter values.
There was a problem hiding this comment.
import stisim
import sciris as sc
pars = sc.objdict()
pars['population'] = sc.objdict()
# Define parameter
pars['population']['n_agents'] = stisim.Parameter(name="n_agents", default_value=10000, dtype=int, description="Number of agents, that is, the size of the agent population.")
# Access value
pars['population']['n_agents']()
# or
pars['population']['n_agents'].value There was a problem hiding this comment.
I love this idea and think it could be very useful for calibration
There was a problem hiding this comment.
Yes, strongly in favor of a richer parameter description!
There was a problem hiding this comment.
Have created issue #22 to track this enhancement.
| raise AlreadyRunError('Simulation already complete (call sim.initialize() to re-run)') | ||
|
|
||
| # Update states, modules, partnerships | ||
| self.people.update_states(sim=self) # This runs modules |
There was a problem hiding this comment.
Maybe we could call this method people.update_living_states`?
I don't think it's necessary to pass the whole sim to people update_states, it only uses the time step size dt and the time index ti.
In people it'd would look like:
def update_living_states(self, dt, ti):
""" Perform all state updates at the current timestep """
self.age[~self.dead] += dt
self.dead[self.ti_dead <= ti] = True
returnI would also suggest to call the module updates/stepping function directly in sim.step, not inside of people.update ... I think this would help with the goal of keeping code modular. And it may be useful to have all the entities that need updating at each time step, and in which order, directly in the place that orchestrates how everything evolves (the sim class).
In sim it would look something like this:
def step()
self.people.update_living_states(sim.dt, sim.ti) # People age or eventually die
self.update_health_states()
self.update_networks()
self.update_connector()
self.update_modules()
# Tidy up
self.ti += 1
if self.ti == self.npts:
self.complete = True
def update_health_states(self):
for module in self.modules.values(): # Other health states progress (disease states or pregnancy stages)
module.update_states(self)
def update_networks(self): # could also be called step_networks, etc
# Update networks -- currently done inside people.update ,
# Perform network updates
for lkey, layer in self.people.networks.items():
layer.update(self.people)
def update_modules(self): # could also be called step_modules, etc
for module in self.modules.values():
module.make_new_cases(self)
module.update_results(self)If the logic in any of the update_ or step_ becomes more complex, then it'd be contained to its local scope rather than obscuring the logic of the main simulation step() method.
There was a problem hiding this comment.
agree with this logic, but think we can do this as a next step after merging?
There was a problem hiding this comment.
Agree, but instead of update_living_states could be update_vital_dynamics or similar.
There was a problem hiding this comment.
I also think we should consider make Demographics a separate module to facilitate additional flexibility.
Also removed an unnecessary dependency and renamed a variable in test_base.py.
| State('uid', sss.default_int), | ||
| State('age', sss.default_float), | ||
| State('female', bool, False), | ||
| State('debut', sss.default_float), |
There was a problem hiding this comment.
Should "debut" be a base State or part of the network? If only used for partnerships / network participation, consider moving to the Network class. But if used more broadly, e.g. for exogeneous fertility, okay to keep here.
There was a problem hiding this comment.
One solution would be for "debut" to be determined by the network itself. For example by the first time index at which an individual forms any network edge. So it would be a State propoerty of the base that would be set by the Network layer(s).
There was a problem hiding this comment.
well this begs the question whether debut is an attribute of a person or of the network. would we expect people to have different ages of debut for different networks. and then what does that mean for a network which is not a sexual network?
| sim.people.add_module(self) | ||
|
|
||
| # Pick some initially infected agents | ||
| self.set_prognoses(sim, np.random.choice(sim.people.uid, sim.pars[self.name]['initial'])) |
There was a problem hiding this comment.
This "bakes in" the np.random.choice. Why not allow the module to choose the people?
| """ Add new cases of module, through transmission, incidence, etc. """ | ||
| pars = sim.pars[self.name] | ||
| for k, layer in sim.people.networks.items(): | ||
| if k in pars['beta']: |
There was a problem hiding this comment.
I don't follow this pars['beta'] thing, what's going on here?
| if k in pars['beta']: | ||
| rel_trans = (sim.people[self.name].infected & ~sim.people.dead).astype(float) | ||
| rel_sus = (sim.people[self.name].susceptible & ~sim.people.dead).astype(float) | ||
| for a, b, beta in [[layer['p1'], layer['p2'], pars['beta'][k][0]], [layer['p2'], layer['p1'], pars['beta'][k][1]]]: |
There was a problem hiding this comment.
The asymmetry in beta is nice!
There was a problem hiding this comment.
Perhaps more efficient to first find module-discordant edges and then only sample transmission one-way?
| raise AlreadyRunError('Simulation already complete (call sim.initialize() to re-run)') | ||
|
|
||
| # Update states, modules, partnerships | ||
| self.people.update_states(sim=self) # This runs modules |
There was a problem hiding this comment.
Agree, but instead of update_living_states could be update_vital_dynamics or similar.
| # >>> sim0.initialize() | ||
| # >>> sim1.initialize() | ||
| # >>> sim0.run() | ||
| # >>> sim1.run() |
There was a problem hiding this comment.
But since these are different instances, why would the results not be the same in either case? (Sorry, feels like a dumb question!)
| raise AlreadyRunError('Simulation already complete (call sim.initialize() to re-run)') | ||
|
|
||
| # Update states, modules, partnerships | ||
| self.people.update_states(sim=self) # This runs modules |
There was a problem hiding this comment.
I also think we should consider make Demographics a separate module to facilitate additional flexibility.
|
|
||
| return | ||
|
|
||
| def update_connectors(self): |
There was a problem hiding this comment.
Feels like this belongs at the sim level.
There was a problem hiding this comment.
isn't it currently at the sim level?
There was a problem hiding this comment.
Yes, but there is a TODO comment on the following line suggesting to move this to people. I like it where it is ;)
| # Set beta. The first entry represents transmission risk from infected p1 -> susceptible p2 | ||
| # Need to be careful to get the ordering right. The set-up here assumes that in the simple | ||
| # sexual network, p1 is male and p2 is female. In the maternal network, p1=mothers, p2=babies. | ||
| hiv.pars['beta'] = {'simple_sexual': [0.0008, 0.0004], 'maternal': [0.2, 0]} |
There was a problem hiding this comment.
Note that in hpv_network, p1 is female.
Adding in basic *sim functionality, with a few changes/exceptions:
people.networksstores instances of theNetworkclass (previously, this was called people.contacts and Networks were called Layers)population.pyis now exclusively home to types of Network. The idea is that we would provide a few defaults but people could easily create their ownmodules.pyoperates as beforeparameters.py-- much of this is still TBClocationfunctionality is not here yet -- whilst this is not particularly tricky to add, it's also data-intensive and maybe not needed just yet, especially since it's still TBD what data we need and what should have responsibility for getting these data -- e.g. if we really just need the hpvsim data then perhaps we could call hpv.data to get it? This would introduce a dependency between stisim and hpvsim that we probably don't want, so perhaps we need to move it all to starsim?Tests are in test_base.py.