Merge pull request #147 from VirtualPlantLab/fix-regression-benchmark

Fixing #146: performance regression

Author: VEZY

src/PlantSimEngine.jl

@@ -43,6 +43,9 @@ include("component_models/RefVector.jl")
 # Simulation table (time-step table, from PlantMeteo):
 include("component_models/TimeStepTable.jl")
 
+# Declaring the dependency graph
+include("dependencies/dependency_graph.jl")
+
 # List of models:
 include("component_models/ModelList.jl")
 include("mtg/MultiScaleModel.jl")
@@ -53,8 +56,7 @@ include("component_models/get_status.jl")
 # Transform into a dataframe:
 include("dataframe.jl")
 
-# Model dependencies:
-include("dependencies/dependency_graph.jl")
+# Computing model dependencies:
 include("dependencies/soft_dependencies.jl")
 include("dependencies/hard_dependencies.jl")
 include("dependencies/traversal.jl")

src/component_models/ModelList.jl

@@ -147,7 +147,8 @@ julia> [typeof(models[i][1]) for i in keys(status(models))]
 struct ModelList{M<:NamedTuple,S}
     models::M
     status::S
-    type_promotion::Union{Nothing, Dict}
+    type_promotion::Union{Nothing,Dict}
+    dependency_graph::DependencyGraph
 end
 
 #=function ModelList(models::M, status::Status) where {M<:NamedTuple{names,T} where {names,T<:NTuple{N,<:AbstractModel} where {N}}}
@@ -160,7 +161,6 @@ function ModelList(
     status=nothing,
     type_promotion::Union{Nothing,Dict}=nothing,
     variables_check::Bool=true,
-    nsteps=nothing,
     kwargs...
 )
 
