new sims, still need to look through them. Go back to sampling households with weights.

CovidAlertVaccinationModel/prem_mixing_approx.jl

+using CovidAlertVaccinationModel:ModelSolution, get_parameters, neighbors, GraphEdge, get_weight
+function approximate_mixing_matricies()
+    sol = abm(get_parameters(),nothing)
+    mean_mixing = zeros(3,3)
+    display(mean_mixing)
+    display(sol.inf_network)
+    for (node_i,demo_i) in enumerate(sol.demographics)
+        for g_list in sol.inf_network.graph_list
+            for g in g_list
+                # display(g.mixing_edges.weights_dict)
+                for node_j in neighbors(g,node_i)
+                    demo_j = sol.demographics[node_j];
+                    mean_mixing[Int(demo_i), Int(demo_j)] += 1#get_weight(g,GraphEdge(node_i,node_j)) /2
+                end 
+             end
+        end
+    end
+    mean_mixing ./= (500 .* length.(sol.index_vectors) * 2)
+    display(mean_mixing)
+    return sol
+end
+approximate_mixing_matricies()
+#rerun with weights
+#
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/model_setup.jl

-
 function get_parameters()
     params = (
         sim_length = 500,

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -12,79 +12,6 @@ function solve_w_parameters(default_p, p_names, new_p_list)
     model = abm(new_params,out)
     return out, model
 end
-#######
-# function fit_pre_inf_behavioural_parameters(p_tuple)
-#     samples = 1
-#     p_names = (:π_base_y,:π_base_m,:π_base_o)
-#     priors = Factored(
-#         Uniform(-10.0,2.0),
-#         Uniform(-10.0,2.0),
-#         Uniform(-10.0,2.0)
-#     )
-    
-#     #simulation begins in july
-#     #60 days for opinion dynamics to stabilize, then immunization begins in september,
-#     #infection is not considered
-    
-#     sim_length = 180
-#     p_tuple_adjust = merge(p_tuple,
-#         (
-#             sim_length = sim_length,
-#             I_0_fraction = 0.000,
-#             immunization_begin_day =60, 
-#             infection_introduction_day = 180,
-#             immunizing = true,
-#         )
-#     )
-#     function cost(p)
-#         output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
-#         target_ymo_vac = ymo_vac .* sum(vaccination_data[1:4]) .* length.(model.index_vectors)
-#         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 = 150, parallel = true)#
-#     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
-# end
-
-# function fit_post_inf_behavioural_parameters(p_tuple)
-#     p_names = (:ω,:β_y,:β_m,:β_o)
-#     priors = Factored(Uniform(0.0,0.1),Uniform(0.0,0.1),Uniform(0.0,0.1),Uniform(0.0,1.0))
-#     #simulation begins in july
-#     #60 days for opinion dynamics to stabilize, then immunization begins in september,
-#     #infection introduced at beginning of december
-#     sim_length = 300
-#     p_tuple_adjust = merge(p_tuple,
-#         (
-#             sim_length = sim_length,
-#             I_0_fraction = 0.005,
-#             immunization_begin_day =60, 
-#             infection_introduction_day = 180,
-#             immunizing = true,
-#         )
-#     )
-#     function cost(p)
-#         output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
-#         target_ymo_vac = ymo_vac .* sum(vaccination_data[1:end]) .* length.(model.index_vectors)
-#         ymo_vaccination_ts = output.daily_immunized_by_age
-#         total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
-
-#         final_size = sum.(eachrow(output.daily_cases_by_age))
-#         target_final_size = ymo_attack_rate .* length.(model.index_vectors)
-#         # display( length.(model.index_vectors))
-#         # display((final_size,target_final_size))
-#         # display((total_postinf_vaccination,target_ymo_vac))
-#         # display((1e-1*sum((total_postinf_vaccination .- target_ymo_vac).^2) , sum((final_size .- target_final_size).^2)))
-#         return 1e-2*sum((total_postinf_vaccination .- target_ymo_vac).^2)   + sum((final_size .- target_final_size).^2)
-#     end
-
-#     # display(cost([0.000,0.001,0.001,1.0]))
-#     out =smc(priors,cost; verbose = true, nparticles = 400, parallel = true)# ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=1000,  parallel = true) #this one has better NaN handling
-#     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
-# end
 
 function fit_all_parameters(p_tuple)
     p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o)
@@ -92,13 +19,13 @@ function fit_all_parameters(p_tuple)
         Uniform(0.0,0.01),
         Uniform(0.0,0.2),
         Uniform(0.0,0.1),
-        Uniform(0.0,1.0),
-        Uniform(-5.0,0.0),
-        Uniform(-5.0,0.0),
-        Uniform(-5.0,5.0),
-        Uniform(0.0,1.0),
-        Uniform(0.0,1.0),
-        Uniform(0.0,1.0),
+        Uniform(0.0,1.0), 
+        Uniform(-5.0,3.0),
+        Uniform(-5.0,3.0),
+        Uniform(-5.0,3.0),
+        TriangularDist(23.0, 44.0, 34.0),
+        TriangularDist(22.0, 54.0, 40.0),
+        TriangularDist(9.0, 59.0, 39.0)   
     )
     #simulation begins in july
     #60 days for opinion dynamics to stabilize, then immunization begins in september,
@@ -108,7 +35,7 @@ function fit_all_parameters(p_tuple)
         (
             sim_length = sim_length,
             I_0_fraction = 0.005,
-            immunization_begin_day =60, 
+            immunization_begin_day = 60, 
             infection_introduction_day = 180,
             immunizing = true,
         )
@@ -130,54 +57,19 @@ function fit_all_parameters(p_tuple)
         total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
         return sum((normalize_by_pop(total_preinfection_vaccination .- target_ymo_vac)).^2)  
         + sum(normalize_by_pop((total_postinf_vaccination .- target_ymo_vac)).^2)  
