spread initial cases out over 36 days

CovidAlertVaccinationModel/src/ABM/model_setup.jl

 function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35,0.35,0.2)
     params = (
-        sim_length = 400,
+        sim_length = 600,
         num_households = 5000,
-        I_0_fraction = 0.005,
-        β_y = 0.00075,
-        β_m = 0.00063,
+        I_0_fraction = 0.003,
+        β_y = 0.00077,
+        β_m = 0.00065,
         β_o = 0.75,
         α_y = 0.4,
         α_m = 0.4,
@@ -27,7 +27,7 @@ function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35
         immunization_delay = 14,
         immunization_begin_day = 60, 
         infection_introduction_day = 180,
-        ζ = 1.25
+        ζ = 1.3
     )
     return params
 end

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -66,17 +66,31 @@ 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_distribution_parameters(default_parameters)
-    serialize(fit_all_parameters_path,output)
+    # output = fit_distribution_parameters(default_parameters)
+    # serialize(fit_all_parameters_path,output)
 
     
     fitted_parameter_tuple = deserialize(fit_all_parameters_path)
     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)))
+    ymo_vaccination_ts = mean.(fitted_sol.daily_immunized_by_age)
+
+    total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
     final_size = sum.(eachrow(mean.(fitted_sol.daily_unvac_cases_by_age)))
-    display(final_vac./avg_populations)
-    display(final_size./avg_populations)
+    total_preinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,1:180]))
+    target_final_size = ymo_attack_rate .*avg_populations
+    target_preinf_vac = ymo_vac .* sum(vaccination_data[1:4]) .* avg_populations
+    target_postinf_vac = ymo_vac .* sum(vaccination_data[5:end]) .*avg_populations
+
+  
+    println("obs final size: $final_size, target: $target_final_size")
+    println("obs preinf vac: $total_preinf_vaccination, target: $target_preinf_vac")
+    println("obs postinf vac: $total_postinf_vaccination,target: $target_postinf_vac")
+
+    total_final_size = sum.(eachrow(mean.(fitted_sol.daily_cases_by_age))) 
+    println("vac + unvac cases proportion: $(total_final_size./avg_populations))")
+    display(sum.(eachrow(ymo_vaccination_ts)) ./avg_populations) 
+
     return fitted_sol
 end 
 
@@ -88,6 +102,7 @@ function plot_fitting_posteriors(fname,particles_tuple,parameters)
     samples = collect(zip(particles_tuple...))
     for p_set in sample(samples,30)
         p_set_named = NamedTuple{names}(p_set)
+        display(p_set_named)
         sol = abm(merge(parameters,p_set_named),output_recorder)
         avg_populations .+= length.(sol.index_vectors)
         fit!(stat_recorder,output_recorder)

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -150,14 +150,14 @@ function weighted_degree(node,network::TimeDepMixingGraph)
     return weighted_degree
 end
 
-function agents_step!(t,modelsol)
+function agents_step!(t,modelsol, init_indices)
 
-    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
-        modelsol.status_totals[Int(Infected)] += length(init_indices)
-    end
     if t>modelsol.params.infection_introduction_day
+        if !isempty(init_indices)
+            inf_index = pop!(init_indices)
+            modelsol.u_inf[inf_index] = Infected
+            modelsol.status_totals[Int(Infected)] += 1
+        end
         update_alert_durations!(t,modelsol)
     end
     update_vaccination_opinion_state!(t,modelsol,modelsol.status_totals[Int(Infected)])
@@ -175,9 +175,10 @@ end
 
 
 function solve!(modelsol,recordings...)
-
+    init_indices = rand(Random.default_rng(Threads.threadid()), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
+    @show length(init_indices)
     for t in 1:modelsol.sim_length
-        agents_step!(t,modelsol) 
+        agents_step!(t,modelsol,init_indices) 
         #advance agent states based on the new network
         for recording in recordings
             record!(t,modelsol,recording)