name pi_base drop parameter zeta

CovidAlertVaccinationModel/src/ABM/model_setup.jl

@@ -3,19 +3,19 @@ function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35
         sim_length = 400,
         num_households = 5000,
         I_0_fraction = 0.005,
-        β_y = 0.00047,
-        β_m = 0.00041,
+        β_y = 0.00075,
+        β_m = 0.00063,
         β_o = 0.75,
         α_y = 0.4,
         α_m = 0.4,
         α_o = 0.4,
         recovery_rate = 1/5,
-        π_base_y = -1.30,
-        π_base_m = -1.31,
+        π_base_y = -1.4,
+        π_base_m = -1.4,
         π_base_o = -0.95,
         η = 0.0,
         κ = 0.0,
-        ω = 0.0005,
+        ω = 0.005,
         ω_en = 0.00,
         γ = 0.0,
         ξ = 5.0,
@@ -27,6 +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
     )
     return params
 end

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -15,7 +15,7 @@ function solve_w_parameters(default_p, p_names, new_p_list)
 end
 
 function fit_distribution_parameters(p_tuple)
-    p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o)
+    p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o)
     priors = Factored(
         Uniform(0.0,0.01),
         Uniform(0.0,0.5),  
@@ -23,15 +23,12 @@ function fit_distribution_parameters(p_tuple)
         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),   
