remove p_en and add recording of age-structured vaccination

CovidAlertVaccinationModel/Manifest.toml

@@ -845,6 +845,12 @@ git-tree-sha1 = "f82a0e71f222199de8e9eb9a09977bd0767d52a0"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
 version = "0.11.0"
 
+[[PairPlots]]
+deps = ["Measures", "NamedTupleTools", "Printf", "RecipesBase", "Statistics", "StatsBase", "Tables"]
+git-tree-sha1 = "b49fc32c74710a8fa62d3ce46cd689c98d01ada1"
+uuid = "43a3c2be-4208-490b-832a-a21dcd55d7da"
+version = "0.1.0"
+
 [[Pandas]]
 deps = ["Compat", "DataValues", "IteratorInterfaceExtensions", "Lazy", "Pkg", "PyCall", "Statistics", "TableTraits", "TableTraitsUtils", "Test"]
 git-tree-sha1 = "40b84de2260989a63792704014b99a86bb8f1425"

CovidAlertVaccinationModel/Project.toml

@@ -30,6 +30,7 @@ LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
 NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
 OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e"
+PairPlots = "43a3c2be-4208-490b-832a-a21dcd55d7da"
 Pandas = "eadc2687-ae89-51f9-a5d9-86b5a6373a9c"
 Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

CovidAlertVaccinationModel/abm_timeseries.jl

@@ -11,4 +11,4 @@ function solve_and_plot_parameters()
 end
 
 out = solve_and_plot_parameters()
-println(sum.(eachrow(mean.(out.new_ymo_immunization))))
\ No newline at end of file
+println(sum.(eachrow(mean.(out.daily_immunized_by_age))))
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/model_setup.jl

@@ -5,8 +5,8 @@ function get_parameters()
         num_households = 5000,
         I_0_fraction = 0.002,
         base_transmission_probability = 0.001,
-        recovery_rate = 1/7,
-        immunization_loss_prob = 0.0001, #mean time of 6 months
+        recovery_rate = 1/5,
+        immunization_loss_prob = 0.0, #no immunization loss
         π_base_y = -5.0,
         π_base_m = -5.0,
         π_base_o = 1.0,
@@ -14,7 +14,6 @@ function get_parameters()
         κ = 0.0,
         ω = 0.00005,
         ω_en = 0.00,
-        ρ_en = [0.0,0.0,0.0],
         γ = 0.0,
         β = 5.0,
         notification_parameter = 0.001,
@@ -69,8 +68,8 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
     app_user_index::Vector{Int}
     status_totals::Vector{Int}
     daily_vaccinators::Int
-    daily_infected::Int
-    daily_immunizations_by_age::Vector{Int}
+    daily_cases_by_age::Vector{Int}
+    daily_immunized_by_age::Vector{Int}
     ws_matrix_tuple::WSMixingDist
     rest_matrix_tuple::RestMixingDist
     immunization_countdown::Vector{Int}
@@ -116,7 +115,7 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
             app_user_index,
             status_totals,
             0,
-            0,
+            [0,0,0],
             [0,0,0],
             ws_matrix_tuple,
             rest_matrix_tuple,

CovidAlertVaccinationModel/src/ABM/output.jl