-        + 2*sum(normalize_by_pop((final_size .- target_final_size)).^2)
+        + sum(normalize_by_pop((final_size .- target_final_size)).^2)
     end
-
-    # display(cost(rand(priors)))
     out = smc(priors,cost; verbose = true, nparticles = 1000, parallel = true)#this one has better NaN handling
     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
 end
 
-# function plot_behavioural_fit(particles,p_tuple)
-#     p_names = (:π_base_y,:π_base_m,:π_base_o)
-#     sim_length = 210
-#     samples = 1
-#     p_tuple_adjust = merge(p_tuple,
-#         (
-#             sim_length = sim_length,
-#             I_0_fraction = 0.000,
-#             immunization_begin_day =60, 
-#             infection_introduction_day = 180,
-#             immunizing = true,
-#         )
-#     )
-#     p = map(e -> mode(e.particles),particles.P)
-#     display(p)
-#     new_params = merge(p_tuple_adjust, NamedTuple{p_names}(ntuple(i -> p[i],length(p_names))))
-#     out = mean_solve(samples, new_params ,DebugRecorder)
-#     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.daily_immunized_by_age)
-#     total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
-#     display(total_preinfection_vaccination)
-#     p = [plot(),plot(),plot()]
-
-#     p = plot_model(nothing,[nothing],[out],new_params.infection_introduction_day,new_params.immunization_begin_day)
-#     savefig(p,"behaviour_fit.pdf")
-#     return out
-# end
-# outbreak_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(default_parameters,0.241) ##outbreak
-# serialize(joinpath(PACKAGE_FOLDER,"abm_parameter_fits","outbreak_inf_parameters.dat"),outbreak_transmission_dist)
-# plot_max_posterior("outbreak", outbreak_transmission_dist,default_parameters)
-
 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")
     fit_all_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","fit_all_parameters.dat")
 
-    # output = fit_all_parameters(default_parameters)
-    # serialize(fit_all_parameters_path,output)
+    output = fit_all_parameters(default_parameters)
+    serialize(fit_all_parameters_path,output)
 
     
     fitted_parameter_tuple = deserialize(fit_all_parameters_path)
@@ -190,33 +82,21 @@ function fit_parameters(default_parameters)
     return fitted_sol
 end 
 
-# function plot_max_posterior(fname,particles,parameters)
-#     samples = 5
-#     base_transmission = mode(particles.base_transmission_probability)
-#     p_tuple_without_vac = merge(parameters,
-#         (
-#             sim_length = 150,
-#             immunization_begin_day = 0,
-#             infection_introduction_day = 1,
-#             immunizing = false,
-#         )
-#     )
-    
-#     new_params = merge(p_tuple_without_vac, (base_transmission_probability = base_transmission,))
-#     out,_ = mean_solve(samples, new_params ,DebugRecorder)
-#     p = plot_model(nothing,[nothing],[out],new_params.infection_introduction_day,new_params.immunization_begin_day)
-#     savefig(p,"$fname.pdf")
-#     hist = StatsBase.fit(Histogram,particles.base_transmission_probability; nbins = 25)
-#     p = plot(hist;legend = false)
-#     savefig(p,"$(fname)_posterior.pdf")
-
-# end
-using PairPlots
 function plot_fitting_posteriors(fname,particles_tuple,parameters)
-    p = merge(parameters,map(mean,particles_tuple))
     # p_adjust = merge(parameters,(β_y = p.β_y*0.1, β_m = p.β_m*0.1 ,β_o = p.β_o)) 
-    out,avg_populations = mean_solve(5, p,DebugRecorder)
-    p = plot_model(nothing,[nothing],[out],parameters.infection_introduction_day,parameters.immunization_begin_day)
+    stat_recorder = DebugRecorder(Variance(), parameters.sim_length)
+    output_recorder = DebugRecorder(parameters.sim_length)
+    avg_populations = [0.0,0.0,0.0]
+    names = keys(particles_tuple)
+    samples = collect(zip(particles_tuple...))
+    for p_set in sample(samples,10)
+        p_set_named = NamedTuple{names}(p_set)
+        sol = abm(merge(parameters,p_set_named),output_recorder)
+        avg_populations .+= length.(sol.index_vectors)
+        fit!(stat_recorder,output_recorder)
+    end
+
+    p = plot_model(nothing,[nothing],[stat_recorder],parameters.infection_introduction_day,parameters.immunization_begin_day)
     savefig(p, "$fname.pdf")
     
     plts = [plot() for i in 1:length(particles_tuple)]
@@ -226,21 +106,5 @@ function plot_fitting_posteriors(fname,particles_tuple,parameters)
     end
     p = plot(plts...; size = (1400,800),bottommargin = 5Plots.mm)
     savefig(p,"$(fname)_posteriors.pdf")
-    return out,avg_populations
+    return stat_recorder,avg_populations
 end
-
-# 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/data.jl

@@ -44,11 +44,11 @@ function read_household_data()
 end
 
 function sample_household_data(n)
-    return sample(Random.default_rng(Threads.threadid()),household_data.households, n)
+    return sample(Random.default_rng(Threads.threadid()),household_data.households,Weights(household_data.weight_vector), n)
 end
 
 function get_household_data_proportions()
-    households_by_demographic_sum = sum.([map(l -> l[i],household_data.households) for i in 1:3])
+    households_by_demographic_sum = sum.([map(((l,w),)-> l[i]*w,zip(household_data.households,household_data.weight_vector)) for i in 1:3])
     return households_by_demographic_sum./sum(households_by_demographic_sum)
     # https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en
 end