From 189c2fa4d2aea281981454c206957a8bf6417abc Mon Sep 17 00:00:00 2001 From: pjentsch <pjentsch@uwaterloo.ca> Date: Wed, 21 Apr 2021 22:16:47 -0400 Subject: [PATCH] improve performance in model loading --- CovidAlertVaccinationModel/src/ABM/abm.jl | 4 +-- CovidAlertVaccinationModel/src/ABM/agents.jl | 37 +++++++++----------- CovidAlertVaccinationModel/src/ABM/solve.jl | 18 ++++++++-- CovidAlertVaccinationModel/src/data.jl | 16 ++++----- 4 files changed, 42 insertions(+), 33 deletions(-) diff --git a/CovidAlertVaccinationModel/src/ABM/abm.jl b/CovidAlertVaccinationModel/src/ABM/abm.jl index 687088d..7b78faf 100644 --- a/CovidAlertVaccinationModel/src/ABM/abm.jl +++ b/CovidAlertVaccinationModel/src/ABM/abm.jl @@ -2,7 +2,7 @@ const population = 14.57e6 #population of ontario const workschool_mixing, rest_mixing = load_mixing_matrices() const age_bins = [(0.0, 25.0),(25.0,65.0),(65.0,Inf)] - +const household_data = read_household_data() default(dpi = 300) default(framestyle = :box) @@ -13,7 +13,7 @@ Runs the model with 5k households for 500 timesteps. """ function bench() Random.seed!(RNG,1) - steps = 500 + steps = 50 model_sol = ModelSolution(steps,get_parameters(),5000) recording = DebugRecorder(steps) output = solve!(model_sol,recording ) diff --git a/CovidAlertVaccinationModel/src/ABM/agents.jl b/CovidAlertVaccinationModel/src/ABM/agents.jl index 5b830aa..81f0937 100644 --- a/CovidAlertVaccinationModel/src/ABM/agents.jl +++ b/CovidAlertVaccinationModel/src/ABM/agents.jl @@ -26,23 +26,29 @@ end Generates a complete graph from a vector of household compositions. Each household composition is a 3 element vectors (one for each demographic group) of integers where each element describes the number of the corresponding demographic group present in that household. This function wires together a graph such that each household is in a complete subgraph. It is much faster than the `mapreduce` solution used previously. + +Also returns a vector of AgentDemographic representing each agent, defined by the household compositions. """ -function complete_graph_from_households_composition(households_composition) +@views function complete_graph_from_households_composition(households_composition) total_household_pop = sum(sum.(households_composition)) + population_list = Vector{AgentDemographic}(undef,total_household_pop) network = SimpleGraph(total_household_pop) vertex_pointer = 1 for household in households_composition num_vertices = sum(household) - for v in vertex_pointer:(vertex_pointer + num_vertices - 1) - for w in vertex_pointer:(vertex_pointer + num_vertices - 1) - if v != w - add_edge!(network,v,w) - end + household_pointer = 0 + for (k,size) in enumerate(household) + @inbounds population_list[vertex_pointer+household_pointer:vertex_pointer+household_pointer+size-1] .= AgentDemographic(k) + household_pointer += size + end + for v in vertex_pointer:(vertex_pointer + num_vertices - 1), w in vertex_pointer:(vertex_pointer + num_vertices - 1) + if v != w + add_edge!(network,v,w) end end vertex_pointer+=num_vertices end - return network + return network,population_list end @@ -59,19 +65,8 @@ index_vectors: Vector of 3 vectors, each of which contains the indexes of all th """ function generate_population(num_households) households_composition = sample_household_data(num_households) - - household_networks = complete_graph_from_households_composition(households_composition) - - #test the complete_graph_from_households_composition function against the mapreduce version to ensure it is correct. - @c_assert all(adjacency_matrix(household_networks) .== adjacency_matrix( mapreduce(LightGraphs.complete_graph,blockdiag, sum.(households_composition); init = SimpleGraph()))) - - households = map( l -> [fill(AgentDemographic(i),n) for (i,n) in enumerate(l)],households_composition) |> - l -> reduce(vcat,l) |> - l -> reduce(vcat,l) - population_list = reduce(vcat,households) + household_networks,population_list = complete_graph_from_households_composition(households_composition) index_vectors = [findall(x -> x == AgentDemographic(i), population_list) for i in 1:(AgentDemographic.size-1)] - - return (; population_list, household_networks, @@ -81,7 +76,7 @@ end """ -Defines pretty printing methods so that `display(AgentStatus)` is more readable. +Defines pretty printing methods so that `display(s::AgentStatus)` is more readable. """ function Base.show(io::IO, status::AgentStatus) if status == Susceptible @@ -96,7 +91,7 @@ function Base.show(io::IO, status::AgentStatus) end """ -Defines pretty printing methods so that `display(AgentDemographic)` is more readable. +Defines pretty printing methods so that `display(s::AgentDemographic)` is more readable. """ function Base.show(io::IO, status::AgentDemographic) if status == Young diff --git a/CovidAlertVaccinationModel/src/ABM/solve.jl b/CovidAlertVaccinationModel/src/ABM/solve.jl index 4470117..737121f 100644 --- a/CovidAlertVaccinationModel/src/ABM/solve.jl +++ b/CovidAlertVaccinationModel/src/ABM/solve.jl @@ -1,10 +1,19 @@ -# using LoopVectorization + + + + function contact_weight(p, contact_time) return 1 - (1-p)^contact_time end + + + Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) +#remove Base.@propagate_inbounds if you get segfaults @unpack notification_parameter = modelsol.params @unpack time_of_last_alert, app_user_index,inf_network,covid_alert_times,app_user = modelsol + + for (i,node) in enumerate(modelsol.app_user_index), mixing_graph in modelsol.inf_network.graph_list[t] for j in 2:14 covid_alert_times[j-1,i] = covid_alert_times[j,i] #shift them all back @@ -15,13 +24,17 @@ Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) covid_alert_times[end,i] += get_weight(mixing_graph,GraphEdge(node,j)) #add the contact times for today to the back end end - if rand(RNG) < 1 - (1- notification_parameter)^sum(covid_alert_times[:,i]) + + covid_alert_total_exposures = 1 - (1 - notification_parameter) ^ sum(covid_alert_times[:,i]) + if rand(RNG) < covid_alert_total_exposures time_of_last_alert[i] = t end end end Base.@propagate_inbounds @views function update_infection_state!(t,modelsol) +#remove Base.@propagate_inbounds if you get segfaults + @unpack base_transmission_probability,immunization_loss_prob,recovery_rate = modelsol.params @unpack u_inf,u_vac,u_next_inf,u_next_vac,demographics,inf_network,status_totals = modelsol @@ -65,6 +78,7 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol) # display(u_next_inf) end Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections) +#remove Base.@propagate_inbounds if you get segfaults @unpack Ï€_base, η,γ, κ, ω, Ï, ω_en,Ï_en,γ,β = modelsol.params @unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list = modelsol diff --git a/CovidAlertVaccinationModel/src/data.jl b/CovidAlertVaccinationModel/src/data.jl index a84afeb..5532743 100644 --- a/CovidAlertVaccinationModel/src/data.jl +++ b/CovidAlertVaccinationModel/src/data.jl @@ -35,20 +35,20 @@ function find_household_composition(df_row) age_distribution[age_resp_to_bin[df_row[:AGERESP]]] += 1 return SVector{3}(age_distribution) end -function sample_household_data(n) +function read_household_data() f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/home_compositions.csv"), ',') df = DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])]) weight_vector::Vector{Float64} = df[!,:WGHT_PER]/sum(df[!,:WGHT_PER]) households = map(find_household_composition,eachrow(df)) - return sample(RNG,households, Weights(weight_vector),n) - # https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en + return (;households,weight_vector) +end + +function sample_household_data(n) + return sample(RNG,household_data.households, Weights(household_data.weight_vector),n) end + function get_household_data_proportions() - f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/home_compositions.csv"), ',') - df = DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])]) - weight_vector::Vector{Float64} = df[!,:WGHT_PER]/sum(df[!,:WGHT_PER]) - households = map(find_household_composition,eachrow(df)) - households_by_demographic_sum = sum.([map(l-> l[i], households) for i in 1:3]) + households_by_demographic_sum = sum.([map(l-> l[i],household_data.households) for i in 1:3]) return households_by_demographic_sum./sum(households_by_demographic_sum) # https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en end -- GitLab