fix bug in mixing matrices generation, add weighted_degree helper function

CovidAlertVaccinationModel/Manifest.toml

@@ -1229,6 +1229,12 @@ git-tree-sha1 = "293aa2c5cbf201e6b98810cb36d9eeafdafdafd1"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 version = "0.19.34"
 
+[[VectorizedRNG]]
+deps = ["Distributed", "Random", "UnPack", "VectorizationBase"]
+git-tree-sha1 = "f60638a5cdd839eba16fd9658909513a8291f862"
+uuid = "33b4df10-0173-11e9-2a0c-851a7edac40e"
+version = "0.2.8"
+
 [[VersionParsing]]
 git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
 uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"

CovidAlertVaccinationModel/Project.toml

@@ -43,6 +43,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+VectorizedRNG = "33b4df10-0173-11e9-2a0c-851a7edac40e"
 ZeroWeightedDistributions = "24733ad3-391a-4e41-8839-c7177de7dea4"
 
 [compat]

CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl

@@ -138,11 +138,11 @@ List of lists of graphs, one list for each day.
     graph_list::Vector{Vector{G}}
 
 """
-struct TimeDepMixingGraph{N,G}
-    resampled_graphs::NTuple{N,G}
+struct TimeDepMixingGraph{N,G,GNT}
+    resampled_graphs::NTuple{N,GNT}
     graph_list::Vector{Vector{G}}
-    function TimeDepMixingGraph(len,resampled_graphs::NTuple{N,G},base_graph_list::Vector{G}) where {G,N}
-        return new{N,G}(
+    function TimeDepMixingGraph(len,resampled_graphs::NTuple{N,GNT},base_graph_list::Vector{G}) where {GNT,G,N}
+        return new{N,G,GNT}(
             resampled_graphs,
             [copy(base_graph_list) for i in 1:len]
         )
@@ -156,17 +156,20 @@ Assumes the simulation begins on Thursday arbitrarily.
 """
 function time_dep_mixing_graphs(len,base_network,demographics,index_vectors,ws_matrix_tuple,rest_matrix_tuple)
     home_static_edges = WeightedGraph(base_network,demographics,index_vectors,contact_time_distributions.hh) #network with households and LTC homes
+   
     ws_static_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.daily,contact_time_distributions.ws)
     ws_weekly_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.twice_a_week,contact_time_distributions.ws)
     ws_daily_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.otherwise,contact_time_distributions.ws)
+    
     rest_static_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.daily,contact_time_distributions.rest)
     rest_weekly_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.twice_a_week,contact_time_distributions.rest)
     rest_daily_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.otherwise,contact_time_distributions.rest)
+   
     inf_network_list = [home_static_edges,rest_static_edges] 
     soc_network_list = [home_static_edges,rest_static_edges,ws_static_edges]
     
-    infected_mixing_graph = TimeDepMixingGraph(len,(ws_daily_edges,rest_daily_edges),inf_network_list)
-    soc_mixing_graph = TimeDepMixingGraph(len,(ws_daily_edges,rest_daily_edges),soc_network_list)
+    infected_mixing_graph = TimeDepMixingGraph(len,((ws_daily_edges,ws_matrix_tuple.daily),(rest_daily_edges,rest_matrix_tuple.daily)),inf_network_list)
+    soc_mixing_graph = TimeDepMixingGraph(len,((ws_daily_edges,ws_matrix_tuple.daily),(rest_daily_edges,rest_matrix_tuple.daily)),soc_network_list)
     # display(infected_mixing_graph.graph_list)
     for (t,l) in enumerate(infected_mixing_graph.graph_list)
         day_of_week = mod(t,7)
@@ -187,8 +190,8 @@ end
 """
 Completely remake all the graphs in `time_dep_mixing_graph.resampled_graphs`.
 """
-function remake!(time_dep_mixing_graph,demographic_index_vectors,mixing_matrix)
-    for weighted_graph in time_dep_mixing_graph.resampled_graphs
+function remake!(time_dep_mixing_graph,demographic_index_vectors)
+    for (weighted_graph,mixing_matrix) in time_dep_mixing_graph.resampled_graphs
         empty!.(weighted_graph.g.fadjlist) #empty all the vector edgelists
         weighted_graph.g.ne = 0
         weighted_graph.mixing_edges = create_mixing_edges(demographic_index_vectors,mixing_matrix,weighted_graph.mixing_edges.sampler_matrix)

CovidAlertVaccinationModel/src/ABM/model_setup.jl

 
 function get_parameters()
     params = (
-        sim_length = 120,
+        sim_length = 500,
         num_households = 5000,
         I_0_fraction = 0.002,
         base_transmission_probability = 0.001,
@@ -21,8 +21,8 @@ function get_parameters()
         notification_threshold = 2,
         immunizing = true,
         immunization_delay = 14,
-        immunization_begin_day = 60,
-        infection_introduction_day = 120
+        immunization_begin_day =60, 
+        infection_introduction_day = 180,
     )
     return params
 end

CovidAlertVaccinationModel/src/ABM/output.jl

