Merge pull request #66 from VirtualPlantLab/61-add-a-way-to-not-use-a…

…n-input-variable-to-build-the-dependency-graph

61 add a way to not use an input variable to build the dependency graph

README.md

@@ -225,15 +225,15 @@ mapping = Dict(
             mapping=[:TT_cu => "Scene"],
         ),
         ToyMaintenanceRespirationModel(1.5, 0.06, 25.0, 0.6, 0.004),
-        Status(biomass=1.0)
+        Status(carbon_biomass=1.0)
     ),
     "Leaf" => (
         MultiScaleModel(
             model=ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
             mapping=[:TT => "Scene",],
         ),
         ToyMaintenanceRespirationModel(2.1, 0.06, 25.0, 1.0, 0.025),
-        Status(biomass=1.0)
+        Status(carbon_biomass=1.0)
     ),
     "Soil" => (
         ToySoilWaterModel(),

docs/make.jl

@@ -16,7 +16,7 @@ makedocs(;
         canonical="https://VirtualPlantLab.github.io/PlantSimEngine.jl",
         edit_link="main",
         assets=String[],
-        size_threshold=250000
+        size_threshold=300000
     ),
     pages=[
         "Home" => "index.md",

docs/src/extending/implement_a_process.md

@@ -68,15 +68,15 @@ Computes the biomass growth of a plant.
 
 # Inputs
 
-- `aPPFD`: the absorbed photosynthetic photon flux density, in mol[PAR] m⁻² d⁻¹
+- `aPPFD`: the absorbed photosynthetic photon flux density, in mol[PAR] m⁻² time-step⁻¹
 
 # Outputs
 
-- `carbon_assimilation`: the assimilation, in gC m⁻² d⁻¹
-- `Rm`: the maintenance respiration, in gC m⁻² d⁻¹
-- `Rg`: the growth respiration, in gC m⁻² d⁻¹
-- `biomass_increment`: the daily biomass increment, in gC m⁻² d⁻¹
-- `biomass`: the plant biomass, in gC m⁻² d⁻¹
+- `carbon_assimilation`: the assimilation, in gC m⁻² time-step⁻¹
+- `Rm`: the maintenance respiration, in gC m⁻² time-step⁻¹
+- `Rg`: the growth respiration, in gC m⁻² time-step⁻¹
+- `biomass_increment`: the daily biomass increment, in gC m⁻² time-step⁻¹
+- `biomass`: the plant biomass, in gC m⁻² time-step⁻¹
 """
 struct ToyAssimGrowthModel{T} <: AbstractGrowthModel
     LUE::T

docs/src/index.md

@@ -217,15 +217,15 @@ mapping = Dict(
             mapping=[:TT_cu => "Scene"],
         ),
         ToyMaintenanceRespirationModel(1.5, 0.06, 25.0, 0.6, 0.004),
-        Status(biomass=1.0)
+        Status(carbon_biomass=1.0)
     ),
     "Leaf" => (
         MultiScaleModel(
             model=ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
             mapping=[:TT => "Scene",],
         ),
         ToyMaintenanceRespirationModel(2.1, 0.06, 25.0, 1.0, 0.025),
-        Status(biomass=1.0)
+        Status(carbon_biomass=1.0)
     ),
     "Soil" => (
         ToySoilWaterModel(),

docs/src/model_coupling/multiscale.md

@@ -132,7 +132,7 @@ mapping = Dict(
             mapping=[:TT => "Scene",],
         ),
         ToyMaintenanceRespirationModel(1.5, 0.06, 25.0, 0.6, 0.004),
-        Status(biomass=1.0),
+        Status(carbon_biomass=1.0),
     ),
     "Leaf" => (
         MultiScaleModel(
@@ -144,7 +144,7 @@ mapping = Dict(
             mapping=[:TT => "Scene",],
         ),
         ToyMaintenanceRespirationModel(2.1, 0.06, 25.0, 1.0, 0.025),
-        Status(biomass=0.5),
+        Status(carbon_biomass=0.5),
     ),
     "Soil" => (
         ToySoilWaterModel(),
@@ -248,8 +248,125 @@ nothing # hide
 
 This is a dictionary with the scale as the key and a vector of `Status` as values, one per node of that scale. So, in this example, the `"Leaf"` scale has two nodes, so the value is a vector of two `Status` objects, and the `"Soil"` scale has only one node, so the value is a vector of one `Status` object.
 
+
+## Avoiding cyclic dependencies
+
+When defining a mapping between models and scales, it is important to avoid cyclic dependencies. A cyclic dependency occurs when a model at a given scale depends on a model at another scale that depends on the first model. Cyclic dependencies are bad because they lead to an infinite loop in the simulation (the dependency graph keeps cycling indefinitely).
+
+PlantSimEngine will detect cyclic dependencies and raise an error if one is found. The error message indicates the models involved in the cycle, and the model that is causing the cycle will be highlighted in red.
+
+For example the following mapping will raise an error:
+
+!!! details
+    <summary>Example mapping</summary>
+
+    ```julia
+    mapping_cyclic = Dict(
+        "Plant" => (
+            MultiScaleModel(
+                model=ToyCAllocationModel(),
+                mapping=[
+                    :carbon_demand => ["Leaf", "Internode"],
+                    :carbon_allocation => ["Leaf", "Internode"]
+                ],
+            ),
+            MultiScaleModel(
+                model=ToyPlantRmModel(),
+                mapping=[:Rm_organs => ["Leaf" => :Rm, "Internode" => :Rm],],
+            ),
+            Status(total_surface=0.001, aPPFD=1300.0, soil_water_content=0.6),
+        ),
+        "Internode" => (
+            ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
+            ToyMaintenanceRespirationModel(1.5, 0.06, 25.0, 0.6, 0.004),
+            Status(TT=10.0, carbon_biomass=1.0),
+        ),
+        "Leaf" => (
+            ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
+            ToyMaintenanceRespirationModel(2.1, 0.06, 25.0, 1.0, 0.025),
+            ToyCBiomassModel(1.2),
+            Status(TT=10.0),
+        )
+    )
+    ```
+
+Let's see what happens when we try to build the dependency graph for this mapping:
+
+```julia
+julia> dep(mapping_cyclic)
+ERROR: Cyclic dependency detected in the graph. Cycle:
+ Plant: ToyPlantRmModel
+ └ Leaf: ToyMaintenanceRespirationModel
+  └ Leaf: ToyCBiomassModel
+   └ Plant: ToyCAllocationModel
+    └ Plant: ToyPlantRmModel
+
+ You can break the cycle using the `PreviousTimeStep` variable in the mapping.
+```
+
+How can we interpret the message? We have a list of five models involved in the cycle. The first model is the one causing the cycle, and the others are the ones that depend on it. In this case, the `ToyPlantRmModel` is the one causing the cycle, and the others are inter-dependent. We can read this as follows:
+
+1. `ToyPlantRmModel` depends on `ToyMaintenanceRespirationModel`, the plant-scale respiration sums up all organs respiration;
+2. `ToyMaintenanceRespirationModel` depends on `ToyCBiomassModel`, the organs respiration depends on the organs biomass;
+3. `ToyCBiomassModel` depends on `ToyCAllocationModel`, the organs biomass depends on the organs carbon allocation;
+4. And finally `ToyCAllocationModel` depends on `ToyPlantRmModel` again, hence the cycle because the carbon allocation depends on the plant scale respiration.
+
+The models can not be ordered in a way that satisfies all dependencies, so the cycle can not be broken. To solve this issue, we need to re-think how models are mapped together, and break the cycle.
+
+There are several ways to break a cyclic dependency:
+
+- **Merge models**: If two models depend on each other because they need *e.g.* recursive computations, they can be merged into a third model that handles the computation and takes the two models as hard dependencies. Hard dependencies are models that are explicitly called by another model and do not participate on the building of the dependency graph.
+- **Change models**: Of course models can be interchanged to avoid cyclic dependencies, but this is not really a solution, it is more a workaround.
+- **PreviousTimeScale**: We can break the dependency graph by defining some variables as taken from the previous time step. A very well known example is the computation of the light interception by a plant that depends on the leaf area, which is usually the result of a model that also depends on the light interception. The cyclic dependency is usually broken by using the leaf area from the previous time step in the interception model, which is a good approximation for most cases.
+
+We can fix our previous mapping by computing the organs respiration using the carbon biomass from the previous time step instead. Let's see how to fix the cyclic dependency in our mapping (look at the leaf and internode scales):
+
+!!! details
+    ```@julia
+    mapping_nocyclic = Dict(
+            "Plant" => (
+                MultiScaleModel(
+                    model=ToyCAllocationModel(),
+                    mapping=[
+                        :carbon_demand => ["Leaf", "Internode"],
+                        :carbon_allocation => ["Leaf", "Internode"]
+                    ],
+                ),
+                MultiScaleModel(
+                    model=ToyPlantRmModel(),
+                    mapping=[:Rm_organs => ["Leaf" => :Rm, "Internode" => :Rm],],
+                ),
+                Status(total_surface=0.001, aPPFD=1300.0, soil_water_content=0.6, carbon_assimilation=5.0),
+            ),
+            "Internode" => (
+                ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
+                MultiScaleModel(
+                    model=ToyMaintenanceRespirationModel(1.5, 0.06, 25.0, 0.6, 0.004),
+                    mapping=[PreviousTimeStep(:carbon_biomass),], #! this is where we break the cyclic dependency (first break)
+                ),
+                Status(TT=10.0, carbon_biomass=1.0),
+            ),
+            "Leaf" => (
+                ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
+                MultiScaleModel(
+                    model=ToyMaintenanceRespirationModel(2.1, 0.06, 25.0, 1.0, 0.025),
+                    mapping=[PreviousTimeStep(:carbon_biomass),], #! this is where we break the cyclic dependency (second break)
+                ),
+                ToyCBiomassModel(1.2),
+                Status(TT=10.0),
+            )
+        );
+    nothing # hide
+    ```
+
+The `ToyMaintenanceRespirationModel` models are now defined as `MultiScaleModel`, and the `carbon_biomass` variable is wrapped in a `PreviousTimeStep` structure. This structure tells PlantSimEngine to take the value of the variable from the previous time step, breaking the cyclic dependency.
+
+!!! note
+    `PreviousTimeStep` tells PlantSimEngine to take the value of the previous time step for the variable it wraps, or the value at initialization for the first time step. The value at initialization is the one provided by default in the models inputs, but is usually provided in the `Status` structure to override this default.
+    A `PreviousTimeStep` is used to wrap the **input** variable of a model, with or without a mapping to another scale *e.g.* `Previous(:carbon_biomass) => "Leaf"`.
+
 ### Wrapping up
 
 In this section, we saw how to define a mapping between models and scales, run a simulation, and access the outputs.
 
-This is just a simple example, but PlantSimEngine can be used to define and combine much more complex models at multiple scales of detail. With its modular architecture and intuitive API, PlantSimEngine is a powerful tool for multi-scale plant growth and development modeling.
+This is just a simple example, but PlantSimEngine can be used to define and combine much more complex models at multiple scales of detail. With its modular architecture and intuitive API, PlantSimEngine is a powerful tool for multi-scale plant growth and development modeling.

examples/ToyAssimGrowthModel.jl

@@ -16,15 +16,15 @@ Computes the biomass growth of a plant.
 
 # Inputs
 
-- `aPPFD`: the absorbed photosynthetic photon flux density, in mol[PAR] m⁻² d⁻¹
+- `aPPFD`: the absorbed photosynthetic photon flux density, in mol[PAR] m⁻² time-step⁻¹
 
 # Outputs
 
-- `carbon_assimilation`: the assimilation, in gC m⁻² d⁻¹
-- `Rm`: the maintenance respiration, in gC m⁻² d⁻¹
-- `Rg`: the growth respiration, in gC m⁻² d⁻¹
-- `biomass_increment`: the daily biomass increment, in gC m⁻² d⁻¹
-- `biomass`: the plant biomass, in gC m⁻² d⁻¹
+- `carbon_assimilation`: the assimilation, in gC m⁻² time-step⁻¹
+- `Rm`: the maintenance respiration, in gC m⁻² time-step⁻¹
+- `Rg`: the growth respiration, in gC m⁻² time-step⁻¹
+- `biomass_increment`: the daily biomass increment, in gC m⁻² time-step⁻¹
+- `biomass`: the plant biomass, in gC m⁻² time-step⁻¹
 """
 struct ToyAssimGrowthModel{T} <: AbstractGrowthModel
     LUE::T

examples/ToyAssimModel.jl

@@ -15,18 +15,18 @@ Computes the assimilation of a plant (= photosynthesis).
 
 # Inputs
 
-- `aPPFD`: the absorbed photosynthetic photon flux density, in mol[PAR] m⁻² d⁻¹
+- `aPPFD`: the absorbed photosynthetic photon flux density, in mol[PAR] m⁻² time-step⁻¹
 - `soil_water_content`: the soil water content, in %
 
 # Outputs
 
-- `carbon_assimilation`: the assimilation or photosynthesis, also sometimes denoted `A`, in gC m⁻² d⁻¹
+- `carbon_assimilation`: the assimilation or photosynthesis, also sometimes denoted `A`, in gC m⁻² time-step⁻¹
 
 # Details
 
 The assimilation is computed as the product of the absorbed photosynthetic photon flux density (aPPFD) and the light use efficiency (LUE),
-so the units of the assimilation usually are in gC m⁻² d⁻¹, but they could be in another spatial or temporal unit depending on the unit of `aPPFD`, *e.g.* 
-if `aPPFD` is in mol[PAR] plant⁻¹ d⁻¹, the assimilation will be in gC plant⁻¹ d⁻¹.
+so the units of the assimilation usually are in gC m⁻² time-step⁻¹, but they could be in another spatial or temporal unit depending on the unit of `aPPFD`, *e.g.* 
+if `aPPFD` is in mol[PAR] plant⁻¹ time-step⁻¹, the assimilation will be in gC plant⁻¹ time-step⁻¹.
 """
 struct ToyAssimModel{T} <: AbstractCarbon_AssimilationModel
     LUE::T
@@ -50,7 +50,7 @@ end
 # Tells Julia what is the type of elements:
 Base.eltype(::ToyAssimModel{T}) where {T} = T
 
-# Implement the growth model:
+# Implement the model:
 function PlantSimEngine.run!(::ToyAssimModel, models, status, meteo, constants, extra)
     # The assimilation is simply the absorbed photosynthetic photon flux density (aPPFD) times the light use efficiency (LUE):
     status.carbon_assimilation = status.aPPFD * models.carbon_assimilation.LUE * status.soil_water_content

examples/ToyCAllocationModel.jl

@@ -14,18 +14,18 @@ to each organ's demand.
 
 # Inputs
 
-- `carbon_assimilation`: a vector of the assimilation of all photosynthetic organs, usually in gC m⁻² d⁻¹
-- `Rm`: the maintenance respiration of the plant, usually in gC m⁻² d⁻¹
-- `carbon_demand`: a vector of the carbon demand of the organs, usually in gC m⁻² d⁻¹
+- `carbon_assimilation`: a vector of the assimilation of all photosynthetic organs, usually in gC m⁻² time-step⁻¹
+- `Rm`: the maintenance respiration of the plant, usually in gC m⁻² time-step⁻¹
+- `carbon_demand`: a vector of the carbon demand of the organs, usually in gC m⁻² time-step⁻¹
 
 # Outputs
 
-- `carbon_assimilation`: the carbon assimilation, usually in gC m⁻² d⁻¹
+- `carbon_assimilation`: the carbon assimilation, usually in gC m⁻² time-step⁻¹
 
 # Details
 
-The units usually are in gC m⁻² d⁻¹, but they could be in another spatial or temporal unit depending on the unit of the inputs, *e.g.*
-in gC plant⁻¹ d⁻¹.
+The units usually are in gC m⁻² time-step⁻¹, but they could be in another spatial or temporal unit depending on the unit of the inputs, *e.g.*
+in gC plant⁻¹ time-step⁻¹.
 """
 struct ToyCAllocationModel <: AbstractCarbon_AllocationModel end
 
@@ -39,7 +39,7 @@ function PlantSimEngine.outputs_(::ToyCAllocationModel)
     (carbon_offer=-Inf, carbon_allocation=[-Inf],)
 end
 
-function PlantSimEngine.run!(::ToyCAllocationModel, models, status, meteo, constants, mtg)
+function PlantSimEngine.run!(::ToyCAllocationModel, models, status, meteo, constants, extra_args)
 
     carbon_demand_tot = sum(status.carbon_demand)
     #Note: this model is multiscale, so status.carbon_demand, status.carbon_allocation, and status.carbon_assimilation are vectors.
@@ -64,5 +64,5 @@ function PlantSimEngine.run!(::ToyCAllocationModel, models, status, meteo, const
     end
 end
 
-# And also over time (time-steps):
+# Can be parallelized over time-steps, but not objects (we have vectors of values coming from other objects as input):
 PlantSimEngine.TimeStepDependencyTrait(::Type{<:ToyCAllocationModel}) = PlantSimEngine.IsTimeStepIndependent()

examples/ToyCBiomassModel.jl

@@ -0,0 +1,46 @@
+# using PlantSimEngine, PlantMeteo # Import the necessary packages, PlantMeteo is used for the meteorology
+
+PlantSimEngine.@process "carbon_biomass" verbose = false
+
+"""
+    ToyCBiomassModel(construction_cost)
+
+Computes the carbon biomass of an organ based on the carbon allocation and construction cost.
+
+# Arguments
+
+- `construction_cost`: the construction cost of the organ, usually in gC gC⁻¹. Should be understood as the amount of carbon needed to build 1g of carbon biomass.
+
+# Inputs
+
+- `carbon_allocation`: the carbon allocation to the organ for the time-step, usually in gC m⁻² time-step⁻¹
+
+# Outputs
+
+- `carbon_biomass_increment`: the increment of carbon biomass, usually in gC time-step⁻¹
+- `carbon_biomass`: the carbon biomass, usually in gC
+- `growth_respiration`: the growth respiration, usually in gC time-step⁻¹
+
+"""
+struct ToyCBiomassModel{T} <: AbstractCarbon_BiomassModel
+    construction_cost::T
+end
+
+# Define inputs:
+function PlantSimEngine.inputs_(::ToyCBiomassModel)
+    (carbon_allocation=-Inf,)
+end
+
+# Define outputs:
+function PlantSimEngine.outputs_(::ToyCBiomassModel)
+    (carbon_biomass_increment=-Inf, carbon_biomass=0.0, growth_respiration=-Inf,)
+end
+
+function PlantSimEngine.run!(m::ToyCBiomassModel, models, status, meteo, constants, extra_args)
+    status.carbon_biomass_increment = status.carbon_allocation / m.construction_cost
+    status.carbon_biomass += status.carbon_biomass_increment
+    status.growth_respiration = status.carbon_allocation - status.carbon_biomass_increment
+end
+
+# Can be parallelized over organs (but not time-steps, as it is incrementally updating the biomass in the status):
+PlantSimEngine.ObjectDependencyTrait(::Type{<:ToyCBiomassModel}) = PlantSimEngine.IsObjectIndependent()