solve.jl 11.7 KB
Newer Older
 Peter Jentsch committed May 11, 2021 1 2 ``````function contact_weight(β, contact_time) return 1 - (1-β)^contact_time `````` 3 ``````end `````` Peter Jentsch committed Apr 21, 2021 4 `````` `````` Peter Jentsch committed May 11, 2021 5 ``````function Φ(payoff,ξ) `````` Peter Jentsch committed May 30, 2021 6 `````` return 1 / (1 + exp(-1*ξ*payoff)) `````` 7 8 ``````end `````` Peter Jentsch committed Jun 30, 2021 9 10 11 `````` @views function update_alert_durations!(t,agent,modelsol) `````` 12 `````` @unpack notification_parameter,notification_threshold = modelsol.params `````` Peter Jentsch committed Jun 30, 2021 13 14 15 16 17 18 19 20 21 `````` @unpack time_of_last_alert,inf_network,covid_alert_notifications,is_app_user, output_data = modelsol for j in 2:14 covid_alert_notifications[j-1,agent] = covid_alert_notifications[j,agent] #shift them all back end total_weight_i = 0 for mixing_graph in inf_network.graph_list[t], neighbor in neighbors(mixing_graph,agent) if is_app_user[neighbor] total_weight_i+= get_weight(mixing_graph,GraphEdge(agent,neighbor)) `````` 22 23 `````` end end `````` Peter Jentsch committed Jun 30, 2021 24 25 26 27 28 29 30 31 32 33 34 `````` coin_flip = 1 - (1 - notification_parameter)^total_weight_i r = rand(Random.default_rng(Threads.threadid())) if r < coin_flip covid_alert_notifications[end,agent] = 1 #add the notifications for today else covid_alert_notifications[end,agent] = 0 end if sum(covid_alert_notifications[:,agent])>=notification_threshold output_data.daily_total_notifications[t] += 1 time_of_last_alert[agent] = t end `````` 35 ``````end `````` Peter Jentsch committed Jul 01, 2021 36 37 ``````function agent_transition!(t,modelsol,agent, from::AgentStatus,to::AgentStatus) @unpack u_next_inf,output_data, immunization_countdown = modelsol `````` Peter Jentsch committed Jul 05, 2021 38 39 40 `````` immunization_countdown[agent] = -1 output_data.recorded_status_totals[Int(from), t+1] -= 1 output_data.recorded_status_totals[Int(to), t+1] += 1 `````` Peter Jentsch committed Jun 27, 2021 41 42 `````` u_next_inf[agent] = to end `````` Peter Jentsch committed Jul 01, 2021 43 ``````function is_susceptible_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec) `````` Peter Jentsch committed Jun 27, 2021 44 45 `````` @unpack u_inf,inf_network, output_data = modelsol `````` Peter Jentsch committed Jun 30, 2021 46 47 48 49 50 51 52 `````` for mixing_graph in inf_network.graph_list[t], neighbor in neighbors(mixing_graph,agent) if u_inf[neighbor] == Infected infection_threshold = contact_weight(β_vec[Int(agent_demo)],get_weight(mixing_graph,GraphEdge(agent,neighbor))) if rand(Random.default_rng(Threads.threadid())) < infection_threshold && agent_inf_status == Susceptible output_data.daily_cases_by_age[Int(agent_demo),t]+=1 output_data.daily_unvac_cases_by_age[Int(agent_demo),t]+=1 `````` Peter Jentsch committed Jul 01, 2021 53 `````` agent_transition!(t,modelsol, agent, Susceptible,Infected) `````` Peter Jentsch committed Jun 30, 2021 54 `````` return true `````` 55 `````` end `````` Peter Jentsch committed Jun 27, 2021 56 57 58 59 `````` end end return false end `````` Peter Jentsch committed Jul 01, 2021 60 ``````function is_immunized_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec) `````` Peter Jentsch committed Jun 30, 2021 61 62 63 64 65 66 67 `````` @unpack u_inf,inf_network, output_data = modelsol for mixing_graph in inf_network.graph_list[t], neighbor in neighbors(mixing_graph,agent) if u_inf[neighbor] == Infected infection_threshold = contact_weight(β_vec[Int(agent_demo)],get_weight(mixing_graph,GraphEdge(agent,neighbor))) if rand(Random.default_rng(Threads.threadid())) < infection_threshold && agent_inf_status == Immunized immunization_ineffective = rand(Random.default_rng(Threads.threadid())) < 1- α_vec[Int(agent_demo)] if immunization_ineffective `````` Peter Jentsch committed Jul 01, 2021 68 `````` agent_transition!(t,modelsol,agent, Immunized,Infected) `````` Peter Jentsch committed Jun 30, 2021 69 70 71 `````` output_data.daily_cases_by_age[Int(agent_demo),t]+=1 return true end `````` 72 73 `````` end end `````` Peter Jentsch committed Jun 30, 2021 74 75 76 77 78 79 `````` end return false end function vaccinate_agent!(t,agent,immunization_countdown) if immunization_countdown[agent] == -1 #not currently in countdown immunization_countdown[agent] = 14 `````` 80 81 `````` end end `````` 82 `````` `````` Peter Jentsch committed Sep 24, 2021 83 ``````function update_vaccination_opinion_state!(t,agent,agent_demo,modelsol,vaccinate_today_flag,total_infections,π_base_vec) `````` Peter Jentsch committed Jul 01, 2021 84 85 `````` @unpack infection_introduction_day,immunization_intervals, η,Γ,ζ, ω, ω_en,ξ = modelsol.params @unpack immunization_countdown,demographics,time_of_last_alert, soc_network,u_vac,u_next_vac,is_app_user,output_data = modelsol `````` Peter Jentsch committed Jun 30, 2021 86 87 88 89 90 91 `````` random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t]) if !isempty(neighbors(random_soc_network,agent)) random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,agent)) app_vac_payoff = 0.0 if is_app_user[agent] && time_of_last_alert[agent]>=0 output_data.daily_total_notified_agents[t] += 1 `````` Peter Jentsch committed Jul 09, 2021 92 93 `````` app_vac_payoff = Γ^((t - time_of_last_alert[agent])) * (η + total_infections*ω_en) `````` Peter Jentsch committed Jun 30, 2021 94 `````` end `````` Peter Jentsch committed Jul 09, 2021 95 `````` output_data.total_app_weight[t] += app_vac_payoff `````` Peter Jentsch committed Jun 30, 2021 96 97 98 99 100 `````` if u_vac[random_neighbour] == u_vac[agent] vac_payoff = π_base_vec[Int(demographics[agent])] + total_infections*ω + app_vac_payoff if u_vac[agent] if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ) u_next_vac[agent] = false `````` Peter Jentsch committed Jun 20, 2021 101 `````` end `````` Peter Jentsch committed Jun 30, 2021 102 103 104 `````` else if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ) u_next_vac[agent] = true `````` Peter Jentsch committed Jul 05, 2021 105 `````` `````` Peter Jentsch committed Jul 08, 2021 106 `````` if vaccinate_today_flag && t>infection_introduction_day `````` Peter Jentsch committed Jul 01, 2021 107 `````` vaccinate_agent!(t, agent, immunization_countdown) `````` Peter Jentsch committed Sep 24, 2021 108 `````` # fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i) `````` Peter Jentsch committed Jul 08, 2021 109 `````` if !is_app_user[agent] `````` Peter Jentsch committed Sep 24, 2021 110 `````` # fit!(output_data.avg_weighted_degree_of_vaccinators_no_EN[Int(agent_demo)],weighted_degree_of_i) `````` Peter Jentsch committed Jul 08, 2021 111 `````` end `````` Peter Jentsch committed Jul 01, 2021 112 `````` end `````` Peter Jentsch committed May 06, 2021 113 `````` end `````` Peter Jentsch committed Jun 30, 2021 114 `````` end `````` Peter Jentsch committed Jul 01, 2021 115 116 `````` return vac_payoff `````` 117 118 `````` end end `````` Peter Jentsch committed Jul 01, 2021 119 `````` return 0 `````` 120 ``````end `````` Peter Jentsch committed Sep 24, 2021 121 122 123 124 125 126 ``````using WeightedOnlineStats function assortativity(t,network::TimeDepMixingGraph,u_inf) total_weight = 0.0 x_stats = WeightedVariance() y_stats = WeightedVariance() `````` 127 `````` `````` Peter Jentsch committed Jun 16, 2021 128 `````` for g in network.graph_list[t] `````` Peter Jentsch committed Sep 24, 2021 129 130 131 132 133 134 135 `````` for e in edges(g.g) v = src(e) w = dst(e) weight = get_weight(g,GraphEdge(v,w)) total_weight += weight fit!(x_stats,(u_inf[v] == Immunized),weight) fit!(y_stats,(u_inf[v] != Immunized),weight) `````` Peter Jentsch committed May 10, 2021 136 137 `````` end end `````` Peter Jentsch committed Sep 24, 2021 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 `````` σ_x = std(x_stats) σ_y = std(y_stats) x_bar = mean(x_stats) y_bar = mean(y_stats) numerator = 0.0 join_count_11 = 0.0 join_count_01 = 0.0 join_count_00 = 0.0 for g in network.graph_list[t] for e in edges(g.g) v = src(e) w = dst(e) weight = get_weight(g,GraphEdge(v,w)) v_inf = u_inf[v] == Immunized w_inf = u_inf[w] != Immunized numerator += weight*(v_inf - x_bar) * (w_inf - y_bar) #ughhhh join_count_11 += ifelse( u_inf[v] == Immunized && u_inf[w] == Immunized,weight,0) join_count_01 += ifelse(xor( u_inf[v] == Immunized, u_inf[w] == Immunized),weight,0) join_count_00 += ifelse(!(u_inf[v] == Immunized) && !(u_inf[w] == Immunized),weight,0) end end return numerator/(total_weight*σ_x*σ_y),(join_count_11,join_count_01,join_count_00) `````` Peter Jentsch committed May 10, 2021 162 ``````end `````` 163 `````` `````` Peter Jentsch committed Jun 14, 2021 164 165 166 167 168 169 170 171 172 173 174 175 176 177 ``````function sample_initial_nodes(nodes,graphs,I_0_fraction) weighted_degrees = zeros(nodes) for v in 1:nodes for g in graphs for w in neighbors(g,v) weighted_degrees[v] += get_weight(g,GraphEdge(v,w)) end end end wv = Weights(weighted_degrees ./sum(weighted_degrees)) num = round(Int,nodes*I_0_fraction) init_indices = sample(Random.default_rng(Threads.threadid()), 1:nodes,wv, num; replace = false) return init_indices end `````` 178 `````` `````` Peter Jentsch committed Jul 08, 2021 179 180 181 `````` function agents_step!(t,modelsol) @unpack EN_intervals,immunizing,infection_introduction_day,recovery_rate,β_y,β_m,β_o,α_y,α_m,α_o,π_base_y,π_base_m,π_base_o,ζ,immunization_intervals = modelsol.params `````` Peter Jentsch committed Sep 24, 2021 182 `````` @unpack u_inf,u_vac,u_next_vac,u_next_inf,demographics,inf_network,soc_network, is_app_user, immunization_countdown, output_data = modelsol `````` Peter Jentsch committed Jun 30, 2021 183 `````` `````` Peter Jentsch committed Jul 08, 2021 184 185 186 187 188 189 190 `````` vaccinate_today_flag = !immunizing ? false : (any(t in I for I in immunization_intervals) ? true : false) EN_today_flag = t in EN_intervals ? true : false `````` Peter Jentsch committed Jun 30, 2021 191 192 `````` u_next_inf .= u_inf u_next_vac .= u_vac `````` Peter Jentsch committed Jul 05, 2021 193 `````` modelsol.output_data.recorded_status_totals[:,t+1] .= modelsol.output_data.recorded_status_totals[:,t] `````` Peter Jentsch committed Jul 09, 2021 194 `````` total_infected = modelsol.output_data.recorded_status_totals[Int(Infected), t] `````` Peter Jentsch committed Jun 30, 2021 195 196 197 198 199 200 201 `````` β_vec = @SVector [β_y,β_m,β_o] α_vec = @SVector [α_y,α_m,α_o] if t0 immunization_countdown[agent] -= 1 end `````` Peter Jentsch committed Jul 11, 2021 236 `````` if agent_is_vaccinator `````` Peter Jentsch committed Sep 24, 2021 237 `````` # fit!(output_data.avg_weighted_degree_of_vaccinator_belief[t],weighted_degree_of_i) `````` Peter Jentsch committed Jul 11, 2021 238 `````` end `````` Peter Jentsch committed Sep 24, 2021 239 `````` update_vaccination_opinion_state!(t, agent,agent_demo, modelsol,vaccinate_today_flag, total_infected,π_base_vec) `````` Peter Jentsch committed Jun 30, 2021 240 241 242 243 `````` end output_data.total_vaccinators[t] = count(==(true),u_vac) u_vac .= u_next_vac u_inf .= u_next_inf `````` Peter Jentsch committed Jul 01, 2021 244 `````` # return total_weight `````` Peter Jentsch committed Jul 01, 2021 245 `````` `````` Peter Jentsch committed Jun 30, 2021 246 247 ``````end `````` Peter Jentsch committed Jun 10, 2021 248 `````` `````` Peter Jentsch committed Jun 23, 2021 249 ``````function solve!(modelsol) `````` Peter Jentsch committed Jun 14, 2021 250 `````` init_indices = sample_initial_nodes(modelsol.nodes, modelsol.inf_network.graph_list[begin], modelsol.params.I_0_fraction) `````` Peter Jentsch committed Jul 01, 2021 251 `````` # total_weight = 0.0 `````` Peter Jentsch committed Jul 01, 2021 252 `````` `````` Peter Jentsch committed Jun 14, 2021 253 `````` for t in 1:modelsol.sim_length `````` Peter Jentsch committed Jun 13, 2021 254 `````` #this also resamples the soc network weights since they point to the same objects, but those are never used `````` Peter Jentsch committed Jun 27, 2021 255 `````` `````` Peter Jentsch committed Jun 13, 2021 256 `````` if t>1 `````` Peter Jentsch committed Jun 23, 2021 257 `````` remake_all!(t,modelsol.inf_network,modelsol.index_vectors,modelsol.demographics) `````` Peter Jentsch committed Jun 13, 2021 258 259 260 261 262 263 `````` end if t>modelsol.params.infection_introduction_day if !isempty(init_indices) inf_index = pop!(init_indices) modelsol.u_inf[inf_index] = Infected `````` Peter Jentsch committed Jul 09, 2021 264 `````` modelsol.output_data.recorded_status_totals[Int(Infected),t] += 1 `````` Peter Jentsch committed Jun 13, 2021 265 `````` end `````` Peter Jentsch committed Jun 17, 2021 266 267 `````` end `````` Peter Jentsch committed Jun 30, 2021 268 `````` `````` Peter Jentsch committed Jun 27, 2021 269 `````` `````` Peter Jentsch committed Jul 08, 2021 270 `````` agents_step!(t,modelsol) `````` 271 `````` #advance agent states based on the new network `````` Peter Jentsch committed Jun 13, 2021 272 `````` `````` Peter Jentsch committed Jun 23, 2021 273 `````` record!(t,modelsol) `````` Peter Jentsch committed Jun 27, 2021 274 `````` # display(mean.(modelsol.output_data.avg_weighted_degree)) `````` Peter Jentsch committed Jun 26, 2021 275 `````` `````` 276 `````` end `````` Peter Jentsch committed Jul 01, 2021 277 `````` # display(total_weight) `````` Peter Jentsch committed Jul 08, 2021 278 ``end``