solve.jl 7.31 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
11
function Φ(payoff,ξ)
    return 1 / (exp(-1*ξ*payoff))
12
13
end

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

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

75
                            elseif agent_status == Susceptible
76
                                modelsol.daily_cases_by_age[Int(agent_demo)]+=1
Peter Jentsch's avatar
Peter Jentsch committed
77
                                modelsol.daily_unvac_cases_by_age[Int(agent_demo)]+=1
78
79
                                agent_transition!(i, Susceptible,Infected)
                            end
80
81
82
83
84
                        end
                    end
                end
            end
        elseif agent_status == Infected
Peter Jentsch's avatar
Peter Jentsch committed
85
            if rand(Random.default_rng(Threads.threadid())) < recovery_rate
86
               agent_transition!(i, Infected,Recovered)
87
88
            end
        end
89
90
        
        if immunization_countdown[i] == 0 
91
            modelsol.daily_immunized_by_age[Int(agent_demo)] += 1
92
93
94
95
            agent_transition!(i, Susceptible,Immunized)
        elseif immunization_countdown[i]>0
            immunization_countdown[i] -= 1
        end
96
97
    end
end
98
99


Peter Jentsch's avatar
Peter Jentsch committed
100
Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
101
    @unpack π_base_y,π_base_m,π_base_o, η,γ, κ, ω, ω_en,γ,ξ = modelsol.params
102
103
104
    @unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list = modelsol
    app_user_pointer = 0
    for i in 1:nodes
Peter Jentsch's avatar
Peter Jentsch committed
105
        π_base = @SVector [π_base_y,π_base_m,π_base_o]
106
107
        vac_payoff = 0
        num_soc_nbrs = 0
Peter Jentsch's avatar
Peter Jentsch committed
108
        random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
109

Peter Jentsch's avatar
Peter Jentsch committed
110
        if !isempty(neighbors(random_soc_network,i))
Peter Jentsch's avatar
Peter Jentsch committed
111
            random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,i))
Peter Jentsch's avatar
Peter Jentsch committed
112
            if u_vac[random_neighbour] == u_vac[i]
113
                vac_payoff += π_base[Int(demographics[i])] + total_infections*ω             
Peter Jentsch's avatar
Peter Jentsch committed
114
                if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
115
                    vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]])) * (η + total_infections*ω_en) 
Peter Jentsch's avatar
Peter Jentsch committed
116
                end
117
        
Peter Jentsch's avatar
Peter Jentsch committed
118
                if u_vac[i]
Peter Jentsch's avatar
Peter Jentsch committed
119
                    if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ)
Peter Jentsch's avatar
Peter Jentsch committed
120
121
122
                        u_next_vac[i] = false
                    end
                else
Peter Jentsch's avatar
Peter Jentsch committed
123
                    if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
Peter Jentsch's avatar
Peter Jentsch committed
124
125
126
                        u_next_vac[i] = true
                    end
                end
127
128
129
            end
        end
    end
130

Peter Jentsch's avatar
Peter Jentsch committed
131
    modelsol.daily_vaccinators = count(==(true),u_vac) #could maybe make this more efficient
132

133
134
end

135
136
137
138
function weighted_degree(node,network::TimeDepMixingGraph)
    weighted_degree = 0
    for g_list in network.graph_list
        for g in g_list
Peter Jentsch's avatar
Peter Jentsch committed
139
            for j in neighbors(g,node)
140
141
142
143
144
145
                weighted_degree += get_weight(g,GraphEdge(node,j))
            end
        end
    end
    return weighted_degree
end
146
147

function agents_step!(t,modelsol)
148

149
150
    remake!(modelsol.inf_network,modelsol.index_vectors)
    for network in modelsol.inf_network.graph_list[t] #this also resamples the soc network weights since they point to the same objects, but those are never used
151
        sample_mixing_graph!(network) #get new contact weights
152
    end
Peter Jentsch's avatar
Peter Jentsch committed
153
    if t == modelsol.params.infection_introduction_day
Peter Jentsch's avatar
Peter Jentsch committed
154
        init_indices = rand(Random.default_rng(Threads.threadid()), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
155
156
157
        modelsol.u_inf[init_indices] .= Infected
        modelsol.status_totals[Int(Infected)] += length(init_indices)
    end
Peter Jentsch's avatar
Peter Jentsch committed
158
    if t>modelsol.params.infection_introduction_day
159
160
        update_alert_durations!(t,modelsol)
    end
161
    update_vaccination_opinion_state!(t,modelsol,modelsol.status_totals[Int(Infected)])
162
    update_infection_state!(t,modelsol)
163

164
165
166
167
168
169
    modelsol.u_vac .= modelsol.u_next_vac
    modelsol.u_inf .= modelsol.u_next_inf
end



170
171
function solve!(modelsol,recordings...)

172
    for t in 1:modelsol.sim_length
Peter Jentsch's avatar
Peter Jentsch committed
173
        agents_step!(t,modelsol) 
174
        #advance agent states based on the new network
175
176
177
        for recording in recordings
            record!(t,modelsol,recording)
        end
178
179
180
181
    end
end