Introduce Topology to check biomass graph disconnections.

[iago-lito]

Start work on addressing #151.

[STATUS] The feature is complete and has landed on dev. Before merging into main, though:

• Try it out, making sure bikeshedding is over.
• Write basic documentation (not much technical about code, but rather high-level about science) [comment].
• Write a first CHANGELOG.md entry to land this into main as our v0.2.1 :)

[iago-lito]

So, this PR addresses #151 in a very deep way by creating the new Topology type.
Values of this type represent the ecological network under a strict topological perspective (neighbourhood and connections) according to the model sketched in #151 (comment).

The reason to put so much effort into Topology semantics and testing is that I don't intend to only use it to solve #151. I expect this type to become a solid candidate to tackle #141 when it'll be time: all the internal representation of what users mean by an "EcologicalNetwork" can be based on such a Topology.

Review welcome @alaindanet @ismael-lajaaiti ;) There is a lot of code so I won't mind if you don't thoroughly scan all of it of course. But if you want to help: please test the feature to check whether it's comfy and whether it corresponds to something that would at least solve #151 :)

Entrypoint into the feature:

• (check out the disconnections branch)
• You can try model.topology to retrieve a topology value, then test the methods in src/topology.jl.
• After simulation you can also try get_topology(solution) to retrieve a similar-looking topology, but with the extinct species removed. If simulation has resulted in the introduction of isolated producers / starving consumers / disconnected components, then you should have seen an @info message about that and you can query it with the same set of methods ;)

Feedback welcome! (here or on the channel ;)

[ismael-lajaaiti]

I started to tried your feature but I quickly encountered an error. Here is my code

using Distributions
using EcologicalNetworksDynamics
using Random
Random.seed!(1234) # For reproducibility.

foodweb = Foodweb(:niche, S=50, C=0.2)
m = default_model(foodweb)
g = get_topology(m)
g == m.topology # Sanity check. Expected to be true.

simulate(m, rand(Uniform(0.1, 1), m.S), 1_000)

which gives me the following error

ERROR: KeyError: key EcologicalNetworksDynamics.Topologies.Abs(3) not found
Stacktrace:
 [1] pop!
   @ ./dict.jl:604 [inlined]
 [2] pop!(s::Set{EcologicalNetworksDynamics.Topologies.Abs}, x::EcologicalNetworksDynamics.Topologies.Abs)
   @ Base ./set.jl:104
 [3] starving_consumers(m::<inner parms>, g::EcologicalNetworksDynamics.Topologies.Topology)
   @ EcologicalNetworksDynamics ~/Documents/Projects/Julia/edn/src/topology.jl:117
 [4] show_degenerated_biomass_graph_properties(model::<inner parms>, biomass::Vector{Float64}, arg::Symbol)
   @ EcologicalNetworksDynamics ~/Documents/Projects/Julia/edn/src/simulate.jl:141
 [5] _simulate(model::<inner parms>, u0::Vector{Float64}, tmax::Int64; kwargs::@Kwargs{model::Model})
   @ EcologicalNetworksDynamics ~/Documents/Projects/Julia/edn/src/simulate.jl:44
 [6] _simulate(::Model, ::Vector{Float64}, ::Vararg{Any}; kwargs::@Kwargs{model::Model})
   @ EcologicalNetworksDynamics ~/Documents/Projects/Julia/edn/src/Framework/method_macro.jl:294
 [7] _simulate
   @ ~/Documents/Projects/Julia/edn/src/Framework/method_macro.jl:288 [inlined]
 [8] simulate(model::Model, u0::Vector{Float64}, tmax::Int64)
   @ EcologicalNetworksDynamics ~/Documents/Projects/Julia/edn/src/simulate.jl:74
 [9] top-level scope
   @ REPL[59]:1

Did I do something wrong?
I wanted to run a simulation and remove some (disconnected) species or to see how I could create a new model from a single disconnected components, or simply from the network of alive species.

Otherwise I was wondering if it would be possible, to have a compact view of the topology for large networks? Because currently it is a bit overhelming (see screenshot below). But this is a very cosmetic feature and in any way shouldn't be (y)our priority.