@@ -121,12 +121,12 @@ function plot_model(varname,univariate_series, output_list::Vector{T},infection_
     sim_length = length(output_list[1].recorded_status_totals.S)
     ts_list(data) = [
         (data.recorded_status_totals.S, "Susceptible over time"),
-        (data.recorded_status_totals.I, "Infected over time"),
         (data.recorded_status_totals.R, "Recovered over time"),
-        (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"),
-        (sum.(eachcol(mean.(data.daily_cases_by_age))),"Daily (incident) cases"),
+        (data.daily_cases_by_age.Y,"Daily (incident) Y cases"),
+        (data.daily_cases_by_age.M,"Daily (incident) M cases"),
+        (data.daily_cases_by_age.O,"Daily (incident) O 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"),

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -3,10 +3,11 @@ using KissABC
 const target_cumulative_vac_proportion = 0.33
 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 = []
+const ymo_attack_rate = [10.376,5.636,7.2]
 
 
 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))))
     out = DebugRecorder(0,default_p.sim_length)
     model = abm(new_params,out)
@@ -36,17 +37,17 @@ function fit_pre_inf_behavioural_parameters(p_tuple)
             immunizing = true,
         )
     )
-
     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[1:4]) .* target_cumulative_vaccinations
-        ymo_vaccination_ts = output.new_ymo_immunization
+        ymo_vaccination_ts = output.daily_immunized_by_age
         total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
         # display(target_ymo_vac)
         # display(total_preinfection_vaccination)
         return sum((total_preinfection_vaccination .- target_ymo_vac).^2)    
     end
+
     out = smc(priors,cost; verbose = true, nparticles = 800, parallel = true)#
     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
 end
@@ -68,22 +69,25 @@ function fit_post_inf_behavioural_parameters(p_tuple)
             immunizing = true,
         )
     )
-    target_cumulative_vac_proportion = 0.33
-    vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
-    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
+        ymo_vaccination_ts = output.daily_immunized_by_age
         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)    
+        final_size = sum.(eachrow(output.daily_cases_by_age))
+        target_final_size = ymo_attack_rate .* model.nodes
+
+        display(final_size)
+        display((sum((total_postinf_vaccination .- target_ymo_vac).^2) , 1e-3*sum((final_size .- target_final_size).^2)))
+        return sum((total_postinf_vaccination .- target_ymo_vac).^2)   + sum((final_size .- target_final_size).^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
+
+    display(cost([0.0001,0.0005]))
+    # 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,))
 end
 
@@ -96,6 +100,7 @@ function plot_behavioural_fit(particles,p_tuple)
             sim_length = sim_length,
             I_0_fraction = 0.000,
             immunization_begin_day =60, 
+            infection_introduction_day = 180,
             immunizing = true,
         )
     )
@@ -106,7 +111,7 @@ function plot_behavioural_fit(particles,p_tuple)
     target_cumulative_vac_proportion = 0.33
     vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
     ymo_vac = @SVector [0.255,0.278,0.602]
-    ymo_vaccination_ts = mean.(out.new_ymo_immunization)
+    ymo_vaccination_ts = mean.(out.daily_immunized_by_age)
     total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
     display(total_preinfection_vaccination)
     p = [plot(),plot(),plot()]
@@ -120,34 +125,25 @@ end
 # plot_max_posterior("outbreak", outbreak_transmission_dist,default_parameters)
 
 function fit_parameters(default_parameters)
-    seasonal_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","outbreak_inf_parameters.dat")
     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)
 
-    # serialize(seasonal_parameters_path,seasonal_transmission_dist)
+    # pre_inf_behaviour_parameters = fit_pre_inf_behavioural_parameters(default_parameters) 
     # 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)
-    
+
+    pre_inf_behaviour_parameters = deserialize(pre_inf_behaviour_parameters_path)
+    post_inf_behaviour_parameters = fit_post_inf_behavioural_parameters(merge(default_parameters,map(mode,pre_inf_behaviour_parameters))) 
+    serialize(post_inf_behaviour_parameters_path, post_inf_behaviour_parameters)
     
     # visualize_π_base(deserialize(pre_inf_behaviour_parameters_path))
-    
-
     fitted_parameters_with_post_inf_behaviour = (;
-        deserialize(seasonal_parameters_path)...,
         deserialize(pre_inf_behaviour_parameters_path)...,
         deserialize(post_inf_behaviour_parameters_path)...
     )
     display(map(mode,fitted_parameters_with_post_inf_behaviour))
 
-
     fitted_sol = plot_fitting_posteriors("post_inf_fitting",fitted_parameters_with_post_inf_behaviour,default_parameters)
     return fitted_sol
 end 

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -122,14 +122,33 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
 
 end
 
+function weighted_degree(node,network::TimeDepMixingGraph)
+    weighted_degree = 0
+    for g_list in network.graph_list
+        for g in g_list
+            for j in neighbors(g.g,node)
+                weighted_degree += get_weight(g,GraphEdge(node,j))
+            end
+        end
+    end
+    return weighted_degree
+end
 
 function agents_step!(t,modelsol)
 
-    remake!(modelsol.inf_network,modelsol.index_vectors,modelsol.ws_matrix_tuple.daily)
-    remake!(modelsol.soc_network,modelsol.index_vectors,modelsol.rest_matrix_tuple.daily)
-    for network in modelsol.inf_network.graph_list[t]
+    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
+    # @show modelsol.inf_network.graph_list[t]
+    # @show weighted_degree(1,modelsol.inf_network)
+    # @show weighted_degree(1,modelsol.soc_network)
+
+
+
+
+
+
     if t == modelsol.params.infection_introduction_day
         init_indices = rand(Random.default_rng(Threads.threadid()), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
         modelsol.u_inf[init_indices] .= Infected