Introduce Topology to check biomass graph disconnections.

docs/make.jl

@@ -45,4 +45,4 @@ makedocs(;
     ],
 )
 
-deploydocs(; repo = "github.com/BecksLab/EcologicalNetworksDynamics.jl", devbranch = "doc")
+deploydocs(; repo = "github.com/BecksLab/EcologicalNetworksDynamics.jl", devbranch = "dev")

src/EcologicalNetworksDynamics.jl

@@ -4,6 +4,7 @@ using Crayons
 using MacroTools
 using OrderedCollections
 using SparseArrays
+using Graphs
 
 #-------------------------------------------------------------------------------------------
 # Shared API internals.
@@ -21,12 +22,29 @@ argerr(mess) = throw(ArgumentError(mess))
 const Option{T} = Union{Nothing,T}
 const SparseMatrix{T} = SparseMatrixCSC{T,Int64}
 
+# Basic equivalence relation for recursive use.
+function equal_fields(a::T, b::T; ignore = Set{Symbol}()) where {T}
+    for name in fieldnames(T)
+        if name in ignore
+            continue
+        end
+        u, v = getfield.((a, b), name)
+        u == v || return false
+    end
+    true
+end
+
 include("./AliasingDicts/AliasingDicts.jl")
 using .AliasingDicts
 
 include("./multiplex_api.jl")
 using .MultiplexApi
 
+# Types to represent the model under a pure topological perspective.
+include("./Topologies/Topologies.jl")
+using .Topologies
+# (will be part of the internals after their refactoring)
+
 #-------------------------------------------------------------------------------------------
 # "Inner" parts: legacy internals.
 
@@ -55,14 +73,6 @@ export add!, properties, blueprints, components
 
 include("./dedicate_framework_to_model.jl")
 
-#-------------------------------------------------------------------------------------------
-# Analysis tools working on the output of the simulation.
-include("output-analysis.jl")
-export richness
-export persistence
-export shannon_diversity
-export total_biomass
-
 #-------------------------------------------------------------------------------------------
 # "Outer" parts: develop user-facing stuff here.
 
@@ -84,11 +94,13 @@ include("./expose_data.jl")
 # The actual user-facing components of the package are defined there,
 # connecting them to the internals via the framework.
 include("./components/main.jl")
-include("./methods/main.jl")
 
 # Additional exposed utils built on top of components and methods.
 include("./default_model.jl")
 include("./nontrophic_layers.jl")
+include("./simulate.jl")
+include("./topology.jl")
+include("./diversity.jl")
 
 # Avoid Revise interruptions when redefining methods and properties.
 Framework.REVISING = true

src/Framework/Framework.jl

@@ -96,6 +96,8 @@ module Framework
 #   - [ ] `depends(other_method_name)` to inherit all dependent components.
 #   - [ ] Recurring pattern: various blueprints types provide 'the same component': reify.
 #   - [ ] Namespace properties into like system.namespace.namespace.property.
+#   - [ ] Hooks need to trigger when special components combination become available.
+#         See for instance the expansion of `Nutrients.Nodes`.
 
 using Crayons
 using MacroTools

src/Framework/method_macro.jl

@@ -293,6 +293,12 @@ macro method(input...)
                 end
                 fn(s._value, args...; kwargs...)
             end
+            # TODO: if unexistent, then
+            # defining a fallback `fn(value, args...; _system=s, kwargs...)` method here
+            # would enable that implementors of `fn` take decisions
+            # depending on the whole system value, and components currently available.
+            # In particular, it would avoid the need
+            # to define both `_simulate` and `simulate` in the exposed lib.
         end,
     )
 

src/Framework/system.jl

@@ -173,6 +173,20 @@ function properties(s::System{V}) where {V}
 end
 export properties
 
+#-------------------------------------------------------------------------------------------
+
+# Basic recursive equivalence relation.
+function equal_fields(a::T, b::T) where {T}
+    for name in fieldnames(T)
+        u, v = getfield.((a, b), name)
+        u == v || return false
+    end
+    true
+end
+Base.:(==)(a::System{U}, b::System{V}) where {U,V} = U == V && equal_fields(a, b)
+Base.:(==)(a::Blueprint{U}, b::Blueprint{V}) where {U,V} =
+    U == V && typeof(a) == typeof(b) && equal_fields(a, b)
+
 #-------------------------------------------------------------------------------------------
 # Display.
 function Base.show(io::IO, sys::System)

src/GraphDataInputs/GraphDataInputs.jl

@@ -250,6 +250,7 @@ using OrderedCollections
 import ..SparseMatrix
 import ..argerr
 import ..checkfails
+import ..join_elided
 
 # Unhygienically define `loc` variable in macros to point to invocation line.
 # Assumes __source__ is available.

src/GraphDataInputs/check.jl

@@ -267,7 +267,7 @@ function outspace((i, j), (n, m))
 end
 either(symbols) =
     length(symbols) == 1 ? "$(repr(first(symbols)))" :
-    "either " * join(repr.(sort(collect(symbols))), ", ", " or ")
+    "either " * join_elided(sort(collect(symbols)), ", ", " or "; max = 12)
 
 #-------------------------------------------------------------------------------------------
 # Assuming the above check passed, check references against a template.