[iago-lito]

Wops! I can reproduce your error with a simple graph indeed:

m = default_model(Foodweb([:a => [:b, :c], :b => :c]))
g = m.topology
starving_consumers(m, g)

This is a bug in the graph visit implementation. I'm fixing it right now. Thank you for reporting.

Otherwise I was wondering if it would be possible, to have a compact view of the topology for large networks? Because currently it is a bit overhelming
Wow, of course! I already have a few utils to elide long vectors.. and maybe I won't display the whole adjacency list if it's long.

BTW are you happy with the default :s1, :s2, .., s:2048 naming, or should we produce shorter identifiers like :a, :b, .. :z, :aa, :ab etc.? Or random fixed-width identifiers like :ux, :bw, :pa, :oj.. ?

[ismael-lajaaiti]

This is a bug in the graph visit implementation. I'm fixing it right now. Thank you for reporting.

Great, thank you for fixing it. 😉
Let me know when it's done, so I can continue experimenting the new feature. 🤓

BTW are you happy with the default :s1, :s2, .., s:2048 naming, or should we produce shorter identifiers like :a, :b, .. :z, :aa, :ab etc.? Or random fixed-width identifiers like :ux, :bw, :pa, :oj.. ?

Yes, I like it. In particular, because with this default it is easy to distinguish between species and nutrients.

[alaindanet]

I have just tried the feature too, I have an error checking topology with a solution object:

using EcologicalNetworksDynamics
using Random
Random.seed!(1234) # For reproducibility.

foodweb = Foodweb(:niche, S=50, C=0.2);
m = default_model(foodweb);
g = get_topology(m);
g == m.topology # Sanity check. Expected to be true.

sol = simulate(m, rand(m.S), 1_000);

julia> get_topology(sol)
ERROR: MethodError: no method matching get_topology(::SciMLBase.ODESolution{…})

Closest candidates are:
  get_topology(::<inner parms>)
   @ EcologicalNetworksDynamics ~/Documents/post-these/sheffield/befwm2_dev/BEFWM2/src/expose_data.jl:602
  get_topology(::Model, Any...; kwargs...)
   @ EcologicalNetworksDynamics ~/Documents/post-these/sheffield/befwm2_dev/BEFWM2/src/Framework/method_macro.jl:288

Stacktrace:
 [1] top-level scope
   @ REPL[35]:1
Some type information was truncated. Use `show(err)` to see complete types.

[ismael-lajaaiti]

I've run a few tests (see the end of this message) and the methods seem to work well. But here are a few remarks:

• Related to @alaindanet's issue it appears that the function to retrieve the community topology is now named topology whereas it was get_topology before. Note that when the argument is a model, both names appear to work. That might be confusing for the users.
• For the method of topology to retrieve the topology from a solution at a given date, I would be in favour of enforcing the user to type topology(solution; date = 42) instead of the current version topology(solution, 42). By doing so, I think it would be easier for a user to understand the difference between the method topology(s::Solution) and topology(s::Solution; date::Number)
• Is there a method to convert a topology object into an adjacency matrix. This should be possible, right?

My tests

using CairoMakie
using DataFrames
using Distributions
using EcologicalNetworksDynamics
using Random
Random.seed!(1234) # For reproducibility.

foodweb = Foodweb(:niche, S=50, C=0.2)
m = default_model(foodweb)
g = topology(m)
g == m.topology # Sanity check. Expected to be true.
sol = simulate(m, rand(Uniform(0.1, 1), m.S), 1_000)

extinctions(sol)
# 49 => 58.8748
# 35 => 511.178
# 16 => 640.776

# Some test to ensure that everything is working as expected.
# Before the first extinction, the topology should be the same as the original one.
all([topology(m) == topology(sol, t) for t in 0:58])
# After extinction, the edges of species 49 should disappear.
topology(sol, 58).outgoing[49]
topology(sol, 59).outgoing[49]

[iago-lito]

Hi @alaindanet @ismael-lajaaiti, and thank you for feedback :)

