solve.jl 7.64 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
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 infection_introduction_day, π_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
106
        π_base = t<infection_introduction_day ? 
        (π_base_y,π_base_m,π_base_o) : 
107
        (π_base_y*ζ,π_base_m*ζ,π_base_o*ζ)
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))
111
            for _ = 1:1
Peter Jentsch's avatar
Peter Jentsch committed
112
113
114
115
116
117
                random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,i))

                if u_vac[random_neighbour] == u_vac[i]
                    vac_payoff = π_base[Int(demographics[i])] + total_infections*ω             
                    if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
                        vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]])) * (η + total_infections*ω_en) 
Peter Jentsch's avatar
Peter Jentsch committed
118
                    end
Peter Jentsch's avatar
Peter Jentsch committed
119
120
121
122
123
124
125
126
127
128
129
            
                    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
130
                    end
Peter Jentsch's avatar
Peter Jentsch committed
131
                    break
Peter Jentsch's avatar
Peter Jentsch committed
132
                end
Peter Jentsch's avatar
Peter Jentsch committed
133

134
135
136
            end
        end
    end
137

Peter Jentsch's avatar
Peter Jentsch committed
138
    modelsol.daily_vaccinators = count(==(true),u_vac) 
139
140
end

141
142
143
144
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
145
            for j in neighbors(g,node)
146
147
148
149
150
151
                weighted_degree += get_weight(g,GraphEdge(node,j))
            end
        end
    end
    return weighted_degree
end
152

153
function agents_step!(t,modelsol, init_indices)
154

Peter Jentsch's avatar
Peter Jentsch committed
155
    if t>modelsol.params.infection_introduction_day
156
157
158
159
160
        if !isempty(init_indices)
            inf_index = pop!(init_indices)
            modelsol.u_inf[inf_index] = Infected
            modelsol.status_totals[Int(Infected)] += 1
        end
161
162
        update_alert_durations!(t,modelsol)
    end
163
    update_vaccination_opinion_state!(t,modelsol,modelsol.status_totals[Int(Infected)])
164
    update_infection_state!(t,modelsol)
165

166
167
    modelsol.u_vac .= modelsol.u_next_vac
    modelsol.u_inf .= modelsol.u_next_inf
Peter Jentsch's avatar
Peter Jentsch committed
168
169
170
171
172

    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
        sample_mixing_graph!(network) #get new contact weights
    end
173
174
175
176
end



177
function solve!(modelsol,recordings...)
178
179
    init_indices = rand(Random.default_rng(Threads.threadid()), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
    @show length(init_indices)
180
    for t in 1:modelsol.sim_length
181
        agents_step!(t,modelsol,init_indices) 
182
        #advance agent states based on the new network
183
184
185
        for recording in recordings
            record!(t,modelsol,recording)
        end
186
187
188
189
    end
end