Merge fd0a662 into 5acdbae

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.2"
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
+BritishNationalGrid = "07503c7f-9676-5157-a8e2-6c0d0a9caee7"
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -67,16 +68,13 @@ IndexedTables = "0.12.6, 0.13, 1"
 Interpolations = "0.13.1"
 JLD = "0.10, 0.11, 0.12"
 JLSO = "2.3"
-MPI = "0.17.2"
 Measures = "0.3.1"
 Missings = "0.4.5, 1.0"
 NetCDF = "0.11.3"
 ORCA = "0.5"
 OnlineStats = "1.5.8"
 OnlineStatsBase = "1.4"
 Phylo = "0.4.18"
-Plots = "0.28.4, 0.29, 1"
-RCall = "0.13.10"
 RecipesBase = "0.7, 0.8, 1"
 Requires = "1"
 Revise = "3.1"

examples/Epidemiology/Scotland_run.jl

@@ -15,9 +15,8 @@ using SQLite
 
 function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, timestep::Unitful.Time; do_plot::Bool = false, do_download::Bool = true, save::Bool = false, savepath::String = pwd())
     # Download and read in population sizes for Scotland
-    DataRegistryUtils.load_array!(db, "human/demographics/population/scotland", "/grid1km/age/persons"; sql_alias="human_demographics_population_scotland_grid1km_age_persons_arr")
-    scotpop = get_3d_km_grid_axis_array(db, ["grid_x", "grid_y", "age_aggr"], "val", "scottish_population_view")
-
+    DataRegistryUtils.load_array!(db, "human/demographics/population/scotland", "/grid area/age/persons"; sql_alias="km_age_persons_arr")
+    scotpop = get_3d_km_grid_axis_array(db, ["grid_area", "age_aggr"], "val", "scottish_population_view")
     # Read number of age categories
     age_categories = size(scotpop, 3)
 
@@ -197,6 +196,7 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
 
     # Populate susceptibles according to actual population spread
     cat_idx = reshape(1:(numclasses * age_categories), age_categories, numclasses)
+    scotpop = shrink_to_active(scotpop)
     reshaped_pop =
         reshape(scotpop[1:size(epienv.active, 1), 1:size(epienv.active, 2), :],
                 size(epienv.active, 1) * size(epienv.active, 2), size(scotpop, 3))'

examples/Epidemiology/Scotland_run_view.sql

@@ -5,10 +5,9 @@ CREATE VIEW scottish_population_view AS
 SELECT age_groups
 , CAST(substr(age_groups, 4, 5) AS INT) AS age
 , (CAST(substr(age_groups, 4, 5) AS INT) / 10) * 10 AS age_aggr
-, CAST(substr(grid_area, 1, instr(grid_area, '-') - 1) AS INT) AS grid_x
-, CAST(substr(grid_area, instr(grid_area, '-') + 1) AS INT) AS grid_y
+, grid_area
 , val
-FROM human_demographics_population_scotland_grid1km_age_persons_arr;
+FROM km_age_persons_arr;
 
 CREATE VIEW pollution_grid_view AS
 SELECT pollutant

examples/Epidemiology/data_config.yaml

@@ -18,7 +18,7 @@ read:
       component: symptom-probability
   - where:
       data_product: human/demographics/population/scotland
-      version: 1.0.1
+      version: 1.0.2
   - where:
       data_product: geography/lookup_table/gridcell_admin_area/scotland
   - where:

examples/env_transitions.jl