• Related to @alaindanet's issue it appears that the function to retrieve the community topology is now named topology whereas it was get_topology before. Note that when the argument is a model, both names appear to work. That might be confusing for the users.

Agreed: that is a source of confusion. That method intersects two naming logics:

• The model properties names: get_A(m) == m.A, get_topology(m) == m.topology etc.
• The solution information: extinctions(sol), model(sol), topology(sol) etc.

Maybe I should switch the whole solution methods suite back to get_extinctions(sol), get_model(sol), get_topology(sol) etc. ? Or would it be too annoying?

---

I would be in favour of enforcing the user to type topology(solution; date = 42)

Yupe, why not, this is just a matter of making it a keyword parameter instead :)

---

Is there a method to convert a topology object into an adjacency matrix. This should be possible, right?

Possible, although not exactly straightforward: the topology abstracts over all possible nodes/edges compartments, so you would have to specify which "adjacency" you actually need: species -> trophic -> species vs. producers -> trophic -> nutrients vs. species -> facilitation -> species etc.

That being said, I can definitely expose something like:

get_adjacency(g::Topology, source_compartment, edge_compartment, target_compartment)::SparseMatrix{Bool}

and then maybe wrap the few typical ones.. ?

get_foodweb_adjacency(g::Topology) = get_adjacency(g, :species, :trophic, :species)
get_nutrients_adjacency(g::Topology) = get_adjacency(g, :species, :trophic, :nutrients) # (damn, this would include consumers)
get_facilitation_adjacency(g::Topology) = get_adjacency(g, :species, :facilitation, :species)
# ...

.. unless maybe just leave get_adjacency exposed and let user do their own business ?

Note that, if you know the extinct species, then removing the corresponding lines + columns from model.A would yield the same result as get_adjacency(g, :species, :trophic, :species).

---

# After extinction, the edges of species 49 should disappear.
topology(sol, 58).outgoing[49]
topology(sol, 59).outgoing[49]

Wops! The Topology fields are unexposed (=private) ^ ^" This works but it may not or it may break in the future without warning. Maybe you are requesting some more basic query primitives like is_edge(g::Topology, :species, 58, :trophic, :species, 49) ? Eww, that looks not nice X) But all this information is needed to query one single edge :
Maybe get_adjacency(top, :species, :trophic, :species)[58, 49] would do the trick instead?

[ismael-lajaaiti]

Maybe I should switch the whole solution methods suite back to get_extinctions(sol), get_model(sol), get_topology(sol) etc. ? Or would it be too annoying?

Sounds good to me.

Possible, although not exactly straightforward: the topology abstracts over all possible nodes/edges compartments, so you would have to specify which "adjacency" you actually need: species -> trophic -> species vs. producers -> trophic -> nutrients vs. species -> facilitation -> species etc.
That being said, I can definitely expose something like:

get_adjacency(g::Topology, source_compartment, edge_compartment, target_compartment)::SparseMatrix{Bool}

and then maybe wrap the few typical ones.. ?

get_foodweb_adjacency(g::Topology) = get_adjacency(g, :species, :trophic, :species)
get_nutrients_adjacency(g::Topology) = get_adjacency(g, :species, :trophic, :nutrients) # (damn, this would include >consumers)
get_facilitation_adjacency(g::Topology) = get_adjacency(g, :species, :facilitation, :species)
# ...

.. unless maybe just leave get_adjacency exposed and let user do their own business ?

Mmh.. I was suspecting that this was not as easy as I thought. IMO, it would be great to leave both options, that is, the 'raw' function and its wrappers.

Wops! The Topology fields are unexposed (=private) ^ ^" This works but it may not or it may break in the future without warning. Maybe you are requesting some more basic query primitives like is_edge(g::Topology, :species, 58, :trophic, :species, 49) ? Eww, that looks not nice X) But all this information is needed to query one single edge :
Maybe get_adjacency(top, :species, :trophic, :species)[58, 49] would do the trick instead?

OK... Then how can we manipulate a Topology object? Ideally, it would be nice to manipulate it like an adjacency list, as it is printed like so. Would it be possible, for example if a topology is stored in a variable a, to allow something like a[49] to access all outgoing link from species 49?

