Merge pull request #122 from VirtualPlantLab/add-modellist-to-mapping-tests-+-automated-model-generation-from-vector

Working version of a modellist to mapping conversion function, add a way to pass in vectors to a status in MultiScale, doc page for tips and workarounds

Author: Samuel-amap

Project.toml

@@ -12,6 +12,7 @@ FLoops = "cc61a311-1640-44b5-9fba-1b764f453329"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MultiScaleTreeGraph = "dd4a991b-8a45-4075-bede-262ee62d5583"
 PlantMeteo = "4630fe09-e0fb-4da5-a846-781cb73437b6"
+SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Term = "22787eb5-b846-44ae-b979-8e399b8463ab"
@@ -25,6 +26,7 @@ FLoops = "0.2"
 Markdown = "1.9"
 MultiScaleTreeGraph = "0.13, 0.14"
 PlantMeteo = "0.6"
+SHA = "0.7.0"
 Statistics = "1.9"
 Tables = "1"
 Term = "1, 2"

docs/src/model_coupling/tips_and_workarounds.md

@@ -0,0 +1,58 @@
+# Tips and workarounds
+
+## PlantSimEngine is actively being developed
+
+PlantSimEngine, despite the somewhat abstract codebase and generic simulation ambitions, is quite grounded in reality. There IS a desire to accomodate for a wide range of possible simulations, without constraining the user too much, but most features are developed on an as-needed basis, and grow out of necessity, partly from the requirements of an increasingly complex and refined implementation of an oil palm model, [XPalm](https://github.com/PalmStudio/XPalm.jl).
+
+Since the oil palm model is actively being developed, and some features aren't ready in PlantSimEngine, or require a lot of rewriting that we're not certain would be worth it (especially if it ends up constraining the codebase or what the user can do), some workarounds and shortcuts are occasionally used to circumvent a limitation. 
+
+There are also a couple of features that are quick hacks or that are meant for quick and dirty prototyping, not for production. 
+
+We'll list a few of them here, and will likely add some entry in the future listing some built-in limitations or implicit expectations of the package.
+
+## Making use of past states in multi-scale simulations
+
+It is possible to make use of the value of a variable in the past simulation timestep via the `PreviousTimeStep` mechanism in the mapping API (In fact, as mentioned elsewhere, it is the default way to break undesirable cyclic dependencies that can come up when coupling models, see : [Avoiding cyclic dependencies](@ref)).
+
+However, it is not possible to go beyond that through the mapping API. Something like `PreviousTimeStep(PreviousTimeStep(PreviousTimeStep(:carbon_biomass)))` is not supported. Don't do that.
+
+One way to access prior variable states is simply to write an ad hoc model that stores a few values into an array or however many variables you might need, which you can then feed into other models that might need it.
+
+## Having a variable simultaneously as input and output of a model 
+
+One current limitation of `PlantSimEngine` that can be occasionally awkward is that using the same variable name as input and output in a single model is unsupported. 
+
+(On a related note : it is not possible to have two variables with the same name *in the same scale*. They are considered as the same variable.)
+
+The reason being that it is usually impossible to automatically determine how the coupling is supposed to work out, when other dependencies latch onto such a model. The user would have to explicitely declare some order of simulation between several models, and some amount of programmer work would also be necessary to implement that extra API feature into `PlantSimEngine`.
+
+We haven't found an approach that was fully satisfactory from both a code simplicity and an API convenience POV. Especially when prototyping and adding in new models, as that might require redeclaring the simulation order for those specific variables. Ideas that came to mind felt like they might constrain the codebase for a more complex API, without enough benefit to justify it.
+
+The current workaround, while a little annoying, is very simple : *rename one of the variables*. It is not ideal, of course, as it means you might not be able to use a predefined model 'out of the box', but it does not have any of the tradeoffs and constrained mentioned above.
+
+## Passing in a vector in a mapping status at a specific scale
+
+You may have noticed that sometimes a vector (1-dimensional array) variable is passed into the `status` component of a `ModelList` in documentation examples (An example here with cumulative thermal time : [Model switching](@ref)).
+
+This is practical for simple simulations, or when quickly prototyping, to avoid having to write a model specifically for it. Whatever models make use of that variable are provided with one element corresponding to the current timestep every iteration.
+
+In multi-scale simulations, this feature is also supported, though not part of the main API. The way outputs and statuses work is a little different, so that little convenience feature is not as straightforward. 
+
+It is more brittle, makes use of not-recommended Julia metaprogramming features (`eval()`), fiddles with global variables, might not work outside of a REPL environment and is not tested for more complex interactions, so it may interact badly with variables that are mapped to different scales or in bizarre dependency couplings.
+
+Due to, uh, implementation quirks, the way to use this is as follows : 
+
+Call the function `replace_mapping_status_vectors_with_generated_models(mapping_with_vectors_in_status, timestep_model_organ_level, nsteps)`on your mapping.
+
+It will parse your mapping, generate custom models to store and feed the vector values each timestep, and return the new mapping you can then use for your simulation. It also slips in a couple of internal models that provide the timestep index to these models (so note that symbols `:current_timestep` and `:next_timestep` will be declared for that mapping). You can decide which scale/organ level you want those models to be in via the `timestep_model_organ_level`parameter. `nsteps`is used as a sanity check, and expects you to provide the amount of simulation timesteps.
+
+!!! note
+    Only subtypes of AbstractVector present in statuses will be affected. In some cases, meteo values might need a small conversion. For instance :
+    ```
+    meteo_day = CSV.read(joinpath(pkgdir(PlantSimEngine), "examples/meteo_day.csv"), DataFrame, header=18)
+    status(TT_cu=cumsum(meteo_day.TT),)```
+
+    cumsum(meteo_day.TT) actually returns a CSV.SentinelArray.ChainedVectors{T, Vector{T}}, which is not a subtype of AbstractVector. 
+    Replacing it with Vector(cumsum(meteo_day.TT)) will provide an adequate type.
+
+This feature is likely to break in simulations that make use of planned future features (such as mixing models with different timesteps), without guarantee of a fix on a short notice. Again, bear in mind it is mostly a convenient shortcut for prototyping, when doing multi-scale simulations.

src/PlantSimEngine.jl

@@ -22,6 +22,9 @@ import MultiScaleTreeGraph: symbol, node_id
 # To compute mean:
 import Statistics
 
+# For avoiding name conflicts when generating models from status vectors
+import SHA: sha1
+
 using PlantMeteo
 
 # UninitializedVar + PreviousTimeStep:
@@ -89,6 +92,9 @@ include("run.jl")
 # Fitting
 include("evaluation/fit.jl")
 
+# Utilities for mapping initialisation
+include("mtg/mapping/model_generation_from_status_vectors.jl")
+
 # Examples
 include("examples_import.jl")
 

src/mtg/initialisation.jl

@@ -307,6 +307,13 @@ automatically passed as is.
 """
 function init_simulation(mtg, mapping; nsteps=1, outputs=nothing, type_promotion=nothing, check=true, verbose=false)
 
+    # Ensure the user called the model generation function to handle vectors passed into a status
+    # before we keep going
+    (organ_with_vector, no_vectors_found) = (check_statuses_contain_no_remaining_vectors(mapping))
+    if !no_vectors_found
+        @assert false "Error : Mapping status at $organ_with_vector level contains a vector. If this was intentional, call the function generate_models_from_status_vectors on your mapping before calling run!. And bear in mind this is not meant for production. If this wasn't intentional, then it's likely an issue on the mapping definition, or an unusual model."
+    end
+
     # Get the status of each node by node type, pre-initialised considering multi-scale variables:
     statuses, status_templates, reverse_multiscale_mapping, vars_need_init =
         init_statuses(mtg, mapping, first(hard_dependencies(mapping; verbose=false)); type_promotion=type_promotion, verbose=verbose, check=check)