🚧 ⌛Start designing Topology to represent whole Model.

src/EcologicalNetworksDynamics.jl

@@ -40,6 +40,9 @@ using .AliasingDicts
 include("./multiplex_api.jl")
 using .MultiplexApi
 
+# Types to represent the model under a pure topological perspective.
+include("./topology.jl") # Will be part of the internals after their refactoring.
+
 #-------------------------------------------------------------------------------------------
 # "Inner" parts: legacy internals.
 

src/Internals/Internals.jl

@@ -58,6 +58,7 @@ const Option{T} = Union{Nothing,T}
 # Since parts of the API is being extracted out of this module to survive,
 # authorize using it here.
 using ..EcologicalNetworksDynamics
+const Topology = EcologicalNetworksDynamics.Topology # Part of future refactoring here.
 const equal_fields = EcologicalNetworksDynamics.equal_fields
 
 include("./macros.jl")

src/Internals/model/model_parameters.jl

@@ -13,6 +13,7 @@ mutable struct ModelParameters
     biorates::BioRates # (this one does but all values are initially 'nothing' inside)
     functional_response::Option{FunctionalResponse}
     producer_growth::Option{ProducerGrowth}
+    topology::Topology # This will actually be part of the future refactoring.
     # Since 'foodweb' is still a mandatory input to construct interaction layers,
     # keep this artificial reference to it to make the system/components API work
     # until we refactor all the internals.

src/components/foodweb_graph.jl

@@ -1,5 +1,7 @@
 # Wrap a SimpleDiGraph into a custom type indexed by stable labels
 # to represent either a foodweb or its biomass-foodweb restriction.
+# HERE: rather a whole model graph with all nodes/edges layers,
+#       and make this only a view into it.
 
 """
 Representation of a restriction of the graph model

src/components/nontrophic_layers/nontrophic_components_utils.jl

@@ -139,10 +139,13 @@ end
 # ==========================================================================================
 # Check/expand full layer components.
 
+has_nontrophic_layers(model) = model.network isa MultiplexNetwork
+export has_nontrophic_layers
+
 # The application procedure differs
 # whether the NTI layer is the first to be set or not.
 function set_layer!(model, interaction, layer)
-    if model.network isa Internals.FoodWeb
+    if !has_nontrophic_layers(model)
         # First NTI component to be added.
         # Switch from plain foodweb to a multiplex network.
         S = model.richness

src/methods/graphs.jl

@@ -1,3 +1,10 @@
+# Retrieve underlying model topology.
+@expose_data graph begin
+    property(topology, model_graph)
+    ref(m -> m.topology)
+    get(m -> deepcopy(m.topology))
+end
+
 """
     restrict_to_live!(g::TrophicGraph, biomasses; threshold = 0)
 
@@ -118,3 +125,34 @@ function starving_consumers(graph::TrophicGraph)
     consumers
 end
 export starving_consumers
+
+# ==========================================================================================
+# Expose as checked model methods.
+
+function biomass_foodweb(
+    m::InnerParms,
+    biomasses;
+    threshold = 0,
+    warn_nontrophic_links = true,
+)
+    g = m.trophic_graph
+    if warn_nontrophic_links && has_nontrophic_layers(m)
+        @warn "Restricting ecological model to only a trophic graph \
+               with extinct species nodes removed \
+               may yield inconsistent expectations regarding \
+               the meaning of remaining non-trophic interactions."
+    end
+    restrict_to_live!(g, biomasses; threshold)
+end
+@method biomass_foodweb depends(Foodweb)
+
+function disconnected_components(m::InnerParms, biomasses; threshold = 0)
+    g = biomass(m, biomasses; threshold)
+    # HERE: this raises too many complications
+    # just because non-trophic interactions
+    # are not represented by the trophic graph.
+    # Make the trophic graph a simple view into the whole model topology,
+    # reified as a true `ModelGraph` type,
+    # with proper layers for nodes and edges.
+    # Drop Graphs.jl in this respect: only locally useful for very particular algorithms.
+end

src/topology.jl