[iago-lito]

Mmh.. I was suspecting that this was not as easy as I thought. IMO, it would be great to leave both options, that is, the 'raw' function and its wrappers.

Trying to list the "wrappers": I can think of get_foodweb_adjacency. The problem with get_nutrients_adjacency is that the topology alone does not know the "producers" so it would either yield a matrix that is larger than expected or it would need an extra information like get_nutrient_adjacency(g, producers_indices) or get_nutrients_adjacency(g, model).. which again lets the user free to input inconsistent arguments. get_<nti>_adjacency() seem okay as far as I can tell, but they would be redundant with each other.. so maybe just get_nti_adjacency(g, :facilitation) etc.?

In the end, maybe this calls for too much API surface X) What about reducing all that with some simpler:

# Lib.
get_adjacency(g, source, edge, target) = ...
get_species_adjacency(g, edge_type) = get_adjacency(g, :species, edge_type, :species)
get_foodweb_adjacency(g) = get_species_adjacency(g, :trophic)

# User.
A = get_foodweb_adjacency(g)
Af = get_species_adjacency(g, :facilitation)
Ar = get_species_adjacency(g, :refuge)
...
An = get_adjacency(g, :species, :trophic, :nutrients)[model.producers_mask, :]

.. for now instead ?

---

OK... Then how can we manipulate a Topology object? Ideally, it would be nice to manipulate it like an adjacency list, as it is printed like so. Would it be possible, for example if a topology is stored in a variable a, to allow something like a[49] to access all outgoing link from species 49?

Yay, well it depends on what you mean by "manipulate" ^ ^" Maybe you can get a sense of the inherent complexity of the
value by noticing that a simple expression like a[49] does not specify:

• that 49 refers to a species node index (and not a nutrient or a spatial patch)
• that you expect :trophic links (not :facilitation or :migrations).
• that you expect "outgoing" links (and not incoming ones).
• that you expect only links to :species targets (and not nutrients targets).

If you want to mean all of this implicitly, then maybe you need to work with A = get_adjacency(a, :species, :trophic, :species) instead of a directly. a contains much more information, and is open to even more in the future.

As a reminder: the introduction of these Topology values serves two purposes:

1. Feature Sanity checks for simulations: disconnected species, extinct species & disconnected graphs #151: the Topology values correctly represent the network. They can be:
   • Created consistently with a model.
   • Pruned from extinct species while remembering their indices.
   • Queried for disconnected components / starving consumers / isolated producers.
2. Become a candidate solution to Internals lack flexibility and could perform better. #141. If accepted, topology values may become completely integrated into the internals and become the fundamental basis of Model itself. Unfortunately, Topology and Model are two distinct types for now :\

So, yeah, in the context of this PR (1), there is not much you can do with topologies yet except for detecting degenerated networks. In the future though, I would be glad to make them easy to query/split/export/whatever, but this is maybe out of scope for #152?..

.. unless you can write down the things you would like to do with them, and they would be easy to feature with small primitives like get_adjacency(topology, s, e, t) in the context of this PR?

[ismael-lajaaiti]

In the end, maybe this calls for too much API surface X) What about reducing all that with some simpler:

# Lib.
get_adjacency(g, source, edge, target) = ...
get_species_adjacency(g, edge_type) = get_adjacency(g, :species, edge_type, :species)
get_foodweb_adjacency(g) = get_species_adjacency(g, :trophic)

# User.
A = get_foodweb_adjacency(g)
Af = get_species_adjacency(g, :facilitation)
Ar = get_species_adjacency(g, :refuge)
...
An = get_adjacency(g, :species, :trophic, :nutrients)[model.producers_mask, :]

.. for now instead ?

Yes, let's go for this for now.

So, yeah, in the context of this PR (1), there is not much you can do with topologies yet except for detecting degenerated networks. In the future though, I would be glad to make them easy to query/split/export/whatever, but this is maybe out of scope for #152?..

.. unless you can write down the things you would like to do with them, and they would be easy to feature with small primitives like get_adjacency(topology, s, e, t) in the context of this PR?

OK, I forgot that you primarily designed the Topology object to solve #151. So, I agree that this is out of the scope of this PR, and for now, it is enough because what I wanted to with a Topology object can be done with the (non-)trophic adjacency matrices. So forgot what I wrote above ahah.

[iago-lito]

Hi @ismael-lajaaiti @alaindanet. I think the last few commits address everything discussed so far here. Can you check whether it's okay? What do you think ? :)

[ismael-lajaaiti]

Everything seems to work good on my side. Thank you @iago-lito !
My only remark would be that we can improve the error message when requesting an adjacency matrix of interaction type that is not present within the model. For example, if I build a model with trophic and refuge links, and that I request the adjacency matrix of facilitation links here is what I get

ERROR: ArgumentError: Invalid edge type label: :facilitation. Valid labels are :refuge and :trophic.

I'm thinking of two things here:

1. Specify in the error message that the model, or the community network, doesn't contain facilitation links and therefore the corresponding adjacency matrix cannot be given. Otherwise, the message is a bit obscure.
2. Why not return an empty adjacency matrix, instead of raising an error?

For next, I also think we should discuss these points, let me know if they deserve their own issue:

1. We could update the method when generating a network with the niche or cascade model so it produces a matrix with no starving consumers or isolated producers. These could be arguments that could be turned on or off, as for the no_cycle argument. Because I think it makes little sense to generate an adjacency matrix where a consumer is starving. Even though one could argue that those species are going to be filtered during the assembly process anyway (that is, they are going extinct during the simulation), but as @alaindanet pointed this can take some time, so maybe better to prune them as soon as possible. Anyway, I think it's good to have this option open.
2. We (I) should also work on expanding the documentation (and especially the online one) as we add new feature as this one. I think this is very important as it is the surface of the package that is exposed to users.

[iago-lito]

Thank you for feedback @ismael-lajaaiti :)

we can improve the error message when requesting an adjacency matrix of interaction type that is not present within the model

ERROR: ArgumentError: Invalid edge type label: :facilitation. Valid labels are :refuge and :trophic.

I'm not sure we can do much better because the list of possible "interaction types" is open. Today, neither the topology or the model values know about all possible interaction types, and so there is no (simple) way they could make a difference between:

• :facilitation: a valid interaction name that just happen not to exist within this model.
• :rfeuge: incorrect spelling of an interaction that happens to exist within this model.

Thus the error message you are seing now. Maybe I can still clarify a little though:

Invalid edge type label: :facilitation. Valid labels in this model are :refuge and :trophic. ..?

Why not return an empty adjacency matrix, instead of raising an error?

For the same reason: in principle, anyone could add the interaction type they want to the system by creating their own custom components (although it's (very) not ready yet). So given input like :rfeuge, there is no way the system could decide whether to error (incorrect name) or to return an empty matrix (valid name), because the list of names is open. So it needs to always take the same decision in both cases. My opinion now is that issuing an error with :facilitation if there is no facilitation in the model is ok, but returning an empty matrix with :rfeuge if it is a typo is absolutely unacceptable ^ ^"

I'm curious though: what would you use this empty matrix for? How come your code needs to request facilitation links in a model without a Facilitation component?

---

update the method when generating a network with the niche or cascade model so it produces a matrix with no starving consumers or isolated producers

Oh, I thought these were guaranteed not to generate such degenerated networks. If they aren't, then of course we could feature that. The algorithm would just be bruteforce like the one to avoid cycles, right? Draw random networks until we find one without this kind of species.
Would you be okay to leave this for another github issue? I am going to screen through every component again in the context of #139. I could maybe pick this new feature up when I get to rewriting the Foodweb blueprints :)

---

We (I) should also work on expanding the documentation

Of course :) My first (modest) contribution was this comment.
It was good not to dive into the documentation too soon because the design is still being discussed here, but as you seem to be mostly satisfied with it now, I'll be happy to review a related doc-labeled PR :)
Again, although the documentation should cover the basics of using topology values to address #151, I would advise against putting too much effort into documenting very sophisticated code snippets using it, because possibly #141 could lead to Topology becoming the internal structure of Model, and this would possibly surface into breaking changes to the API introduced here in #152.
So, yeah, focus on the science between these new features and their meaning (that should not change). I hope this comment should be helpful in this respect. Then don't be overly precise about the code details ;)

