Sanity checks for simulations: disconnected species, extinct species & disconnected graphs #151
Description
The goal would be to warn users that the ecological network simulations can have properties that can create problems for interpreting the results or computing output metrics.
Those problematic features are:
-
Disconnected species:
- a consumer whose resources/preys are all dead (i.e. their biomass is 0)
- a producer whose the consumers/predators are all dead (i.e. their biomass is 0)
-
Disconnected graphs:
- Two disconnected trophic chains
-
Tasks:
- Functions for disconnected species
- Function for disconnected graphs
- Message/Warning implemented in
simulate - Write a workaround doc for user to simulate dynamics until there are no more extinct or disconnected species.
- Expose processing metrics (CV, trophic levels, ...)
A set of function for disconnected species:
"""
kill_disconnected_species(A; alive_species = nothing, bm = nothing)
Returns a vector of species biomass, all the disconnected and extinct species being set to a
biomass of 0.
# Arguments:
- `A`: Adjacency matrix
- `alive_species`: vector of alive species indices. See [`living_species`](@ref)
- `bm`: vector of species biomass
# Examples
ti = [
0 0 0 0;
0 0 0 0;
1 0 0 0;
0 1 0 0
]
bm = [1, 1, 3, 4]
alive = [1, 3, 4]
kill_disconnected_species(ti, alive_species = [1, 2, 3, 4], bm = bm) == [1, 1, 3, 4]
kill_disconnected_species(ti, alive_species = [], bm = bm) == [0, 0, 0, 0]
kill_disconnected_species(ti, alive_species = [1, 3, 4], bm = bm) == [1, 0, 3, 0]
kill_disconnected_species(ti, alive_species = [1, 2, 4], bm = bm) == [0, 1, 0, 4]
kill_disconnected_species(ti, alive_species = [1, 2, 4], bm = bm) == [0, 1, 0, 4]
# Filter species biomass and then adjacency matrix
## Case 1
to = kill_disconnected_species(ti, alive_species = [1, 2, 3, 4], bm = bm)
mask = to .== 0
A = ti[mask, mask]
new_bm = bm[mask]
dim(A) == length(new_bm)
## Case 2
to = kill_disconnected_species(ti, alive_species = [1, 2, 4], bm = bm)
mask = to .== 0
A = ti[mask, mask]
new_bm = bm[mask]
dim(A) == length(new_bm)
"""
function kill_disconnected_species(A; alive_species = nothing, bm = nothing)
bm = deepcopy(bm)
disconnected_alive_species = check_disconnected_species(A, alive_species)
in_disconnected = in(disconnected_alive_species)
mask_alive_disconnected = in_disconnected.(alive_species)
alive_to_keep = alive_species[.!mask_alive_disconnected]
bm_to_set_to_zero = 1:length(bm) .∉ [alive_to_keep]
bm[bm_to_set_to_zero] .= 0
bm
end
"""
remove_disconnected_species(A, alive_species)
# Examples
ti = [
0 0 0 0;
0 0 0 0;
1 0 0 0;
0 1 0 0
]
remove_disconnected_species(ti, [1, 2, 3])
"""
function remove_disconnected_species(A, alive_species)
disconnected_alive_species = check_disconnected_species(A, alive_species)
in_disconnected = in(disconnected_alive_species)
mask_alive_disconnected = in_disconnected.(alive_species)
to_keep = .!(mask_alive_disconnected)
living_A = A[alive_species, alive_species]
living_A[to_keep, to_keep]
end
"""
check_disconnected_species(A, alive_species)
# Examples
ti = [
0 0 0 0;
0 0 0 0;
1 0 0 0;
0 1 0 0
]
check_disconnected_species(ti, [1, 2, 3]) == [2]
check_disconnected_species(ti, [1, 2, 4]) == [1]
check_disconnected_species(ti, [2, 3, 4]) == [3]
check_disconnected_species(ti, [1, 3, 4]) == [4]
check_disconnected_species(ti, [1, 2, 3, 4]) == []
"""
function check_disconnected_species(A, alive_species; verbose = false)
# Get binary matrix
A = A .> 0
living_A = A[alive_species, alive_species]
cons = sum.(eachrow(A)) .!= 0
prod = sum.(eachrow(A)) .== 0
alive_cons = cons[alive_species]
alive_prod = prod[alive_species]
species_with_no_pred = sum.(eachcol(living_A)) .== 0
species_with_no_prey = sum.(eachrow(living_A)) .== 0
disconnected_prod = alive_prod .&& species_with_no_pred
disconnected_cons = alive_cons .&& species_with_no_prey
if sum([disconnected_prod; disconnected_cons]) > 0 & verbose
println("There are $(sum(disconnected_prod)) disconnected producers
and $(sum(disconnected_cons)) consumers.")
end
alive_species[ disconnected_prod .|| disconnected_cons ]
endThe goal is that simulate produces a useful message when there is a disconnected species. For example, in that case when there is a disconnected consumer:
# Disconnected consumer that goes extinct, all fine
A = [0 0 0 0; 0 0 0 0; 1 0 0 0; 0 1 0 0]
foodweb = FoodWeb(A);
params = ModelParameters(foodweb);
B0 = [0.5, 0, .5, .5];
julia> sol = simulate(params, B0);
┌ Info: Species [4] went extinct at time t = 38.37898034897611.
└ 2 out of 4 species are extinct.
julia> check_disconnected_species(get_parameters(sol).network.A, living_species(sol).idxs,
verbose = true
)
Int64[]
julia> kill_disconnected_species(get_parameters(sol).network.A;
alive_species = living_species(sol).idxs,
bm = biomass(sol).species
)
4-element Vector{Float64}:
0.19199402785494718
0.0
0.33605893664276615
0.0
# Disconnected producers
B0 = [0.5, 0.5, 0, .5];
sol = simulate(params, B0);
check_disconnected_species(get_parameters(sol).network.A, living_species(sol).idxs,
verbose = true
)
kill_disconnected_species(get_parameters(sol).network.A;
alive_species = living_species(sol).idxs,
bm = biomass(sol).species
)
julia> check_disconnected_species(get_parameters(sol).network.A, living_species(sol).idxs,
verbose = true
)
There are 1 disconnected producers
and 0 consumers.
1-element Vector{Int64}:
1
julia> kill_disconnected_species(get_parameters(sol).network.A;
alive_species = living_species(sol).idxs,
bm = biomass(sol).species
)
There are 1 disconnected producers
and 0 consumers.
4-element Vector{Float64}:
0.0
0.19194243058968705
0.0
0.36856095125952454
On top of sanity checks, we can also provide a snippet in the doc to enable users to run the simulations until there are no more extinct species or disconnected species:
global ti = false
global i = 1
global remove_disconnected = false
global starting_bm = [0.5, 0.5, 0, .5];
#starting_bm = rand(richness(foodweb))
while ti != true
println("Iteration = $i.")
output = simulate(params, starting_bm)
# Retrieve network and alive species
old_A = get_parameters(output).network.A
alive_species = living_species(output).idxs
# Check
disconnected_species = check_disconnected_species(old_A, alive_species)
ti = length(disconnected_species) == 0
if ti == true
println("No disconnected species, all fine.")
return output
end
# Build the vector of biomasses
biomass_vector = biomass(output).species
killed_species = kill_disconnected_species(old_A;
alive_species = alive_species,
bm = biomass_vector)
mask_sp_to_keep = killed_species .!= 0.0
idxs = 1:size(old_A, 1)
idxs_to_keep = idxs[mask_sp_to_keep]
# Rebuilding network or set disconnected species to 0
if remove_disconnected == true
A = old_A[mask_sp_to_keep, mask_sp_to_keep]
starting_bm = biomass_vector[mask_sp_to_keep]
foodweb = FoodWeb(A);
params = ModelParameters(foodweb);
println("Rebuilding model without disconnected species.")
else
starting_bm = killed_species
A = A
println("Keep the same model without disconnected species.")
end
i = i + 1
end
Metadata
Metadata
Assignees
Labels
No labels