solve.jl 8.38 KB
Newer Older
1

2
3
4
#remove Base.@propagate_inbounds during tests if you get segfaults or change anything,
# as that might mean there is an indexing bug that is not being caught,
#and therefore accessing the wrong memory is causing background errors
5

6
7
function contact_weight(β, contact_time) 
    return 1 - (1-β)^contact_time
8
end
9

10
function Φ(payoff,ξ)
Peter Jentsch's avatar
Peter Jentsch committed
11
    return 1 / (1 + exp(-1*ξ*payoff))
12
13
end

14
15
Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # Base.@propagate_inbounds 
    @unpack notification_parameter,notification_threshold = modelsol.params
Peter Jentsch's avatar
Peter Jentsch committed
16
    @unpack time_of_last_alert, app_user_index,inf_network,covid_alert_notifications,app_user, output_data = modelsol
Peter Jentsch's avatar
Peter Jentsch committed
17
    
Peter Jentsch's avatar
Peter Jentsch committed
18
    for (i,node) in enumerate(app_user_index)
19
        for j in 2:14
20
            covid_alert_notifications[j-1,i] = covid_alert_notifications[j,i] #shift them all back  
21
        end
22
23
        total_weight_i = 0
        for mixing_graph in inf_network.graph_list[t]
Peter Jentsch's avatar
Peter Jentsch committed
24
            for (j,w) in neighbors_and_weights(mixing_graph,node)
25
                if app_user[j]
Peter Jentsch's avatar
Peter Jentsch committed
26
                    total_weight_i += w[]
27
                end
28
29
            end
        end
30
        coin_flip = 1 - (1 - notification_parameter)^total_weight_i
Peter Jentsch's avatar
Peter Jentsch committed
31
        r = rand(Random.default_rng(Threads.threadid()))
32
33
34
35
36
37
        if r < coin_flip
            covid_alert_notifications[end,i] = 1  #add the notifications for today
        else
            covid_alert_notifications[end,i] = 0
        end
        if sum(covid_alert_notifications[:,i])>=notification_threshold
Peter Jentsch's avatar
Peter Jentsch committed
38
            output_data.daily_total_notifications[t] += 1
39
40
41
42
            time_of_last_alert[i] = t
        end
    end
end
43

Peter Jentsch's avatar
Peter Jentsch committed
44
Base.@propagate_inbounds @views function update_infection_state!(t,modelsol; record_degrees = false)
45
    @unpack β_y,β_m,β_o,α_y,α_m,α_o,recovery_rate,immunizing,immunization_begin_day = modelsol.params
Peter Jentsch's avatar
Peter Jentsch committed
46
    @unpack u_inf,u_vac,u_next_inf,demographics,inf_network,status_totals, immunization_countdown, output_data = modelsol
Peter Jentsch's avatar
Peter Jentsch committed
47
    
48
    function agent_transition!(node, from::AgentStatus,to::AgentStatus) 
49
        immunization_countdown[node] = -1 
50
51
        status_totals[Int(from)] -= 1
        status_totals[Int(to)] += 1
52
53
54
        u_next_inf[node] = to
    end
    u_next_inf .= u_inf
55
56
57
58
    for i in 1:modelsol.nodes
        agent_status = u_inf[i]
        is_vaccinator = u_vac[i]
        agent_demo = demographics[i]
59
60
61
62
63
        if agent_status == Susceptible && is_vaccinator && immunizing && immunization_countdown[i] == -1 && t> immunization_begin_day
            immunization_countdown[i] = 14
        end
        if  agent_status == Susceptible || agent_status == Immunized 
            for mixing_graph in inf_network.graph_list[t]
Peter Jentsch's avatar
Peter Jentsch committed
64
                for (j,w) in neighbors_and_weights(mixing_graph,i)   
65
66
67
                    if u_inf[j] == Infected && u_next_inf[i] != Infected
                        β_vec = (β_y,β_m,β_o)
                        α_vec = (α_y,α_m,α_o)