What I can do on my side is to start keeping track of API changes in a toplevel CHANGELOG.md file. This disconnections branch may become our first version bump to EcologicalNetworksDynamics v0.2.1 ;)

Unfortunately, I am leaving tonight for three weeks off. What I can do is merge this PR into dev right now so that we can sync up on this version which also happens to bundle a few bugfixes like f1b740d. But I'll leave it open and won't merge it into main until we are confident with the design and we have written the corresponding doc and the CHANGELOG :)

See you about end of August <3

[ismael-lajaaiti]

Invalid edge type label: :facilitation. Valid labels in this model are :refuge and :trophic. ..?

Yes, IMO this is already better!

For the same reason: in principle, anyone could add the interaction type they want to the system by creating their own custom components (although it's (very) not ready yet). So given input like :rfeuge, there is no way the system could decide whether to error (incorrect name) or to return an empty matrix (valid name), because the list of names is open. So it needs to always take the same decision in both cases. My opinion now is that issuing an error with :facilitation if there is no facilitation in the model is ok, but returning an empty matrix with :rfeuge if it is a typo is absolutely unacceptable ^ ^"

OK, I see.

I'm curious though: what would you use this empty matrix for? How come your code needs to request facilitation links in a model without a Facilitation component?

First, my reasoning was that a model without a Facilitation component and a model with an empty Facilitation component is essentially the same model. Second, for instance, let's say you have simulated lots of models with and without different non-trophic components, and you want to investigate the impact of the number of non-trophic links on a given property (e.g., species diversity). We may want to loop on this set of simulated models, extract their number of non-trophic links at the end of the simulation to plot that against the final diversity. Or something along those lines.

Oh, I thought these were guaranteed not to generate such degenerated networks. If they aren't, then of course we could feature that. The algorithm would just be bruteforce like the one to avoid cycles, right? Draw random networks until we find one without this kind of species. Would you be okay to leave this for another github issue? I am going to screen through every component again in the context of #139. I could maybe pick this new feature up when I get to rewriting the Foodweb blueprints :)

Yes, and yes. ;)
Issue opened #159.

Unfortunately, I am leaving tonight for three weeks off. What I can do is merge this PR into dev right now so that we can sync up on this version which also happens to bundle a few bugfixes like f1b740d. But I'll leave it open and won't merge it into main until we are confident with the design and we have written the corresponding doc and the CHANGELOG :)

See you by the end of August <3

Enjoy your (well deserved) holidays! 🌴 ☀️

[iago-lito]

I'm curious though: what would you use this empty matrix for? How come your code needs to request facilitation links in a model without a Facilitation component?

First, my reasoning was that a model without a Facilitation component and a model with an empty Facilitation component is essentially the same model.

That is.. very true. They are the same in principle, but they differ in the code because the former does not store extra facilitation data whereas the latter does. This is an unfortunate quirk.. :\ But I'm not sure what to do about it without also stating that "every model is a model with an empty Brakbaorgas layer" 8\ Luckily this is rather a theoretical/philosophical quirk, and we haven't yet come across a situation where it's really annoying. But yeah, m1 == m2 would fail to capture the identity between two such models :(

say you have simulated lots of models with and without different non-trophic components, and you want to investigate the impact of the number of non-trophic links on a given property (e.g., species diversity). We may want to loop on this set of simulated models, extract their number of non-trophic links at the end of the simulation to plot that against the final diversity.

Okay, I see two approaches in this situations:

1. Simulate all your models with a FacilitationLayer, some of which happen to be empty. Then species_adjacency_matrix(g, :facilitation) should never error and return null matrices when empty.
2. OR guard your final extraction with something like nl = has_component(model, FacilitationLayer) ? model.n_facilitation_links : 0.

Would either be okay ? :)