Commit e1f51783 authored by Peter Jentsch's avatar Peter Jentsch
Browse files

new parameter fit.. still not working, manually tuned parameter fit

parent 5779702d
using CovidAlertVaccinationModel using CovidAlertVaccinationModel
using OnlineStats using OnlineStats
using Plots using Plots
const samples = 1 const samples = 5
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]
const ymo_attack_rate = [10.376,5.636,7.2]./100
function solve_and_plot_parameters() function solve_and_plot_parameters()
p = CovidAlertVaccinationModel.get_parameters() p = CovidAlertVaccinationModel.get_parameters()
out = mean_solve(samples, p,DebugRecorder) out,avg_populations = mean_solve(samples, p,DebugRecorder)
p = plot_model(nothing,[nothing],[out],p.infection_introduction_day,p.immunization_begin_day) p = plot_model(nothing,[nothing],[out],p.infection_introduction_day,p.immunization_begin_day)
savefig(p,"timeseries.pdf")
return out
end
out = solve_and_plot_parameters() ymo_vaccination_ts = mean.(out.daily_immunized_by_age)
total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,1180:end]))
final_size = sum.(eachrow(mean.(out.daily_unvac_cases_by_age)))
total_preinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,1:1180]))
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(sum.(eachrow(mean.(out.daily_immunized_by_age))))
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)
savefig(p,"timeseries.pdf")
return out
end
final_size = sum.(eachrow(mean.(out.daily_cases_by_age))) out = solve_and_plot_parameters()
display(final_size)
\ No newline at end of file
using CovidAlertVaccinationModel:abm,ModelSolution, get_parameters, neighbors, GraphEdge, get_weight using CovidAlertVaccinationModel:abm,ModelSolution, get_parameters, neighbors, GraphEdge, get_weight
function approximate_mixing_matricies() using LightGraphs
p = get_parameters() function approximate_mixing_matricies_1(sol)
sol = abm(p,nothing) mean_mixing_degree = zeros(3,3)
mean_mixing = zeros(3,3) mean_mixing_weighted_degree = zeros(3,3)
display(p.O_distribution_shift) display(p.O_distribution_shift)
te = 0
for g_list in sol.inf_network.graph_list for g_list in sol.inf_network.graph_list
for g in g_list for g in g_list
# display(g.mixing_edges.weights_dict) # display(g.mixing_edges.weights_dict)
for e in keys(g.mixing_edges.weights_dict) for e in keys(g.mixing_edges.weights_dict)
demo_i = sol.demographics[e.a] demo_i = sol.demographics[e.a]
demo_j = sol.demographics[e.b] demo_j = sol.demographics[e.b]
mean_mixing[Int(demo_i), Int(demo_j)] += 1#get_weight(g,GraphEdge(node_i,node_j)) /2 if(has_edge(g.g,e.a,e.b))
mean_mixing[Int(demo_j), Int(demo_i)] += 1#get_weight(g,GraphEdge(node_i,node_j)) /2 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(e.a,e.b))
mean_mixing_weighted_degree[Int(demo_j), Int(demo_i)] += get_weight(g,GraphEdge(e.a,e.b))
else
if(demo_i == demo_j)
te += 1
end
end
end end
end end
end end
mean_mixing ./= (500 .* length.(sol.index_vectors) * 2) mean_mixing_degree ./= (500 .* length.(sol.index_vectors) * 2)
display(mean_mixing) mean_mixing_weighted_degree ./= (500 .* length.(sol.index_vectors) * 2)
display(mean_mixing_degree)
display(mean_mixing_weighted_degree)
display(te)
return sol return sol
end end
approximate_mixing_matricies(); function approximate_mixing_matricies_2(sol)
mean_mixing_degree = zeros(3,3)
mean_mixing_weighted_degree = zeros(3,3)
display(p.O_distribution_shift)
for g_list in sol.inf_network.graph_list
for g in g_list
# display(g.mixing_edges.weights_dict)
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 ./= (500 .* length.(sol.index_vectors) * 2)
mean_mixing_weighted_degree ./= (500 .* length.(sol.index_vectors) * 2)
display(mean_mixing_degree)
display(mean_mixing_weighted_degree)
return sol
end
p = get_parameters()
sol = abm(p,nothing)
approximate_mixing_matricies_1(sol);
approximate_mixing_matricies_2(sol);
println("dasdlkjasdas") println("dasdlkjasdas")
#rerun with weights #rerun with weights
# #
\ No newline at end of file
...@@ -37,6 +37,7 @@ function sample_mixing_graph!(mixing_graph) ...@@ -37,6 +37,7 @@ function sample_mixing_graph!(mixing_graph)
rand!(Random.default_rng(Threads.threadid()), mixing_edges.sampler_matrix[j,i],mixing_edges.sample_cache[j,i]) rand!(Random.default_rng(Threads.threadid()), mixing_edges.sampler_matrix[j,i],mixing_edges.sample_cache[j,i])
for k in 1:length(mixing_edges.contact_array[j,i]) for k in 1:length(mixing_edges.contact_array[j,i])
edge_weight_k = mixing_edges.sample_cache[j,i][k] edge_weight_k = mixing_edges.sample_cache[j,i][k]
set!(mixing_edges.weights_dict, mixing_edges.contact_array[j,i][k], edge_weight_k) set!(mixing_edges.weights_dict, mixing_edges.contact_array[j,i][k], edge_weight_k)
end end
end end
...@@ -90,7 +91,6 @@ function create_mixing_edges(demographic_index_vectors,mixing_matrix,weights_dis ...@@ -90,7 +91,6 @@ function create_mixing_edges(demographic_index_vectors,mixing_matrix,weights_dis
num_degrees_ij = similar(demographic_index_vectors[i]) num_degrees_ij = similar(demographic_index_vectors[i])
num_degrees_ji = similar(demographic_index_vectors[j]) num_degrees_ji = similar(demographic_index_vectors[j])
generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],num_degrees_ij,num_degrees_ji) generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],num_degrees_ij,num_degrees_ji)
m = sum(num_degrees_ij) m = sum(num_degrees_ij)
stubs_i = Vector{Int}(undef,m) stubs_i = Vector{Int}(undef,m)
stubs_j = Vector{Int}(undef,m) stubs_j = Vector{Int}(undef,m)
...@@ -101,8 +101,8 @@ function create_mixing_edges(demographic_index_vectors,mixing_matrix,weights_dis ...@@ -101,8 +101,8 @@ function create_mixing_edges(demographic_index_vectors,mixing_matrix,weights_dis
end end
contact_array[j,i] = GraphEdge.(stubs_i,stubs_j) contact_array[j,i] = GraphEdge.(stubs_i,stubs_j)
end end
return MixingEdges(tot,contact_array,weights_distribution_matrix) return MixingEdges(tot,contact_array,weights_distribution_matrix)
end end
...@@ -203,7 +203,10 @@ function graph_from_mixing_edges(g,mixing_edges) ...@@ -203,7 +203,10 @@ function graph_from_mixing_edges(g,mixing_edges)
for i in 1:size(mixing_edges.contact_array)[1], j in 1:i #diagonal for i in 1:size(mixing_edges.contact_array)[1], j in 1:i #diagonal
for k in 1:length(mixing_edges.contact_array[i,j]) for k in 1:length(mixing_edges.contact_array[i,j])
e = mixing_edges.contact_array[i,j][k] e = mixing_edges.contact_array[i,j][k]
add_edge!(g,e.a,e.b) success = add_edge!(g,e.a,e.b)
# if !success
# display((has_edge(g,e.a,e.b)))
# end
end end
end end
end end
......
function get_parameters() 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 = ( params = (
sim_length = 500, sim_length = 1300,
num_households = 5000, num_households = 5000,
I_0_fraction = 0.005, I_0_fraction = 0.005,
β_y = 0.001, β_y =0.00051,
β_m = 0.001, β_m = 0.00047,
β_o = 0.001, β_o = 0.07,
α_y = 1.0, α_y = 0.35,
α_m = 1.0, α_m = 0.35,
α_o = 1.0, α_o = 0.35,
recovery_rate = 1/5, recovery_rate = 1/5,
π_base_y = -5.0, π_base_y = -1.35,
π_base_m = -5.0, π_base_m = -1.35,
π_base_o = 0.0, π_base_o = -0.9,
η = 0.0, η = 0.0,
κ = 0.0, κ = 0.0,
ω = 0.00005, ω = 0.00,
ω_en = 0.00, ω_en = 0.00,
γ = 0.0, γ = 0.0,
ξ = 5.0, ξ = 5.0,
notification_parameter = 0.001, notification_parameter = 0.001,
vaccinator_prob = 0.2, vaccinator_prob = 0.2,
app_user_fraction = 0.0, app_user_fraction = 0.0,
notification_threshold = 2, notification_threshold = 20,
immunizing = true, immunizing = true,
immunization_delay = 14, immunization_delay = 14,
immunization_begin_day =60, immunization_begin_day = 1060,
infection_introduction_day = 180, infection_introduction_day = 1180,
O_distribution_shift = 1.0, O_distribution_shift = 0.5,
) )
return params return params
end end
......
...@@ -8,66 +8,12 @@ const ymo_attack_rate = [10.376,5.636,7.2]./100 ...@@ -8,66 +8,12 @@ const ymo_attack_rate = [10.376,5.636,7.2]./100
function solve_w_parameters(default_p, p_names, new_p_list) function solve_w_parameters(default_p, p_names, new_p_list)
new_params = merge(default_p, NamedTuple{p_names}(ntuple(i -> new_p_list[i],length(p_names)))) new_params = merge(default_p, NamedTuple{p_names}(ntuple(i -> new_p_list[i],length(p_names))))
out = DebugRecorder(0,default_p.sim_length) out = DebugRecorder(0,default_p.sim_length)
model = abm(new_params,out) model = abm(new_params,out)
return out, model return out, model
end end
function fit_all_parameters(p_tuple)
p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o)
priors = Factored(
Uniform(0.0,0.01),
Uniform(0.0,0.2),
Uniform(0.0,0.1),
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,
#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
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))
target_final_size = ymo_attack_rate .* length.(model.index_vectors)
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))
# return sum((normalize_by_pop(total_preinfection_vaccination .- target_ymo_vac)).^2)
# + sum(normalize_by_pop((total_postinf_vaccination .- target_ymo_vac)).^2)
# + sum(normalize_by_pop((final_size .- target_final_size)).^2)
display(final_size)
display( sum.(eachrow(ymo_vaccination_ts[:,1:end])))
display( length.(model.index_vectors))
return sum(normalize_by_pop((final_size .- target_final_size)).^2)
end
display(cost((0.0,0.0005,0.0005,1.0,-3.0,-3.0,-0.2,0.25,0.25,0.0)))
# 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 fit_distribution_parameters(p_tuple) function fit_distribution_parameters(p_tuple)
p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o,:O_distribution_shift) p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o,:O_distribution_shift)
priors = Factored( priors = Factored(
...@@ -117,7 +63,7 @@ function fit_distribution_parameters(p_tuple) ...@@ -117,7 +63,7 @@ function fit_distribution_parameters(p_tuple)
+ sum((total_postinf_vaccination .- total_postinf_vaccination).^2) + sum((total_postinf_vaccination .- total_postinf_vaccination).^2)
+ 5*sum((final_size .- target_final_size).^2) + 5*sum((final_size .- target_final_size).^2)
end end
display(cost((0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35,0.35,0.2))) 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) # out = smc(priors,cost; verbose = true, nparticles = 600, parallel = true)
return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names))) return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
...@@ -126,9 +72,6 @@ function fit_parameters(default_parameters) ...@@ -126,9 +72,6 @@ function fit_parameters(default_parameters)
# pre_inf_behaviour_parameters_path =joinpath(PACKAGE_FOLDER,"abm_parameter_fits","pre_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") # 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") fit_all_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","fit_all_parameters.dat")
fit_dist_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","fit_dist_parameters.dat")
output = fit_distribution_parameters(default_parameters) output = fit_distribution_parameters(default_parameters)
serialize(fit_all_parameters_path,output) serialize(fit_all_parameters_path,output)
...@@ -144,7 +87,6 @@ function fit_parameters(default_parameters) ...@@ -144,7 +87,6 @@ function fit_parameters(default_parameters)
end end
function plot_fitting_posteriors(fname,particles_tuple,parameters) function plot_fitting_posteriors(fname,particles_tuple,parameters)
# p_adjust = merge(parameters,(β_y = p.β_y*0.1, β_m = p.β_m*0.1 ,β_o = p.β_o))
stat_recorder = DebugRecorder(Variance(), parameters.sim_length) stat_recorder = DebugRecorder(Variance(), parameters.sim_length)
output_recorder = DebugRecorder(parameters.sim_length) output_recorder = DebugRecorder(parameters.sim_length)
avg_populations = [0.0,0.0,0.0] avg_populations = [0.0,0.0,0.0]
......
...@@ -150,7 +150,6 @@ function agents_step!(t,modelsol) ...@@ -150,7 +150,6 @@ function agents_step!(t,modelsol)
for network in modelsol.inf_network.graph_list[t] #this also resamples the soc network weights since they point to the same objects, but those are never used for network in modelsol.inf_network.graph_list[t] #this also resamples the soc network weights since they point to the same objects, but those are never used
sample_mixing_graph!(network) #get new contact weights sample_mixing_graph!(network) #get new contact weights
end end
if t == modelsol.params.infection_introduction_day 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)) 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.u_inf[init_indices] .= Infected
......
No preview for this file type
No preview for this file type
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment