Extend extinct_species()

[ismael-lajaaiti]

TODO

• rename get_extinct_species to get_extinction_events as the function return events (which species has gone extinct and when) and not only extinct species list
• create an get_extinct_species(extinction_events, t) that given the dictionary of extinction_events return the list of extinct species at the time step t
• create get_alive_species, the sister of get_extinction_species that return alive species instead of dead species, as in some cases in can be more convenient in this way

[alaindanet]

Thanks a lot!
Do you think that it is possible that those functions returns the timestep number in plus that t? It would be more convenient to extract the solution object.
And I do not know if it is possible to have both the index of species position in the solution object (which I call idxs like in the DifferentialEquations.jl package) and the species names. It is also more handy also to extract the solution object.

[ismael-lajaaiti]

Do you think that it is possible that those functions returns the timestep number in plus that t? It would be more convenient to extract the solution object.

You mean the index of the corresponding time step?

And I do not know if it is possible to have both the index of species position in the solution object (which I call idxs like in the DifferentialEquations.jl package) and the species names. It is also more handy also to extract the solution object.

Mmmh... this will be redundant information, because then one element of the dictionary would be something like (name = :turtle, idx = 3) => (t = 120, idx = 30). If this is what you're suggesting, this would be too heavy IMO. But we could write a functions to easily access species names from their index e.g. species_name(idx, network)

[alaindanet]

You mean the index of the corresponding time step?

Yes 🙏

Mmmh... this will be redundant information, because then one element of the dictionary would be something like (name = :turtle, idx = 3) => (t = 120, idx = 30). If this is what you're suggesting, this would be too heavy IMO. But we could write a functions to easily access species names from their index e.g. species_name(idx, network)

Excellent! species_name(idx, network) would be perfect for me as well as species_idxs(name, network)

[iago-lito]

I sense that this pair of functions will be increasingly useful as the simulated variables will be soon extended to contain not only species biomasses (see #113). Althought I would prefer they be named species_names(indexes, network) and species_indexes(names, network) instead, according to Julia's naming style 0:)

• conciseness is valued, but avoid abbreviation (indexin rather than indxin) as it becomes difficult to remember whether and how particular words are abbreviated.

[ismael-lajaaiti]

About the index of the time step, I am also thinking of creating a sweet function timestep_index(t, solution), that you could use to extract extinct species, but I am not sure this is necessary to always include the time step index in the dictionary of extinction events. Would it also be good for you this way?

[iago-lito]

I'm sorry I don't understand what you mean here. Do you expect that we could extract extinct/alive species from a solution at a given time dy doing things like:

alive_then = alive_species(timestep_index(100.0, solution))

or something?

[iago-lito]

Don't forget to address this weakness here within the scope of this PR.

[alaindanet]

Currently I use this framework to determine alive species:

function living_species(
    solution::SciMLBase.ODESolution;
    threshold = 0,
    idxs = nothing,
    kwargs...,
)

    measure_on = extract_last_timesteps(solution; idxs, kwargs...)

    alive_sp = living_species(measure_on; threshold)

    sp = get_parameters(solution).network.species

    tmp_idxs = process_idxs(solution, idxs)
    idxs = tmp_idxs[alive_sp]
    species = sp[idxs]

    (; species, idxs)
end

living_species(mat::AbstractMatrix; threshold = 0) =
    findall(>(threshold), biomass(mat).species)

living_species(n::AbstractVector; threshold = 0) = findall(>(threshold), n)

A living species at time t is a species that has a biomass > 0.
A living species in a solution object is a species that has a biomass >0 over the last timesteps.

On top of this, I warn user if their simulation extraction over the last timesteps contains an extinction events:

function extract_last_timesteps(solution; idxs = nothing, quiet = false, kwargs...)
   ....
    quiet || check_last_extinction(solution; idxs = idxs, last = last)
    ...
end

To do this, I need to get the timestep indexes of extinction events. To do this, I rely on the fact that a species became extinct when its biomass reached exactly 0:

function get_extinction_timesteps(solution; idxs = nothing)
    idxs = process_idxs(solution, idxs;)
    sp = get_parameters(solution).network.species[idxs]

    ext_t = findfirst.(isequal(0), eachrow(solution[idxs, :]))
    extinct = ext_t .!== nothing

    (
        species = sp[extinct],
        idxs = idxs[extinct],
        extinction_timestep = something.(ext_t[extinct]),
    )
end

function get_extinction_timesteps(m::AbstractVector; threshold = 0)
    findfirst(x -> x <= threshold, m)
end

It is there that I would need something like timestep_index(t, solution) to extract directly the extinction timesteps @ismael-lajaaiti 🙏

[iago-lito]

Blocked on #137.