@@ -18,10 +18,10 @@ DebugRecorder should store everything we might want to know about the model outp
 """
 struct DebugRecorder{ElType,ArrT1<:AbstractArray{ElType},ArrT2<:AbstractArray{ElType},ArrT3<:AbstractArray{ElType}} <: AbstractRecorder{ElType}
     recorded_status_totals::ArrT1
-    daily_cases::ArrT2
+    daily_cases_by_age::ArrT3
     total_vaccinators::ArrT2
     mean_time_since_last_notification::ArrT2
-    new_ymo_immunization::ArrT3
+    daily_immunized_by_age::ArrT3
 end
 # struct AgeSpecificVaccination{ElType,ArrType1<:AbstractArray{ElType}} <: AbstractRecorder{ElType}
 #     new_ymo_vaccinations::ArrType1
@@ -36,14 +36,16 @@ function DebugRecorder(sim_length)
     total_vaccinators = Vector{Int}(undef,sim_length)
     mean_time_since_last_notification = Vector{Int}(undef,sim_length)
     daily_cases= Vector{Int}(undef,sim_length)
-    new_ymo_immunization = @LArray Array{Int}(undef,3,sim_length) (Y = (1,:),M = (2,:),O = (3,:))
+    daily_ymo_by_age = @LArray Array{Int}(undef,3,sim_length) (Y = (1,:),M = (2,:),O = (3,:))
     
+    daily_cases_by_age = @LArray Array{Int}(undef,3,sim_length) (Y = (1,:),M = (2,:),O = (3,:))
+    daily_immunized_by_age = @LArray Array{Int}(undef,3,sim_length) (Y = (1,:),M = (2,:),O = (3,:))
     return DebugRecorder(
         state_totals,
-        daily_cases,
+        daily_cases_by_age,
         total_vaccinators,
         mean_time_since_last_notification,
-        new_ymo_immunization
+        daily_immunized_by_age
     )
 
 end
@@ -55,23 +57,24 @@ function DebugRecorder(val::T, sim_length) where T
     totals = [copy(val) for i in 1:4, j in 1:sim_length]
     state_totals = @LArray totals (S = (1,:),I = (2,:),R = (3,:), V = (4,:))
     total_vaccinators = [copy(val) for j in 1:sim_length]
-
     mean_time_since_last_notification = [copy(val) for j in 1:sim_length]
-    daily_cases =  [copy(val) for j in 1:sim_length]
-    totals_immunization = [copy(val) for i in 1:3, j in 1:sim_length]
-    new_ymo_immunization = @LArray totals_immunization (Y = (1,:),M = (2,:),O = (3,:))
+    ymo = [copy(val) for i in 1:3, j in 1:sim_length]
+
+    totals_ymo = [copy(val) for i in 1:3, j in 1:sim_length]
+    daily_cases_by_age = @LArray deepcopy(totals_ymo) (Y = (1,:),M = (2,:),O = (3,:))
+    daily_immunized_by_age = @LArray deepcopy(totals_ymo) (Y = (1,:),M = (2,:),O = (3,:))
     return DebugRecorder(
         state_totals,
-        daily_cases,
+        daily_cases_by_age,
         total_vaccinators,
         mean_time_since_last_notification,
-        new_ymo_immunization
+        daily_immunized_by_age
     )
 end
 
 function record!(t,modelsol, recorder::DebugRecorder)
-    recorder.total_vaccinators[t] = modelsol.daily_vaccinators
-    recorder.daily_cases[t] = modelsol.daily_infected
+    recorder.total_vaccinators[t] =  modelsol.daily_vaccinators
+    recorder.daily_cases_by_age[:,t] .= modelsol.daily_cases_by_age
     recorder.recorded_status_totals[:,t] .= modelsol.status_totals
     alerts = filter(>(0),modelsol.time_of_last_alert)
     if !isempty(alerts)
@@ -81,7 +84,7 @@ function record!(t,modelsol, recorder::DebugRecorder)
         recorder.mean_time_since_last_notification[t] = 0
     end
     # display(modelsol.daily_immunizations_by_age)
-    recorder.new_ymo_immunization[:,t].=modelsol.daily_immunizations_by_age
+    recorder.daily_immunized_by_age[:,t] .= modelsol.daily_immunized_by_age
 end
 
 function record!(t,modelsol, recorder::Nothing)
@@ -123,10 +126,10 @@ function plot_model(varname,univariate_series, output_list::Vector{T},infection_
         (data.recorded_status_totals.V, "Vaccinated over time"),
         (data.total_vaccinators,  "No. vaccinators over time"),
         (data.mean_time_since_last_notification, "Mean time since last notification"),
-        (data.daily_cases,"Daily (incident) cases"),
-        (data.new_ymo_immunization.Y, "new Y vaccinations each day"),
-        (data.new_ymo_immunization.M, "new M vaccinations each day"),
-        (data.new_ymo_immunization.O, "new O vaccinations each day"),
+        (sum.(eachcol(mean.(data.daily_cases_by_age))),"Daily (incident) cases"),
+        (data.daily_immunized_by_age.Y, "new Y vaccinations each day"),
+        (data.daily_immunized_by_age.M, "new M vaccinations each day"),
+        (data.daily_immunized_by_age.O, "new O vaccinations each day"),
     ]
 
     l = length(ts_list(output_list[1]))

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

 using KissABC
-using CurveFit
-using Loess
 
 const target_cumulative_vac_proportion = 0.33
-const vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
-const ymo_vac = @SVector [0.255,0.278,0.602]
+const vaccination_data = [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
+const ymo_vac = [0.255,0.278,0.602]
+const ymo_attack_rate = []
+
 
 function solve_w_parameters(default_p, p_names, new_p_list)
     new_params = merge(default_p, NamedTuple{p_names}(ntuple(i -> new_p_list[i],length(p_names))))
@@ -12,38 +12,6 @@ function solve_w_parameters(default_p, p_names, new_p_list)
     model = abm(new_params,out)
     return out, model
 end
-@views function exp_growth_rate_estimation(incident_cases)
-    diff(l) = [l[i] - l[i-1] for i=2:length(l)]
-    l = length(incident_cases)
-    loess_interp = loess(1:l,incident_cases)
-    interp_incident_cases = Loess.predict(loess_interp, collect(1.0:l))
-    second_deriative =diff(diff(interp_incident_cases))
-    inflection_pt = any(second_deriative.<=-1) ? findfirst(<=(-1),second_deriative) : l  #completely heuristic
-    ic,exp_growth_rate = exp_fit(1:inflection_pt, incident_cases[1:inflection_pt]) #fit y = a exp(b x), b is exp growth rate
-    return ic,exp_growth_rate
-end
-
-function fit_epi_parameters(p_tuple, target_growth_rate)
-    p_tuple_without_vac = merge(p_tuple,
-        (
-            sim_length = 30,
-            immunization_begin_day = 0,
-            infection_introduction_day = 1,
-            immunizing = false,
-        )
-    )
-    p_names = (:base_transmission_probability,)
-    priors = Factored(Uniform(0.0001,0.005))    
-    function cost(p)
-        output,model = solve_w_parameters(p_tuple_without_vac,p_names ,p)
-        incident_cases = output.daily_cases
-        _,exp_growth_rate = exp_growth_rate_estimation(incident_cases)
-        return (target_growth_rate - exp_growth_rate)^2
-    end
-
-    out = ABCDE(priors,cost,0.01; verbose=true, nparticles= 200,generations=20,  parallel = true) #this one has better NaN handling
-    return NamedTuple{p_names}((out.P.particles,))
-end
 
 function fit_pre_inf_behavioural_parameters(p_tuple)
     samples = 1
@@ -79,15 +47,14 @@ function fit_pre_inf_behavioural_parameters(p_tuple)
         # display(total_preinfection_vaccination)
         return sum((total_preinfection_vaccination .- target_ymo_vac).^2)    
     end
-    #smc(priors,cost; verbose = true, nparticles = 100, parallel = true)# smc(priors,cost; verbose = true, nparticles =300, parallel = true)#
     out = smc(priors,cost; verbose = true, nparticles = 800, parallel = true)#
     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
 end
 
 function fit_post_inf_behavioural_parameters(p_tuple)
     samples = 1
-    p_names = (:ω,)
-    priors = Factored(Uniform(0.0,0.1))
+    p_names = (:ω,:base_transmission_probability)
+    priors = Factored(Uniform(0.0,0.1),Uniform(0.0,0.1))
     #simulation begins in july
     #60 days for opinion dynamics to stabilize, then immunization begins in september,
     #infection introduced at beginning of december
@@ -106,11 +73,15 @@ function fit_post_inf_behavioural_parameters(p_tuple)
     ymo_vac = @SVector [0.255,0.278,0.602]
     function cost(p)
         output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
+
+
         target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
         target_ymo_vac = ymo_vac .* sum(vaccination_data[5:end]) .* target_cumulative_vaccinations
         ymo_vaccination_ts = output.new_ymo_immunization
-        total_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
-        return sum((total_vaccination .- target_ymo_vac).^2)    
+        total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
+
+        final_size = output.recorded_status_totals.R[end]
+        return sum((total_postinf_vaccination .- target_ymo_vac).^2)    
     end
     out =smc(priors,cost; verbose = true, nparticles = 200, parallel = true)# ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=1000,  parallel = true) #this one has better NaN handling
     return NamedTuple{p_names}((out.P.particles,))
@@ -153,18 +124,20 @@ function fit_parameters(default_parameters)
     pre_inf_behaviour_parameters_path =joinpath(PACKAGE_FOLDER,"abm_parameter_fits","pre_inf_behaviour_parameters.dat")
     post_inf_behaviour_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","post_inf_behaviour_parameters.dat")
 
-    seasonal_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(default_parameters,0.073) ##seasonal
-    plot_max_posterior("seasonal", seasonal_transmission_dist,default_parameters)
-    pre_inf_behaviour_parameters = CovidAlertVaccinationModel.fit_pre_inf_behavioural_parameters(default_parameters)
-    fitted_parameters = map(mode,((;seasonal_transmission_dist...,pre_inf_behaviour_parameters...)))
-    fitted_parameter_tuple = merge(default_parameters,fitted_parameters)
+    # seasonal_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(default_parameters,0.073) ##seasonal
+    # plot_max_posterior("seasonal", seasonal_transmission_dist,default_parameters)
+    # pre_inf_behaviour_parameters = CovidAlertVaccinationModel.fit_pre_inf_behavioural_parameters(default_parameters)
+    # fitted_parameters = map(mode,((;seasonal_transmission_dist...,pre_inf_behaviour_parameters...)))
+    # fitted_parameter_tuple = merge(default_parameters,fitted_parameters)
 
-    serialize(seasonal_parameters_path,seasonal_transmission_dist)
-    serialize(pre_inf_behaviour_parameters_path, pre_inf_behaviour_parameters)
+    # serialize(seasonal_parameters_path,seasonal_transmission_dist)
+    # serialize(pre_inf_behaviour_parameters_path, pre_inf_behaviour_parameters)
 
-    post_inf_behaviour_parameters = fit_post_inf_behavioural_parameters(fitted_parameter_tuple) 
-    serialize(post_inf_behaviour_parameters_path, post_inf_behaviour_parameters)
+    # post_inf_behaviour_parameters = fit_post_inf_behavioural_parameters(fitted_parameter_tuple) 
+    # serialize(post_inf_behaviour_parameters_path, post_inf_behaviour_parameters)
+    
     
+    # visualize_π_base(deserialize(pre_inf_behaviour_parameters_path))
     
 
     fitted_parameters_with_post_inf_behaviour = (;
@@ -173,7 +146,10 @@ function fit_parameters(default_parameters)
         deserialize(post_inf_behaviour_parameters_path)...
     )
     display(map(mode,fitted_parameters_with_post_inf_behaviour))
-    return plot_fitting_posteriors("post_inf_fitting",fitted_parameters_with_post_inf_behaviour,default_parameters)
+
+
+    fitted_sol = plot_fitting_posteriors("post_inf_fitting",fitted_parameters_with_post_inf_behaviour,default_parameters)
+    return fitted_sol
 end 
 
 function plot_max_posterior(fname,particles,parameters)
@@ -195,7 +171,7 @@ function plot_max_posterior(fname,particles,parameters)
     p = plot(hist;legend = false)
     savefig(p,"$(fname)_posterior.pdf")
 end
-
+using PairPlots
 function plot_fitting_posteriors(fname,particles_tuple,parameters)
     p_tuple_adjust = merge(parameters,
         (
@@ -212,10 +188,28 @@ function plot_fitting_posteriors(fname,particles_tuple,parameters)
     
     plts = [plot() for i in 1:length(particles_tuple)]
     for (plt,(k,v)) in zip(plts,pairs(particles_tuple))
-        hist = StatsBase.fit(Histogram,v; nbins = 40)
+        hist = StatsBase.fit(Histogram,v; nbins = 30)
         plot!(plt,hist;legend = false,xlabel = k)            
     end
-    p = plot(plts...; size = (1000,600))
+    p = plot(plts...; size = (1400,800),bottommargin = 5Plots.mm)
     savefig(p,"$(fname)_posteriors.pdf")
     return out
 end
+
+#error function for inf
+#
+function visualize_π_base(particles_tuple)
+
+    param_keys = [:π_base_y,:π_base_m,:π_base_o]
+    # π_bases_array = map(f -> getproperty(particles_tuple,f),param_keys) |> l -> mapreduce(t-> [t...],hcat,zip(l...))
+    # display(π_bases_array)
+
+    # p = cornerplot(π_bases_array'; labels = string.(param_keys))
+    params = NamedTuple{(param_keys...,)}(particles_tuple)
+    # display(params)
+    p = corner(params)
+    display(p)
+
+
+    # display(scatter(map(f -> getproperty(particles_tuple,f),param_keys)...))
+end
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -43,8 +43,8 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
     @unpack base_transmission_probability,immunization_loss_prob,recovery_rate,immunizing,immunization_begin_day = modelsol.params
     @unpack u_inf,u_vac,u_next_inf,u_next_vac,demographics,inf_network,status_totals, immunization_countdown = modelsol
     
-    modelsol.daily_infected = 0
-    modelsol.daily_immunizations_by_age.= 0
+    modelsol.daily_cases_by_age .= 0
+    modelsol.daily_immunized_by_age .= 0
     function agent_transition!(node, from::AgentStatus,to::AgentStatus) 
         immunization_countdown[node] = -1 
         status_totals[Int(from)] -= 1
@@ -64,7 +64,7 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
                     for j in neighbors(mixing_graph.g,i)    
                         if u_inf[j] == Infected && u_next_inf[i] != Infected 
                             if rand(Random.default_rng(Threads.threadid())) < contact_weight(base_transmission_probability,get_weight(mixing_graph,GraphEdge(i,j)))
-                                modelsol.daily_infected+=1
+                                modelsol.daily_cases_by_age[Int(agent_demo)]+=1
                                 agent_transition!(i, Susceptible,Infected)
                             end
                         end
@@ -82,63 +82,17 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
         end
         
         if immunization_countdown[i] == 0 
-            modelsol.daily_immunizations_by_age[Int(agent_demo)] += 1
+            modelsol.daily_immunized_by_age[Int(agent_demo)] += 1
             agent_transition!(i, Susceptible,Immunized)
         elseif immunization_countdown[i]>0
             immunization_countdown[i] -= 1
         end
     end
 end
-# Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
-#     @unpack π_base, η,γ, κ, ω, ρ_y,ρ_m,ρ_o, ω_en,ρ_en,γ,β = modelsol.params
-#     @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
-#         vac_payoff = 0
-#         soc_nbrs_vac = @MArray [0,0,0]
-#         soc_nbrs_nonvac = 0
-#         num_soc_nbrs = 0
-#         for sc_g in soc_network.graph_list[t]
-#             soc_nbrs = neighbors(sc_g.g,i)
-#             num_soc_nbrs += length(soc_nbrs)
-#             for nbr in soc_nbrs
-#                 if u_vac[nbr]
-#                     soc_nbrs_vac[Int(demographics[nbr])] += 1
-#                 else
-#                     soc_nbrs_nonvac += 1
-#                 end
-#             end
-#         end
-#         ρ = @SVector [ρ_y,ρ_m,ρ_o]
-
-#         vac_payoff += π_base + dot(ρ,soc_nbrs_vac) + total_infections*ω +
-#             ifelse(num_soc_nbrs> 0, κ * ((sum(soc_nbrs_vac) - soc_nbrs_nonvac)/num_soc_nbrs),0)
-        
-#         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]]))* (η + dot(ρ_en,soc_nbrs_vac) + total_infections*ω_en) 
-#         end
-        
-#         if u_vac[i]
-#             if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,β)
-#                 u_next_vac[i] = !u_vac[i]
-#             else
-#                 u_next_vac[i] = u_vac[i]
-#             end
-#         else
-#             if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,β)
-#                 u_next_vac[i] = !u_vac[i]
-#             else
-#                 u_next_vac[i] = u_vac[i]
-#             end
-#         end
-#     end
-
-#     modelsol.daily_vaccinators = count(==(true),u_vac) #could maybe make this more efficient
-# end
 
 
 Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
-    @unpack π_base_y,π_base_m,π_base_o, η,γ, κ, ω, ω_en,ρ_en,γ,β = modelsol.params
+    @unpack π_base_y,π_base_m,π_base_o, η,γ, κ, ω, ω_en,γ,β = modelsol.params
     @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
@@ -152,7 +106,7 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
             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) 
+                    vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]])) * (η + total_infections*ω_en) 
                 end
         
                 if u_vac[i]