Skip to content

Commit

Permalink
some modifications for epidemics
Browse files Browse the repository at this point in the history
  • Loading branch information
@benmaier
benmaier committed May 19, 2020
1 parent 4ab7bd9 commit e51b358
Show file tree
Hide file tree
Showing 8 changed files with 305 additions and 33 deletions.
72 changes: 67 additions & 5 deletions _tacoma/FlockworkPModel.cpp
View file
Expand Up @@ -90,6 +90,41 @@ void
}
}

void FlockworkPModel::update_aggregated_network(
double const &t,
vector < pair < size_t, size_t > > &e_in,
vector < pair < size_t, size_t > > &e_out
)
{
if (verbose)
cout << "===================== updating at time t " << endl;

for(auto const &edge: e_out)
{
if (verbose)
cout << "edge out: " << edge.first << " " << edge.second << endl;
double ti = edge_activation_time[edge];

if (verbose)
cout << "edge was initialized at t = " << ti << endl;
double duration = t - ti;

aggregated_network[edge].value += duration;
if (verbose)
cout << "added to aggregated network" << endl;
cout << "aggregated_network[" << edge.first << " " << edge.second << "].value =" << aggregated_network[edge].value << endl;
}

for(auto const &edge: e_in)
{
if (verbose)
cout << "edge in: " << edge.first << " " << edge.second << endl;
edge_activation_time[edge] = t;
if (verbose)
cout << "added activation time" << endl;
}
}


void FlockworkPModel::get_rates_and_Lambda(
vector < double > &_rates,
Expand All @@ -109,13 +144,33 @@ void FlockworkPModel::make_event(
)
{

if (verbose)
cout << "attempting to make event " << event << endl;

e_in.clear();
e_out.clear();

if (event==0)
rewire(1.0,e_in,e_out);
else if (event == 1)
rewire(0.0,e_in,e_out);

if (verbose)
cout << "rewired." << endl;

if ( (e_in.size() > 0) or (e_out.size() > 0) )
{
if (save_aggregated_network)
update_aggregated_network(t, e_in, e_out);
if (save_temporal_network)
{
edg_chg.t.push_back(t);
edg_chg.edges_in.push_back(e_in);
edg_chg.edges_out.push_back(e_out);
}

}

}

void FlockworkPModel::simulate(
Expand All @@ -129,7 +184,8 @@ void FlockworkPModel::simulate(

double t = t0;

edg_chg.t0 = t0;
if (save_network)
edg_chg.t0 = t0;



Expand All @@ -148,11 +204,17 @@ void FlockworkPModel::simulate(

rewire(P,e_in,e_out);

if (save_network)
if ( (e_in.size() > 0) or (e_out.size() > 0) )
{
edg_chg.t.push_back(t);
edg_chg.edges_in.push_back(e_in);
edg_chg.edges_out.push_back(e_out);
if (save_network)
{
edg_chg.t.push_back(t);
edg_chg.edges_in.push_back(e_in);
edg_chg.edges_out.push_back(e_out);
}

if (save_aggregated_network)
update_aggregated_network(t, e_in, e_out);
}
}
}
Expand Down
95 changes: 70 additions & 25 deletions _tacoma/FlockworkPModel.h
View file
Expand Up @@ -58,10 +58,14 @@ class FlockworkPModel
size_t seed;
bool verbose;
bool save_temporal_network;
bool save_aggregated_network;
vector < pair < size_t, size_t > > E;
vector < double > rates;
double Lambda;

map < pair < size_t, size_t > , double > edge_activation_time;
map < pair < size_t, size_t > , edge_weight > aggregated_network;


vector < set < size_t > > G;

Expand All @@ -77,6 +81,7 @@ class FlockworkPModel
double _P,
double _t0 = 0.0,
bool _save_temporal_network = false,
bool _save_aggregated_network = false,
size_t _seed = 0,
bool _verbose = false
)
Expand All @@ -89,6 +94,7 @@ class FlockworkPModel
verbose = _verbose;
seed = _seed;
save_temporal_network = _save_temporal_network;
save_aggregated_network = _save_aggregated_network;

alpha = gamma*P;
beta = gamma*(1-P);
Expand All @@ -109,59 +115,73 @@ class FlockworkPModel
cout << "beta = " << beta << endl;
}

//reset();

}

void set_initial_configuration( double _t0, vector < set < size_t > > &_G)
{
vector < pair < size_t, size_t > > _E;
edgelist_from_graph(_E, _G);
reset(t0, E);
}

void reset(double _t0, vector < pair < size_t, size_t > > &_E)
{
t0 = _t0;

if (verbose)
cout << "resetting edge_changes..." << endl;

// reset observables
edg_chg.N = N;

edg_chg.t0 = t0;
edg_chg.t.clear();
edg_chg.edges_in.clear();
edg_chg.edges_out.clear();

if (verbose)
cout << "seeding engine and constructing graph..." << endl;

// seed engine
if (seed == 0)
randomly_seed_engine(*generator);
else
generator->seed(seed);

G = vector < set < size_t > >(N);
graph_from_edgelist(G, _E);

G.clear();
for(size_t node=0; node<N; node++)
{
set < size_t > neigh = _G[node];
G.push_back(neigh);
}
if (verbose)
cout << "initializing edge list..." << endl;

if (save_temporal_network)
edgelist_from_graph(edg_chg.edges_initial, G);

}

void reset()
{
// reset observables
edg_chg.N = N;
if (save_temporal_network)
edg_chg.edges_initial = _E;

edg_chg.t.clear();
edg_chg.edges_in.clear();
edg_chg.edges_out.clear();
if (verbose)
cout << "initializing aggrgated network..." << endl;

// seed engine
if (seed == 0)
randomly_seed_engine(*generator);
else
generator->seed(seed);
if (save_aggregated_network)
{
edge_activation_time.clear();
aggregated_network.clear();
for(auto &edge: _E)
{
edge_activation_time[edge] = t0;
}

G = vector < set < size_t > >(N);
graph_from_edgelist(G, E);
}
if (verbose)
cout << "done resetting..." << endl;

}

if (save_temporal_network)
edg_chg.edges_initial = E;
void reset()
{
reset(t0, E);
}

void get_rates_and_Lambda(vector < double > &rates,
Expand Down Expand Up @@ -204,6 +224,25 @@ class FlockworkPModel
delete generator;
}

map < pair <size_t, size_t>, double > finish_and_get_aggregated_network(double tmax)
{
// finish the contributions of edges at the end of the process
vector < pair < size_t, size_t > > _E;
vector < pair < size_t, size_t > > _dummy_e_in;
edgelist_from_graph(_E,G);
update_aggregated_network(tmax, _dummy_e_in, _E);

map < pair <size_t, size_t>, double > result;
for(auto const & entry: aggregated_network)
{
result[entry.first] = entry.second.value;
}

aggregated_network.clear();

return result;
}

private:

bool has_external_generator = false;
Expand All @@ -219,6 +258,12 @@ class FlockworkPModel
{
}

void update_aggregated_network(
double const &t,
vector < pair < size_t, size_t > > &e_in,
vector < pair < size_t, size_t > > &e_out
);



};
Expand Down
6 changes: 6 additions & 0 deletions _tacoma/SIR.cpp
View file
Expand Up @@ -180,6 +180,12 @@ void SIR::infection_event()
// change node status of this node
node_status[this_susceptible] = EPI::I;

// save the infection edge in observable
if(save_infection_events)
{
infection_events.push_back(SI_edges.at(this_susceptible_index));
}

// erase all edges in the SI set where this susceptible is part of
SI_edges.erase(
remove_if(
Expand Down
25 changes: 24 additions & 1 deletion _tacoma/SIR.h
View file
Expand Up @@ -54,7 +54,9 @@ class SIR
size_t number_of_initially_infected;
size_t number_of_initially_vaccinated;
size_t seed;
bool save_infection_events;
bool verbose;
bool stop_simulation_when_all_initially_infected_recovered;
double sampling_dt;

mt19937_64 generator;
Expand All @@ -65,6 +67,8 @@ class SIR
vector < size_t > SI;
vector < size_t > I;
vector < size_t > R;
vector < pair <size_t, size_t>> infection_events;
vector < size_t > initially_infected;

SIR(
size_t _N,
Expand All @@ -75,6 +79,8 @@ class SIR
size_t _number_of_initially_vaccinated = 0,
double _sampling_dt = 0.0,
size_t _seed = 0,
bool _save_infection_events = false,
bool _stop_simulation_when_all_initially_infected_recovered = false,
bool _verbose = false
)
{
Expand All @@ -84,6 +90,8 @@ class SIR
recovery_rate = _recovery_rate;
number_of_initially_vaccinated = _number_of_initially_vaccinated;
number_of_initially_infected = _number_of_initially_infected;
save_infection_events = _save_infection_events || _stop_simulation_when_all_initially_infected_recovered;
stop_simulation_when_all_initially_infected_recovered = _stop_simulation_when_all_initially_infected_recovered;
verbose = _verbose;
seed = _seed;
sampling_dt = _sampling_dt;
Expand All @@ -104,6 +112,8 @@ class SIR
SI.clear();
I.clear();
R.clear();
infection_events.clear();
initially_infected.clear();

// seed engine
if (seed == 0)
Expand Down Expand Up @@ -152,6 +162,11 @@ class SIR
infected.push_back( node_ints[node] );
}

for(auto const &inf: infected)
{
initially_infected.push_back(inf);
}

if (verbose)
{
cout << "infected set has size = " << infected.size() << endl;
Expand All @@ -166,7 +181,15 @@ class SIR

bool simulation_ended()
{
return (infected.size() == 0);
if (stop_simulation_when_all_initially_infected_recovered)
{
bool stop = true;
for(auto const &i: initially_infected)
stop = (stop) and (node_status[i] != EPI::I);
return stop;
}
else
return (infected.size() == 0);
};

void get_rates_and_Lambda(vector < double > &rates,
Expand Down

0 comments on commit e51b358

Please sign in to comment.