@@ -0,0 +1,286 @@
+module Topologies
+
+argerr(mess) = throw(ArgumentError(mess))
+
+# Mark removed nodes.
+struct Tombstone end
+
+"""
+Values of this type are constructed from a model value,
+to represent its pure topology:
+
+  - Nodes identity and types: species, nutrients, patches..
+  - Edges types: trophic interaction, non-trophic interactions, migrations..
+
+They are supposed not to be mutated,
+as they carry faithful topological information reflecting the model
+at the moment it has been extracted from it.
+Nodes and edge information can be queried from it using labels or indices.
+They *may* be removed to represent *e.g.* species extinction,
+and study topological consequences or removing them.
+But the indices and labels remain stable,
+always consistent with their indices from the model value when extracted.
+As a consequence, every node (at least) in the topology
+can be queried for having been 'removed' or not: tombstones remain.
+No tombstone remain for edges: once removed, there is no trace left of them.
+"""
+struct Topology
+    # List/index possible types for nodes and edges.
+    node_types_labels::Vector{Symbol}
+    node_types_index::Dict{Symbol,Int}
+    edge_types_labels::Vector{Symbol}
+    edge_types_index::Dict{Symbol,Int}
+    # List nodes and their associated types.
+    # Sorted are *sorted by type*:
+    # so all nodes of the same type are stored contiguously in this array.
+    nodes_labels::Vector{Symbol}
+    nodes_index::Dict{Symbol,Int}
+    nodes_types::Vector{UnitRange{Int}} # {type -> (start, end) of nodes with this type}
+    # Topological information: paired, layered adjacency lists.
+    outgoing::Vector{Union{Tombstone,Vector{Vector{Int}}}}
+    incoming::Vector{Union{Tombstone,Vector{Vector{Int}}}}
+    # [node: [edgetype: [nodeid]]]
+    # ^--------------------------^- : Adjacency list: one entry per node 'N'.
+    #        ^-------------------^- : Tombstone (removed node) or one entry per edge type.
+    #                   ^--------^- : One entry per neighbour of 'N': its index.
+    # Cached redundant information.
+    n_edges::Vector{Int} # Per edge type.
+    n_nodes::Vector{Int} # Per node type, not counting tombstones.
+    Topology() = new([], Dict(), [], Dict(), [], Dict(), [], [], [], [])
+end
+export Topology
+
+#-------------------------------------------------------------------------------------------
+# Construction primitives.
+
+# Only push whole slices of nodes of a new type at once.
+function push_nodes!(top::Topology, labels, type::Symbol)
+    haskey(top.node_types_index, type) &&
+        argerr("Node type :$type has already been pushed.")
+    haskey(top.edge_types_index, type) &&
+        argerr("Node type :$type could be confused with edge type :$type.")
+    push!(top.node_types_labels, type)
+    top.node_types_index[type] = length(top.node_types_labels)
+    n_before = length(top.nodes_labels)
+    nl = top.nodes_labels
+    for lab::Symbol in labels
+        haskey(nl, lab) && argerr("Label :$lab was already given \
+                                   to a node of type :$(node_type(top, lab)).")
+        push!(nl, lab)
+        top.nodes_index[lab] = length(nl)
+        for adj in (top.outgoing, top.incoming)
+            push!(adj, Vector{Vector{Int}}())
+        end
+    end
+    n_after = length(nl)
+    push!(top.nodes_types, n_before+1:n_after)
+    top
+end
+export push_nodes!
+
+function new_edge_type!(top::Topology, type::Symbol)
+    haskey(top.edge_types_index, type) && argerr("Edge type :$type has already been added.")
+    haskey(top.node_types_index, type) &&
+        argerr("Edge type :$type would be confused with node type :$type.")
+    push!(top.edge_types_labels, type)
+    top.edge_types_index[type] = length(top.edge_types_labels)
+    for adj in (top.outgoing, top.incoming)
+        for node in adj
+            node isa Tombstone && continue
+            push!(node, [])
+        end
+    end
+    push!(top.n_edges, 0)
+    top
+end
+export new_edge_type!
+
+#-------------------------------------------------------------------------------------------
+# Basic unchecked queries.
+# (ref-leaking methods are still _protected)
+
+const imap = Iterators.map
+const ifilter = Iterators.filter
+idmap(x) = imap(identity, x) # Useful to not leak refs to private collections.
+
+# Information about types.
+n_node_types(top::Topology) = length(top.node_types_labels)
+n_edge_types(top::Topology) = length(top.edge_types_labels)
+node_type_label(top::Topology, i::Int) = top.node_types_labels[i]
+node_type_index(top::Topology, lab::Symbol) = top.node_types_index[lab]
+node_type_label(::Topology, lab::Symbol) = lab
+node_type_index(::Topology, i::Int) = i
+edge_type_label(top::Topology, i::Int) = top.edge_types_labels[i]
+edge_type_index(top::Topology, lab::Symbol) = top.edge_types_index[lab]
+edge_type_label(::Topology, lab::Symbol) = lab
+edge_type_index(::Topology, i::Int) = i
+node_types(top::Topology) = idmap(identity, top.node_types_labels)
+edge_types(top::Topology) = idmap(identity, top.edge_types_labels)
+
+# General information about nodes.
+n_nodes(top::Topology, type) = top.n_nodes[node_type_index(type)]
+n_nodes_and_removed(top::Topology, type) = length(top.nodes_types[node_type_index(type)])
+nodes_indices(top::Topology, type) = top.nodes_types[node_type_index(type)] # Okay to leak: immutable.
+_nodes_labels(top::Topology, type) = top.nodes_labels[nodes_indices(top, type)]
+node_labels(top::Topology, type) = idmap(identity, _nodes_labels(top, type))
+
+# Particular information about nodes.
+node_label(top::Topology, i::Int) = top.nodes_labels[i]
+node_index(top::Topology, label::Symbol) = top.nodes_index[label]
+node_label(::Topology, lab::Symbol) = lab
+node_index(::Topology, i::Int) = i
+node_type_index(top::Topology, id) =
+    findfirst(range -> node_index(top, id) in range, top.nodes_types)
+node_type(top::Topology, id) = node_type_label(node_type_index(top, id))
+
+# Information about edges.
+n_edges(top::Topology, type) = top.n_edges[edge_type_index(type)]
+
+# Direct neighbourhood when querying particular edge type.
+# (assuming target is not a tombstone)
+function _outgoing_indices(top::Topology, node, edge_type)
+    i_type = edge_type_index(top, edge_type)
+    _outgoing_indices(top, node)[i_type]
+end
+function _incoming_indices(top::Topology, node, edge_type)
+    i_type = edge_type_index(top, edge_type)
+    _incoming_indices(top, node)[i_type]
+end
+outgoing_indices(top::Topology, node, type) = idmap(_outgoing_indices(top, node, type))
+incoming_indices(top::Topology, node, type) = idmap(_incoming_indices(top, node, type))
+outgoing_labels(top::Topology, node, type) =
+    imap(i -> top.nodes_labels[i], _outgoing_indices(top, node, type))
+incoming_labels(top::Topology, node, type) =
+    imap(i -> top.nodes_labels[i], _incoming_indices(top, node, type))
+
+# Direct neighbourhood: return twolevel iterator:
+# first iterate on edge types, then neighbours.
+# (assuming target is not a tombstone)
+function _outgoing_indices(top::Topology, node)
+    i_node = node_index(top, node)
+    top.outgoing[i_node]
+end
+function _incoming_indices(top::Topology, node)
+    i_node = node_index(top, node)
+    top.incoming[i_node]
+end
+outgoing_indices(top::Topology, node) = imap(
+    (i_edge_type, _neighbours) -> (i_edge_type, idmap(_neighbours)),
+    enumerate(_outgoing_indices(top, node)),
+)
+incoming_indices(top::Topology, node) = imap(
+    (i_edge_type, _neighbours) -> (i_edge_type, idmap(_neighbours)),
+    enumerate(_incoming_indices(top, node)),
+)
+outgoing_labels(top::Topology, node) = imap(
+    (i_edge, _neighbours) -> (
+        top.edge_types_labels[i_edge],
+        imap(i_node -> top.nodes_labels[i_node], _neighbours),
+    ),
+    enumerate(_outgoing_indices(top, node)),
+)
+incoming_labels(top::Topology, node) = imap(
+    (i_edge, _neighbours) -> (
+        top.edge_types_labels[i_edge],
+        imap(i_node -> top.nodes_labels[i_node], _neighbours),
+    ),
+    enumerate(_incoming_indices(top, node)),
+)
+
+# Filter adjacency iterators given one particular edge type.
+# Also return twolevel iterators: focal node, then its neighbours.
+function _outgoing_edges_indices(top::Topology, edge_type)
+    i_type = edge_type_index(top, edge_type)
+    imap(ifilter(enumerate(top.outgoing)) do (_, node)
+        !(node isa Tombstone)
+    end) do (i, _neighbours)
+        (i, _neighbours[i_type])
+    end
+end
+function _incoming_edges_indices(top::Topology, edge_type)
+    i_type = edge_type_index(top, edge_type)
+    imap(ifilter(enumerate(top.incoming)) do (_, node)
+        !(node isa Tombstone)
+    end) do (i, _neighbours)
+        (i, _neighbours[i_type])
+    end
+end
+outgoing_edges_indices(top::Topology, edge_type) = imap(
+    (i_node, _neighbours) -> (i_node, idmap(_neighbours)),
+    _outgoing_edges_indices(top, edge_type),
+)
+incoming_edges_indices(top::Topology, edge_type) = imap(
+    (i_node, _neighbours) -> (i_node, idmap(_neighbours)),
+    _incoming_edges_indices(top, edge_type),
+)
+outgoing_edges_labels(top::Topology, edge_type) = imap(
+    (i_node, _neighbours) ->
+        (node_label(top, i_node), idmap(i -> node_label(top, i), _neighbours)),
+    _outgoing_edges_indices(top, edge_type),
+)
+incoming_edges_labels(top::Topology, edge_type) = imap(
+    (i_node, _neighbours) ->
+        (node_label(top, i_node), idmap(i -> node_label(top, i), _neighbours)),
+    _incoming_edges_indices(top, edge_type),
+)
+
+
+#-------------------------------------------------------------------------------------------
+# Display.
+
+function Base.show(io::IO, top::Topology)
+    nnt = n_node_types(top)
+    net = new_edge_type!(top)
+    s = n -> n > 1 ? "s" : ""
+    print(
+        io,
+        "Topology for $nnt node type$(s(nnt)) \
+         and $net edge type$s(net):\n",
+    )
+    tombs = [] # Collect removed nodes to display at the end.
+    for (i_type, type) in enumerate(node_types(top))
+        tomb = []
+        push!(tombs, tomb)
+        print(io, "  Nodes '$type':")
+        empty = true
+        for node in _nodes_labels(i_type)
+            if node isa Tombstone
+                push!(tomb, node)
+                continue
+            end
+            print(io, " :$node")
+            empty = false
+        end
+        if empty
+            print(io, " <none>")
+        end
+        println(io)
+    end
+    for (i_type, type) in enumerate(edge_types(top))
+        print(io, "  Edges of type '$type':")
+        empty = true
+        for (i_source, _neighbours) in _outgoing_edges_indices(top, i_type)
+            isempty(_neighbours) && continue
+            source = node_label(top, i_source)
+            print(io, "    $source:")
+            for target in _neighbours
+                print(io, " :$target")
+            end
+            empty = false
+        end
+        if empty
+            print(io, " <none>")
+        end
+        println(io)
+    end
+end
+
+end
+
+# ==========================================================================================
+# DEBUG
+using .Topologies
+top = Topology()
+push_nodes!(top, Symbol.(collect("abc")), :species)
+new_edge_type!(top, :trophic)
+top #  HERE: debug all the above at least until displaying works.