Commit 4d549016 authored by Peter Jentsch's avatar Peter Jentsch
Browse files

first heatmaps

parent 2bda668f
...@@ -7,6 +7,7 @@ version = "0.1.0" ...@@ -7,6 +7,7 @@ version = "0.1.0"
AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5" AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
CurveFit = "5a033b19-8c74-5913-a970-47c3779ef25c" CurveFit = "5a033b19-8c74-5913-a970-47c3779ef25c"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
......
...@@ -6,7 +6,7 @@ using CovidAlertVaccinationModel:vaccination_data,ymo_vac,ymo_attack_rate ...@@ -6,7 +6,7 @@ using CovidAlertVaccinationModel:vaccination_data,ymo_vac,ymo_attack_rate
const samples = 10 const samples = 10
function solve_and_plot_parameters() function solve_and_plot_parameters()
p = CovidAlertVaccinationModel.get_parameters() p = CovidAlertVaccinationModel.get_app_parameters()
display(p) display(p)
out,avg_populations = 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)
...@@ -14,7 +14,7 @@ function solve_and_plot_parameters() ...@@ -14,7 +14,7 @@ function solve_and_plot_parameters()
ymo_vaccination_ts = mean.(out.daily_immunized_by_age) ymo_vaccination_ts = mean.(out.daily_immunized_by_age)
total_postinf_vaccination = mean.(out.total_postinf_vaccination)#sum.(eachrow(ymo_vaccination_ts[:,180:end])) total_postinf_vaccination = mean.(out.total_postinf_vaccination)#sum.(eachrow(ymo_vaccination_ts[:,180:end]))
final_size = mean.(out.final_size_by_age)#sum.(eachrow(mean.(out.daily_unvac_cases_by_age))) final_size = mean.(out.unvac_final_size_by_age)#sum.(eachrow(mean.(out.daily_unvac_cases_by_age)))
total_preinf_vaccination = mean.(out.total_preinf_vaccination)#sum.(eachrow(ymo_vaccination_ts[:,1:180])) total_preinf_vaccination = mean.(out.total_preinf_vaccination)#sum.(eachrow(ymo_vaccination_ts[:,1:180]))
target_final_size = ymo_attack_rate .*avg_populations target_final_size = ymo_attack_rate .*avg_populations
target_preinf_vac = ymo_vac .* sum(vaccination_data[1:4]) .* avg_populations target_preinf_vac = ymo_vac .* sum(vaccination_data[1:4]) .* avg_populations
......
using CovidAlertVaccinationModel:abm,ModelSolution, get_parameters, neighbors, GraphEdge, get_weight, contact_time_distributions
using LightGraphs
function approximate_mixing_matricies_1(sol)
mean_mixing_degree = zeros(3,3)
mean_mixing_weighted_degree = zeros(3,3)
display(p.O_distribution_shift)
te_missed = zeros(8)
te = zeros(8)
edgelist = []
for g_list in sol.inf_network.graph_list
for (k,g) in enumerate(g_list)
# display(g.mixing_edges.weights_dict)
for e in keys(g.mixing_edges.weights_dict)
demo_i = sol.demographics[e.a]
demo_j = sol.demographics[e.b]
if(has_edge(g.g,e.a,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))
te[k] += 1
else
push!(edgelist,e)
te_missed[k] += 1
end
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)
display(te_missed)
return sol
end
function approximate_mixing_matricies_2(sol)
mean_mixing_degree = zeros(3,3)
mean_mixing_weighted_degree = zeros(3,3)
te_missed = zeros(8)
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)
return sol
end
p = get_parameters()
sol = abm(p,nothing)
approximate_mixing_matricies_1(sol);
approximate_mixing_matricies_2(sol);
println("dasdlkjasdas")
#rerun with weights
#
\ No newline at end of file
...@@ -15,8 +15,8 @@ function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35 ...@@ -15,8 +15,8 @@ function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35
π_base_o = -0.95, π_base_o = -0.95,
η = 0.0, η = 0.0,
κ = 0.0, κ = 0.0,
ω = 0.0055, ω = 0.0,
ω_en = 0.00, ω_en = 0.0005,
Γ = 1/7, Γ = 1/7,
ξ = 5.0, ξ = 5.0,
notification_parameter = 0.001, notification_parameter = 0.001,
......
...@@ -17,7 +17,7 @@ if !ispath(univariate_path) ...@@ -17,7 +17,7 @@ if !ispath(univariate_path)
mkdir(univariate_path) mkdir(univariate_path)
end end
function univariate_simulations() function univariate_simulations()
len = 10 len = 20
univarate_test_list = ( univarate_test_list = (
# (:I_0_fraction, range(0.0, 0.05; length = len)), # (:I_0_fraction, range(0.0, 0.05; length = len)),
# (:base_transmission_probability, range(0.0002, 0.002; length = len)), # (:base_transmission_probability, range(0.0002, 0.002; length = len)),
...@@ -48,14 +48,14 @@ end ...@@ -48,14 +48,14 @@ end
using AxisKeys using AxisKeys
function multivariate_simulations() function multivariate_simulations()
len = 10 len = 15
samples = 10 samples = 10
app_simulations = ( app_simulations = (
(:η, range(0.0, 0.01; length = len)), (:η, range(0.0, 0.01; length = len)),
(:ω_en, range(0.0, 0.0005; length = len)), (:ω_en, range(0.0, 0.05; length = len)),
# (:notification_threshold, (1:len)), # (:notification_threshold, (1:len)),
) )
run_multivariate_sims(app_simulations,1) run_multivariate_sims(app_simulations,samples)
# for ((varname,_),p) in zip(univarate_test_list,plt_list) # for ((varname,_),p) in zip(univarate_test_list,plt_list)
...@@ -65,26 +65,52 @@ end ...@@ -65,26 +65,52 @@ end
using ProgressMeter using ProgressMeter
function run_multivariate_sims(sims,samples) function run_multivariate_sims(sims,samples)
varnames, sim_ranges = zip(sims...) varnames, sim_ranges = zip(sims...)
without_app_future = @spawn mean_solve(samples,get_parameters(),HeatmapRecorder)
default_parameters = get_app_parameters() default_parameters = get_app_parameters()
simvars = Iterators.product(sim_ranges...) simvars = Iterators.product(sim_ranges...)
progmeter = Progress(length(simvars)) progmeter = Progress(length(simvars))
output = ThreadsX.map(simvars) do vars app_output = ThreadsX.map(simvars) do vars
vars_with_names = NamedTuple{varnames}(vars) vars_with_names = NamedTuple{varnames}(vars)
parameters = merge(default_parameters,vars_with_names) parameters = merge(default_parameters,vars_with_names)
out,_ = mean_solve(samples, parameters,DebugRecorder) out,_ = mean_solve(1, parameters,HeatmapRecorder)
next!(progmeter) next!(progmeter)
return out return out
end end
display(length(simvars)) display(length(simvars))
fname = join(string.(varnames),"_") fname = join(string.(varnames),"_")
keyed_output = KeyedArray(output;NamedTuple{varnames}(sim_ranges)...) keyed_output = KeyedArray(app_output;NamedTuple{varnames}(sim_ranges)...)
path = joinpath(PACKAGE_FOLDER,"abm_output","$fname.dat") path = joinpath(PACKAGE_FOLDER,"abm_output","$fname.dat")
serialize(path,keyed_output) serialize(path,(fetch(without_app_future)[1],keyed_output))
return fname return fname
end end
using ColorSchemes
using LaTeXStrings
function plot_parameter_plane(input_fname) function plot_parameter_plane(input_fname)
map(1: length(axiskeys(output)[end])) do i path = joinpath(PACKAGE_FOLDER,"abm_output","$input_fname.dat")
output_no_app, output = deserialize(path)
var_ranges = axiskeys(output)
vars = (L"\eta",L"\omega_{en}")
mean_final_size(p) = mean(reduce(merge!,p.final_size_by_age))
base_outcome = mean_final_size(output_no_app)
final_size_map = map(x-> (mean_final_size(x) - base_outcome),output)
mean_weighted_degree_change(p) = mean(p.avg_weighted_degree_of_vaccinators)-mean(p.avg_weighted_degree)
weighted_degree_map = map(mean_weighted_degree_change,output)
cs = cgrad([:blue,:orange])
datamaps = [weighted_degree_map, final_size_map]
fnames = ["wdg_change.pdf", "final_size_change.pdf"]
titles = [
"Average w. deg. of vaccinators minus average w. deg.",
"Change in average final size of infection outbreak",
]
for (fname,title,datamap) in zip(fnames,titles,datamaps)
p = heatmap(var_ranges[1],var_ranges[2],datamap; title, xlabel = vars[1], ylabel = vars[2], seriescolor=cs, size = (800,600))
savefig(p,joinpath(PACKAGE_FOLDER,"plots","app_heatmaps","$fname"))
end end
end end
\ No newline at end of file
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