diff --git a/_tacoma/Events.cpp b/_tacoma/Events.cpp index cd61b0b..a8a8470 100644 --- a/_tacoma/Events.cpp +++ b/_tacoma/Events.cpp @@ -174,6 +174,70 @@ pair < vector < pair < size_t, size_t > >, vector < pair < size_t, size_t > > > return make_pair(edges_out,edges_in); } +pair < vector < pair < size_t, size_t > >, vector < pair < size_t, size_t > > > + rewire_P_without_SI_checking( + vector < set < size_t > > & G, //Adjacency matrix + double P, //probability to connect with neighbors of neighbor + mt19937_64 & generator, + uniform_real_distribution & distribution + ) +{ + //choose two nodes + size_t N = G.size(); + double r1 = distribution(generator); + double r2 = distribution(generator); + size_t i,j; + choose(N,i,j,r1,r2); + + bool do_rewiring = distribution(generator) < P; + + vector < pair < size_t, size_t > > edges_out; + vector < pair < size_t, size_t > > edges_in; + + //check if new neighbor is actually an old neighbor + //and if this is the case return an empty event + if ( do_rewiring and (G[i].find(j) != G[i].end()) ) + { + return make_pair(edges_out,edges_in); + } + + //loop through the neighbors of i + for(auto neigh_i : G[i] ) + { + //and erase the link to i + G[neigh_i].erase(i); + + pair < size_t, size_t > current_edge = get_sorted_pair(i,neigh_i); + edges_out.push_back( current_edge ); + } + + //erase the links from the perspective of i + G[i].clear(); + + if ( do_rewiring ) + { + //loop through the neighbors of j + for(auto neigh_j : G[j] ) + { + G[ neigh_j ].insert( i ); + G[ i ].insert( neigh_j ); + + pair < size_t, size_t > current_edge = get_sorted_pair(i,neigh_j); + edges_in.push_back( current_edge ); + } + + //add j as neighbor and count additional edge + G[ j ].insert( i ); + G[ i ].insert( j ); + + pair < size_t, size_t > current_edge = get_sorted_pair(i,j); + edges_in.push_back( current_edge ); + } + + return make_pair(edges_out,edges_in); +} + + pair < vector < pair < size_t, size_t > >, vector < pair < size_t, size_t > > > rewire_P_without_SI_checking_single_node( diff --git a/_tacoma/Events.h b/_tacoma/Events.h index 1c5783b..cf31964 100644 --- a/_tacoma/Events.h +++ b/_tacoma/Events.h @@ -72,6 +72,14 @@ pair < vector < pair < size_t, size_t > >, vector < pair < size_t, size_t > > > uniform_real_distribution & distribution ); +pair < vector < pair < size_t, size_t > >, vector < pair < size_t, size_t > > > + rewire_P_without_SI_checking( + vector < set < size_t > > & G, //Adjacency matrix + double P, //probability to connect with neighbors of neighbor + mt19937_64 & generator, + uniform_real_distribution & distribution + ); + pair < vector < pair < size_t, size_t > >, vector < pair < size_t, size_t > > > rewire_P_without_SI_checking_single_node( size_t i, diff --git a/_tacoma/FlockworkPModel.cpp b/_tacoma/FlockworkPModel.cpp new file mode 100644 index 0000000..3d24de9 --- /dev/null +++ b/_tacoma/FlockworkPModel.cpp @@ -0,0 +1,110 @@ +/* + * The MIT License (MIT) + * Copyright (c) 2018, Benjamin Maier + * + * Permission is hereby granted, free of charge, to any person + * obtaining a copy of this software and associated documentation + * files (the "Software"), to deal in the Software without + * restriction, including without limitation the rights to use, + * copy, modify, merge, publish, distribute, sublicense, and/or + * sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall + * be included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NON- + * INFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS + * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN + * AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF + * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS + * IN THE SOFTWARE. + */ + +#include "FlockworkPModel.h" + +using namespace std; +//namespace osis = SIS; +void + FlockworkPModel::rewire( + double const &_P, //probability to connect with neighbors of neighbor + vector < pair < size_t, size_t > > &edges_in, + vector < pair < size_t, size_t > > &edges_out + ) +{ + //choose two nodes + double r1 = randuni(*generator); + double r2 = randuni(*generator); + size_t i,j; + choose(N,i,j,r1,r2); + + bool do_rewiring = randuni(*generator) < _P; + + //check if new neighbor is actually an old neighbor + //and if this is the case return an empty event + if (not ( do_rewiring and (G[i].find(j) != G[i].end()) )) + { + //loop through the neighbors of i + for(auto neigh_i : G[i] ) + { + //and erase the link to i + G[neigh_i].erase(i); + + pair < size_t, size_t > current_edge = get_sorted_pair(i,neigh_i); + edges_out.push_back( current_edge ); + } + + //erase the links from the perspective of i + G[i].clear(); + + if ( do_rewiring ) + { + //loop through the neighbors of j + for(auto neigh_j : G[j] ) + { + G[ neigh_j ].insert( i ); + G[ i ].insert( neigh_j ); + + pair < size_t, size_t > current_edge = get_sorted_pair(i,neigh_j); + edges_in.push_back( current_edge ); + } + + //add j as neighbor and count additional edge + G[ j ].insert( i ); + G[ i ].insert( j ); + + pair < size_t, size_t > current_edge = get_sorted_pair(i,j); + edges_in.push_back( current_edge ); + } + } +} + + +void FlockworkPModel::get_rates_and_Lambda( + vector < double > &_rates, + double &_Lambda + ) +{ + // delete current rates + _rates = rates; + _Lambda = Lambda; +} + +void FlockworkPModel::make_event( + size_t const &event, + double t, + vector < pair < size_t, size_t > > &e_in, + vector < pair < size_t, size_t > > &e_out + ) +{ + + 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); +} diff --git a/_tacoma/FlockworkPModel.h b/_tacoma/FlockworkPModel.h new file mode 100644 index 0000000..2abd497 --- /dev/null +++ b/_tacoma/FlockworkPModel.h @@ -0,0 +1,213 @@ +/* + * The MIT License (MIT) + * Copyright (c) 2018, Benjamin Maier + * + * Permission is hereby granted, free of charge, to any person + * obtaining a copy of this software and associated documentation + * files (the "Software"), to deal in the Software without + * restriction, including without limitation the rights to use, + * copy, modify, merge, publish, distribute, sublicense, and/or + * sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall + * be included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NON- + * INFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS + * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN + * AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF + * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS + * IN THE SOFTWARE. + */ + +#ifndef __FLOCKWORKP_MODEL_CLASS_H__ +#define __FLOCKWORKP_MODEL_CLASS_H__ + +#include "Utilities.h" +#include "ResultClasses.h" +#include "activity_model.h" + +//#include +//#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +class FlockworkPModel +{ + public: + size_t N; + double gamma; + double P; + double alpha; + double beta; + double t0; + size_t seed; + bool verbose; + bool save_temporal_network; + vector < pair < size_t, size_t > > E; + vector < double > rates; + double Lambda; + + + vector < set < size_t > > G; + + mt19937_64 * generator; + uniform_real_distribution randuni; + + edge_changes edg_chg; + + FlockworkPModel( + vector < pair < size_t, size_t > > _E, + size_t _N, + double _gamma, + double _P, + double _t0 = 0.0, + bool _save_temporal_network = false, + size_t _seed = 0, + bool _verbose = false + ) + { + E = _E; + N = _N; + P = _P; + gamma = _gamma; + t0 = _t0; + verbose = _verbose; + seed = _seed; + save_temporal_network = _save_temporal_network; + + alpha = gamma*P; + beta = gamma*(1-P); + + rates.push_back(N*alpha); + rates.push_back(N*beta); + Lambda = N*alpha + N*beta; + + + generator = new mt19937_64; + randuni = uniform_real_distribution < double > (0.0, 1.0); + + if (verbose) + { + cout << "gamma = " << gamma << endl; + cout << "P = " << P << endl; + cout << "alpha = " << alpha << endl; + cout << "beta = " << beta << endl; + } + + } + + void set_initial_configuration( double _t0, vector < set < size_t > > &_G) + { + t0 = _t0; + + // reset observables + edg_chg.N = N; + + edg_chg.t.clear(); + edg_chg.edges_in.clear(); + edg_chg.edges_out.clear(); + + // seed engine + if (seed == 0) + randomly_seed_engine(*generator); + else + generator->seed(seed); + + + G.clear(); + for(size_t node=0; node neigh = _G[node]; + G.push_back(neigh); + } + + if (save_temporal_network) + edgelist_from_graph(edg_chg.edges_initial, G); + + } + + void reset() + { + // reset observables + edg_chg.N = N; + + edg_chg.t.clear(); + edg_chg.edges_in.clear(); + edg_chg.edges_out.clear(); + + // seed engine + if (seed == 0) + randomly_seed_engine(*generator); + else + generator->seed(seed); + + G = vector < set < size_t > >(N); + graph_from_edgelist(G, E); + + + if (save_temporal_network) + edg_chg.edges_initial = E; + } + + void get_rates_and_Lambda(vector < double > &rates, + double &Lambda + ); + + void make_event(size_t const &event, + double t, + vector < pair < size_t, size_t > > &e_in, + vector < pair < size_t, size_t > > &e_out + ); + + void print() + { + } + + void set_generator(mt19937_64 &_generator) + { + delete generator; + generator = &_generator; + has_external_generator = true; + } + + ~FlockworkPModel() + { + if (not has_external_generator) + delete generator; + } + + private: + + bool has_external_generator = false; + + void + rewire( + double const &_P, //probability to connect with neighbors of neighbor + vector < pair < size_t, size_t > > &edges_in, + vector < pair < size_t, size_t > > &edges_out + ); + + void print_internal() + { + } + + + +}; + +#endif diff --git a/_tacoma/_tacoma.cpp b/_tacoma/_tacoma.cpp index 3057146..c3b4b79 100644 --- a/_tacoma/_tacoma.cpp +++ b/_tacoma/_tacoma.cpp @@ -65,6 +65,7 @@ #include "flockwork_parameter_estimation.h" #include "activity_model.h" #include "EdgeActivityModel.h" +#include "FlockworkPModel.h" #include "slice.h" using namespace std; @@ -817,6 +818,63 @@ PYBIND11_MODULE(_tacoma, m) py::arg("reset_simulation_objects") = true, py::arg("verbose") = false); + m.def("gillespie_SIS_on_FlockworkPModel", &gillespie_on_model, + "Perform a Gillespie SIS simulation on the Flockwork P-model.", + py::arg("model"), + py::arg("SIS"), + py::arg("reset_simulation_objects") = true, + py::arg("verbose") = false); + + m.def("gillespie_QS_SIS_on_FlockworkPModel", &gillespie_on_model, + "Perform a quasi-stationary Gillespie SIS simulation on the Flockwork P-model.", + py::arg("model"), + py::arg("QS_SIS"), + py::arg("reset_simulation_objects") = false, + py::arg("verbose") = false); + + m.def("gillespie_eSIS_on_FlockworkPModel", &gillespie_on_model, + R"pbdoc(Perform a Gillespie :math:`\varepsilon`-SIS simulation on the Flockwork P-model.)pbdoc", + py::arg("model"), + py::arg("eSIS"), + py::arg("reset_simulation_objects") = true, + py::arg("verbose") = false); + + m.def("gillespie_SIR_on_FlockworkPModel", &gillespie_on_model, + "Perform a Gillespie SIR simulation on the Flockwork P-model.", + py::arg("model"), + py::arg("SIR"), + py::arg("reset_simulation_objects") = true, + py::arg("verbose") = false); + + m.def("gillespie_SIRS_on_FlockworkPModel", &gillespie_on_model, + "Perform a Gillespie SIRS simulation on the Flockwork P-model.", + py::arg("model"), + py::arg("SIRS"), + py::arg("reset_simulation_objects") = true, + py::arg("verbose") = false); + + m.def("gillespie_SI_on_FlockworkPModel", &gillespie_on_model, + "Perform a Gillespie SI simulation on the Flockwork P-model.", + py::arg("model"), + py::arg("SI"), + py::arg("reset_simulation_objects") = true, + py::arg("verbose") = false); + + m.def("gillespie_node_based_SIS_on_FlockworkPModel", &gillespie_on_model, + "Perform a node-based Gillespie SIS simulation on the Flockwork P-model.", + py::arg("model"), + py::arg("SI"), + py::arg("reset_simulation_objects") = true, + py::arg("verbose") = false); + + m.def("markov_SIS_on_FlockworkPModel", &markov_on_model, + "Perform a mixed Markov-Gillespie SIS simulation on the Flockwork P-model.", + py::arg("model"), + py::arg("MARKOV_SIS"), + py::arg("max_dt"), + py::arg("reset_simulation_objects") = true, + py::arg("verbose") = false); + py::class_(m, "edge_changes", R"pbdoc(Description of a temporal network by listing the changes of edges at a certain time.)pbdoc") .def(py::init<>()) @@ -1534,4 +1592,45 @@ PYBIND11_MODULE(_tacoma, m) .def_readwrite("N", &EdgeActivityModel::N, R"pbdoc(Number of nodes.)pbdoc"); + py::class_(m, "FlockworkPModel", "Base class for the simulation of a simple Flockwork-P model. Pass this to :func:`tacoma.api.gillespie_epidemics`") + .def(py::init< vector < pair < size_t, size_t > > , size_t, double, double, double, bool, size_t, bool>(), + py::arg("E"), + py::arg("N"), + py::arg("gamma"), + py::arg("P"), + py::arg("t0") = 0.0, + py::arg("save_temporal_network") = false, + py::arg("seed") = 0, + py::arg("verbose") = false, + R"pbdoc( + Parameters + ---------- + E : list of pair of int + Initial edge list. + N : int + Number of nodes in the temporal network. + gamma : float + The probability per unit time per node that any event happens. + P : float + The probability to reconnect when an event happened. + t0 : float, default = 0.0 + initial time + save_temporal_network : bool, default: False + If this is `True`, the changes are saved in an instance of + :func:`_tacoma.edge_changes` (in the attribute `edge_changes`. + seed : int, default = 0 + Seed for RNG initialization. If this is 0, the seed will be initialized randomly. + However, the generator will be rewritten + in :func:`tacoma.api.gillespie_SIS_EdgeActivityModel` anyway. + verbose : bool, default = False + Be talkative. + )pbdoc") + .def("set_initial_configuration",&FlockworkPModel::set_initial_configuration, + R"pbdoc(Reset the state of the network to a certain graph (:obj:`list` of :obj:`set` of :obj:`int`))pbdoc") + .def_readwrite("edge_changes", &FlockworkPModel::edg_chg, + R"pbdoc(An instance of :class:`_tacoma.edge_changes` with the saved temporal network (only if + `save_temporal_network` is `True`).)pbdoc") + .def_readwrite("N", &FlockworkPModel::N, + R"pbdoc(Number of nodes.)pbdoc"); + } diff --git a/sandbox/test_markov_SIS.py b/sandbox/test_markov_SIS.py index 9aa2d63..2ebd54a 100644 --- a/sandbox/test_markov_SIS.py +++ b/sandbox/test_markov_SIS.py @@ -12,7 +12,7 @@ t_run_total = 500.0 tmax = 1000.0 seed = 7925 -infection_rate = 0.5 +infection_rate = 0.1 recovery_rate = 0.1 E = flockwork_P_equilibrium_configuration(N,P[0]) @@ -55,6 +55,17 @@ end = time.time() print("Markov on edge_lists took", end-start,"seconds") +pl.plot(mv_SIS.time, mv_SIS.I) + +FW = tc.FlockworkPModel(E, N, rewiring_rate[0][1], P[0]) +mv_SIS = tc.MARKOV_SIS(N,t_run_total*2,infection_rate,recovery_rate,0.01,N//2,seed = seed) + +start = time.time() +_tc.markov_SIS_on_FlockworkPModel(FW, mv_SIS,max_dt = 0.001) +end = time.time() +print("Markov on Model took", end-start,"seconds") + + pl.plot(mv_SIS.time, mv_SIS.I) pl.ylim([0,N]) diff --git a/setup.py b/setup.py index b1f61a6..e57865a 100644 --- a/setup.py +++ b/setup.py @@ -27,6 +27,7 @@ def __str__(self): [ '_tacoma/Utilities.cpp', '_tacoma/Events.cpp', + '_tacoma/FlockworkPModel.cpp', '_tacoma/slice.cpp', '_tacoma/eSIS.cpp', '_tacoma/QS_SIS.cpp', diff --git a/tacoma/api.py b/tacoma/api.py index 0365ec0..17dc8dd 100644 --- a/tacoma/api.py +++ b/tacoma/api.py @@ -31,6 +31,8 @@ eSIS = _tc.eSIS EdgeActivityModel = _tc.EdgeActivityModel QS_SIS = _tc.QS_SIS +FlockworkPModel = _tc.FlockworkPModel +MARKOV_SIS = _tc.MARKOV_SIS def _get_raw_temporal_network(temporal_network): @@ -334,7 +336,7 @@ def gillespie_epidemics(temporal_network_or_model, epidemic_object, is_static=Fa Parameters ---------- - temporal_network_or_model : :class:`_tacoma.edge_changes`, :class:`_tacoma.edge_lists`, :class:`_tacoma.edge_changes_with_histograms`, :class:`_tacoma.edge_lists_with_histograms`, or :class:`_tacoma.EdgeActivityModel`. + temporal_network_or_model : :class:`_tacoma.edge_changes`, :class:`_tacoma.edge_lists`, :class:`_tacoma.edge_changes_with_histograms`, :class:`_tacoma.edge_lists_with_histograms`, :class:`_tacoma.EdgeActivityModel`, or :class:`_tacoma.FlockworkPModel`. An instance of a temporal network or network model. epidemic_object : :class:`_tacoma.SI`, :class:`_tacoma.SIS`, :class:`_tacoma.SIR`, :class:`_tacoma.SIRS`, :class:`_tacoma.node_based_SIS`, :class:`_tacoma.eSIS` An initialized epidemic object. @@ -400,6 +402,22 @@ def gillespie_epidemics(temporal_network_or_model, epidemic_object, is_static=Fa _tc.gillespie_eSIS_on_EdgeActivityModel(temporal_network, epidemic_object, verbose) else: raise ValueError('Invalid epidemic object type: ' + str(type(tn_or_mdl))) + elif type(tn_or_mdl) in [FlockworkPModel]: + temporal_network = tn_or_mdl + if type(epidemic_object) == SI: + _tc.gillespie_SI_on_FlockworkPModel(temporal_network, epidemic_object, verbose) + elif type(epidemic_object) == SIS: + _tc.gillespie_SIS_on_FlockworkPModel(temporal_network, epidemic_object, verbose) + elif type(epidemic_object) == SIR: + _tc.gillespie_SIR_on_FlockworkPModel(temporal_network, epidemic_object, verbose) + elif type(epidemic_object) == SIRS: + _tc.gillespie_SIRS_on_FlockworkPModel(temporal_network, epidemic_object, verbose) + elif type(epidemic_object) == node_based_SIS: + _tc.gillespie_node_based_SIS_on_FlockworkPModel(temporal_network, epidemic_object, verbose) + elif type(epidemic_object) == eSIS: + _tc.gillespie_eSIS_on_FlockworkPModel(temporal_network, epidemic_object, verbose) + else: + raise ValueError('Invalid epidemic object type: ' + str(type(tn_or_mdl))) else: raise ValueError('Invalid temporal network/model format: ' + str(type(tn_or_mdl))) diff --git a/tacoma/metadata.py b/tacoma/metadata.py index cf6686a..f4c087d 100644 --- a/tacoma/metadata.py +++ b/tacoma/metadata.py @@ -3,7 +3,7 @@ Contains a bunch of information about this package """ -__version__ = "0.1.11" +__version__ = "0.1.12" __author__ = "Benjamin F. Maier" __copyright__ = "Copyright 2018, Benjamin F. Maier"