solve.jl 11.7 KB
Newer Older
1
2
function contact_weight(β, contact_time) 
    return 1 - (1-β)^contact_time
3
end
4

5
function Φ(payoff,ξ)
Peter Jentsch's avatar
Peter Jentsch committed
6
    return 1 / (1 + exp(-1*ξ*payoff))
7
8
end

Peter Jentsch's avatar
Peter Jentsch committed
9
10
11


@views function update_alert_durations!(t,agent,modelsol) 
12
    @unpack notification_parameter,notification_threshold = modelsol.params
Peter Jentsch's avatar
Peter Jentsch committed
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's avatar
Peter Jentsch committed
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's avatar
Peter Jentsch committed
36
37
function agent_transition!(t,modelsol,agent, from::AgentStatus,to::AgentStatus) 
    @unpack u_next_inf,output_data, immunization_countdown = modelsol
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
41
42
    u_next_inf[agent] = to
end
Peter Jentsch's avatar
Peter Jentsch committed
43
function is_susceptible_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
44
45
    @unpack u_inf,inf_network, output_data = modelsol

Peter Jentsch's avatar
Peter Jentsch committed
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's avatar
Peter Jentsch committed
53
                agent_transition!(t,modelsol, agent, Susceptible,Infected)
Peter Jentsch's avatar
Peter Jentsch committed
54
                return true
55
            end
56
57
58
59
        end
    end
    return false
end
Peter Jentsch's avatar
Peter Jentsch committed
60
function is_immunized_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
Peter Jentsch's avatar
Peter Jentsch committed
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's avatar
Peter Jentsch committed
68
                    agent_transition!(t,modelsol,agent, Immunized,Infected)
Peter Jentsch's avatar
Peter Jentsch committed
69
70
71
                    output_data.daily_cases_by_age[Int(agent_demo),t]+=1
                    return true
                end
72
73
            end
        end
Peter Jentsch's avatar
Peter Jentsch committed
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

83
function update_vaccination_opinion_state!(t,agent,agent_demo,modelsol,vaccinate_today_flag,total_infections,π_base_vec)
Peter Jentsch's avatar
Peter Jentsch committed
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's avatar
Peter Jentsch committed
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
92
93
            app_vac_payoff = Γ^((t - time_of_last_alert[agent])) * (η + total_infections*ω_en)

Peter Jentsch's avatar
Peter Jentsch committed
94
        end
95
        output_data.total_app_weight[t] += app_vac_payoff
Peter Jentsch's avatar
Peter Jentsch committed
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's avatar
Peter Jentsch committed
101
                end
Peter Jentsch's avatar
Peter Jentsch committed
102
103
104
            else
                if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
                    u_next_vac[agent] = true
105

Peter Jentsch's avatar
Peter Jentsch committed
106
                    if vaccinate_today_flag && t>infection_introduction_day
Peter Jentsch's avatar
Peter Jentsch committed
107
                        vaccinate_agent!(t, agent, immunization_countdown)
108
                        # fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
Peter Jentsch's avatar
Peter Jentsch committed
109
                        if !is_app_user[agent]
110
                            # fit!(output_data.avg_weighted_degree_of_vaccinators_no_EN[Int(agent_demo)],weighted_degree_of_i)
Peter Jentsch's avatar
Peter Jentsch committed
111
                        end
Peter Jentsch's avatar
Peter Jentsch committed
112
                    end
Peter Jentsch's avatar
Peter Jentsch committed
113
                end
Peter Jentsch's avatar
Peter Jentsch committed
114
            end
Peter Jentsch's avatar
Peter Jentsch committed
115
116
            return vac_payoff

117
118
        end
    end
Peter Jentsch's avatar
Peter Jentsch committed
119
    return 0
120
end
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's avatar
Peter Jentsch committed
128
    for g in network.graph_list[t]
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)
136
137
        end
    end
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)
162
end
163

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's avatar
Peter Jentsch committed
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
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's avatar
Peter Jentsch committed
183

Peter Jentsch's avatar
Peter Jentsch committed
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's avatar
Peter Jentsch committed
191
192
    u_next_inf .= u_inf
    u_next_vac .= u_vac
193
    modelsol.output_data.recorded_status_totals[:,t+1] .= modelsol.output_data.recorded_status_totals[:,t]
194
    total_infected = modelsol.output_data.recorded_status_totals[Int(Infected), t]
Peter Jentsch's avatar
Peter Jentsch committed
195
196
197
198
199
200
201
    β_vec = @SVector [β_y,β_m,β_o]
    α_vec = @SVector [α_y,α_m,α_o]
    if t<last(minimum(immunization_intervals))
        π_base_vec = @SVector [π_base_y,π_base_m,π_base_o]
    else
        π_base_vec = @SVector [π_base_y*ζ,π_base_m*ζ,π_base_o*ζ]
    end
202
    if t == last(minimum(immunization_intervals)) + 14