@@ -45,7 +45,8 @@ virus_decay = 1.0/2days
 age_mixing = fill(1.0, ages, ages)
 param = (birth = birth, death = death, virus_growth = virus_growth,
         virus_decay = virus_decay, beta_env = beta_env,
-        beta_force = beta_force, age_mixing = age_mixing, env_exposure = 2.5)
+        beta_force = beta_force, age_mixing = age_mixing,
+        env_exposure = 2.0/mean(epienv.habitat.matrix))
 
 cat_idx = reshape(1:(nrow(abun_h) * ages), ages, nrow(abun_h))
 transitiondat = DataFrame([
@@ -81,9 +82,9 @@ for loc in eachindex(epienv.habitat.matrix)
     for age in 1:ages
         addtransition!(transitions, ForceProduce(cat_idx[age, 3], loc, param.virus_growth[age]))
         addtransition!(transitions, ForceDisperse(cat_idx[age, 3], loc))
-        addtransition!(transitions, EnvExposure(transitiondat[1, :from_id][age], loc,
-            transitiondat[1, :to_id][age], transitiondat[1, :prob].force[age], transitiondat[1, :prob].env[age],
-            param.env_exposure, x -> x.matrix))
+        addtransition!(transitions, EnvTransition(Exposure(transitiondat[1, :from_id][age], loc,
+            transitiondat[1, :to_id][age], transitiondat[1, :prob].force[age], transitiondat[1, :prob].env[age]),
+            param.env_exposure))
         addtransition!(transitions, Infection(transitiondat[2, :from_id][age], loc,
             transitiondat[2, :to_id][age], transitiondat[2, :prob][age]))
         addtransition!(transitions, Recovery(transitiondat[3, :from_id][age], loc,

examples/pollution_transitions.jl

@@ -13,11 +13,33 @@ using DataFrames
 using Plots
 using SQLite
 
-function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, timestep::Unitful.Time; do_plot::Bool = false, do_download::Bool = true, save::Bool = false, savepath::String = pwd())
+module PollutionRun
+    using EcoSISTEM
+    using Statistics
+    import EcoSISTEM.getprob
+    # Define how probabilities for transitions should be altered by pollution
+    getprob(eco::Ecosystem, rule::Exposure) = (rule.force_prob * 5.0 * get_env(eco.abenv.habitat, rule.location), rule.virus_prob)
+    getprob(eco::Ecosystem, rule::DevelopSymptoms) = rule.prob * 2.0 * get_env(eco.abenv.habitat, rule.location)
+    getprob(eco::Ecosystem, rule::Hospitalise) = rule.prob * 2.0 * get_env(eco.abenv.habitat, rule.location)
+    getprob(eco::Ecosystem, rule::DeathFromInfection) = rule.prob * 5.0 * get_env(eco.abenv.habitat, rule.location)
+    export getprob
+end
+using Main.PollutionRun
+
+const stochasticmode = false
 
+function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, timestep::Unitful.Time; do_plot::Bool = false, save::Bool = false, savepath::String = pwd(), death_boost::Float64 = 2.0)
     # Download and read in population sizes for Scotland
-    DataRegistryUtils.load_array!(db, "human/demographics/population/scotland", "/grid1km/age/persons"; sql_alias="human_demographics_population_scotland_grid1km_age_persons_arr")
-    scotpop = get_3d_km_grid_axis_array(db, ["grid_x", "grid_y", "age_aggr"], "val", "scottish_population_view")
+    DataRegistryUtils.load_array!(db, "human/demographics/population/scotland", "/grid area/age/persons"; sql_alias="km_age_persons_arr")
+    scotpop = get_3d_km_grid_axis_array(db, ["grid_area", "age_aggr"], "val", "scottish_population_view")
+    scotpop = shrink_to_active(scotpop)
+
+    DataRegistryUtils.load_array!(db, "records/pollution", "/array"; sql_alias="records_pollution_array_arr")
+    pollution = get_3d_km_grid_axis_array(db, ["grid_x", "grid_y"], "val", "pollution_grid_view", 1.0μg * m^-3)
+    pollution = pollution[531500m .. 1215000m, 54500m .. 465000m]
+
+    rngtype = stochasticmode ? Random.MersenneTwister : MedianGenerator
+    seed_fun = stochasticmode ? seedinfected! : deterministic_seed!
 
     # Read number of age categories
     age_categories = size(scotpop, 3)
@@ -49,14 +71,6 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
         (AxisArrays.axes(scotpop, 2)[end] + AxisArrays.axes(scotpop, 2)[2] -
         2 * AxisArrays.axes(scotpop, 2)[1]) * 1.0
 
-    # Sum up age categories and turn into simple matrix
-    total_pop = dropdims(sum(Float64.(scotpop), dims=3), dims=3)
-    total_pop = AxisArray(total_pop, AxisArrays.axes(scotpop)[1], AxisArrays.axes(scotpop)[2])
-    total_pop.data[total_pop .≈ 0.0] .= NaN
-
-    # Shrink to smallest bounding box. The NaNs are inactive.
-    total_pop = shrink_to_active(total_pop);
-
     # Prob of developing symptoms
     p_s = fill(read_estimate(db,
         "human/infection/SARS-CoV-2/%",
@@ -98,7 +112,6 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
     )[1] * Unitful.hr)
     @show T_asym
 
-
     # Time pre-symptomatic
     T_presym = 1.5days
     # Time symptomatic
@@ -172,9 +185,12 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
         (from="Hospitalised", from_id=cat_idx[:, 6], to="Dead", to_id=cat_idx[:, 8], prob=death_hospital)
     ])
 
