add mixing check to tests, some new graphs

CovidAlertVaccinationModel/abm_sensitivity_analysis.jl

 using CovidAlertVaccinationModel
+using CovidAlertVaccinationModel:get_parameters
 using OnlineStats
 using ThreadsX
 using Plots
-const samples = 5
+const samples = 10
 ##Univariate tests
-const len = 4 #number of points to evaluate
+const len = 10 #number of points to evaluate
 gr()
 
-const default_parameters = (
-    sim_length = 600,
-    num_households = 5000,
-    I_0_fraction = 0.005,
-    base_transmission_probability = 0.001,
-    recovery_rate = 1/7,
-    immunization_loss_prob = 0.0055, #mean time of 6 months
-    π_base = -4.0,
-    η = 0.0,
-    κ = 1.0,
-    ω = 0.001,
-    ρ = [0.0,0.0,0.0],
-    ω_en = 0.00,
-    ρ_en =  [0.0,0.0,0.0],
-    γ = 0.0,
-    β = 5.0,
-    notification_parameter = 0.001,
-    vaccinator_prob = 0.2,
-    app_user_fraction = 0.4,
-    notification_threshold = 2,
-    immunizing = true,
-    immunization_delay = 14,
-    immunization_begin_delay = 50,
-    infection_introduction_delay = 100
-)
-
 #add beta
 #run model without vaccination
 #use derivative of log of no infections to calibrate p
@@ -42,31 +17,32 @@ const default_parameters = (
 #
 
 const univarate_test_list = (
-    (:I_0_fraction, range(0.0, 0.05; length = len)), 
-    (:base_transmission_probability, range(0.0002, 0.002; length = len)),
-    (:recovery_rate, range(0.1, 0.5; length = len)),
-    (:immunization_loss_prob, range(0.00, 0.05; length = len)),
-    (:π_base, range(-4.5, -3.5;  length = len)),
-    (:η, range(0.0, 0.01; length = len)),
-    (:κ, range(0.5, 1.5; length = len)),
-    (:ω, range(0.0, 0.01; length = len)),
-    (:ω_en, range(0.0, 0.5; length = len)),
-    (:γ, range(0.0, 0.5; length = len)),
-    (:β, range(2, 30; length = len)),
-    (:notification_parameter, range(0.00, 0.05; length = len)),
-    (:app_user_fraction, range(0.05, 0.25; length = len)),
-    (:notification_threshold, (1:len)),
-    (:immunization_delay, [7,10,14,20]),
+    # (:I_0_fraction, range(0.0, 0.05; length = len)), 
+    # (:base_transmission_probability, range(0.0002, 0.002; length = len)),
+    # (:recovery_rate, range(0.1, 0.5; length = len)),
+    # (:immunization_loss_prob, range(0.00, 0.05; length = len)),
+    # (:π_base, range(-4.5, -3.5;  length = len)),
+    # (:η, range(0.0, 0.01; length = len)),
+    # (:κ, range(0.5, 1.5; length = len)),
+    # (:ω, range(0.0, 0.01; length = len)),
+    # (:ω_en, range(0.0, 0.5; length = len)),
+    # (:γ, range(0.0, 0.5; length = len)),
+    (:ξ, range(1, 15; length = len)),
+    # (:notification_parameter, range(0.00, 0.05; length = len)),
+    # (:app_user_fraction, range(0.05, 0.25; length = len)),
+    # (:notification_threshold, (1:len)),
+    # (:immunization_delay, [7,10,14,20]),
 )
 
 const univariate_path = "CovidAlertVaccinationModel/plots/univariate/" 
 function univarate_test(variable, variable_range)
+    default_parameters = get_parameters() 
     parameter_range_list = [merge(default_parameters,NamedTuple{(variable,)}((value,))) for value in variable_range]
-    solve_fn(p) = mean_solve(samples, p,DebugRecorder)
+    solve_fn(p) = mean_solve(samples, p,DebugRecorder)[1]
 
     univariate_outlist = ThreadsX.map(solve_fn, parameter_range_list)
     
-    p = plot_model(variable,parameter_range_list,univariate_outlist,default_parameters.infection_introduction_delay,default_parameters.immunization_begin_delay)
+    p = plot_model(variable,parameter_range_list,univariate_outlist,default_parameters.infection_introduction_day,default_parameters.immunization_begin_day)
     return p
 end
 

CovidAlertVaccinationModel/abm_timeseries.jl

 using CovidAlertVaccinationModel
 using OnlineStats
 using Plots
-const samples = 1
+const samples = 10
 
 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]