@@ -187,10 +187,13 @@ function ModelList(
     ts_kwargs = homogeneous_ts_kwargs(status)
     ts_kwargs = add_model_vars(ts_kwargs, mods, type_promotion)
 
+
+
     model_list = ModelList(
         mods,
         ts_kwargs,
-        type_promotion
+        type_promotion,
+        dep(; verbose=true, mods...)
     )
     variables_check && !is_initialized(model_list)
 
@@ -219,8 +222,8 @@ function add_model_vars(x, models, type_promotion)
 
     # If the user gave a status, we check if all the variables are already initialized:
     vars_in_x = status_keys(x)
-    status_x = 
-    all([k in vars_in_x for k in keys(ref_vars)]) && return isa(x, Status) ? x : Status(x)  # If so, we return the input
+    status_x =
+        all([k in vars_in_x for k in keys(ref_vars)]) && return isa(x, Status) ? x : Status(x)  # If so, we return the input
 
     # Else, we add the variables by making a new object (carefull, this is a copy so it takes more time):
 
@@ -229,15 +232,15 @@ function add_model_vars(x, models, type_promotion)
 
     # If the user gave an empty status, we initialize all variables to their default values:
     if x === nothing
-       return Status(ref_vars)
+        return Status(ref_vars)
     end
-    
+
     if Tables.istable(x)
         # This situation only occurs if the user provided a table instead of a status
         # Meaning we have a status of vector values, all initialized up to a certain point
         # Unsure this is desirable, as that means run! does nothing or overwrites everything
         # Anyway, we wish to create a NamedTuple() of Vectors here
-        x_full = (;zip(propertynames(x), Tables.columns(x))...)
+        x_full = (; zip(propertynames(x), Tables.columns(x))...)
         x_full = merge(ref_vars, x_full)
 
     else
@@ -286,7 +289,7 @@ PlantSimEngine.homogeneous_ts_kwargs((Tₗ=[25.0, 26.0], aPPFD=1000.0))
 function homogeneous_ts_kwargs(kwargs::NamedTuple{N,T}) where {N,T}
     length(kwargs) == 0 && return kwargs
     vars_vals = collect(Any, values(kwargs))
- 
+
     vars_array = NamedTuple{keys(kwargs)}(j for j in vars_vals)
 
     return vars_array
@@ -325,15 +328,17 @@ function Base.copy(m::T) where {T<:ModelList}
     ModelList(
         m.models,
         deepcopy(m.status),
-        deepcopy(m.type_promotion)
+        deepcopy(m.type_promotion),
+        deepcopy(m.dependency_graph)
     )
 end
 
 function Base.copy(m::T, status) where {T<:ModelList}
     ModelList(
         m.models,
         status,
-        deepcopy(m.type_promotion)
+        deepcopy(m.type_promotion),
+        deepcopy(m.dependency_graph)
     )
 end
 
@@ -465,7 +470,7 @@ function convert_vars!(mapped_vars::Dict{String,Dict{Symbol,Any}}, type_promotio
 end
 
 function Base.show(io::IO, ::MIME"text/plain", t::ModelList)
-    print(io, dep(t, verbose=false))
+    print(io, dep(t))
     print(io, status(t))
 end
 

src/component_models/Status.jl

@@ -135,26 +135,28 @@ end
 
 # Returns a status with all vector variables replaced with their first value (ie a Status ready for simulation)
 # also returns a tuple of symbols corresponding to the vector variables
-function flatten_status(s::Status)
-    status_values_flattened = NamedTuple()
-    vector_variables = NamedTuple()
-
-    for (var, value) in zip(keys(s), s)
-        if length(value) > 1
-            vector_variables = (vector_variables..., var)
-            status_values_flattened = (status_values_flattened..., value[1])
-        else
-            status_values_flattened = (status_values_flattened..., value)
-        end
+function flatten_status(s::Status{T}) where {T}
+    n_vars_several_values = findall(x -> length(x) > 1, s)
+    if length(n_vars_several_values) == 0
+        return s, n_vars_several_values
+    else
+        return Status{keys(s)}(first.(values(s))), n_vars_several_values
     end
-
-    return Status(; zip(keys(s), status_values_flattened)...), vector_variables
 end
 
-# Update to the next timestep the variables that were passed in as vectors by the user
-function update_vector_variables(s::Status, sf::Status, vector_variables, i)
-    for vec in vector_variables
-        sf[vec] = s[vec][i]
+"""
+    set_variables_at_timestep!(status_timestep::Status, user_status::Status, variables_to_update, timestep)
+
+Update `status_timestep` to the current values at the `timestep` for all `variables_to_update` in the status provided by the user (`user_status`).
+The variables to update are given in `variables_to_update`, which is a vector of symbols.
+
+`status_timestep` is a status representing a single time-step. `user_status` is the status provided that gives values for variables that are not computed by any model.
+It may give constant values or vectors of values, in which case the `timestep` is used to select the value to use for the current time step.
+
+"""
+function set_variables_at_timestep!(status_timestep::Status, user_status::Status, variables_to_update, timestep)
+    for vec in variables_to_update
+        status_timestep[vec] = user_status[vec][timestep]
     end
 end
 

src/dependencies/dependencies.jl

@@ -1,8 +1,9 @@
 dep(::T, nsteps=1) where {T<:AbstractModel} = NamedTuple()
 
 """
-    dep(m::ModelList, nsteps=1; verbose::Bool=true)
+    dep(m::ModelList)
     dep(mapping::Dict{String,T}; verbose=true)
+    dep!(m::ModelList, nsteps=1)
 
 Get the model dependency graph given a ModelList or a multiscale model mapping. If one graph is returned, 
 then all models are coupled. If several graphs are returned, then only the models inside each graph are coupled, and
@@ -34,7 +35,12 @@ to other scales if needed. Then we transform all these nodes into soft dependenc
 Then we traverse all these and we set nodes that need outputs from other nodes as inputs as children/parents.
 If a node has no dependency, it is set as a root node and pushed into a new Dict (independant_process_root). This Dict is the returned dependency graph. And 
 it presents root nodes as independent starting points for the sub-graphs, which are the models that are coupled together. We can then traverse each of 
-these graphs independently to retrieve the models that are coupled together, in the right order of execution.
+these graphs independently to r
+
+# Notes
+
+The difference between `dep(m::ModelList)` and `dep!(m::ModelList, nsteps)` is that the first one returns the dependency graph found in the model list, while the 
+second one returns the dependency graph with the specified number of steps, modifying the simulation IDs of each node in the graph (`simulation_id=fill(0, nsteps)`).
 
 # Examples
 
@@ -75,8 +81,18 @@ function dep(nsteps=1; verbose::Bool=true, vars...)
     return deps
 end
 
-function dep(m::ModelList, nsteps=1; verbose::Bool=true)
-    dep(nsteps; verbose=verbose, m.models...)
+function dep(m::ModelList)
+    m.dependency_graph
+end
+
+function dep!(m::ModelList, nsteps=1)
+    traverse_dependency_graph!(m.dependency_graph; visit_hard_dep=false) do node
+        if length(node.simulation_id) != nsteps
+            node.simulation_id = fill(0, nsteps)
+        end
+    end
+
+    return m.dependency_graph
 end
 
 

src/mtg/save_results.jl

@@ -112,7 +112,7 @@ julia> collect(keys(preallocated_vars["Leaf"]))
 
 function pre_allocate_outputs(statuses, statuses_template, reverse_multiscale_mapping, vars_need_init, outs, nsteps; type_promotion=nothing, check=true)
     outs_ = Dict{String,Vector{Symbol}}()
-    
+
     # default behaviour : track everything
     if isnothing(outs)
         for organ in keys(statuses)
@@ -130,7 +130,7 @@ function pre_allocate_outputs(statuses, statuses_template, reverse_multiscale_ma
         end
     end
 
-    len = Dict{String, Int}()
+    len = Dict{String,Int}()
     for (organ, vals) in outs_
         len[organ] = length(outs_[organ])
         unique!(outs_[organ])
@@ -210,15 +210,15 @@ function pre_allocate_outputs(statuses, statuses_template, reverse_multiscale_ma
     node_type = only(node_type)
 
     # I don't know if this function barrier is necessary
-    preallocated_outputs = Dict{String, Vector}()
+    preallocated_outputs = Dict{String,Vector}()
     complete_preallocation_from_types!(preallocated_outputs, nsteps, outs_, node_type, statuses_template)
     return preallocated_outputs
 end
 
 function complete_preallocation_from_types!(preallocated_outputs, nsteps, outs_, node_type, statuses_template)
     types = Vector{DataType}()
     for organ in keys(outs_)
-        
+
         outs_no_node = filter(x -> x != :node, outs_[organ])
 
         #types = [typeof(status_from_template(statuses_template[organ])[var]) for var in outs[organ]]
@@ -230,14 +230,14 @@ function complete_preallocation_from_types!(preallocated_outputs, nsteps, outs_,
         symbols_tuple = (:timestep, :node, outs_no_node...,)
         # using node_type.parameters[1] is clunky, but covers both NodeMTG and AbstractNodeMTG types
         values_tuple = (1, MultiScaleTreeGraph.Node((node_type.parameters[1])("/", "Uninitialized", 0, 0),), values...,)
-    
+
         # Dummy value to make accessing the type easier 
         # (empty arrays don't have references to an instance, so their types can't be inspected and manipulated as easily)
-        dummy_status = (;zip(symbols_tuple, values_tuple)...)
+        dummy_status = (; zip(symbols_tuple, values_tuple)...)
         data = typeof(Status(dummy_status))[]
         resize!(data, nsteps)
-        
-        for ii in 1:nsteps 
+
+        for ii in 1:nsteps
             data[ii] = Status(dummy_status)
         end
         preallocated_outputs[organ] = data
@@ -278,22 +278,22 @@ function save_results!(object::GraphSimulation, i)
         # So there may be possible simplifications (maybe no need for a function barrier, perhaps the resizing could be made a one-liner...)
         # But this should work without causing visible performance regressions on XPalm
         len = length(outs[organ])
-        if length(statuses[organ]) + index - 1 > len            
+        if length(statuses[organ]) + index - 1 > len
             min_required = max(length(statuses[organ]) + index - len, index)
-            
-            extra_length = 2*min_required - len
+
+            extra_length = 2 * min_required - len
             data = eltype(outs[organ])[]
             resize!(data, extra_length)
             dummy_value = NamedTuple(outs[organ][1])
             # TODO set timestep to 0 for clarity ?
-            
+
             # Using fill! caused Ref issues, so call a Status constructor here instead of passing a prebuilt value
             # This will avoid having all array entries point to the same ref but keep construction cost at a minimum
             for new_entry in 1:extra_length
                 data[new_entry] = Status(dummy_value)
             end
 
-            outs[organ] = cat(outs[organ], data, dims=1)           
+            outs[organ] = cat(outs[organ], data, dims=1)
             #println("len : ", len, " statuses #", length(statuses[organ]), " index ", index)
             #println("min_required : ", min_required, " extra_length ", extra_length, " new len ", length(outs[organ]))
         end
@@ -316,15 +316,9 @@ function copy_tracked_outputs_into_vector!(outs_organ, i, statuses_organ, tracke
     return j
 end
 
-
 function pre_allocate_outputs(m::ModelList, outs, nsteps; type_promotion=nothing, check=true)
-    
-    # NOTE : init_variables recreates a DependencyGraph, it's not great
-    # TODO : copy ?
-    out_vars_pre_type_promotion = merge(init_variables(m; verbose=false)...)
-    
-    # bit hacky, could be cleaned up
-    out_vars_all = convert_vars(out_vars_pre_type_promotion, m.type_promotion)
+    st, = flatten_status(status(m))
+    out_vars_all = convert_vars(st, type_promotion)
 
     out_keys_requested = Symbol[]
     if !isnothing(outs)
@@ -337,9 +331,10 @@ function pre_allocate_outputs(m::ModelList, outs, nsteps; type_promotion=nothing
 
     # default implicit behaviour, track everything
     if isempty(out_keys_requested)
-        out_vars_requested = out_vars_all
+        # We already have the status here, just repeating its value:
+        out_vars_requested = NamedTuple(out_vars_all)
     else
-        unexpected_outputs = setdiff(out_keys_requested, status_keys(status(m)))
+        unexpected_outputs = setdiff(out_keys_requested, keys(st))
 
         if !isempty(unexpected_outputs)
             e = string(
@@ -354,22 +349,21 @@ function pre_allocate_outputs(m::ModelList, outs, nsteps; type_promotion=nothing
                 @info e
                 [delete!(unexpected_outputs, i) for i in unexpected_outputs]
             end
-        end    
+        end
 
         out_defaults_requested = (out_vars_all[i] for i in out_keys_requested)
-        out_vars_requested = (;zip(out_keys_requested, out_defaults_requested)...)
+        out_vars_requested = (; zip(out_keys_requested, out_defaults_requested)...)
     end
 
-    outputs_timestep = fill(out_vars_requested, nsteps)
-    return TimeStepTable([Status(i) for i in outputs_timestep])
+    return TimeStepTable([Status(out_vars_requested) for i in Base.OneTo(nsteps)])
 end
 
 function save_results!(status_flattened::Status, outputs, i)
-    if length(outputs) == 0 
-        return 
+    if length(outputs) == 0
+        return
     end
     outs = outputs[i]
-    
+
     for var in keys(outs)
         outs[var] = status_flattened[var]
     end

src/processes/model_initialisation.jl

@@ -65,8 +65,8 @@ mapping = Dict(
 to_initialize(mapping)
 ```
 """
-function to_initialize(m::ModelList; verbose::Bool=true)
-    needed_variables = to_initialize(dep(m; verbose=verbose))
+function to_initialize(m::ModelList)
+    needed_variables = to_initialize(dep(m))
     to_init = Dict{Symbol,Tuple}()
     for (process, vars) in needed_variables
         # default_values = needed_variables[:process1]
@@ -245,7 +245,7 @@ function init_variables(model::T; verbose::Bool=true) where {T<:AbstractModel}
 end
 
 function init_variables(m::ModelList; verbose::Bool=true)
-    init_variables(dep(m; verbose=verbose))
+    init_variables(dep(m))
 end
 
 function init_variables(m::DependencyGraph)
@@ -300,7 +300,7 @@ is_initialized(models)
 ```
 """
 function is_initialized(m::T; verbose=true) where {T<:ModelList}
-    var_names = to_initialize(m; verbose=verbose)
+    var_names = to_initialize(m)
 
     if any([length(to_init) > 0 for (process, to_init) in pairs(var_names)])
         verbose && @info "Some variables must be initialized before simulation: $var_names (see `to_initialize()`)" maxlog = 1