-    epienv = simplehabitatAE(298.0K, size(total_pop), area, Lockdown(20days))
-
-    movement_balance = (home = fill(0.5, numclasses * age_categories), work = fill(0.5, numclasses * age_categories))
+    total_pop = dropdims(sum(Float64.(scotpop), dims=3), dims=3)
+    total_pop[total_pop .≈ 0.0] .= NaN
+    axis1 = pollution.axes[1]; axis2 = pollution.axes[2]
+    poll = pollution ./ mean(pollution)
+    pollution = AxisArray(poll, axis1, axis2)
+    epienv = ukclimateAE(pollution, area, Lockdown(20days))
 
     # Dispersal kernels for virus and disease classes
     dispersal_dists = fill(1.0km, length(total_pop))
@@ -185,25 +201,27 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
     # Import commuter data (for now, fake table)
     active_cells = findall(.!isnan.(total_pop[1:end]))
     from = active_cells
-    to = sample(active_cells, weights(total_pop[active_cells]), length(active_cells))
+    to = sample(rngtype(), active_cells, weights(total_pop[active_cells]), length(active_cells))
     count = round.(total_pop[to]/10)
     home_to_work = DataFrame(from=from, to=to, count=count)
     work = Commuting(home_to_work)
     movement = EpiMovement(home, work)
 
     # Traits for match to environment (turned off currently through param choice, i.e. virus matches environment perfectly)
     traits = GaussTrait(fill(298.0K, numvirus), fill(0.1K, numvirus))
+    movement_balance = (home = fill(0.5, numclasses * age_categories), work = fill(0.5, numclasses * age_categories))
     epilist = SpeciesList(traits, abun_v, abun_h, movement, transitiondat, param, age_categories, movement_balance)
     rel = Gauss{eltype(epienv.habitat)}()
 
     transitions = create_transition_list()
     addtransition!(transitions, UpdateEpiEnvironment(update_epi_environment!))
-    addtransition!(transitions, SeedInfection(seedinfected!))
+    addtransition!(transitions, SeedInfection(seed_fun))
+
+    initial_infecteds = 10_000
 
-    initial_infecteds = 100
     # Create epi system with all information
-    @time epi = Ecosystem(epilist, epienv, rel, total_pop, UInt32(1),
-    initial_infected = initial_infecteds, transitions = transitions)
+    @time epi = Ecosystem(epilist, epienv, rel, scotpop, UInt32(1),
+    initial_infected = initial_infecteds, transitions = transitions, rngtype = rngtype)
 
     for loc in epi.cache.ordered_active
         addtransition!(epi.transitions, ViralLoad(loc, param.virus_decay))
@@ -216,8 +234,8 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
             addtransition!(epi.transitions, ForceDisperse(cat_idx[age, 4], loc))
             addtransition!(epi.transitions, ForceDisperse(cat_idx[age, 5], loc))
             # Exposure