203
204
205
206
207
208
209
        weighted_assortativity, join_counts = assortativity(t,soc_network,u_inf)
        fit!(modelsol.output_data.assortativity_at_end,weighted_assortativity)
        fit!(modelsol.output_data.vac_join_11_counts_at_end,join_counts[1])
        fit!(modelsol.output_data.vac_join_01_counts_at_end,join_counts[2])
        fit!(modelsol.output_data.vac_join_00_counts_at_end,join_counts[3])
    end        
    for (agent,(agent_inf_status,agent_is_vaccinator,agent_demo)) in enumerate(zip(u_inf,u_vac,demographics))       
Peter Jentsch's avatar
Peter Jentsch committed
210
        if agent_inf_status == Susceptible 
Peter Jentsch's avatar
Peter Jentsch committed
211
            is_susceptible_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
Peter Jentsch's avatar
Peter Jentsch committed
212
            if vaccinate_today_flag && agent_is_vaccinator && t<(infection_introduction_day-1)
213
                # fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
Peter Jentsch's avatar
Peter Jentsch committed
214
                if !is_app_user[agent]
215
                    # fit!(output_data.avg_weighted_degree_of_vaccinators_no_EN[Int(agent_demo)],weighted_degree_of_i)
Peter Jentsch's avatar
Peter Jentsch committed
216
                end
Peter Jentsch's avatar
Peter Jentsch committed
217
218
219
220
                vaccinate_agent!(t, agent, immunization_countdown)
            end

        elseif agent_inf_status == Immunized
Peter Jentsch's avatar
Peter Jentsch committed
221
            is_immunized_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
Peter Jentsch's avatar
Peter Jentsch committed
222
223
        elseif agent_inf_status == Infected
            if rand(Random.default_rng(Threads.threadid())) < recovery_rate
Peter Jentsch's avatar
Peter Jentsch committed
224
                agent_transition!(t,modelsol, agent, Infected,Recovered)
Peter Jentsch's avatar
Peter Jentsch committed
225
226
            end
        end
Peter Jentsch's avatar
Peter Jentsch committed
227
228
229
        if EN_today_flag && is_app_user[agent]
            update_alert_durations!(t, agent,modelsol)
        end
Peter Jentsch's avatar
Peter Jentsch committed
230
231
        if immunization_countdown[agent] == 0 
            output_data.daily_immunized_by_age[Int(agent_demo),t] += 1
Peter Jentsch's avatar
Peter Jentsch committed
232
            agent_transition!(t,modelsol, agent, Susceptible,Immunized)
Peter Jentsch's avatar
Peter Jentsch committed
233
234
235
        elseif immunization_countdown[agent]>0
            immunization_countdown[agent] -= 1
        end
Peter Jentsch's avatar
Peter Jentsch committed
236
        if agent_is_vaccinator
237
            # fit!(output_data.avg_weighted_degree_of_vaccinator_belief[t],weighted_degree_of_i)
Peter Jentsch's avatar
Peter Jentsch committed
238
        end
239
        update_vaccination_opinion_state!(t, agent,agent_demo, modelsol,vaccinate_today_flag, total_infected,π_base_vec)
Peter Jentsch's avatar
Peter Jentsch committed
240
241
242
243
    end
    output_data.total_vaccinators[t] = count(==(true),u_vac) 
    u_vac .= u_next_vac
    u_inf .= u_next_inf
244
    # return total_weight
Peter Jentsch's avatar
Peter Jentsch committed
245

Peter Jentsch's avatar
Peter Jentsch committed
246
247
end

Peter Jentsch's avatar
Peter Jentsch committed
248

Peter Jentsch's avatar
Peter Jentsch committed
249
function solve!(modelsol)
250
    init_indices = sample_initial_nodes(modelsol.nodes, modelsol.inf_network.graph_list[begin], modelsol.params.I_0_fraction)
251
    # total_weight = 0.0
Peter Jentsch's avatar
Peter Jentsch committed
252
    
253
    for t in 1:modelsol.sim_length  
Peter Jentsch's avatar
Peter Jentsch committed
254
        #this also resamples the soc network weights since they point to the same objects, but those are never used
255

Peter Jentsch's avatar
Peter Jentsch committed
256
        if t>1
Peter Jentsch's avatar
Peter Jentsch committed
257
            remake_all!(t,modelsol.inf_network,modelsol.index_vectors,modelsol.demographics)
Peter Jentsch's avatar
Peter Jentsch committed
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
264
                modelsol.output_data.recorded_status_totals[Int(Infected),t] += 1
Peter Jentsch's avatar
Peter Jentsch committed
265
            end
Peter Jentsch's avatar
Peter Jentsch committed
266
267
           
        end
Peter Jentsch's avatar
Peter Jentsch committed
268

269

Peter Jentsch's avatar
Peter Jentsch committed
270
        agents_step!(t,modelsol)
271
        #advance agent states based on the new network
Peter Jentsch's avatar
Peter Jentsch committed
272
        
Peter Jentsch's avatar
Peter Jentsch committed
273
        record!(t,modelsol)
274
        # display(mean.(modelsol.output_data.avg_weighted_degree))
Peter Jentsch's avatar
Peter Jentsch committed
275

276
    end
277
    # display(total_weight)
Peter Jentsch's avatar
Peter Jentsch committed
278
end