+        Uniform(-5.0,3.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
+    sim_length = 400
     p_tuple_adjust = merge(p_tuple,
         (
             sim_length = sim_length,
@@ -58,21 +55,21 @@ function fit_distribution_parameters(p_tuple)
         # display((total_postinf_vaccination,target_postinf_vac))
         # display(sum.(eachrow(ymo_vaccination_ts)) ./length.(model.index_vectors))
 
-        # return sum((total_preinf_vaccination .- target_preinf_vac).^2)  
-        # + sum((total_postinf_vaccination .- total_postinf_vaccination).^2)  
-        return sum((final_size .- target_final_size).^2)
+        return sum((total_preinf_vaccination .- target_preinf_vac).^2)  
+        + sum((total_postinf_vaccination .- total_postinf_vaccination).^2)  
+        + sum((final_size .- target_final_size).^2)
     end 
-    #    display(cost((0.0000,0.00051,0.00046,0.18,-1.30,-1.24,-0.8,0.35,0.35,0.35,0.2)))
-    #
-    out = smc(priors,cost; verbose = true, nparticles = 600, parallel = true)
+    #    display(cost((0.005,0.00075, 0.00063,0.75,-1.4,-1.4,-0.95)))
+    
+    out = smc(priors,cost; verbose = true, nparticles = 1000, parallel = true)
     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
 end
 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)
@@ -91,13 +88,13 @@ function plot_fitting_posteriors(fname,particles_tuple,parameters)
     avg_populations = [0.0,0.0,0.0]
     names = keys(particles_tuple)
     samples = collect(zip(particles_tuple...))
-    for p_set in sample(samples,20)
+    for p_set in sample(samples,1)
         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
-    avg_populations ./= 20
+    avg_populations ./= 1
     p = plot_model(nothing,[nothing],[stat_recorder],parameters.infection_introduction_day,parameters.immunization_begin_day)
     savefig(p, "$fname.pdf")
     

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -98,13 +98,13 @@ end
 
 
 Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
-    @unpack infection_introduction_day, π_base_y,π_base_m,π_base_o, η,γ, κ, ω, ω_en,γ,ξ = modelsol.params
+    @unpack infection_introduction_day, π_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
         π_base = t<infection_introduction_day ? 
         (π_base_y,π_base_m,π_base_o) : 
-        (π_base_y*1.25,π_base_m*1.25,π_base_o)
+        (π_base_y*ζ,π_base_m*ζ,π_base_o*ζ)
         random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
 
         if !isempty(neighbors(random_soc_network,i))

CovidAlertVaccinationModel/test/ABM/mixing_test.jl

@@ -15,24 +15,24 @@ const model_sizes = [1000,5000]
 
 #WRITE MORE TESTS
 
-# @testset "mixing matrices, size: $sz" for sz in model_sizes
-#     for rep = 1:samples
-#         m = ModelSolution(100,get_parameters(),sz)
-#         ws_dist = m.ws_matrix_tuple
-#         r_dist = m.rest_matrix_tuple 
-#         index_vec =m.index_vectors
-#         @testset "workschool" for i in dem_cat, j in dem_cat
-#             for t in 1:length(ws_dist)
-#                 @test mean(ws_dist[t][i,j])*length(index_vec[i]) == mean(ws_dist[t][j,i])*length(index_vec[j])
-#             end
-#         end
-#         @testset "rest" for i in dem_cat, j in dem_cat
-#             for t in 1:length(ws_dist)
-#                 @test mean(r_dist[t][i,j])*length(index_vec[i]) == mean(r_dist[t][j,i])*length(index_vec[j])
-#             end
-#         end
-#     end
-# end
+@testset "mixing matrices, size: $sz" for sz in model_sizes
+    for rep = 1:samples
+        m = ModelSolution(100,get_parameters(),sz)
+        ws_dist = m.ws_matrix_tuple
+        r_dist = m.rest_matrix_tuple 
+        index_vec =m.index_vectors
+        @testset "workschool" for i in dem_cat, j in dem_cat
+            for t in 1:length(ws_dist)
+                @test mean(ws_dist[t][i,j])*length(index_vec[i]) == mean(ws_dist[t][j,i])*length(index_vec[j])
+            end
+        end
+        @testset "rest" for i in dem_cat, j in dem_cat
+            for t in 1:length(ws_dist)
+                @test mean(r_dist[t][i,j])*length(index_vec[i]) == mean(r_dist[t][j,i])*length(index_vec[j])
+            end
+        end
+    end
+end
 
 function approximate_mixing_matricies_weights_dict(p,sol)
     mean_mixing_degree = zeros(3,3)
@@ -73,46 +73,47 @@ function approximate_mixing_matricies_graph(p,sol)
         mean_mixing_weighted_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
     return mean_mixing_degree, mean_mixing_weighted_degree
 end
-# @testset "weights dict and graph match" for _ in samples
-#     p = get_parameters()
-#     sol = abm(p,nothing)
-#     graph_deg_matrix,graph_wdeg_matrix = approximate_mixing_matricies_graph(p,sol)
-#     dict_deg_matrix,dict_wdeg_matrix = approximate_mixing_matricies_weights_dict(p,sol)
+@testset "weights dict and graph match" for _ in samples
+    p = get_parameters()
+    sol = abm(p,nothing)
+    graph_deg_matrix,graph_wdeg_matrix = approximate_mixing_matricies_graph(p,sol)
+    dict_deg_matrix,dict_wdeg_matrix = approximate_mixing_matricies_weights_dict(p,sol)
 
-#     @test all(graph_deg_matrix .≈ dict_deg_matrix)
-#     @test all(graph_wdeg_matrix .≈ dict_wdeg_matrix)
-# end
+    @test all(graph_deg_matrix .≈ dict_deg_matrix)
+    @test all(graph_wdeg_matrix .≈ dict_wdeg_matrix)
+end
 
 using OnlineStats
 
 
-# @testset "contact vector generation" for _ in samples
-#     model = ModelSolution(1, get_parameters(), 1_000_000) #big
-#     @unpack demographics, index_vectors, rest_matrix_tuple,ws_matrix_tuple = model
-#     ThreadsX.map((ws_matrix_tuple...,rest_matrix_tuple...)) do dist
-#         for i in 1:3, j in 1:i  #diagonal
-#             if i != j 
-#                 num_degrees_ij = similar(index_vectors[i])
-#                 num_degrees_ji = similar(index_vectors[j])
-#                 generate_contact_vectors!(dist[i,j],dist[j,i],num_degrees_ij,num_degrees_ji) 
-#                 @test sum(num_degrees_ji) == sum(num_degrees_ij)
-#                 @test mean(num_degrees_ij) ≈ mean(dist[i,j]) atol = 1e-2
-#                 @test mean(num_degrees_ji) ≈ mean(dist[j,i]) atol = 1e-2
-#                 @test var(num_degrees_ij) ≈ var(dist[i,j]) atol = 1e-1
-#                 @test var(num_degrees_ji) ≈ var(dist[j,i]) atol = 1e-1
-#                 @test skewness(num_degrees_ij) ≈ skewness(dist[i,j]) atol = 5e-1
-#                 @test skewness(num_degrees_ji) ≈ skewness(dist[j,i]) atol = 5e-1
-#             else
-#                 num_degrees_ii = similar(index_vectors[i])
-#                 generate_contact_vectors!(dist[i,j],num_degrees_ii) 
-#                 @test iseven(sum(num_degrees_ii))
-#                 @test mean(num_degrees_ii) ≈ mean(dist[i,j]) atol = 1e-2
-#                 @test var(num_degrees_ii) ≈ var(dist[i,j]) atol = 1e-1
-#                 @test skewness(num_degrees_ii) ≈ skewness(dist[i,j]) atol = 5e-1
-#             end
-#         end
-#     end
-# end
+@testset "contact vector generation" for _ in samples
+    model = ModelSolution(1, get_parameters(), 1_000_000) #big
+    @unpack demographics, index_vectors, rest_matrix_tuple,ws_matrix_tuple = model
+    
+    for dist in (ws_matrix_tuple...,rest_matrix_tuple...)
+        for i in 1:3, j in 1:i  #diagonal
+            if i != j 
+                num_degrees_ij = similar(index_vectors[i])
+                num_degrees_ji = similar(index_vectors[j])
+                generate_contact_vectors!(dist[i,j],dist[j,i],num_degrees_ij,num_degrees_ji) 
+                @test sum(num_degrees_ji) == sum(num_degrees_ij)
+                @test mean(num_degrees_ij) ≈ mean(dist[i,j]) atol = 1e-2
+                @test mean(num_degrees_ji) ≈ mean(dist[j,i]) atol = 1e-2
+                @test var(num_degrees_ij) ≈ var(dist[i,j]) atol = 1e-1
+                @test var(num_degrees_ji) ≈ var(dist[j,i]) atol = 1e-1
+                @test skewness(num_degrees_ij) ≈ skewness(dist[i,j]) atol = 5e-1
+                @test skewness(num_degrees_ji) ≈ skewness(dist[j,i]) atol = 5e-1
+            else
+                num_degrees_ii = similar(index_vectors[i])
+                generate_contact_vectors!(dist[i,j],num_degrees_ii) 
+                @test iseven(sum(num_degrees_ii))
+                @test mean(num_degrees_ii) ≈ mean(dist[i,j]) atol = 1e-2
+                @test var(num_degrees_ii) ≈ var(dist[i,j]) atol = 1e-1
+                @test skewness(num_degrees_ii) ≈ skewness(dist[i,j]) atol = 5e-1
+            end
+        end
+    end
+end
 
 @testset "degree distribution generation"  for _ in samples
     model = ModelSolution(1, get_parameters(), 100_000)