-            addtransition!(epi.transitions, Exposure(transitiondat[1, :from_id][age], loc,
-                transitiondat[1, :to_id][age], transitiondat[1, :prob].force[age], transitiondat[1, :prob].env[age]))
+            addtransition!(epi.transitions, Exposure(transitiondat[1, :from_id][age], loc, transitiondat[1, :to_id][age],
+                transitiondat[1, :prob].force[age], transitiondat[1, :prob].env[age]))
             # Infected but asymptomatic
             addtransition!(epi.transitions, Infection(transitiondat[2, :from_id][age], loc,
                 transitiondat[2, :to_id][age], transitiondat[2, :prob][age]))
@@ -247,22 +265,21 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
                 transitiondat[10, :to_id][age], transitiondat[10, :prob][age]))
         end
     end
-
-    # Populate susceptibles according to actual population spread
-    scotpop = shrink_to_active(scotpop)
-    reshaped_pop =
-        reshape(scotpop[1:size(epienv.active, 1), 1:size(epienv.active, 2), :],
-                size(epienv.active, 1) * size(epienv.active, 2), size(scotpop, 3))'
-    epi.abundances.matrix[cat_idx[:, 1], :] = reshaped_pop
-
-    N_cells = size(epi.abundances.matrix, 2)
+
+    # Store BNG names in abenv
+    norths = ustrip.(AxisArrays.axes(scotpop)[1].val)
+    easts = ustrip.(AxisArrays.axes(scotpop)[2].val)
+    n_e = collect(Base.product(norths, easts))[1:end]
+    epi.abenv.names = [get_bng(n[2], n[1]) for n in n_e]
 
     # Turn off work moves for <20s and >70s
     epi.spplist.species.home_balance[cat_idx[1:2, :]] .= 1.0
     epi.spplist.species.home_balance[cat_idx[7:10, :]] .= 1.0
     epi.spplist.species.work_balance[cat_idx[1:2, :]] .= 0.0
     epi.spplist.species.work_balance[cat_idx[7:10, :]] .= 0.0
 
+    N_cells = size(epi.abundances.matrix, 2)
+    println("Ecosystem initiated")
     # Run simulation
     abuns = zeros(UInt32, size(epi.abundances.matrix, 1), N_cells, floor(Int, times/timestep) + 1)
     @time simulate_record!(abuns, epi, times, interval, timestep, save = save, save_path = savepath)
@@ -283,20 +300,17 @@ function run_model(db::SQLite.DB, times::Unitful.Time, interval::Unitful.Time, t
             "Deaths" => cat_idx[:, 8],
         )
         display(plot_epidynamics(epi, abuns, category_map = category_map))
-        display(plot_epiheatmaps(epi, abuns, steps = [30]))
+        display(plot_epiheatmaps(epi, abuns, steps = [30], match_dimensions = false))
     end
     return abuns
 end
 
+cd("examples")
 data_dir= "Epidemiology/data/"
 config = "Epidemiology/data_config.yaml"
 view_sql = "Epidemiology/Scotland_run_view.sql"
 db = initialise_local_registry(data_dir, data_config = config, sql_file = view_sql)
 
 times = 1month; interval = 1day; timestep = 1day
-run_model(db, times, interval, timestep, do_plot = true)
-
-#
-# pollution = parse_pollution(api)
-# pollution = pollution[5513m .. 470513m, 531500m .. 1221500m, "pm2-5"]
-# ukclimateAE(pollution, size(total_pop), area, fill(true, total_pop), Lockdown(20days))
+run_model(db, 1day, interval, timestep)
+run_model(db, times, interval, timestep, save = true, do_plot = true)

src/EcoSISTEM.jl

@@ -64,13 +64,13 @@ include("Biodiversity/Landscape.jl")
 export GridLandscape, CachedGridLandscape
 
 include("Transitions/Transitions.jl")
-export TransitionList, create_transition_list, addtransition!
+export TransitionList, create_transition_list, addtransition!, getprob
 
 include("Epidemiology/MedianGenerator.jl")
 export MedianGenerator
 
 include("Epidemiology/data_utils.jl")
-export parse_hdf5, get_3d_km_grid_axis_array
+export parse_hdf5, get_3d_km_grid_axis_array, get_bng, get_en
 
 include("Epidemiology/EpiControl.jl")
 export NoControl, Lockdown