@@ -10,23 +10,27 @@ const ymo_attack_rate = [10.376,5.636,7.2]./100
 
 function solve_and_plot_parameters()
     p =  CovidAlertVaccinationModel.get_parameters()
+    display(p)
     out,avg_populations = mean_solve(samples, p,DebugRecorder)
     p = plot_model(nothing,[nothing],[out],p.infection_introduction_day,p.immunization_begin_day)
 
     ymo_vaccination_ts = mean.(out.daily_immunized_by_age)
 
-    total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,1180:end]))
+    total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
     final_size = sum.(eachrow(mean.(out.daily_unvac_cases_by_age)))
-    total_preinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,1:1180]))
+    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
 
   
-    display((final_size,target_final_size))
-    display((total_preinf_vaccination,target_preinf_vac))
-    display((total_postinf_vaccination,target_postinf_vac))
-    display(sum.(eachrow(ymo_vaccination_ts)) ./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.(out.daily_cases_by_age))) 
+    println("vac + unvac cases proportion: $(total_final_size./avg_populations))")
+    display(sum.(eachrow(ymo_vaccination_ts)) ./avg_populations) 
 
     savefig(p,"timeseries.pdf")
     return out

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 = 500,
+        sim_length = 400,
         num_households = 5000,
         I_0_fraction = 0.005,
-        β_y = 0.00042,
-        β_m = 0.00042,
-        β_o = 0.1,
-        α_y = 0.35,
-        α_m = 0.35,
-        α_o = 0.35,
+        β_y = 0.00046,
+        β_m = 0.00044,
+        β_o = 0.6,
+        α_y = 0.4,
+        α_m = 0.4,
+        α_o = 0.4,
         recovery_rate = 1/5,
-        π_base_y = -1.0,
-        π_base_m = -1.0,
-        π_base_o = -1.0,
+        π_base_y = -1.30,
+        π_base_m = -1.31,
+        π_base_o = -0.95,
         η = 0.0,
         κ = 0.0,
         ω = 0.00,
@@ -20,14 +20,13 @@ function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35
         γ = 0.0,
         ξ = 5.0,
         notification_parameter = 0.001,
-        vaccinator_prob = 0.2,
+        vaccinator_prob = 0.6,
         app_user_fraction = 0.0,
         notification_threshold = 20,
         immunizing = true,
         immunization_delay = 14,
         immunization_begin_day = 60, 
         infection_introduction_day = 180,
-        O_distribution_shift = 0.0,
     )
     return params
 end
@@ -86,14 +85,14 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
         ws_mixing_tuple_preshift = deepcopy(workschool_mixing)
         rest_mixing_tuple_preshift = deepcopy(rest_mixing)
 
-        for md in ws_mixing_tuple_preshift
-            adjust_distributions_mean!(md[1:3,3],params.O_distribution_shift)
-            adjust_distributions_mean!(md[3,1:2],params.O_distribution_shift)#dont shift OO twice
-        end
-        for md in rest_mixing_tuple_preshift
-            adjust_distributions_mean!( md[1:3,3],params.O_distribution_shift)
-            adjust_distributions_mean!( md[3,1:2],params.O_distribution_shift) #dont shift OO twice
-        end
+        # for md in ws_mixing_tuple_preshift
+        #     adjust_distributions_mean!(md[1:3,3],params.O_distribution_shift)
+        #     adjust_distributions_mean!(md[3,1:2],params.O_distribution_shift)#dont shift OO twice
+        # end
+        # for md in rest_mixing_tuple_preshift
+        #     adjust_distributions_mean!( md[1:3,3],params.O_distribution_shift)
+        #     adjust_distributions_mean!( md[3,1:2],params.O_distribution_shift) #dont shift OO twice
+        # end
         map_symmetrize(m_tuple) = map(md -> symmetrize_means(pop_sizes,md), m_tuple)
 
         ws_matrix_tuple = map_symmetrize(ws_mixing_tuple_preshift) 