Peter Jentsch's avatar
Peter Jentsch committed
68
                        if rand(Random.default_rng(Threads.threadid())) < contact_weight(β_vec[Int(agent_demo)],w[])
Peter Jentsch's avatar
Peter Jentsch committed
69
                            if agent_status == Immunized && rand(Random.default_rng(Threads.threadid())) < 1- α_vec[Int(agent_demo)]
70
                                agent_transition!(i, Immunized,Infected)
Peter Jentsch's avatar
Peter Jentsch committed
71
                                output_data.daily_cases_by_age[Int(agent_demo),t]+=1
Peter Jentsch's avatar
Peter Jentsch committed
72

73
                            elseif agent_status == Susceptible
Peter Jentsch's avatar
Peter Jentsch committed
74
75
                                output_data.daily_cases_by_age[Int(agent_demo),t]+=1
                                output_data.daily_unvac_cases_by_age[Int(agent_demo),t]+=1
76
77
                                agent_transition!(i, Susceptible,Infected)
                            end
78
79
80
81
82
                        end
                    end
                end
            end
        elseif agent_status == Infected
Peter Jentsch's avatar
Peter Jentsch committed
83
            if rand(Random.default_rng(Threads.threadid())) < recovery_rate
84
               agent_transition!(i, Infected,Recovered)
85
86
            end
        end
Peter Jentsch's avatar
Peter Jentsch committed
87
88

        weighted_degree_of_i::Int = output_data.record_degrees_flag ? weighted_degree(t,i,inf_network) : 0
89
        if immunization_countdown[i] == 0 
Peter Jentsch's avatar
Peter Jentsch committed
90
91
            output_data.daily_immunized_by_age[Int(agent_demo),t] += 1
            fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
92
            agent_transition!(i, Susceptible,Immunized)
Peter Jentsch's avatar
Peter Jentsch committed
93

94
        elseif immunization_countdown[i]>0
Peter Jentsch's avatar
Peter Jentsch committed
95
            fit!(output_data.avg_weighted_degree[Int(agent_demo)],weighted_degree_of_i)
96
            immunization_countdown[i] -= 1
Peter Jentsch's avatar
Peter Jentsch committed
97
        else
Peter Jentsch's avatar
Peter Jentsch committed
98
            fit!(output_data.avg_weighted_degree[Int(agent_demo)],weighted_degree_of_i)
99
        end
100
101
    end
end
102
103


Peter Jentsch's avatar
Peter Jentsch committed
104
Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
Peter Jentsch's avatar
Peter Jentsch committed
105
    @unpack infection_introduction_day, π_base_y,π_base_m,π_base_o, η,Γ,ζ, ω, ω_en,ξ = modelsol.params
Peter Jentsch's avatar
Peter Jentsch committed
106
    @unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list,output_data = modelsol
Peter Jentsch's avatar
Peter Jentsch committed
107
   
108
    for i in 1:nodes
Peter Jentsch's avatar
Peter Jentsch committed
109
110
        π_base = t<infection_introduction_day ? 
        (π_base_y,π_base_m,π_base_o) : 
111
        (π_base_y*ζ,π_base_m*ζ,π_base_o*ζ)
Peter Jentsch's avatar
Peter Jentsch committed
112
        random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
Peter Jentsch's avatar
Peter Jentsch committed
113
        if !isempty(neighbors(random_soc_network,i))
Peter Jentsch's avatar
Peter Jentsch committed
114
                random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network,i))
Peter Jentsch's avatar
Peter Jentsch committed
115
116
                app_vac_payoff = 0.0    
                if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
Peter Jentsch's avatar
Peter Jentsch committed
117
                    output_data.daily_total_notified_agents[t] += 1
Peter Jentsch's avatar
Peter Jentsch committed
118
119
120
                    app_vac_payoff = Γ^((t - time_of_last_alert[app_user_list[i]])) * (η + total_infections*ω_en) 
                    # display(t - time_of_last_alert[app_user_list[i]])
                end