@@ -152,10 +152,10 @@ export BirthProcess, DeathProcess, GenerateSeed, AllDisperse, SeedDisperse,
 UpdateEnergy, UpdateEnvironment, update_environment!
 
 include("Transitions/EpiTransitions.jl")
-export ForceProduce, ForceDisperse, ViralLoad, Exposure, EnvExposure, Infection,
+export ForceProduce, ForceDisperse, ViralLoad, EnvViralLoad, Exposure, Infection,
 DevelopSymptoms, Hospitalise, DeathFromInfection,
 Recovery, SeedInfection, UpdateEpiEnvironment, update_epi_environment!,
-get_env
+get_env, deterministic_seed!
 
 include("Transitions/BiodiversityRun.jl")
 include("Transitions/EpiRun.jl")

src/Epidemiology/EpiEnv.jl

@@ -34,7 +34,7 @@ mutable struct GridEpiEnv{H, C, A} <: AbstractEpiEnv{H, C}
     ) where {H, C, A <: AbstractMatrix{Bool}}
         countsubcommunities(habitat) == length(names) ||
             error("Number of subcommunities must match subcommunity names")
-        if (size(habitat.matrix, 1), size(habitat.matrix, 2)) != size(active)
+        if (size(habitat, 1), size(habitat, 2)) != size(active)
             throw(DimensionMismatch(
                 "size(habitat.matrix)=$(size(habitat.matrix)) != " *
                 "size(active)=$(size(active))"
@@ -121,12 +121,34 @@ else one is created with all grid cells active.
     The simulation grid will be shrunk so that it tightly wraps the active values
 """
 function ukclimateAE(
-    climatearray::AxisArray,
+    climatearray::AxisArray{T, 2},
     dimension::Tuple{Int64, Int64},
     area::Unitful.Area{Float64},
     active::AbstractMatrix{Bool},
     control::C,
-) where C <: AbstractControl
+) where {C <: AbstractControl, T}
+    if typeof(first(climatearray)) <: Unitful.Temperature
+        climatearray .= uconvert.(K, climatearray)
+    end
+    area = uconvert(km^2, area)
+    gridsquaresize = sqrt(area / (dimension[1] * dimension[2]))
+
+    # Shrink to active region
+    # This doesn't change the gridsquaresize
+    climatearray = _shrink_to_active(climatearray, active)
+    active = _shrink_to_active(active, active)
+
+    hab = ContinuousHab(Array(climatearray), gridsquaresize, HabitatUpdate(NoChange, 0.0/s, Unitful.Dimensions{()}))
+    return GridEpiEnv{typeof(hab), typeof(control)}(hab, active, control)
+end
+
+function ukclimateAE(
+    climatearray::AxisArray{T, 3},
+    dimension::Tuple{Int64, Int64},
+    area::Unitful.Area{Float64},
+    active::AbstractMatrix{Bool},
+    control::C,
+) where {C <: AbstractControl, T}
     if typeof(first(climatearray)) <: Unitful.Temperature
         climatearray .= uconvert.(K, climatearray)
     end

src/Epidemiology/EpiPlots.jl

@@ -64,7 +64,7 @@ end
     seriescolor --> :heat
 
     subplot = 1
-    gridsize = (size(epi.abenv.habitat.matrix, 1), size(epi.abenv.habitat.matrix, 2))
+    gridsize = (size(epi.abenv.habitat, 1), size(epi.abenv.habitat, 2))
     for step in steps
         data = Float64.(reshape(abuns[idx, :, step], gridsize...))
         data[.!epi.abenv.active] .= NaN

src/Epidemiology/MedianGenerator.jl

@@ -33,3 +33,13 @@ function rand!(::MedianGenerator, dist::Distribution{Multivariate,S},
     container .= median(dist)
     return nothing
 end
+
+import StatsBase.sample
+function sample(::MedianGenerator, a::AbstractArray, wv::AbstractWeights, n::Int64)
+  order = sortperm(wv.values, rev =true)
+  return a[order[1:n]]
+end
+
+function sample(::MedianGenerator, a::AbstractArray, n::Int64)
+  return a[1:n]
+end