mixing_graphs.jl 8.4 KB
 1 ``````using Dictionaries `````` Peter Jentsch committed Apr 21, 2021 2 3 4 5 6 ``````""" A type that stores a pair of nodes representing an edge in the graph. We need to define a custom type for these so we can define a hash function on graph edges, in order to more efficiently use them in hashmaps (dictionaries) """ `````` 7 8 9 10 ``````struct GraphEdge a::Int b::Int end `````` Peter Jentsch committed Apr 21, 2021 11 12 13 14 15 16 `````` """ Define a hash on GraphEdge such that ``hash(a,b) = hash(b,a)`` (hash is commutative). This is helpful because then we only need to store (a,b) in the graph edges weights dictionary, rather than both (a,b) and (b,a). """ `````` 17 18 19 ``````function Base.hash(e::GraphEdge) return hash(minmax(e.a,e.b)) end `````` Peter Jentsch committed Apr 21, 2021 20 21 `````` """ `````` Peter Jentsch committed Apr 23, 2021 22 ``````Define symmetric edge equality, matches the hash function. `````` Peter Jentsch committed Apr 21, 2021 23 ``````""" `````` 24 25 26 ``````function Base.isequal(e1::GraphEdge,e2::GraphEdge) return isequal(minmax(e1.a,e1.b),minmax(e2.a,e2.b)) end `````` Peter Jentsch committed Apr 21, 2021 27 28 29 30 31 32 33 `````` """ sample_mixing_graph!(mixing_graph::Graph) Resample all the weights in `mixing_graph` """ `````` Peter Jentsch committed Jun 14, 2021 34 35 36 37 38 39 40 41 ``````function sample_mixing_edges!(weights_dict,sampler_matrix,demographics) indices = keys(weights_dict) for ind in indices e = ind i = e.a j = e.b weight = rand(Random.default_rng(Threads.threadid()),sampler_matrix[Int(demographics[j]),Int(demographics[i])]) setindex!(weights_dict,weight,ind) `````` 42 43 44 `````` end end `````` Peter Jentsch committed Apr 21, 2021 45 46 47 48 49 ``````""" Stores the full time dependent mixing graph for the model. I think this might be a weird abstraction for this idea but it works fine. #Fields `````` Peter Jentsch committed Jun 13, 2021 50 51 `````` remade_graphs::NTuple{N,G} `````` Peter Jentsch committed Apr 21, 2021 52 53 ``````These are references to the graphs that get resampled everyday. `````` Peter Jentsch committed Jun 13, 2021 54 55 56 `````` resampled_graphs::NTuple{N,G} These are references to the graphs that get resampled everyday. `````` Peter Jentsch committed Apr 21, 2021 57 58 59 60 `````` graph_list::Vector{Vector{G}} `````` Peter Jentsch committed Jun 13, 2021 61 62 ``````List of lists of graphs, one list for each day. `````` Peter Jentsch committed Apr 21, 2021 63 ``````""" `````` Peter Jentsch committed Jun 13, 2021 64 65 66 ``````struct TimeDepMixingGraph{G,T1,T2} remade_graphs::T1 resampled_graphs::T2 `````` 67 `````` graph_list::Vector{Vector{G}} `````` Peter Jentsch committed Jun 13, 2021 68 69 70 `````` function TimeDepMixingGraph(len,remade_graphs::T1,resampled_graphs::T2,base_graph_list::Vector{G}) where {G,T1,T2} return new{G,T1,T2}( remade_graphs, `````` 71 72 73 74 75 `````` resampled_graphs, [copy(base_graph_list) for i in 1:len] ) end end `````` Peter Jentsch committed Apr 21, 2021 76 77 78 79 80 81 `````` """ Creates the `TimeDepMixingGraph` for our specific model. Assumes the simulation begins on Thursday arbitrarily. """ `````` 82 ``````function time_dep_mixing_graphs(len,base_network,demographics,index_vectors,ws_matrix_tuple,rest_matrix_tuple) `````` Peter Jentsch committed Jun 13, 2021 83 `````` `````` Peter Jentsch committed Jun 14, 2021 84 `````` home_static_edges = WeightedGraph(base_network,demographics,contact_time_distributions.hh) #network with households and LTC homes `````` Peter Jentsch committed Jun 13, 2021 85 `````` `````` 86 87 `````` ws_static_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.daily,contact_time_distributions.ws) ws_weekly_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.twice_a_week,contact_time_distributions.ws) `````` Peter Jentsch committed Jun 14, 2021 88 `````` ws_justonce_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.otherwise,contact_time_distributions.ws) `````` Peter Jentsch committed Jun 13, 2021 89 `````` `````` 90 91 `````` rest_static_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.daily,contact_time_distributions.rest) rest_weekly_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.twice_a_week,contact_time_distributions.rest) `````` Peter Jentsch committed Jun 14, 2021 92 `````` rest_justonce_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.otherwise,contact_time_distributions.rest) `````` Peter Jentsch committed May 10, 2021 93 `````` `````` Peter Jentsch committed Jun 14, 2021 94 `````` inf_network_list = [home_static_edges,rest_static_edges,ws_justonce_edges,rest_justonce_edges] `````` 95 `````` soc_network_list = [home_static_edges,rest_static_edges,ws_static_edges] `````` Peter Jentsch committed Jun 13, 2021 96 `````` `````` Peter Jentsch committed Jun 14, 2021 97 `````` remade_graphs = (ws_justonce_edges,rest_justonce_edges) `````` Peter Jentsch committed Jun 13, 2021 98 99 100 101 `````` resampled_graphs = (home_static_edges,rest_static_edges,ws_static_edges,rest_weekly_edges,ws_weekly_edges) infected_mixing_graph = TimeDepMixingGraph(len,remade_graphs,resampled_graphs,inf_network_list) `````` 102 103 104 105 106 `````` for (t,l) in enumerate(infected_mixing_graph.graph_list) day_of_week = mod(t,7) if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess push!(l, ws_static_edges) end `````` Peter Jentsch committed May 18, 2021 107 `````` if rand(Random.default_rng(Threads.threadid()))<5/7 `````` 108 109 110 111 112 `````` push!(l, ws_weekly_edges) push!(l, rest_weekly_edges) end end `````` Peter Jentsch committed Jun 13, 2021 113 `````` soc_mixing_graph = TimeDepMixingGraph(len,remade_graphs,resampled_graphs,soc_network_list) `````` 114 115 116 `````` return infected_mixing_graph,soc_mixing_graph end `````` Peter Jentsch committed Apr 21, 2021 117 118 119 ``````""" Completely remake all the graphs in `time_dep_mixing_graph.resampled_graphs`. """ `````` Peter Jentsch committed Jun 14, 2021 120 ``````function remake!(t,time_dep_mixing_graph,index_vectors,demographics) `````` Peter Jentsch committed Jun 13, 2021 121 122 123 `````` for wg in time_dep_mixing_graph.remade_graphs remake!(wg,index_vectors) end `````` Peter Jentsch committed Jun 14, 2021 124 `````` # display_degree(time_dep_mixing_graph.resampled_graphs[1]) `````` Peter Jentsch committed Jun 13, 2021 125 126 `````` for wg in time_dep_mixing_graph.resampled_graphs if wg in time_dep_mixing_graph.graph_list[t] `````` Peter Jentsch committed Jun 14, 2021 127 `````` sample_mixing_edges!(wg.weights_dict,wg.sampler_matrix,demographics) `````` Peter Jentsch committed Jun 13, 2021 128 129 `````` end end `````` Peter Jentsch committed Jun 14, 2021 130 `````` # display_degree(time_dep_mixing_graph.resampled_graphs[1]) `````` Peter Jentsch committed Apr 20, 2021 131 ``````end `````` Peter Jentsch committed Apr 21, 2021 132 `````` `````` Peter Jentsch committed Jun 14, 2021 133 `````` `````` Peter Jentsch committed Apr 21, 2021 134 ``````""" `````` Peter Jentsch committed Jun 13, 2021 135 ``````Weighted graph type. Stores the graph in `g`, and the weights and edges in `mixing_edges`. `````` Peter Jentsch committed Jun 13, 2021 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 ``````Fields g::SimpleGraph Stores the actual graph structure weights_dict::Dictionary{GraphEdge,UInt8} Stores the weights used in the graph, so they can be easily resampled. mixing_matrix::M1 Matrix of distributions determining node degrees sampler_matrix::M Matrix of distributions determining the edge weights `````` Peter Jentsch committed Apr 21, 2021 155 ``````""" `````` Peter Jentsch committed Jun 14, 2021 156 ``````struct WeightedGraph{G,M1,M2} `````` Peter Jentsch committed Apr 08, 2021 157 `````` g::G `````` Peter Jentsch committed Jun 13, 2021 158 159 160 `````` weights_dict::Dictionary{GraphEdge,UInt8} mixing_matrix::M1 sampler_matrix::M2 `````` Peter Jentsch committed Jun 14, 2021 161 `````` function WeightedGraph(demographics::AbstractVector,index_vectors,mixing_matrix::M1, sampler_matrix::M2) where {M1,M2} `````` 162 `````` g = Graph(length(demographics)) `````` Peter Jentsch committed Jun 14, 2021 163 `````` weights_dict = Dictionary{GraphEdge,UInt8}() `````` Peter Jentsch committed Jun 14, 2021 164 `````` assemble_graph!(g,weights_dict,index_vectors,mixing_matrix,sampler_matrix) `````` Peter Jentsch committed Jun 13, 2021 165 `````` return new{typeof(g),M1,M2}( `````` Peter Jentsch committed Apr 08, 2021 166 `````` g, `````` Peter Jentsch committed Jun 13, 2021 167 168 `````` weights_dict, mixing_matrix, `````` Peter Jentsch committed Jun 14, 2021 169 `````` sampler_matrix, `````` Peter Jentsch committed Apr 08, 2021 170 171 `````` ) end `````` Peter Jentsch committed Jun 14, 2021 172 `````` function WeightedGraph(g::G,demographics,sampler_matrix::M2) where {G<:SimpleGraph,M2} `````` Peter Jentsch committed Jun 14, 2021 173 `````` weights_dict = Dictionary{GraphEdge,UInt8}(;sizehint = ne(g)) `````` Peter Jentsch committed Jun 13, 2021 174 `````` for e in edges(g) `````` Peter Jentsch committed Jun 14, 2021 175 176 177 178 `````` j = src(e) i = dst(e) weight = rand(Random.default_rng(Threads.threadid()),sampler_matrix[Int(demographics[j]),Int(demographics[i])]) set!(weights_dict,GraphEdge(j,i),weight) `````` Peter Jentsch committed Jun 13, 2021 179 180 `````` end return new{typeof(g),Nothing,M2}( `````` Peter Jentsch committed Apr 08, 2021 181 `````` g, `````` Peter Jentsch committed Jun 13, 2021 182 183 `````` weights_dict, nothing, `````` Peter Jentsch committed Jun 14, 2021 184 `````` sampler_matrix, `````` Peter Jentsch committed Apr 08, 2021 185 `````` ) `````` Peter Jentsch committed Jan 30, 2021 186 187 `````` end end `````` Peter Jentsch committed Jun 13, 2021 188 189 190 191 192 `````` function remake!(wg::WeightedGraph,index_vectors) empty!.(wg.g.fadjlist) #empty all the vector edgelists wg.g.ne = 0 empty!(wg.weights_dict) `````` Peter Jentsch committed Jun 14, 2021 193 194 195 196 197 `````` assemble_graph!(wg.g,wg.weights_dict,index_vectors,wg.mixing_matrix,wg.sampler_matrix) end function assemble_graph!(g,weights_dict,index_vectors,mixing_matrix,sampler_matrix) `````` Peter Jentsch committed Jun 14, 2021 198 199 `````` for i in 1:3, j in 1:i #diagonal if i != j `````` Peter Jentsch committed Jun 14, 2021 200 `````` edges = fast_chung_lu_bipartite(index_vectors[i],index_vectors[j],mixing_matrix[i,j],mixing_matrix[j,i]) `````` Peter Jentsch committed Jun 14, 2021 201 `````` else #from one group to itself we need another algorithm `````` Peter Jentsch committed Jun 14, 2021 202 `````` edges = fast_chung_lu(index_vectors[i],mixing_matrix[i,i]) `````` Peter Jentsch committed Jun 14, 2021 203 204 `````` end for e in edges `````` Peter Jentsch committed Jun 14, 2021 205 206 `````` edge_weight_k = rand(Random.default_rng(Threads.threadid()),sampler_matrix[j,i]) set!(weights_dict, e, edge_weight_k) `````` Peter Jentsch committed Jun 14, 2021 207 208 `````` end end `````` Peter Jentsch committed Jun 14, 2021 209 210 211 `````` for e in keys(weights_dict) add_edge!(g,e.a,e.b) end `````` Peter Jentsch committed Jun 13, 2021 212 ``````end `````` Peter Jentsch committed Jun 14, 2021 213 214 215 216 217 218 219 220 221 222 223 224 225 ``````function fast_chung_lu_bipartite(pop_i,pop_j,mixing_dist_ij,mixing_dist_ji) num_degrees_ij = similar(pop_i) num_degrees_ji = similar(pop_j) generate_contact_vectors!(mixing_dist_ij,mixing_dist_ji,num_degrees_ij,num_degrees_ji) num_edges = sum(num_degrees_ij) stubs_i = Vector{Int}(undef,num_edges) stubs_j = similar(stubs_i) if num_edges>0 sample!(Random.default_rng(Threads.threadid()),pop_i,Weights(num_degrees_ij./num_edges),stubs_i) sample!(Random.default_rng(Threads.threadid()),pop_j,Weights(num_degrees_ji./num_edges),stubs_j) end return GraphEdge.(stubs_i,stubs_j) end `````` Peter Jentsch committed Jun 14, 2021 226 `````` `````` Peter Jentsch committed Jun 14, 2021 227 228 229 230 231 232 233 234 235 236 ``````function fast_chung_lu(pop_i,mixing_dist) num_degrees_ii = similar(pop_i) generate_contact_vectors!(mixing_dist,num_degrees_ii) m = sum(num_degrees_ii) num_edges= div(m,2) stubs_i = Vector{Int}(undef,num_edges) stubs_j = similar(stubs_i) if m>0 sample!(Random.default_rng(Threads.threadid()),pop_i,Weights(num_degrees_ii./m),stubs_i) sample!(Random.default_rng(Threads.threadid()),pop_i,Weights(num_degrees_ii./m),stubs_j) `````` Peter Jentsch committed Jun 13, 2021 237 `````` end `````` Peter Jentsch committed Jun 14, 2021 238 `````` return GraphEdge.(stubs_i,stubs_j) `````` Peter Jentsch committed Jun 14, 2021 239 240 241 242 243 244 245 246 247 ``````end neighbors(g::WeightedGraph,i) = neighbors(g.g,i) get_weight(g::WeightedGraph,e) = g.weights_dict[e] function Base.show(io::IO, g::WeightedGraph) print(io, "WG \$(ne(g.g))") end ``````