Skip to content

Commit

Permalink
Merge pull request #65 from boydorr/claireh93/birth-update
Browse files Browse the repository at this point in the history 
Plant demographics
  • Loading branch information
@claireh93
claireh93 committed Jul 14, 2021
2 parents 95a3a7c + 11fc0bf commit 0fd854e
Show file tree
Hide file tree
Showing 6 changed files with 53 additions and 47 deletions.
Binary file added examples/Africa.tif
View file
Binary file not shown.
51 changes: 13 additions & 38 deletions examples/Africa_run.jl
View file
Expand Up @@ -8,7 +8,7 @@ using Unitful.DefaultSymbols
using Distances
using StatsBase
using Plots
file = "../Chapter5/data/Africa.tif"
file = "Africa.tif"
africa = readfile(file, -25°, 50°, -35°, 40°)
active = Array{Bool, 2}(.!isnan.(africa'))

Expand Down Expand Up @@ -171,7 +171,7 @@ using Unitful
using Unitful.DefaultSymbols
using JLD
using Printf
file = "Documents/Chapter5/data/Africa.tif"
file = "Africa.tif"
africa = readfile(file, -25°, 50°, -35°, 40°)
active = Array{Bool, 2}(.!isnan.(africa'))
# Set up initial parameters for ecosystem
Expand Down Expand Up @@ -218,18 +218,6 @@ rel = Gauss{typeof(1.0K)}()
eco = Ecosystem(sppl, abenv, rel)
eco.abundances.matrix[50_000, :] .= 0

import EcoSISTEM.simulate!
function simulate!(eco::Ecosystem, times::Unitful.Time, timestep::Unitful.Time, cacheInterval::Unitful.Time, cacheFolder::String, scenario_name::String)
time_seq = 0s:timestep:times
counting = 0
for i in 1:length(time_seq)
update!(eco, timestep);
# Save cache of abundances
if mod(time_seq[i], cacheInterval) == 0year
JLD.save(joinpath(cacheFolder, scenario_name * (@sprintf "%02d.jld" uconvert(NoUnits,time_seq[i]/cacheInterval))), "abun", eco.abundances.matrix)
end
end
end

# Simulation Parameters
burnin = 100years; times = 100years; timestep = 1month; record_interval = 12months;
Expand All @@ -250,7 +238,7 @@ using Unitful.DefaultSymbols
using JLD
using Printf

file = "Documents/Chapter5/data/Africa.tif"
file = "Africa.tif"
africa = readfile(file, -25°, 50°, -35°, 40°)
active = Array{Bool, 2}(.!isnan.(africa'))
# Set up initial parameters for ecosystem
Expand Down Expand Up @@ -295,19 +283,6 @@ rel = Gauss{typeof(1.0K)}()
#Create ecosystem
eco = Ecosystem(sppl, abenv, rel)

import EcoSISTEM.simulate!
function simulate!(eco::Ecosystem, times::Unitful.Time, timestep::Unitful.Time, cacheInterval::Unitful.Time, cacheFolder::String, scenario_name::String)
time_seq = 0s:timestep:times
counting = 0
for i in 1:length(time_seq)
update!(eco, timestep);
# Save cache of abundances
if mod(time_seq[i], cacheInterval) == 0year
JLD.save(joinpath(cacheFolder, scenario_name * (@sprintf "%02d.jld" uconvert(NoUnits,time_seq[i]/cacheInterval))), "abun", eco.abundances.matrix)
end
end
end

# Simulation Parameters
burnin = 10years; times = 100years; timestep = 1month; record_interval = 12months;
lensim = length(0years:record_interval:times)
Expand All @@ -317,7 +292,7 @@ lensim = length(0years:record_interval:times)
using JLD
using Plots
using Diversity
abuns = load("examples/Biodiversity/Africa_run_coexist100.jld", "abun")
abuns = load("examples/Africa_run_coexist100.jld", "abun")
meta = Metacommunity(abuns)
div = norm_sub_alpha(meta, 0)
sumabuns = reshape(div[!, :diversity], 100, 100)
Expand All @@ -329,7 +304,7 @@ heatmap(sumabuns,
clim = (0, 50_000), margin = 0.5 * Plots.mm,
title = "A", titleloc = :left)

abuns = load("examples/Biodiversity/Africa_run50.jld", "abun")
abuns = load("examples/Africa_run50.jld", "abun")
meta = Metacommunity(abuns)
div = norm_sub_alpha(meta, 0)
sumabuns = reshape(div[!, :diversity], 100, 100)
Expand All @@ -341,7 +316,7 @@ heatmap!(sumabuns,
clim = (0, 50_000), right_margin = 2.0 * Plots.mm,
title = "B", titleloc = :left)

abuns = load("examples/Biodiversity/Africa_run100.jld", "abun")
abuns = load("examples/Africa_run100.jld", "abun")
meta = Metacommunity(abuns)
div = norm_sub_alpha(meta, 0)
sumabuns = reshape(div[!, :diversity], 100, 100)
Expand All @@ -353,16 +328,16 @@ heatmap!(sumabuns,
clim = (0, 50_000), right_margin = 2.0 * Plots.mm,
title = "C", titleloc = :left)


abuns = load("examples/Biodiversity/Africa_run50.jld", "abun")
using Diversity.Ecology
abuns = load("examples/Africa_run50.jld", "abun")
meta = Metacommunity(abuns)
div = norm_sub_rho(meta, 1.0)
sumabuns = reshape(div[!, :diversity], 100, 100)
diver = shannon(meta)
sumabuns = reshape(diver[!, :diversity], 100, 100)
heatmap!(sumabuns,
background_color = :lightblue,
background_color_outside=:white,
grid = false, color = :algae,
aspect_ratio = 1, subplot = 4,
right_margin = 2.0 * Plots.mm,
title = "D", titleloc = :left, clim = (0, 1))
Plots.pdf("examples/Biodiversity/Africa.pdf")
right_margin = 2.0 * Plots.mm,
title = "D", titleloc = :left, clim = (0, 10))
Plots.pdf("examples/Africa.pdf")
16 changes: 8 additions & 8 deletions src/Generate.jl
View file
Expand Up @@ -67,17 +67,17 @@ function update!(eco::Ecosystem, timestep::Unitful.Time)
# Check if grid cell currently active
if eco.abenv.active[x, y] && (eco.cache.totalE[i, 1] > 0)
# Calculate effective rates
birthprob = params.birth[j] * timestep * adjusted_birth
deathprob = params.death[j] * timestep * adjusted_death
birthrate = params.birth[j] * timestep * adjusted_birth
deathparam = params.death[j] * timestep * adjusted_death

# Put probabilities into 0 - 1
newbirthprob = 1.0 - exp(-birthprob)
newdeathprob = 1.0 - exp(-deathprob)
# Turn deathparam into probability and cancel units of birthrate
birthrate += 0.0
deathprob = 1.0 - exp(-deathparam)

(newbirthprob >= 0) & (newdeathprob >= 0) || error("Birth: $newbirthprob \n Death: $newdeathprob \n \n i: $i \n j: $j")
(birthrate >= 0) & (deathprob >= 0) || error("Birth: $birthrate \n Death: $deathprob \n \n i: $i \n j: $j")
# Calculate how many births and deaths
births = rand(rng, Poisson(eco.abundances.matrix[j, i] * newbirthprob))
deaths = rand(rng, Binomial(eco.abundances.matrix[j, i], newdeathprob))
births = rand(rng, Poisson(eco.abundances.matrix[j, i] * birthrate))
deaths = rand(rng, Binomial(eco.abundances.matrix[j, i], deathprob))

# Update population
eco.abundances.matrix[j, i] += (births - deaths)
Expand Down
28 changes: 28 additions & 0 deletions src/Helper.jl
View file
Expand Up @@ -2,6 +2,7 @@ using Unitful
using Unitful.DefaultSymbols
using Diversity
using EcoSISTEM.Units
using Printf
import Diversity.Gamma

"""
Expand All @@ -17,6 +18,13 @@ function simulate!(eco::AbstractEcosystem, duration::Unitful.Time, timestep::Uni
update!(eco, timestep)
end
end

"""
simulate!(cache::CachedEcosystem, srt::Unitful.Time, timestep::Unitful.Time)
Function to run a cached ecosystem, `cache` at a specified timepoint, `srt`,
for a particular timestep, 'timestep'.
"""
function simulate!(cache::CachedEcosystem, srt::Unitful.Time, timestep::Unitful.Time)
eco = Ecosystem{typeof(cache.abenv), typeof(cache.spplist),
typeof(cache.relationship)}(copy(cache.abundances.matrix[srt]),
Expand All @@ -26,6 +34,26 @@ function simulate!(cache::CachedEcosystem, srt::Unitful.Time, timestep::Unitful.
cache.abundances.matrix[srt + timestep] = eco.abundances
end

"""
simulate!(eco::Ecosystem, times::Unitful.Time, timestep::Unitful.Time, cacheInterval::Unitful.Time,
cacheFolder::String, scenario_name::String)
Function to run an ecosystem, `eco` for specified length of times, `duration`,
for a particular timestep, 'timestep'. A cache interval and folder/file name
are specified for saving output.
"""
function simulate!(eco::Ecosystem, times::Unitful.Time, timestep::Unitful.Time, cacheInterval::Unitful.Time,
cacheFolder::String, scenario_name::String)
time_seq = zero(times):timestep:times
for i in 1:length(time_seq)
update!(eco, timestep);
# Save cache of abundances
if mod(time_seq[i], cacheInterval) == zero(time_seq[i])
JLD.save(joinpath(cacheFolder, scenario_name * (@sprintf "%02d.jld" uconvert(NoUnits,time_seq[i]/cacheInterval))), "abun", eco.abundances.matrix)
end
end
end

function generate_storage(eco::Ecosystem, times::Int64, reps::Int64)
numSpecies = length(eco.spplist.abun)
gridSize = _countsubcommunities(eco.abenv.habitat)
Expand Down
2 changes: 1 addition & 1 deletion test/canonical/test_biodiversity.jl
View file
Expand Up @@ -378,7 +378,7 @@ for i in eachindex(species)
local survival = 0.1
local boost = 1.0

local size_mean = rand(Normal(1.0, 0.5), numSpecies) .* m^2
local size_mean = rand(Normal(1.0, 0.05), numSpecies) .* m^2
# Set up how much energy each species consumes
local energy_vec1 = SolarRequirement(abs.(req[1] .* size_mean))
local energy_vec2 = WaterRequirement(abs.(req[2] .* size_mean))
Expand Down
3 changes: 3 additions & 0 deletions test/test_Helper.jl
View file
Expand Up @@ -24,6 +24,9 @@ end
@test size(abun) == (size(eco.abundances.matrix, 1), size(eco.abundances.matrix, 2), lensim, 1)
@test_nowarn simulate!(eco, burnin, timestep)
@test_nowarn simulate_record!(abun, eco, times, interval, timestep)
isdir("data") || mkdir("data")
@test_nowarn simulate!(eco, times, interval, timestep, "data", "testrun")
rm("data", recursive = true)
end
@testset "scenarios" begin
# Run with scenarios
Expand Down

0 comments on commit 0fd854e

Please sign in to comment.