new simulations

CovidAlertVaccinationModel/src/ABM/output.jl

@@ -93,11 +93,14 @@ end
 function mean_solve(samples,parameter_tuple,recorder)
     stat_recorder = recorder(Variance(), parameter_tuple.sim_length)
     output_recorder = recorder(parameter_tuple.sim_length)
-    for i in 1:samples
-        global sol = abm(parameter_tuple,output_recorder)
+    avg_populations = [0.0,0.0,0.0]
+    for _ in 1:samples
+        sol = abm(parameter_tuple,output_recorder)
+        avg_populations .+= length.(sol.index_vectors)
         fit!(stat_recorder,output_recorder)
     end
-    return stat_recorder
+    avg_populations ./= samples
+    return stat_recorder,avg_populations
 end
 
 function OnlineStats.fit!(stat_type::R,pt) where {T,OS<:OnlineStat{T}, R<:AbstractRecorder{OS}}

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -12,79 +12,79 @@ 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)
-    )
+#######
+# 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
+#     #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
+#     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)
@@ -117,6 +117,8 @@ function fit_all_parameters(p_tuple)
         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
+
+        normalize_by_pop(v) = v./length.(model.index_vectors)
         total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
 
         final_size = sum.(eachrow(output.daily_unvac_cases_by_age))
@@ -126,112 +128,94 @@ function fit_all_parameters(p_tuple)
 
         ymo_vaccination_ts = output.daily_immunized_by_age
         total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
-        return sum((total_preinfection_vaccination .- target_ymo_vac).^2)  + sum((total_postinf_vaccination .- target_ymo_vac).^2)   + sum((final_size .- target_final_size).^2)
+        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)
     end
 
-    # display(cost([0.000,0.001,0.001,1.0]))
-    out =smc(priors,cost; verbose = true, nparticles = 1000, parallel = true)# ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=1000,  parallel = true) #this one has better NaN handling
+    # 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
+# 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")
+    # 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")
 
-
-    # pre_inf_behaviour_parameters = fit_pre_inf_behavioural_parameters(default_parameters) 
-    # serialize(pre_inf_behaviour_parameters_path, pre_inf_behaviour_parameters)
-
-
-    # pre_inf_behaviour_parameters = deserialize(pre_inf_behaviour_parameters_path)
-    # display(map(mode,pre_inf_behaviour_parameters))
-    # 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)
-    
-
-
-    # fitted_parameter_tuple = (;
-    #     deserialize(pre_inf_behaviour_parameters_path)...,
-    #     deserialize(post_inf_behaviour_parameters_path)...
-    # )
-    # display(map(mode,fitted_parameters_with_post_inf_behaviour))
-
     # output = fit_all_parameters(default_parameters)
     # serialize(fit_all_parameters_path,output)
 
     
     fitted_parameter_tuple = deserialize(fit_all_parameters_path)
-    fitted_sol = plot_fitting_posteriors("post_inf_fitting",fitted_parameter_tuple,default_parameters)
+    fitted_sol,avg_populations = plot_fitting_posteriors("post_inf_fitting",fitted_parameter_tuple,default_parameters)
+
+    final_vac = sum.(eachrow(mean.(fitted_sol.daily_immunized_by_age)))
+    final_size = sum.(eachrow(mean.(fitted_sol.daily_unvac_cases_by_age)))
+    display(final_vac./avg_populations)
+    display(final_size./avg_populations)
     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,
-        )
-    )
+# 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
+#     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_tuple_adjust = merge(parameters,
-        (
-            sim_length = 500,
-            I_0_fraction = 0.002,
-            immunization_begin_day =60, 
-            infection_introduction_day = 180,
-            immunizing = true,
-        )
-    )
-    out = mean_solve(5, merge(p_tuple_adjust,map(mode,particles_tuple)) ,DebugRecorder)
+    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)
     savefig(p, "$fname.pdf")
     
@@ -242,7 +226,7 @@ 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
+    return out,avg_populations
 end
 
 # function visualize_π_base(particles_tuple)