src/GraphDataInputs/convert.jl

@@ -125,6 +125,31 @@ function graphdataconvert(::Type{BinMap{<:Any}}, input; expected_I = nothing)
     res
 end
 
+# The binary case *can* accept boolean masks.
+function graphdataconvert(
+    ::Type{BinMap{<:Any}},
+    input::AbstractVector{Bool};
+    expected_I = Int64,
+)
+    res = BinMap{expected_I}()
+    for (i, val) in enumerate(input)
+        val && push!(res, i)
+    end
+    res
+end
+
+function graphdataconvert(
+    ::Type{BinMap{<:Any}},
+    input::AbstractSparseVector{Bool,I};
+    expected_I = I,
+) where {I}
+    res = BinMap{expected_I}()
+    for i in findnz(input)[1]
+        push!(res, i)
+    end
+    res
+end
+
 #-------------------------------------------------------------------------------------------
 # Similar, nested logic for adjacency maps.
 
@@ -193,6 +218,38 @@ function graphdataconvert(::Type{BinAdjacency{<:Any}}, input; expected_I = nothi
     res
 end
 
+# The binary case *can* accept boolean matrices.
+function graphdataconvert(
+    ::Type{BinAdjacency{<:Any}},
+    input::AbstractMatrix{Bool},
+    expected_I = Int64,
+)
+    res = BinAdjacency{expected_I}()
+    for (i, row) in enumerate(eachrow(input))
+        adj_line = BinMap(j for (j, val) in enumerate(row) if val)
+        isempty(adj_line) && continue
+        res[i] = adj_line
+    end
+    res
+end
+
+function graphdataconvert(
+    ::Type{BinAdjacency{<:Any}},
+    input::AbstractSparseMatrix{Bool,I},
+    expected_I = I,
+) where {I}
+    res = BinAdjacency{expected_I}()
+    nzi, nzj, _ = findnz(input)
+    for (i, j) in zip(nzi, nzj)
+        if haskey(res, i)
+            push!(res[i], j)
+        else
+            res[i] = BinMap([j])
+        end
+    end
+    res
+end
+
 # Alias if types matches exactly.
 graphdataconvert(::Type{Map{<:Any,T}}, input::Map{Symbol,T}) where {T} = input
 graphdataconvert(::Type{Map{<:Any,T}}, input::Map{Int64,T}) where {T} = input
@@ -273,9 +330,11 @@ end
 # Example usage:
 #   @tographdata var {Sym, Scal, SpVec}{Float64}
 #   @tographdata var YSN{Float64}
-macro tographdata(var, input)
+macro tographdata(var::Symbol, input)
     @defloc
-    var isa Symbol || argerr("Not a variable: $(repr(var)) at $loc.")
+    tographdata(loc, var, input)
+end
+function tographdata(loc, var, input)
     @capture(input, types_{Target_} | types_{})
     isnothing(types) && argerr("Invalid @tographdata target types at $loc.\n\
                                 Expected @tographdata var {aliases...}{Target}. \
@@ -317,3 +376,14 @@ function _tographdata(vsym, var, targets)
             The value received is $(repr(var)) ::$(typeof(var)).")
 end
 export @tographdata
+
+# Convenience to re-bind in local scope, avoiding the akward following pattern:
+#   long_var_name = @tographdata long_var_name <...>
+# In favour of:
+#   @tographdata! long_var_name <...>
+macro tographdata!(var::Symbol, input)
+    @defloc
+    evar = esc(var)
+    :($evar = $(tographdata(loc, var, input)))
+end
+export @tographdata!

src/GraphDataInputs/types.jl

@@ -217,7 +217,7 @@ inspace((a, b), (x, y)) = inspace(a, x) && inspace(b, y)
 # Pretty display for maps and adjacency lists.
 
 display_short(map::Map) = "{$(join(("$(repr(k)): $v" for (k, v) in map), ", "))}"
-function display_long(map::Map, level = 0)
+function display_long(map::Map; level = 0)
     res = "{"
     ind(n) = "\n" * repeat("  ", level + n)
     for (k, v) in map
@@ -227,7 +227,7 @@ function display_long(map::Map, level = 0)
 end
 
 display_short(map::BinMap) = "{$(join(("$(repr(k))" for k in map), ", "))}"
-function display_long(map::BinMap, level = 0)
+function display_long(map::BinMap; level = 0)
     res = "{"
     ind(n) = "\n" * repeat("  ", level + n)
     for k in map
@@ -238,11 +238,11 @@ end
 
 display_short(adj::Union{Adjacency,BinAdjacency}) =
     "{$(join(("$(repr(k)): $(display_short(list))" for (k, list) in adj), ", "))}"
-function display_long(adj::Union{Adjacency,BinAdjacency}, level = 0)
+function display_long(adj::Union{Adjacency,BinAdjacency}; level = 0)
     res = "{"
     ind(n) = "\n" * repeat("  ", level + n)
     for (k, list) in adj
-        res *= ind(1) * "$(repr(k)) => $(display_long(list, level + 1)),"
+        res *= ind(1) * "$(repr(k)) => $(display_long(list; level = level + 1)),"
     end
     res * ind(0) * "}"
 end

src/Internals/Internals.jl

@@ -58,6 +58,10 @@ 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 equal_fields = EcologicalNetworksDynamics.equal_fields
+
+# Part of future refactoring here.
+const Topology = EcologicalNetworksDynamics.Topologies.Topology
 
 include("./macros.jl")
 include("./inputs/foodwebs.jl")

src/Internals/inputs/biological_rates.jl

@@ -16,6 +16,8 @@ mutable struct BioRates
 end
 #### end ####
 
+Base.:(==)(a::BioRates, b::BioRates) = equal_fields(a, b)
+
 #### Type display ####
 """
 One line [`BioRates`](@ref) display.

src/Internals/inputs/environment.jl

@@ -2,6 +2,8 @@ mutable struct Environment
     T::Float64
 end
 
+Base.:(==)(a::Environment, b::Environment) = equal_fields(a, b)
+
 """
 One line Environment display.
 """

src/Internals/inputs/foodwebs.jl

@@ -34,6 +34,8 @@ mutable struct FoodWeb <: EcologicalNetwork
         new(A, sp, M, mc, mth, Dict(Symbol(s) => i for (i, s) in enumerate(sp)))
 end
 
+Base.:(==)(a::FoodWeb, b::FoodWeb) = equal_fields(a, b)
+
 """
     FoodWeb(
         A::AbstractMatrix;

src/Internals/inputs/functional_response.jl

@@ -27,6 +27,9 @@ struct LinearResponse <: FunctionalResponse
 end
 #### end ####
 
+Base.:(==)(a::U, b::V) where {U<:FunctionalResponse,V<:FunctionalResponse} =
+    U == V && equal_fields(a, b)
+
 #### Type display ####
 """
 One line display FunctionalResponse

src/Internals/inputs/producer_growth.jl

@@ -30,6 +30,9 @@ mutable struct NutrientIntake <: ProducerGrowth
     NutrientIntake(args...) = new(args...)
 end
 
+Base.:(==)(a::U, b::V) where {U<:ProducerGrowth,V<:ProducerGrowth} =
+    U == V && equal_fields(a, b)
+
 """
     length(n::NutrientIntake)
 

src/Internals/model/model_parameters.jl

@@ -11,6 +11,7 @@ mutable struct ModelParameters
     # These don't exactly have an 'empty' variant.
     network::Option{EcologicalNetwork}
     biorates::BioRates # (this one does but all values are initially 'nothing' inside)
+    _topology::Topology # This will actually be part of the future refactoring.
     functional_response::Option{FunctionalResponse}
     producer_growth::Option{ProducerGrowth}
     # Since 'foodweb' is still a mandatory input to construct interaction layers,
@@ -38,6 +39,7 @@ mutable struct ModelParameters
             NoTemperatureResponse(),
             nothing,
             BioRates(),
+            Topology(),
             repeat([nothing], 3)...,
             Dict(),
             Dict(),
@@ -46,10 +48,16 @@ mutable struct ModelParameters
     ModelParameters(args...) = new(args...)
 end
 # Required to fork the system.
-# Ok as long as the value above contains no critical self-references.
+# Self-references (very) fortunately are not a problem:
+# https://discourse.julialang.org/t/how-is-deepcopy-so-clever-regarding-aliasing-and-self-reference/113235
 Base.copy(m::ModelParameters) = deepcopy(m)
 #### end ####
 
+# Ignore cached values for equivalence relation.
+const model_equivalence_ignored_fields = Set([:_cache])
+Base.:(==)(a::ModelParameters, b::ModelParameters) =
+    equal_fields(a, b; ignore = model_equivalence_ignored_fields)
+
 #### Type display ####
 """
 One line ModelParameters display.
@@ -175,6 +183,7 @@ function ModelParameters(
         temperature_response,
         network,
         biorates,
+        Topology(),
         functional_response,
         producer_growth,
         nothing,

src/Internals/model/simulate.jl

@@ -120,8 +120,16 @@ function simulate(
         ExtinctionCallback(extinction_threshold, params, verbose),
     ),
     diff_code_data = (dudt!, params),
+    # FROM THE FUTURE - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+    # Record original user component-based model within the solution.
+    # The design should change during refactoring of the internals.
+    model = nothing,
     kwargs...,
 )
+    isnothing(model) ||
+        model._value === params ||
+        throw("Inconsistent input to `simulate`: this is a bug in the package.")
+
     # Interpret parameters and check them for consistency.
     S = richness(params)
     all(B0 .>= 0) ||
@@ -178,7 +186,13 @@ function simulate(
         u0 = B0
     end
 
-    p = (params = data, extinct_sp = extinct_sp, original_params = params)
+    p = (
+        params = data,
+        extinct_sp = extinct_sp,
+        original_params = params,
+        # Own the copy to not allow post-simulation modifications.
+        model = isnothing(model) ? nothing : copy(model),
+    )
     timespan = (t0, tmax)
     problem = ODEProblem(fun, u0, timespan, p)
     sol = solve(