tapios / risk-networks Goto Github PK
View Code? Open in Web Editor NEWCode for risk networks: a blend of compartmental models, graphs, data assimilation and semi-supervised learning
License: Other
Code for risk networks: a blend of compartmental models, graphs, data assimilation and semi-supervised learning
License: Other
I think it'd be forward-looking to move the code that is meant to produce data and figures for the paper to move to its own directory.
The purpose of examples/
is to demonstrate how the code works. Its effectiveness is dependent on the smallness and simplicity of its contents. It is currently a bit bloated.
Consider that the code in paper/
will be frozen to a certain commit at the time of submission, even while the code continues to evolve. We will want to keep examples/
working for the current version of the code, but we will not want to change the code in paper/
, as it is meant to reproduce results reported in the paper.
It could even make sense at some point to move paper/
to a separate repository (this will make reproducing the results a bit easier, since we can record the relevant commit hash of epiforecast
). It doesn't really make sense to mix a code (which is a living, evolving thing) with the scripts that are needed to reproduce paper results (which are fixed in time and depend on a specific version of the code / git commit). We do not have to decide on this point now, however; this could be the last thing that's done before publication and opening up of the repo.
I was doing some work with the DA
and the kinetic model
. When I merge into master I found out that the class DiurnalContactInceptionRate
is no longer available. Is this reflected on the Readme.md
?
The piece of the code that needs it is:
epidemic_simulator = EpidemicSimulator(
contact_network = contact_network,
mean_contact_lifetime = mean_contact_lifetime,
contact_inception_rate = DiurnalContactInceptionRate(minimum = 2, maximum = 22),
transition_rates = transition_rates,
community_transmission_rate = community_transmission_rate,
hospital_transmission_reduction = hospital_transmission_reduction,
static_contact_interval = static_contact_interval,
health_service = health_service,
start_time = start_time
)
In a recent meeting it was discussed to include Exogenous infections, where users of the app can be infected by those without it.
This could be treated by
UserBase
(this is already in a function there)MasterEquationModelEnsemble
Can we use the rate
parameter here:
to set beta
and beta_prime
? This is a better solution than multiplying the weights of the contact network by the transmission rate --- it's simpler if the contact network weights represent only the average contact duration.
cc @odunbar
Future PRs must be created to add:
This would crash if it ever goes inside that if
:
because ContactSimulator
doesn't have a method mean_contact_duration
Many of the functions in this file:
appear to replicate functionality implemented elsewhere; for example, in
https://github.com/dburov190/risk-networks/blob/master/epiforecast/samplers.py
We should ensure that the same functions are used for the same functionality throughout the code.
After merging the new "ContactSimulator", I realized that the reproducibility issues appeared again. I performed some simulations using "simulate_idealized_NYC_epidemic_nointervention.py" and always obtain different trajectories.
@glwagner Any ideas what causes these issues (dictionaries, random access elements, ...)?
I couldn't find any additional libraries in the new contact simulator function that could cause these reproducibility issues. I also checked different "seeding" protocols, but it didn't solve the issue.
It should be possible to initialize ContactGenerator
with parameters such as:
edge_list
a list that denotes which w_ji
are non-zeroWe probably also want to have another class --- say DiurnalContactModulation
which contains information about how mean rates of contact are modulated throughout the day. This will allow us to change this function, for example, from (1 - cos(t)^4)^4 to (1-cos(t)^2)^2. This class should have a function of the form mean_contact_duration(time_of_day)
.
The ContactGenerator
should have a function --- something like contact_generator.simulate_contacts
. This method should return a time-averaged contact network in some form, whether that be a dictionary corresponding to w_ji
, or a sparse matrix. The inputs should be something like
start_time
: the time of day to start, referenced to midnightaveraging_interval
: the time-interval over which to run the contact networkThis code:
import numpy as np
import os, sys; sys.path.append(os.path.join(".."))
from epiforecast.contacts import ContactGenerator,StaticNetworkTimeSeries
np.random.seed(12345)
edge_list_file = os.path.join('..', 'data', 'networks', 'edge_list_SBM_1e3.txt')
active_edge_list_frac = 0.034
mean_contact_duration = 1.0 / 1920.0 # unit: days
# construct the generator
contact_generator= ContactGenerator(edge_list_file,
active_edge_list_frac,
mean_contact_duration)
# create the list of static networks (for 1 day)
static_intervals_per_day = int(2) #must be int
static_network_interval = 1.0/static_intervals_per_day
day_of_contact_networks = StaticNetworkTimeSeries()
contact_generator.generate_static_networks(static_network_list=day_of_contact_networks,
dt_averaging=static_network_interval)
print("Start time:", contact_generator.t0)
print("Stop time:", contact_generator.t1)
print("Birth/death times:", contact_generator.times)
produces:
Start time: 0.0
Stop time: 1.0
Birth/death times: [0. 0.00401336 0.00812667 0.01216957 0.01609839 0.02003514
0.02400519 0.02800519 0.03208136 0.03601827 0.0400094 0.04400589
0.04806824 0.05201493 0.05603721 0.06003129 0.06403446 0.06802471
0.07211888 0.07612433 0.08005874 0.08405014 0.0880175 0.09204642
0.09621028 0.10030389 0.10403115 0.1080373 0.11202207 0.11601673
0.12000081 0.1240108 0.12806288 0.13205449 0.13606364 0.14000924
0.14401498 0.1480657 0.15206046 0.15608879 0.16025658 0.16402694
0.16804135 0.17204249 0.17600942 0.18051468 0.18406326 0.18808714
0.1920246 0.1961438 0.20000846 0.20432648 0.20812067 0.21206776
0.21605961 0.22018659 0.22400221 0.22806346 0.23200533 0.23604351
0.24037212 0.24402917 0.248025 0.25205581 0.25608147 0.26002858
0.26400065 0.26802923 0.2720463 0.27604605 0.28002648 0.28402813
0.28802523 0.29200641 0.29600096 0.30000615 0.30404507 0.30801393
0.31200526 0.31615849 0.32002301 0.32401775 0.32807395 0.33204853
0.33606264 0.34004273 0.34400716 0.34807622 0.35202952 0.35602722
0.36002322 0.36406266 0.36804313 0.37211678 0.37605928 0.38013173
0.38403772 0.38813364 0.39217248 0.39609233 0.4001655 0.40404885
0.40806934 0.41232591 0.41610893 0.42074454 0.42400289 0.42821755
0.43300812 0.43605162 0.44008399 0.44448888 0.44898761 0.45224854
0.456072 0.46009569 0.4647293 0.46880973 0.4721274 0.47695183
0.48005003 0.48533513 0.48837028 0.49304084 0.49931858 0.50176534
0.50654111 0.51156591 0.51236821 0.51631232 0.53107713 0.53140987
0.53188511 0.53345727 0.55193081 0.55922347 0.56213265 0.56228804
0.56306499 0.56772741 0.57091388 0.58156196 0.60149484 0.61111747
0.64823248]
It looks like the output times do not extent very close to the final time, t1=1.0
. Why is this? Is this correct?
I input 2 uniform(0,1) states for 10 ensemble members each. No model parameters need updating, just the state.
We have true data of 1.0, variance of 1e-9 -> essentially certain data.
I expect therefore, that the output state to be an ensemble very close to 1 with small variance. Instead i get an ensemble almost identical to the input.
This snippet is the MWE (github doesn't allow me to attach .py file here):
import os, sys; sys.path.append(os.path.join(".."))
import numpy as np
import copy
from epiforecast.ensemble_adjusted_kalman_filter import EnsembleAdjustedKalmanFilter
np.random.seed(999111999)
ensemble_size = 10
#data
truth = np.array([1.0,1.0])
var = np.array([1e-9,1e-9])
cov=np.diag(var)
#generate ensemble_state
ensemble_state=np.random.random(size=(ensemble_size,2))
#pass in 'no' model parameters
ensemble_transition_rates = np.empty((ensemble_size, 0), dtype=float)
ensemble_transmission_rate = np.empty((ensemble_size, 0), dtype=float)
eakf=EnsembleAdjustedKalmanFilter()
prev_ensemble_state = copy.deepcopy(ensemble_state)
(ensemble_state[:,:],
new_ensemble_transition_rates,
new_ensemble_transmission_rate
) = eakf.update(ensemble_state[:,:],
ensemble_transition_rates,
ensemble_transmission_rate,
truth,
cov)
print(np.mean(prev_ensemble_state,axis=0))
print(" ")
print(np.mean(ensemble_state,axis=0))
Running the same epidemic simulation twice doesn't produce the same output. As an example, I used "super_simple_epidemic.py" and obtained
and (for the same parameters):
The origin of this issue may be the numba parallelization routines in "ContactSimulator", as pointed out by @glwagner:
Whats the largest network we can run this on, @agarbuno . It was killed on my laptop when attempting the 10^4 NYC data network.
contact_network
is being passed into find_observation_states
but is never used; is this a bug (and it should be used in some way) or is it safe to remove contact_network
from the argument list?
@lubo93 said:
Probably we have to go back to my original implementation where I simulated only one diurnal cycle (it's useful for system sizes N > 1e3). The current (continuous) contact simulation seems to be relatively slow.
We need to know the requirements of this system. What is the size of the network that we want to run, and how quickly do we need it to run?
@odunbar, @agarbuno, and I propose to make the size of the contact network equal to the size of the population, omitting placeholder nodes.
The ContactSimulator
will simulate all possible contacts, which includes contacts between healthcare workers and hospital beds.
Given the hospitalization state of a population, some of these edges are "phantom edges", since they represent connections between community members and people who are hospitalized, or connections between heath care workers and empty hospital beds.
In consequence:
len(contact_network) = population
n_contacts = population + hospital_beds
The HealthService
keeps track of who is in the hospital. This mapping determines the structure of the contact network, because when a node is hospitalized, their neighbors change. The HealthService
will also manage setting the weights of the contact network, since the HealthService
manages the mapping of contacts (which includes "phantom edges") to "actual" edges on the contact network.
Both the KineticModel
and the MasterEquationModelEnsemble
will simulate an epidemic on the "actual" contact network, which does not include placeholder nodes.
There are two advantages of this approach. One is that it corresponds more closely to the "true" system --- nodes represent people, and "hospitalization" changes the contact network, not the "identity" of a node.
The second advantage is corollary to the first: clinical information attached to nodes does not need to be "copied" to placeholder nodes when somebody is hospitalized (which requires extra lines of code for the correct simulation of an epidemic with correct rates). "Hospitalization" merely implies a change in the contact network, not a change in the clinical parameters of a node. This reduces the complexity of the code; the hospitalization procedure is the sole domain of the HealthService
.
@lubo93 reports that this script:
produces a KeyError
that comes from Gillespie_simple_contagion
:
In two places, here:
https://github.com/dburov190/risk-networks/blob/8d8bbfc0b61621d9ef48df22b33251dd8f386b64/epiforecast/measurements.py#L325-L330
and here:
https://github.com/dburov190/risk-networks/blob/8d8bbfc0b61621d9ef48df22b33251dd8f386b64/epiforecast/measurements.py#L424-L429
noisy_measurement
is being passed into two observe
methods as an argument; again, is this okay to delete or something actually should take place?
This function
may be misleading. It looks like beta_dict
is not a dictionary of transmission rates, but rather a dictionary of the conflated quantity
beta * w_ij
Is this correct?
In other words, the transmission rate (which we currently assume is constant across the network, and is a parameter that needs to be inferred in data assimilation) is multiplied by the average contact rate. I think it is incorrect to call this product of terms "beta".
The comments here:
are not correct. For example, the transition rate from exposed to infected is not the latent period itself, but is the inverse of the latent period. This has created confusion in the data assimilation, so we should not conflate the two.
We specify tmax
in the call to Gillespie_simple_contagion
:
Note, however, that this is not the final time of the simulation. Similar to ContactSimulator
, the Gillespie_simple_contagion
"overshoots" the end of the interval. In other words, self.time
in
is not equal to tmax
(which I think is assumed in this function).
We need to take this into account when specifying interval-wise simulations. We can probably reuse the same method / algorithm developed for the same purpose in ContactSimulator
.
If not, we don't need this check:
@lubo93 @glwagner It seems we were having some conceptual issues with the how the cycle contacts would work with the full epidemic simulator with DA.
I feel this feature is actually unnecessary as it does not work with the changing contact network such as the HealthService
and thus should be removed rather than accommodated.
Greg's new simulator means that setting the contact durations is going to be very fast, and not the bottleneck in the system so the cycling provides no benefit.
A Possibly incomplete list of tasks from the meeting
Bugs that need fixing/testing:
Imminently required features:
Other features:
The issues discussed in #58 and #60 arise because the node id's generated from 'edge_list_SBM_1e4.txt'
are not a range from 0 to population-1
.
In contrast, many initialization functions --- populate_ages
, sample_distribution
, and TransitionRates
assume that the nodes are a range from 0 to population-1
, and therefore correspond 1-1 with the node id's.
To relax the assumption that node id's are a range, we need to either use dictionaries, or the contact network itself to assign ages and clinical statistics to the population.
I propose that we use the contact network itself. In this abstraction, we will assign ages to the contact network according to a distribution. Thus populate_ages
will be replaced with assign_ages(contact_network, distribution)
. The generation of distributions for the clinical parameters, as well as the constructor for TransitionRates
will have to be modified accordingly.
I don't think this will take long. @lubo93, please advise.
Resolution of issue appears to be setting BOTH numpy and random random number generators:
import numpy as np
import random
np.random.seed(111)
random.seed(222)
We want to use a single random number generator
This is well documented in EoN
springer-math/Mathematics-of-Epidemics-on-Networks#22 (comment)
Creating MasterEquationModelEnsemble
with reduced_system=False
will result in a crash when simulate
method is executed:
Traceback (most recent call last):
File "./DA_interventions.py", line 197, in <module>
states_ensemble = master_eqn_ensemble.simulate(day)
File "../../epiforecast/risk_simulator.py", line 274, in simulate
self.y0[mm] += self.dt * self.do_step_full(t, self.y0[mm], member, closure = closure)
File "../../epiforecast/risk_simulator.py", line 226, in do_step_full
member.y_dot = member.coeffs.dot(member.y0)
File "/Library/Python/3.7/site-packages/scipy/sparse/base.py", line 363, in dot
return self * other
File "/Library/Python/3.7/site-packages/scipy/sparse/base.py", line 502, in __mul__
raise ValueError('dimension mismatch')
ValueError: dimension mismatch
Is this file
https://github.com/dburov190/risk-networks/blob/master/epiforecast/simulation_investigation.py
used anywhere? from search results it seems like it's not used anywhere. If so, we should delete it.
it should be
if observation.obs_state.size > 0:
instead
Let's implement a forward simulation example where we
KineticModel
.KineticModel
.MasterEquationEnsembleModel
.KineticModel
forward.MasterEquationEnsembleModel
forward.We should also update simulate_kinetic_model.py
to do steps 1-4 and 6 separately, without the MasterEquationModelEnsemble
.
Once this is done, we can nuke contacts.py
, ContactGenerator
, and StaticNetworkTimeSeries
.
This function
https://github.com/dburov190/risk-networks/blob/8d8bbfc0b61621d9ef48df22b33251dd8f386b64/epiforecast/populations.py#L8
is only used by a couple old examples (simple_measurements.py
and simulate_master_eqns.py
); the function should be deleted and those examples should be changed to use the new syntax with TranslationRates
directly.
There are currently two classes named Observations
defined in:
observations.py
measurements.py
In the latest simulation examples (like the ones with idealized NYC data) it seems that the latter is used; can we get rid of the former? Or merge? or rename? If we need both, what is the purpose of each one of them?
This code:
import numpy as np
from epiforecast.contacts import ContactGenerator, StaticNetworkTimeSeries
np.random.seed(12345)
second = 1 / 60 / 60 / 24 # in units of days
contact_generator = ContactGenerator(os.path.join('..', 'data', 'networks', 'edge_list_SBM_1e3.txt'), 0.034, 45 * second)
networks_time_series = StaticNetworkTimeSeries()
for dt_averaging in (1, 1/2, 1/4, 1/8):
contact_generator.generate_static_networks(static_network_list = networks_time_series,
dt_averaging = dt_averaging)
print('dt_averaging:', dt_averaging)
print(networks_time_series.community_networks_dict[-1][(0, 1)])
produces
dt_averaging: 1
0.9793103448275862
dt_averaging: 0.5
1.0
/Users/gregorywagner/opt/anaconda3/envs/risknet/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.
out=out, **kwargs)
/Users/gregorywagner/opt/anaconda3/envs/risknet/lib/python3.6/site-packages/numpy/core/_methods.py:154: RuntimeWarning: invalid value encountered in true_divide
ret, rcount, out=ret, casting='unsafe', subok=False)
dt_averaging: 0.25
nan
dt_averaging: 0.125
1.0
Thus when dt_averaging = 1/4
, there are NaNs in the weights. This can also occur for dt_averaging = 1/8
, depending on the seed, I think.
(Running the above requires epiforecast
to be on one's path, eg import os, sys; sys.path.append(os.path.join(".."))
from /examples
.)
cc @odunbar
It produces this error:
(risknet) ~/Projects/risk-networks/examples$ python3 simulate_idealized_NYC_epidemic.py
Traceback (most recent call last):
File "simulate_idealized_NYC_epidemic.py", line 138, in <module>
static_population_network = contact_network)
TypeError: __init__() got an unexpected keyword argument 'patient_capacity'
I suggest that we nuke this example, as there are several similar examples for an NYC epidemic.
As discussed elsewhere, the name epiforecast
may be a misnomer --- the main functionality of our package does not provide forecasting, due to the fact that the contact network cannot be inferred.
It's relative easy to change the name so this is not urgent. But, good ideas as they come up can be left here.
It's been mentioned that the name should evoke the concept of "risk".
Possible starting points: epirisk
, epitrace
, risktracer
, safecontact
...
README and examples need to be updated to reflect new abstractions: ContactSimulator
, KineticModel.set_mean_contact_duration
, etc.
I'm going to put together some tests using pytest
.
It's fairly simple: we'll just have a folder called /test
, and files named test_*
that implement tests according to the pytest
style.
I'll put in basic unit tests that ensure that critical functions runs, and some "regression" style tests that ensure output is reproducible.
However, we will also need tests that ensure behavior is correct. These are more difficult to come up with. We can discuss ideas on this issue. A few ideas:
" # If the node is longer in hospitalized state then the value != 'H'" --> # If the node is NO longer in hospitalized state then the value != 'H'
the use of &
(bitwise 'and' operator) is dubious; I'd use and
for consistency.
@lubo93 has pointed out that the initial fraction of active edges set, for example, here:
should be set to the equilibrium solution
initial_fraction_active_edges = lambda / (1 / mean_contact_duration + lambda)
where lambda
is the current number of mean contacts.
This means that the function outputting the mean contact rate, which is used by the ContactSimulator
here:
needs to be provided during instantiation of ContactSimulator
.
In the README, the function that implements the diurnal modulation of the mean contact rate is an input argument to the contact network generator. However, in TemporalAdjacency
, this function is hard-coded:
This design makes it difficult to experiment with different modulation functions.
This screws up syntax highlighting on GitHub, etc. T
should only be capitalized if it is a class.
Myself and @agarbuno were looking at efficacy of different tests, and noticed that we seem to get very different paced epidemics when running master equations and kinetic model.
I have run the master and kinetic equations with:
transition_rates
(set all ensembles to be the transition_rates chosen forcommunity_transmission_rate
This is achieved by setting
transition_rates_ensemble[i] = transition_rates
transmission_rate_community_ensemble[i] = transmission_rate_community
The lines are from the master equation solution, the 'X's are totals for infected, hospitalized and deceased from the kinetic model.
These two lines:
do not do anything. For example:
Is this a bug? Or unnecessary code?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.