CovidAlertVaccinationModel/src/ABM/output.jl

@@ -157,5 +157,5 @@ function plot_model(varname,univariate_series, output_list::Vector{T},infection_
         vline!(p,[infection_begin]; label = "infection begin", line =:dot)
         vline!(p,[vac_begin]; label = "vaccination begin",line = :dot)
     end
-    return plot(plts...;layout = (l,1),size=(800,400*l),leftmargin = 12Plots.mm)
+    return plot(plts...;layout = (l,1),size=(800,400*l),leftmargin = 12Plots.mm, legend = :outerright)
 end

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -106,7 +106,7 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
         random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
 
         if !isempty(neighbors(random_soc_network,i))
-            for _ = 1:10
+            for _ = 1:1
                 random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,i))
 
                 if u_vac[random_neighbour] == u_vac[i]

CovidAlertVaccinationModel/test/ABM/mixing_test.jl

 
 using CovidAlertVaccinationModel
-using CovidAlertVaccinationModel:ModelSolution,AgentDemographic,mean,AgentStatus,get_u_0,get_parameters,solve!, DebugRecorder
+using CovidAlertVaccinationModel:ModelSolution,AgentDemographic,mean,AgentStatus,get_u_0,get_parameters,solve!, DebugRecorder, WeightedGraph
 using Random
 using Plots
 using ThreadsX
@@ -30,5 +30,57 @@ const model_sizes = [1000,5000]
     end
 end
 
+function approximate_mixing_matricies_weights_dict(p,sol)
+    mean_mixing_degree = zeros(3,3)
+    mean_mixing_weighted_degree = zeros(3,3)
+        for g_list in sol.inf_network.graph_list
+            for (k,g) in enumerate(g_list)
+                for e in keys(g.mixing_edges.weights_dict)
+                    demo_i = sol.demographics[e.a]
+                    demo_j = sol.demographics[e.b]
+                    mean_mixing_degree[Int(demo_i), Int(demo_j)] += 1
+                    mean_mixing_degree[Int(demo_j), Int(demo_i)] += 1
+                    mean_mixing_weighted_degree[Int(demo_i), Int(demo_j)] += get_weight(g,GraphEdge(e.a,e.b))
+                    mean_mixing_weighted_degree[Int(demo_j), Int(demo_i)] += get_weight(g,GraphEdge(e.a,e.b))
+                end 
+        end
+    end
+    mean_mixing_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
+    mean_mixing_weighted_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
+    return mean_mixing_degree, mean_mixing_weighted_degree
+end
+
+function approximate_mixing_matricies_graph(p,sol)
+    mean_mixing_degree = zeros(3,3)
+    mean_mixing_weighted_degree = zeros(3,3)
+        for g_list in sol.inf_network.graph_list
+            for g in g_list
+                for e in edges(g.g)  
+                        demo_i = sol.demographics[src(e)]
+                        demo_j = sol.demographics[dst(e)]
+                        mean_mixing_degree[Int(demo_i), Int(demo_j)] += 1#get_weight(g,GraphEdge(node_i,node_j)) /2
+                        mean_mixing_degree[Int(demo_j), Int(demo_i)] += 1#get_weight(g,GraphEdge(node_i,node_j)) /2
+                        mean_mixing_weighted_degree[Int(demo_i), Int(demo_j)] += get_weight(g,GraphEdge(src(e),dst(e)))
+                        mean_mixing_weighted_degree[Int(demo_j), Int(demo_i)] += get_weight(g,GraphEdge(src(e),dst(e)))
+                end 
+            end
+        end
+        mean_mixing_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
+        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)
+
+    @test all(graph_deg_matrix .≈ dict_deg_matrix)
+    @test all(graph_wdeg_matrix .≈ dict_wdeg_matrix)
+end
+
+
 
+@testset "degree distribution generation" for _ in samples
 
+end
\ No newline at end of file