Peter Jentsch's avatar
Peter Jentsch committed
121
                if u_vac[random_neighbour] == u_vac[i]
Peter Jentsch's avatar
Peter Jentsch committed
122
                    vac_payoff = π_base[Int(demographics[i])] + total_infections*ω + app_vac_payoff
Peter Jentsch's avatar
Peter Jentsch committed
123
124
125
126
127
128
129
130
131
132
                    if u_vac[i]
                        # display(1 - Φ(vac_payoff,ξ))
                        if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ)
                            u_next_vac[i] = false
                        end
                    else
                        # display( Φ(vac_payoff,ξ))
                        if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
                            u_next_vac[i] = true
                        end
Peter Jentsch's avatar
Peter Jentsch committed
133
134
                    end
                end
135
136
        end
    end
Peter Jentsch's avatar
Peter Jentsch committed
137
    output_data.total_vaccinators[t] = count(==(true),u_vac) 
138
139
end

Peter Jentsch's avatar
Peter Jentsch committed
140
function weighted_degree(t,node,network::TimeDepMixingGraph)
141
    weighted_degree = 0
Peter Jentsch's avatar
Peter Jentsch committed
142
    for g in network.graph_list[t]
Peter Jentsch's avatar
Peter Jentsch committed
143
144
        for (j,w) in neighbors_and_weights(mixing_graph,node)
            weighted_degree += w[]
145
146
147
148
        end
    end
    return weighted_degree
end
149

150
151
152
153
function sample_initial_nodes(nodes,graphs,I_0_fraction)
    weighted_degrees = zeros(nodes)
    for v in 1:nodes
        for g in graphs
Peter Jentsch's avatar
Peter Jentsch committed
154
155
            for (j,w) in neighbors_and_weights(g,v)
                weighted_degrees[v] += w[]
156
157
158
159
160
161
162
163
            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    
164

Peter Jentsch's avatar
Peter Jentsch committed
165

Peter Jentsch's avatar
Peter Jentsch committed
166
function solve!(modelsol)
167
168
    init_indices = sample_initial_nodes(modelsol.nodes, modelsol.inf_network.graph_list[begin], modelsol.params.I_0_fraction)
    for t in 1:modelsol.sim_length  
Peter Jentsch's avatar
Peter Jentsch committed
169
170
        #this also resamples the soc network weights since they point to the same objects, but those are never used
        if t>1
Peter Jentsch's avatar
Peter Jentsch committed
171
            remake_all!(t,modelsol.inf_network,modelsol.index_vectors,modelsol.demographics)
Peter Jentsch's avatar
Peter Jentsch committed
172
173
174
175
176
177
178
179
        end

        if t>modelsol.params.infection_introduction_day
            if !isempty(init_indices)
                inf_index = pop!(init_indices)
                modelsol.u_inf[inf_index] = Infected
                modelsol.status_totals[Int(Infected)] += 1
            end
Peter Jentsch's avatar
Peter Jentsch committed
180
181
182
           
        end
        if t>modelsol.params.immunization_begin_day
Peter Jentsch's avatar
Peter Jentsch committed
183
184
            update_alert_durations!(t,modelsol)
        end
Peter Jentsch's avatar
Peter Jentsch committed
185
 
Peter Jentsch's avatar
Peter Jentsch committed
186
        update_infection_state!(t,modelsol; record_degrees = true)
Peter Jentsch's avatar
Peter Jentsch committed
187
        update_vaccination_opinion_state!(t,modelsol,modelsol.status_totals[Int(Infected)])
188
        #advance agent states based on the new network
Peter Jentsch's avatar
Peter Jentsch committed
189
190

        
Peter Jentsch's avatar
Peter Jentsch committed
191
        record!(t,modelsol)
Peter Jentsch's avatar
Peter Jentsch committed
192
193
194

        modelsol.u_vac .= modelsol.u_next_vac
        modelsol.u_inf .= modelsol.u_next_inf
195
196
197
198
    end
end