From 908c168fad98666ebda367006dc799708ea22722 Mon Sep 17 00:00:00 2001 From: Alejandro Morales Sierra Date: Sat, 11 Nov 2023 16:56:50 +0100 Subject: [PATCH] Adapt to Literate workflow --- test/algae.jl | 229 ++++---- test/context.jl | 139 ++--- test/forest.jl | 567 ++++++++++---------- test/growthforest.jl | 790 ++++++++++++++------------- test/raytracedforest.jl | 1058 ++++++++++++++++++------------------- test/relationalqueries.jl | 585 ++++++++++---------- test/runtests.jl | 32 +- test/snowflakes.jl | 314 +++++------ test/tree.jl | 420 ++++++--------- 9 files changed, 2019 insertions(+), 2115 deletions(-) diff --git a/test/algae.jl b/test/algae.jl index edd80f1..56ef2ad 100644 --- a/test/algae.jl +++ b/test/algae.jl @@ -1,7 +1,9 @@ #= -Algae growth +# Algae growth -Alejandro Morales - May 2023 +Alejandro Morales + +Centre for Crop Systems Analysis - Wageningen University In this first example, we learn how to create a `Graph` and update it dynamically with rewriting rules. @@ -10,12 +12,13 @@ The model described here is based on the non-branching model of [algae growth](https://en.wikipedia.org/wiki/L-system#Example_1:_Algae) proposed by Lindermayer as one of the first L-systems. -First, we need to load the VirtualPlantLab metapackage, which will automatically load all -the packages in the VirtualPlantLab ecosystem. +First, we need to load the VPL metapackage, which will automatically load all +the packages in the VPL ecosystem. + =# using VirtualPlantLab - #= + The rewriting rules of the L-system are as follows: **axiom**: A @@ -24,115 +27,127 @@ The rewriting rules of the L-system are as follows: **rule 2**: B $\rightarrow$ A -In VirtualPlantLab, this L-system would be implemented as a graph where the nodes can be of +In VPL, this L-system would be implemented as a graph where the nodes can be of type `A` or `B` and inherit from the abstract type `Node`. It is advised to include type definitions in a module to avoid having to restart the Julia session whenever we want to redefine them. Because each module is an independent -namespace, we need to import `Node` from the VirtualPlantLab package inside the module: +namespace, we need to import `Node` from the VPL package inside the module: + =# module algae -import VirtualPlantLab: Node -struct A <: Node end -struct B <: Node end + import VirtualPlantLab: Node + struct A <: Node end + struct B <: Node end end import .algae +#= -let - #= - Note that in this very example we do not need to store any data or state inside - the nodes, so types `A` and `B` do not require fields. - - The axiom is simply defined as an instance of type of `A`: - =# - axiom = algae.A() - - #= - The rewriting rules are implemented in VirtualPlantLab as objects of type `Rule`. In VirtualPlantLab, a - rewriting rule substitutes a node in a graph with a new node or subgraph and is - therefore composed of two parts: - - 1. A condition that is tested against each node in a graph to choose which nodes - to rewrite. - 2. A subgraph that will replace each node selected by the condition above. - - In VirtualPlantLab, the condition is split into two components: - - 1. The type of node to be selected (in this example that would be `A` or `B`). - 2. A function that is applied to each node in the graph (of the specified type) - to indicate whether the node should be selected or not. This function is - optional (the default is to select every node of the specified type). - - The replacement subgraph is specified by a function that takes as input the node - selected and returns a subgraph defined as a combination of node objects. - Subgraphs (which can also be used as axioms) are created by linearly combining - objects that inherit from `Node`. The operation `+` implies a linear - relationship between two nodes and `[]` indicates branching. - - The implementation of the two rules of algae growth model in VirtualPlantLab is as follows: - =# - rule1 = Rule(algae.A, rhs = x -> algae.A() + algae.B()) - rule2 = Rule(algae.B, rhs = x -> algae.A()) - - #= - Note that in each case, the argument `rhs` is being assigned an anonymous (aka - *lambda*) function. This is a function without a name that is defined directly - in the assigment to the argument. That is, the Julia expression `x -> A() + B()` - is equivalent to the following function definition: - =# - function rule_1(x) - algae.A() + algae.B() - end - - #= - For simple rules (especially if the right hand side is just a line of code) it - is easier to just define the right hand side of the rule with an anonymous - function rather than creating a standalone function with a meaningful name. - However, standalone functions are easier to debug as you can call them directly - from the REPL. - - With the axiom and rules we can now create a `Graph` object that represents the - algae organism. The first argument is the axiom and the second is a tuple with - all the rewriting rules: - =# - organism = Graph(axiom = axiom, rules = (rule1, rule2)) - - #= - If we apply the rewriting rules iteratively, the graph will grow, in this case - representing the growth of the algae organism. The rewriting rules are applied - on the graph with the function `rewrite!()`: - =# - rewrite!(organism) +Note that in this very example we do not need to store any data or state inside +the nodes, so types `A` and `B` do not require fields. + +The axiom is simply defined as an instance of type of `A`: + +=# +axiom = algae.A() +#= + +The rewriting rules are implemented in VPL as objects of type `Rule`. In VPL, a +rewriting rule substitutes a node in a graph with a new node or subgraph and is +therefore composed of two parts: + +1. A condition that is tested against each node in a graph to choose which nodes + to rewrite. +2. A subgraph that will replace each node selected by the condition above. - #= - Since there was only one node of type `A`, the only rule that was applied was - `rule1`, so the graph should now have two nodes of types `A` and `B`, - respectively. We can confirm this by drawing the graph. We do this with the - function `draw()` which will always generate the same representation of the - graph, but different options are available depending on the context where the - code is executed. By default, `draw()` will create a new window where an - interactive version of the graph will be drawn and one can zoom and pan with the - mouse. - =# - import GLMakie - draw(organism) - - #= - Notice that each node in the network representation is labelled with the type of - node (`A` or `B` in this case) and a number in parenthesis. This number is a - unique identifier associated to each node and it is useful for debugging - purposes (this will be explained in more advanced examples). - - Applying multiple iterations of rewriting can be achieved with a simple loop: - =# - for i = 1:4 - rewrite!(organism) - end - - # And we can verify that the graph grew as expected: - draw(organism) - - # The network is rather boring as the system is growing linearly (no branching) - # but it already illustrates how graphs can grow rapidly in just a few - # iterations. Remember that the interactive visualization allows adjusting the - # zoom, which is handy when graphs become large. +In VPL, the condition is split into two components: + +1. The type of node to be selected (in this example that would be `A` or `B`). +2. A function that is applied to each node in the graph (of the specified type) + to indicate whether the node should be selected or not. This function is + optional (the default is to select every node of the specified type). + +The replacement subgraph is specified by a function that takes as input the node +selected and returns a subgraph defined as a combination of node objects. +Subgraphs (which can also be used as axioms) are created by linearly combining +objects that inherit from `Node`. The operation `+` implies a linear +relationship between two nodes and `[]` indicates branching. + +The implementation of the two rules of algae growth model in VPL is as follows: + +=# +rule1 = Rule(algae.A, rhs = x -> algae.A() + algae.B()) +rule2 = Rule(algae.B, rhs = x -> algae.A()) +#= + +Note that in each case, the argument `rhs` is being assigned an anonymous (aka +*lambda*) function. This is a function without a name that is defined directly +in the assigment to the argument. That is, the Julia expression `x -> A() + B()` +is equivalent to the following function definition: + +=# +function rule_1(x) + algae.A() + algae.B() end +#= + +For simple rules (especially if the right hand side is just a line of code) it +is easier to just define the right hand side of the rule with an anonymous +function rather than creating a standalone function with a meaningful name. +However, standalone functions are easier to debug as you can call them directly +from the REPL. + +With the axiom and rules we can now create a `Graph` object that represents the +algae organism. The first argument is the axiom and the second is a tuple with +all the rewriting rules: + +=# +organism = Graph(axiom = axiom, rules = (rule1, rule2)) +#= + +If we apply the rewriting rules iteratively, the graph will grow, in this case +representing the growth of the algae organism. The rewriting rules are applied +on the graph with the function `rewrite!()`: + +=# +rewrite!(organism) +#= + +Since there was only one node of type `A`, the only rule that was applied was +`rule1`, so the graph should now have two nodes of types `A` and `B`, +respectively. We can confirm this by drawing the graph. We do this with the +function `draw()` which will always generate the same representation of the +graph, but different options are available depending on the context where the +code is executed. By default, `draw()` will create a new window where an +interactive version of the graph will be drawn and one can zoom and pan with the +mouse (in this online document a static version is shown, see +[Backends](../manual/Visualization.md) for details): + +=# +import GLMakie +draw(organism) +#= + +Notice that each node in the network representation is labelled with the type of +node (`A` or `B` in this case) and a number in parenthesis. This number is a +unique identifier associated to each node and it is useful for debugging +purposes (this will be explained in more advanced examples). + +Applying multiple iterations of rewriting can be achieved with a simple loop: + +=# +for i in 1:4 + rewrite!(organism) +end +#= + +And we can verify that the graph grew as expected: + +=# +draw(organism) +#= + +The network is rather boring as the system is growing linearly (no branching) +but it already illustrates how graphs can grow rapidly in just a few iterations. +Remember that the interactive visualization allows adjusting the zoom, which is +handy when graphs become large. +=# diff --git a/test/context.jl b/test/context.jl index 06ba461..ad8b060 100644 --- a/test/context.jl +++ b/test/context.jl @@ -1,8 +1,9 @@ #= +# Context sensitive rules -Context sensitive rules +Alejandro Morales -Alejandro Morales - May 2023 +Centre for Crop Systems Analysis - Wageningen University This examples goes back to a very simple situation: a linear sequence of 3 cells. The point of this example is to introduce relational growth rules and @@ -18,7 +19,7 @@ also to use properties of those neighbours in the right hand side component of the rule. This is know as "capturing the context" of the node being updated. This can be done by returning the additional nodes from the `lhs` component (in addition to `true` or `false`) and by accepting these additional nodes in the -`rhs` component. In addition, we tell VirtualPlantLab that this rule is capturing the +`rhs` component. In addition, we tell VPL that this rule is capturing the context with `captures = true`. In the example below, each `Cell` keeps track of a `state` variable (which is @@ -29,75 +30,77 @@ the right hand side, we replace the cell with a new cell with the state of the parent node that was captured. Note that that now, the rhs component gets a new argument, which corresponds to the context of the father node captured in the lhs. + =# using VirtualPlantLab module types -using VirtualPlantLab -struct Cell <: Node - state::Int64 -end + using VirtualPlantLab + struct Cell <: Node + state::Int64 + end end import .types: Cell - -let - function transfer(context) - if has_parent(context) - return (true, (parent(context),)) - else - return (false, ()) - end +function transfer(context) + if has_parent(context) + return (true, (parent(context), )) + else + return (false, ()) end - rule = Rule( - Cell, - lhs = transfer, - rhs = (context, father) -> Cell(data(father).state), - captures = true, - ) - axiom = Cell(1) + Cell(0) + Cell(0) - pop = Graph(axiom = axiom, rules = rule) - - #= - In the original state defined by the axiom, only the first node contains a state - of 1. We can retrieve the state of each node with a query. A `Query` object is a - like a `Rule` but without a right-hand side (i.e., its purpose is to return the - nodes that match a particular condition). In this case, we just want to return - all the `Cell` nodes. A `Query` object is created by passing the type of the - node to be queried as an argument to the `Query` function. Then, to actually - execute the query we need to use the `apply` function on the graph. - =# - getCell = Query(Cell) - apply(pop, getCell) - - # If we rewrite the graph one we will see that a second cell now has a state of 1. - rewrite!(pop) - apply(pop, getCell) - - # And a second iteration results in all cells have a state of 1 - rewrite!(pop) - apply(pop, getCell) - - #= - Note that queries may not return nodes in the same order as they were created - because of how they are internally stored (and because queries are meant to - return collection of nodes rather than reconstruct the topology of a graph). If - we need to process nodes in a particular order, then it is best to use a - traversal algorithm on the graph that follows a particular order (for example - depth-first traversal with `traverse_dfs()`). This algorithm requires a function - that applies to each node in the graph. In this simple example we can just store - the `state` of each node in a vector (unlike Rules and Queries, this function - takes the actual node as argument rather than a `Context` object, see the - documentation for more details): - =# - - pop = Graph(axiom = axiom, rules = rule) - states = Int64[] - traverse_dfs(pop, fun = node -> push!(states, node.state)) - states - - # Now the states of the nodes are in the same order as they were created: - rewrite!(pop) - states = Int64[] - traverse_dfs(pop, fun = node -> push!(states, node.state)) - states - end +rule = Rule(Cell, lhs = transfer, rhs = (context, father) -> Cell(data(father).state), captures = true) +axiom = Cell(1) + Cell(0) + Cell(0) +pop = Graph(axiom = axiom, rules = rule) +#= + +In the original state defined by the axiom, only the first node contains a state +of 1. We can retrieve the state of each node with a query. A `Query` object is a +like a `Rule` but without a right-hand side (i.e., its purpose is to return the +nodes that match a particular condition). In this case, we just want to return +all the `Cell` nodes. A `Query` object is created by passing the type of the +node to be queried as an argument to the `Query` function. Then, to actually +execute the query we need to use the `apply` function on the graph. + +=# +getCell = Query(Cell) +apply(pop, getCell) +#= + +If we rewrite the graph one we will see that a second cell now has a state of 1. + +=# +rewrite!(pop) +apply(pop, getCell) +#= + +And a second iteration results in all cells have a state of 1 + +=# +rewrite!(pop) +apply(pop, getCell) +#= + +Note that queries may not return nodes in the same order as they were created +because of how they are internally stored (and because queries are meant to +return collection of nodes rather than reconstruct the topology of a graph). If +we need to process nodes in a particular order, then it is best to use a +traversal algorithm on the graph that follows a particular order (for example +depth-first traversal with `traverse_dfs()`). This algorithm requires a function +that applies to each node in the graph. In this simple example we can just store +the `state` of each node in a vector (unlike Rules and Queries, this function +takes the actual node as argument rather than a `Context` object, see the +documentation for more details): + +=# +pop = Graph(axiom = axiom, rules = rule) +states = Int64[] +traverse_dfs(pop, fun = node -> push!(states, node.state)) +states +#= + +Now the states of the nodes are in the same order as they were created: + +=# +rewrite!(pop) +states = Int64[] +traverse_dfs(pop, fun = node -> push!(states, node.state)) +states diff --git a/test/forest.jl b/test/forest.jl index 71eed84..c409274 100644 --- a/test/forest.jl +++ b/test/forest.jl @@ -1,324 +1,301 @@ #= - Forest Alejandro Morales Sierra Centre for Crop Systems Analysis - Wageningen -University +# Forest -In this example we extend the tree example into a forest, where each tree is -described by a separate graph object and parameters driving the growth of these -trees vary across individuals following a predefined distribution. This tutorial -requires using the Distributions.jl package: +Alejandro Morales +Centre for Crop Systems Analysis - Wageningen University + +In this example we extend the tree example into a forest, where +each tree is described by a separate graph object and parameters driving the +growth of these trees vary across individuals following a predefined distribution. The data types, rendering methods and growth rules are the same as in the binary tree example: + =# using VirtualPlantLab -using Distributions, ColorTypes +using Distributions, Plots, ColorTypes import GLMakie -import Base.Threads: @threads -# Data types +## Data types module TreeTypes -import VirtualPlantLab -# Meristem -struct Meristem <: VirtualPlantLab.Node end -# Bud -struct Bud <: VirtualPlantLab.Node end -# Node -struct Node <: VirtualPlantLab.Node end -# BudNode -struct BudNode <: VirtualPlantLab.Node end -# Internode (needs to be mutable to allow for changes over time) -Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node - length::Float64 = 0.10 # Internodes start at 10 cm + import VirtualPlantLab + ## Meristem + struct Meristem <: VirtualPlantLab.Node end + ## Bud + struct Bud <: VirtualPlantLab.Node end + ## Node + struct Node <: VirtualPlantLab.Node end + ## BudNode + struct BudNode <: VirtualPlantLab.Node end + ## Internode (needs to be mutable to allow for changes over time) + Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node + length::Float64 = 0.10 # Internodes start at 10 cm + end + ## Leaf + Base.@kwdef struct Leaf <: VirtualPlantLab.Node + length::Float64 = 0.20 # Leaves are 20 cm long + width::Float64 = 0.1 # Leaves are 10 cm wide + end + ## Graph-level variables + Base.@kwdef struct treeparams + growth::Float64 = 0.1 + budbreak::Float64 = 0.25 + phyllotaxis::Float64 = 140.0 + leaf_angle::Float64 = 30.0 + branch_angle::Float64 = 45.0 + end end -# Leaf -Base.@kwdef struct Leaf <: VirtualPlantLab.Node - length::Float64 = 0.20 # Leaves are 20 cm long - width::Float64 = 0.1 # Leaves are 10 cm wide + +import .TreeTypes + +# Create geometry + color for the internodes +function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, data) + ## Rotate turtle around the head to implement elliptical phyllotaxis + rh!(turtle, data.phyllotaxis) + HollowCylinder!(turtle, length = i.length, height = i.length/15, width = i.length/15, + move = true, color = RGB(0.5,0.4,0.0)) + return nothing end -# Graph-level variables -Base.@kwdef struct treeparams - growth::Float64 = 0.1 - budbreak::Float64 = 0.25 - phyllotaxis::Float64 = 140.0 - leaf_angle::Float64 = 30.0 - branch_angle::Float64 = 45.0 + +# Create geometry + color for the leaves +function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, data) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -data.leaf_angle) + ## Generate the leaf + Ellipse!(turtle, length = l.length, width = l.width, move = false, + color = RGB(0.2,0.6,0.2)) + ## Rotate turtle back to original direction + ra!(turtle, data.leaf_angle) + return nothing end + +# Insertion angle for the bud nodes +function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, data) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -data.branch_angle) end -import .TreeTypes -let - # Create geometry + color for the internodes - function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, data) - # Rotate turtle around the head to implement elliptical phyllotaxis - rh!(turtle, data.phyllotaxis) - HollowCylinder!( - turtle, - length = i.length, - height = i.length / 15, - width = i.length / 15, - move = true, - color = RGB(0.5, 0.4, 0.0), - ) - return nothing +# Rules +meristem_rule = Rule(TreeTypes.Meristem, rhs = mer -> TreeTypes.Node() + + (TreeTypes.Bud(), TreeTypes.Leaf()) + + TreeTypes.Internode() + TreeTypes.Meristem()) + +function prob_break(bud) + ## We move to parent node in the branch where the bud was created + node = parent(bud) + ## We count the number of internodes between node and the first Meristem + ## moving down the graph + check, steps = has_descendant(node, condition = n -> data(n) isa TreeTypes.Meristem) + steps = Int(ceil(steps/2)) # Because it will count both the nodes and the internodes + ## Compute probability of bud break and determine whether it happens + if check + prob = min(1.0, steps*graph_data(bud).budbreak) + return rand() < prob + ## If there is no meristem, an error happened since the model does not allow + ## for this + else + error("No meristem found in branch") end +end +branch_rule = Rule(TreeTypes.Bud, + lhs = prob_break, + rhs = bud -> TreeTypes.BudNode() + TreeTypes.Internode() + TreeTypes.Meristem()) +#= - # Create geometry + color for the leaves - function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, data) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -data.leaf_angle) - # Generate the leaf - Ellipse!( - turtle, - length = l.length, - width = l.width, - move = false, - color = RGB(0.2, 0.6, 0.2), - ) - # Rotate turtle back to original direction - ra!(turtle, data.leaf_angle) - return nothing - end +The main difference with respect to the tree is that several of the parameters +will vary per TreeTypes. Also, the location of the tree and initial orientation of +the turtle will also vary. To achieve this we need to: - # Insertion angle for the bud nodes - function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, data) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -data.branch_angle) - end +(i) Add two additional initial nodes that move the turtle to the starting position +of each tree and rotates it. +(ii) Wrap the axiom, rules and the creation of the graph into a function that +takes the required parameters as inputs. - # Rules - meristem_rule = Rule( - TreeTypes.Meristem, - rhs = mer -> - TreeTypes.Node() + - (TreeTypes.Bud(), TreeTypes.Leaf()) + - TreeTypes.Internode() + - TreeTypes.Meristem(), - ) - - function prob_break(bud) - # We move to parent node in the branch where the bud was created - node = parent(bud) - # We count the number of internodes between node and the first Meristem - # moving down the graph - check, steps = has_descendant(node, condition = n -> data(n) isa TreeTypes.Meristem) - steps = Int(ceil(steps / 2)) # Because it will count both the nodes and the internodes - # Compute probability of bud break and determine whether it happens - if check - prob = min(1.0, steps * graph_data(bud).budbreak) - return rand() < prob - # If there is no meristem, an error happened since the model does not allow - # for this - else - error("No meristem found in branch") - end - end - branch_rule = Rule( - TreeTypes.Bud, - lhs = prob_break, - rhs = bud -> TreeTypes.BudNode() + TreeTypes.Internode() + TreeTypes.Meristem(), - ) - - #= - The main difference with respect to the tree is that several of the parameters - will vary per TreeTypes. Also, the location of the tree and initial orientation - of the turtle will also vary. To achieve this we need to: - - (i) Add two additional initial nodes that move the turtle to the starting - position of each tree and rotates it. - - (ii) Wrap the axiom, rules and the creation of the graph into a function that - takes the required parameters as inputs. - =# - function create_tree(origin, growth, budbreak, orientation) - axiom = T(origin) + RH(orientation) + TreeTypes.Internode() + TreeTypes.Meristem() - tree = Graph( - axiom = axiom, - rules = (meristem_rule, branch_rule), - data = TreeTypes.treeparams(growth = growth, budbreak = budbreak), - ) - return tree - end +=# +function create_tree(origin, growth, budbreak, orientation) + axiom = T(origin) + RH(orientation) + TreeTypes.Internode() + TreeTypes.Meristem() + tree = Graph(axiom = axiom, rules = (meristem_rule, branch_rule), + data = TreeTypes.treeparams(growth = growth, budbreak = budbreak)) + return tree +end +#= - #= - The code for elongating the internodes to simulate growth remains the same as - for the binary tree example - =# - getInternode = Query(TreeTypes.Internode) +The code for elongating the internodes to simulate growth remains the same as for +the binary tree example - function elongate!(tree, query) - for x in apply(tree, query) - x.length = x.length * (1.0 + data(tree).growth) - end - end +=# +getInternode = Query(TreeTypes.Internode) - function growth!(tree, query) - elongate!(tree, query) - rewrite!(tree) +function elongate!(tree, query) + for x in apply(tree, query) + x.length = x.length*(1.0 + data(tree).growth) end +end - function simulate(tree, query, nsteps) - new_tree = deepcopy(tree) - for i = 1:nsteps - growth!(new_tree, query) - end - return new_tree - end +function growth!(tree, query) + elongate!(tree, query) + rewrite!(tree) +end - #= - Let's simulate a forest of 10 x 10 trees with a distance between (and within) - rows of 2 meters. First we generate the original positions of the trees. For the - position we just need to pass a `Vec` object with the x, y, and z coordinates of - the location of each TreeTypes. The code below will generate a matrix with the - coordinates: - =# - origins = [Vec(i, j, 0) for i = 1:2.0:20.0, j = 1:2.0:20.0] - - #= - We may assume that the initial orientation is uniformly distributed between 0 - and 360 degrees: - =# - orientations = [rand() * 360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0] - - #= - For the `growth` and `budbreak` parameters we will assumed that they follow a - LogNormal and Beta distribution, respectively. We can generate random values - from these distributions using the `Distributions` package. For the relative - growth rate: - =# - growths = rand(LogNormal(-2, 0.3), 10, 10) - #histogram(vec(growths)) - - #= - And for the budbreak parameter: - =# - budbreaks = rand(Beta(2.0, 10), 10, 10) - #histogram(vec(budbreaks)) - - #= - Now we can create our forest by calling the `create_tree` function we defined - earlier with the correct inputs per tree: - =# - forest = vec(create_tree.(origins, growths, budbreaks, orientations)) - - #= - By vectorizing `create_tree()` over the different arrays, we end up with an - array of trees. Each tree is a different Graph, with its own nodes, rewriting - rules and variables. This avoids having to create a large graphs to include all - the plants in a simulation. Below we will run a simulation, first using a - sequential approach (i.e. using one core) and then using multiple cores in our - computers (please check https://docs.julialang.org/en/v1/manual/multi-threading/ - if the different cores are not being used as you may need to change some - settings in your computer). - - ## Sequential simulation - - We can simulate the growth of each tree by applying the method `simulate` to - each tree, creating a new version of the forest (the code below is an array - comprehension) - =# - newforest = [simulate(tree, getInternode, 2) for tree in forest] - - #= - And we can render the forest with the function `render` as in the binary tree - example but passing the whole forest at once - =# - render(Scene(newforest)) - - #= - If we iterate 4 more iterations we will start seeing the different individuals - diverging in size due to the differences in growth rates - =# - newforest = [simulate(tree, getInternode, 15) for tree in newforest] - render(Scene(newforest)) - - #= - ## Multithreaded simulation - - In the previous section, the simulation of growth was done sequentially, one - tree after another (since the growth of a tree only depends on its own - parameters). However, this can also be executed in multiple threads. In this - case we use an explicit loop and execute the iterations of the loop in multiple - threads using the macro `@threads`. Note that the rendering function can also be - ran in parallel (i.e. the geometry will be generated separately for each plant - and the merge together): - =# - using Base.Threads - newforest = deepcopy(forest) - @threads for i = 1:length(forest) - newforest[i] = simulate(forest[i], getInternode, 6) - end - render(Scene(newforest), parallel = true) - - #= - An alternative way to perform the simulation is to have an outer loop for each - timestep and an internal loop over the different trees. Although this approach - is not required for this simple model, most FSP models will probably need such a - scheme as growth of each individual plant will depend on competition for - resources with neighbouring plants. In this case, this approach would look as - follows: - =# - newforest = deepcopy(forest) - for step = 1:15 - @threads for i = 1:length(newforest) - newforest[i] = simulate(newforest[i], getInternode, 1) - end +function simulate(tree, query, nsteps) + new_tree = deepcopy(tree) + for i in 1:nsteps + growth!(new_tree, query) end - render(Scene(newforest), parallel = true) - - #= - # Customizing the scene - - Here we are going to customize the scene of our simulation by adding a - horizontal tile represting soil and tweaking the 3D representation. When we want - to combine plants generated from graphs with any other geometric element it is - best to combine all these geometries in a `GLScene` object. We can start the - scene with the `newforest` generated in the above: - =# - scene = Scene(newforest) - - #= - We can create the soil tile directly, without having to create a graph. The - simplest approach is two use a special constructor `Rectangle` where one species - a corner of the rectangle and two vectors defining the two sides of the vectors. - Both the sides and the corner need to be specified with `Vec` just like in the - above when we determined the origin of each plant. VirtualPlantLab offers some shortcuts: - `O()` returns the origin (`Vec(0.0, 0.0, 0.0)`), whereas `X`, `Y` and `Z` - returns the corresponding axes and you can scale them by passing the desired - length as input. Below, a rectangle is created on the XY plane with the origin - as a corner and each side being 11 units long: - =# - soil = Rectangle(length = 21.0, width = 21.0) - rotatey!(soil, pi / 2) - VirtualPlantLab.translate!(soil, Vec(0.0, 10.5, 0.0)) - - #= - We can now add the `soil` to the `scene` object with the `add!` function. - =# - VirtualPlantLab.add!(scene, mesh = soil, color = RGB(1, 1, 0)) - - #= - We can now render the scene that combines the random forest of binary trees and - a yellow soil. Notice that in all previous figures, a coordinate system with - grids was being depicted. This is helpful for debugging your code but also to - help setup the scene (e.g. if you are not sure how big the soil tile should be). - Howver, it may be distracting for the visualization. It turns out that we can - turn that off with `show_axes = false`: - =# - render(scene, axes = false) - - #= - We may also want to save a screenshot of the scene. For this, we need to store - the output of the `render` function. We can then resize the window rendering the - scene, move around, zoom, etc. When we have a perspective that we like, we can - run the `save_scene` function on the object returned from `render`. The argument - `resolution` can be adjusted in both `render` to increase the number of pixels - in the final image. A helper function `calculate_resolution` is provided to - compute the resolution from a physical width and height in cm and a dpi (e.g., - useful for publications and posters): - =# - res = calculate_resolution(width = 16.0, height = 16.0, dpi = 1_000) - output = render(scene, axes = false, resolution = res) - #export_scene(scene = output, filename = "nice_trees.png") + return new_tree +end +#= + +Let's simulate a forest of 10 x 10 trees with a distance between (and within) rows +of 2 meters. First we generate the original positions of the trees. For the +position we just need to pass a `Vec` object with the x, y, and z coordinates of +the location of each TreeTypes. The code below will generate a matrix with the coordinates: + +=# +origins = [Vec(i,j,0) for i = 1:2.0:20.0, j = 1:2.0:20.0] +#= + +We may assume that the initial orientation is uniformly distributed between 0 and 360 degrees: + +=# +orientations = [rand()*360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0] +#= + +For the `growth` and `budbreak` parameters we will assumed that they follow a +LogNormal and Beta distribution, respectively. We can generate random +values from these distributions using the `Distributions` package. For the +relative growth rate: + +=# +growths = rand(LogNormal(-2, 0.3), 10, 10) +histogram(vec(growths)) +#= + +And for the budbreak parameter: + +=# +budbreaks = rand(Beta(2.0, 10), 10, 10) +histogram(vec(budbreaks)) +#= + +Now we can create our forest by calling the `create_tree` function we defined earlier +with the correct inputs per tree: + +=# +forest = vec(create_tree.(origins, growths, budbreaks, orientations)); +#= + +By vectorizing `create_tree()` over the different arrays, we end up with an array +of trees. Each tree is a different Graph, with its own nodes, rewriting rules +and variables. This avoids having to create a large graphs to include all the +plants in a simulation. Below we will run a simulation, first using a sequential +approach (i.e. using one core) and then using multiple cores in our computers (please check +https://docs.julialang.org/en/v1/manual/multi-threading/ if the different cores are not being used +as you may need to change some settings in your computer). + +## Sequential simulation + +We can simulate the growth of each tree by applying the method `simulate` to each +tree, creating a new version of the forest (the code below is an array comprehension) + +=# +newforest = [simulate(tree, getInternode, 2) for tree in forest]; +#= + +And we can render the forest with the function `render` as in the binary tree +example but passing the whole forest at once + +=# +render(Scene(newforest)) +#= + +If we iterate 4 more iterations we will start seeing the different individuals +diverging in size due to the differences in growth rates + +=# +newforest = [simulate(tree, getInternode, 15) for tree in newforest]; +render(Scene(newforest)) +#= + +## Multithreaded simulation + +In the previous section, the simulation of growth was done sequentially, one tree +after another (since the growth of a tree only depends on its own parameters). However, +this can also be executed in multiple threads. In this case we use an explicit loop +and execute the iterations of the loop in multiple threads using the macro `@threads`. +Note that the rendering function can also be ran in parallel (i.e. the geometry will be +generated separately for each plant and the merge together): + +=# +using Base.Threads +newforest = deepcopy(forest) +@threads for i in eachindex(forest) + newforest[i] = simulate(forest[i], getInternode, 6) +end +render(Scene(newforest), parallel = true) +#= + +An alternative way to perform the simulation is to have an outer loop for each timestep and an internal loop over the different trees. Although this approach is not required for this simple model, most FSP models will probably need such a scheme as growth of each individual plant will depend on competition for resources with neighbouring plants. In this case, this approach would look as follows: +=# +newforest = deepcopy(forest) +for step in 1:15 + @threads for i in eachindex(newforest) + newforest[i] = simulate(newforest[i], getInternode, 1) + end end +render(Scene(newforest), parallel = true) +#= + +# Customizing the scene + +Here we are going to customize the scene of our simulation by adding a horizontal tile represting soil and +tweaking the 3D representation. When we want to combine plants generated from graphs with any other +geometric element it is best to combine all these geometries in a `GLScene` object. We can start the scene +with the `newforest` generated in the above: + +=# +scene = Scene(newforest); +#= + +We can create the soil tile directly, without having to create a graph. The simplest approach is two use +a special constructor `Rectangle` where one species a corner of the rectangle and two vectors defining the +two sides of the vectors. Both the sides and the corner need to be specified with `Vec` just like in the +above when we determined the origin of each plant. VPL offers some shortcuts: `O()` returns the origin +(`Vec(0.0, 0.0, 0.0)`), whereas `X`, `Y` and `Z` returns the corresponding axes and you can scale them by +passing the desired length as input. Below, a rectangle is created on the XY plane with the origin as a +corner and each side being 11 units long: + +=# +soil = Rectangle(length = 21.0, width = 21.0) +rotatey!(soil, pi/2) +VirtualPlantLab.translate!(soil, Vec(0.0, 10.5, 0.0)) +#= + +We can now add the `soil` to the `scene` object with the `add!` function. + +=# +VirtualPlantLab.add!(scene, mesh = soil, color = RGB(1,1,0)) +#= + +We can now render the scene that combines the random forest of binary trees and a yellow soil. Notice that +in all previous figures, a coordinate system with grids was being depicted. This is helpful for debugging +your code but also to help setup the scene (e.g. if you are not sure how big the soil tile should be). +Howver, it may be distracting for the visualization. It turns out that we can turn that off with +`show_axes = false`: + +=# +render(scene, axes = false) +#= + +We may also want to save a screenshot of the scene. For this, we need to store the output of the `render` function. +We can then resize the window rendering the scene, move around, zoom, etc. When we have a perspective that we like, +we can run the `save_scene` function on the object returned from `render`. The argument `resolution` can be adjusted in both +`render` to increase the number of pixels in the final image. A helper function `calculate_resolution` is provided to +compute the resolution from a physical width and height in cm and a dpi (e.g., useful for publications and posters): + +=# +res = calculate_resolution(width = 16.0, height = 16.0, dpi = 1_000) +output = render(scene, axes = false, resolution = res) +export_scene(scene = output, filename = "nice_trees.png") diff --git a/test/growthforest.jl b/test/growthforest.jl index 3944173..7adffb0 100644 --- a/test/growthforest.jl +++ b/test/growthforest.jl @@ -1,6 +1,8 @@ #= - Growth forest -Alejandro Morales Sierra +# Growth forest + +Alejandro Morales + Centre for Crop Systems Analysis - Wageningen University @@ -13,16 +15,17 @@ biomass of each organ is then translated in to an area or volume and the dimensions of the organs are updated accordingly (assuming a particular shape). The following packages are needed: + =# using VirtualPlantLab, ColorTypes using Base.Threads: @threads +using Plots import Random using FastGaussQuadrature using Distributions -import GLMakie Random.seed!(123456789) - #= + ## Model definition ### Node types @@ -34,439 +37,422 @@ module. The differences with respect to the previous example are: - A relative sink strength approach is used to allocate biomass to organs - The geometry of the organs is updated based on the new biomass - Bud break probability is a function of distance to apical meristem + =# -# Data types +## Data types module TreeTypes -using VirtualPlantLab -using Distributions -# Meristem -Base.@kwdef mutable struct Meristem <: VirtualPlantLab.Node - age::Int64 = 0 # Age of the meristem -end -# Bud -struct Bud <: VirtualPlantLab.Node end -# Node -struct Node <: VirtualPlantLab.Node end -# BudNode -struct BudNode <: VirtualPlantLab.Node end -# Internode (needs to be mutable to allow for changes over time) -Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node - age::Int64 = 0 # Age of the internode - biomass::Float64 = 0.0 # Initial biomass - length::Float64 = 0.0 # Internodes - width::Float64 = 0.0 # Internodes - sink::Exponential{Float64} = Exponential(5) + using VirtualPlantLab + using Distributions + ## Meristem + Base.@kwdef mutable struct Meristem <: VirtualPlantLab.Node + age::Int64 = 0 ## Age of the meristem + end + ## Bud + struct Bud <: VirtualPlantLab.Node end + ## Node + struct Node <: VirtualPlantLab.Node end + ## BudNode + struct BudNode <: VirtualPlantLab.Node end + ## Internode (needs to be mutable to allow for changes over time) + Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node + age::Int64 = 0 ## Age of the internode + biomass::Float64 = 0.0 ## Initial biomass + length::Float64 = 0.0 ## Internodes + width::Float64 = 0.0 ## Internodes + sink::Exponential{Float64} = Exponential(5) + end + ## Leaf + Base.@kwdef mutable struct Leaf <: VirtualPlantLab.Node + age::Int64 = 0 ## Age of the leaf + biomass::Float64 = 0.0 ## Initial biomass + length::Float64 = 0.0 ## Leaves + width::Float64 = 0.0 ## Leaves + sink::Beta{Float64} = Beta(2,5) + end + ## Graph-level variables -> mutable because we need to modify them during growth + Base.@kwdef mutable struct treeparams + ## Variables + biomass::Float64 = 2e-3 ## Current total biomass (g) + ## Parameters + RGR::Float64 = 1.0 ## Relative growth rate (1/d) + IB0::Float64 = 1e-3 ## Initial biomass of an internode (g) + SIW::Float64 = 0.1e6 ## Specific internode weight (g/m3) + IS::Float64 = 15.0 ## Internode shape parameter (length/width) + LB0::Float64 = 1e-3 ## Initial biomass of a leaf + SLW::Float64 = 100.0 ## Specific leaf weight (g/m2) + LS::Float64 = 3.0 ## Leaf shape parameter (length/width) + budbreak::Float64 = 1/0.5 ## Bud break probability coefficient (in 1/m) + plastochron::Int64 = 5 ## Number of days between phytomer production + leaf_expansion::Float64 = 15.0 ## Number of days that a leaf expands + phyllotaxis::Float64 = 140.0 + leaf_angle::Float64 = 30.0 + branch_angle::Float64 = 45.0 + end end -# Leaf -Base.@kwdef mutable struct Leaf <: VirtualPlantLab.Node - age::Int64 = 0 # Age of the leaf - biomass::Float64 = 0.0 # Initial biomass - length::Float64 = 0.0 # Leaves - width::Float64 = 0.0 # Leaves - sink::Beta{Float64} = Beta(2, 5) + +import .TreeTypes +#= + +### Geometry + +The methods for creating the geometry and color of the tree are the same as in +the previous example. + +=# +## Create geometry + color for the internodes +function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, vars) + ## Rotate turtle around the head to implement elliptical phyllotaxis + rh!(turtle, vars.phyllotaxis) + HollowCylinder!(turtle, length = i.length, height = i.width, width = i.width, + move = true, color = RGB(0.5,0.4,0.0)) + return nothing end -# Graph-level variables -> mutable because we need to modify them during growth -Base.@kwdef mutable struct treeparams - # Variables - biomass::Float64 = 2e-3 # Current total biomass (g) - # Parameters - RGR::Float64 = 1.0 # Relative growth rate (1/d) - IB0::Float64 = 1e-3 # Initial biomass of an internode (g) - SIW::Float64 = 0.1e6 # Specific internode weight (g/m3) - IS::Float64 = 15.0 # Internode shape parameter (length/width) - LB0::Float64 = 1e-3 # Initial biomass of a leaf - SLW::Float64 = 100.0 # Specific leaf weight (g/m2) - LS::Float64 = 3.0 # Leaf shape parameter (length/width) - budbreak::Float64 = 1 / 0.5 # Bud break probability coefficient (in 1/m) - plastochron::Int64 = 5 # Number of days between phytomer production - leaf_expansion::Float64 = 15.0 # Number of days that a leaf expands - phyllotaxis::Float64 = 140.0 - leaf_angle::Float64 = 30.0 - branch_angle::Float64 = 45.0 + +## Create geometry + color for the leaves +function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -vars.leaf_angle) + ## Generate the leaf + Ellipse!(turtle, length = l.length, width = l.width, move = false, + color = RGB(0.2,0.6,0.2)) + ## Rotate turtle back to original direction + ra!(turtle, vars.leaf_angle) + return nothing end + +## Insertion angle for the bud nodes +function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -vars.branch_angle) end +#= -import .TreeTypes +### Development -let - #= - ### Geometry - - The methods for creating the geometry and color of the tree are the same as in - the previous example. - =# - # Create geometry + color for the internodes - function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, vars) - # Rotate turtle around the head to implement elliptical phyllotaxis - rh!(turtle, vars.phyllotaxis) - HollowCylinder!( - turtle, - length = i.length, - height = i.width, - width = i.width, - move = true, - color = RGB(0.5, 0.4, 0.0), - ) - return nothing - end +The meristem rule is now parameterized by the initial states of the leaves and +internodes and will only be triggered every X days where X is the plastochron. - # Create geometry + color for the leaves - function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -vars.leaf_angle) - # Generate the leaf - Ellipse!( - turtle, - length = l.length, - width = l.width, - move = false, - color = RGB(0.2, 0.6, 0.2), - ) - # Rotate turtle back to original direction - ra!(turtle, vars.leaf_angle) - return nothing - end +=# +## Create right side of the growth rule (parameterized by the initial states +## of the leaves and internodes) +function create_meristem_rule(vleaf, vint) + meristem_rule = Rule(TreeTypes.Meristem, + lhs = mer -> mod(data(mer).age, graph_data(mer).plastochron) == 0, + rhs = mer -> TreeTypes.Node() + + (TreeTypes.Bud(), + TreeTypes.Leaf(biomass = vleaf.biomass, + length = vleaf.length, + width = vleaf.width)) + + TreeTypes.Internode(biomass = vint.biomass, + length = vint.length, + width = vint.width) + + TreeTypes.Meristem()) +end +#= - # Insertion angle for the bud nodes - function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -vars.branch_angle) - end - #= - - ### Development - - The meristem rule is now parameterized by the initial states of the leaves and - internodes and will only be triggered every X days where X is the plastochron. - =# - # Create right side of the growth rule (parameterized by the initial states - # of the leaves and internodes) - function create_meristem_rule(vleaf, vint) - meristem_rule = Rule( - TreeTypes.Meristem, - lhs = mer -> mod(data(mer).age, graph_data(mer).plastochron) == 0, - rhs = mer -> - TreeTypes.Node() + - ( - TreeTypes.Bud(), - TreeTypes.Leaf( - biomass = vleaf.biomass, - length = vleaf.length, - width = vleaf.width, - ), - ) + - TreeTypes.Internode( - biomass = vint.biomass, - length = vint.length, - width = vint.width, - ) + - TreeTypes.Meristem(), - ) - end +The bud break probability is now a function of distance to the apical meristem +rather than the number of internodes. An adhoc traversal is used to compute this +length of the main branch a bud belongs to (ignoring the lateral branches). - #= - The bud break probability is now a function of distance to the apical meristem - rather than the number of internodes. An adhoc traversal is used to compute this - length of the main branch a bud belongs to (ignoring the lateral branches). - =# - # Compute the probability that a bud breaks as function of distance to the meristem - function prob_break(bud) - # We move to parent node in the branch where the bud was created - node = parent(bud) - # Extract the first internode - child = filter(x -> data(x) isa TreeTypes.Internode, children(node))[1] - data_child = data(child) - # We measure the length of the branch until we find the meristem - distance = 0.0 - while !isa(data_child, TreeTypes.Meristem) - # If we encounter an internode, store the length and move to the next node - if data_child isa TreeTypes.Internode - distance += data_child.length - child = children(child)[1] - data_child = data(child) - # If we encounter a node, extract the next internode - elseif data_child isa TreeTypes.Node +=# +## Compute the probability that a bud breaks as function of distance to the meristem +function prob_break(bud) + ## We move to parent node in the branch where the bud was created + node = parent(bud) + ## Extract the first internode + child = filter(x -> data(x) isa TreeTypes.Internode, children(node))[1] + data_child = data(child) + ## We measure the length of the branch until we find the meristem + distance = 0.0 + while !isa(data_child, TreeTypes.Meristem) + ## If we encounter an internode, store the length and move to the next node + if data_child isa TreeTypes.Internode + distance += data_child.length + child = children(child)[1] + data_child = data(child) + ## If we encounter a node, extract the next internode + elseif data_child isa TreeTypes.Node child = filter(x -> data(x) isa TreeTypes.Internode, children(child))[1] data_child = data(child) - else - error("Should be Internode, Node or Meristem") - end + else + error("Should be Internode, Node or Meristem") end - # Compute the probability of bud break as function of distance and - # make stochastic decision - prob = min(1.0, distance * graph_data(bud).budbreak) - return rand() < prob end + ## Compute the probability of bud break as function of distance and + ## make stochastic decision + prob = min(1.0, distance*graph_data(bud).budbreak) + return rand() < prob +end - # Branch rule parameterized by initial states of internodes - function create_branch_rule(vint) - branch_rule = Rule( - TreeTypes.Bud, +## Branch rule parameterized by initial states of internodes +function create_branch_rule(vint) + branch_rule = Rule(TreeTypes.Bud, lhs = prob_break, - rhs = bud -> - TreeTypes.BudNode() + - TreeTypes.Internode( - biomass = vint.biomass, - length = vint.length, - width = vint.width, - ) + - TreeTypes.Meristem(), - ) - end + rhs = bud -> TreeTypes.BudNode() + + TreeTypes.Internode(biomass = vint.biomass, + length = vint.length, + width = vint.width) + + TreeTypes.Meristem()) +end +#= - #= - ### Growth - - We need some functions to compute the length and width of a leaf or internode - from its biomass - =# - function leaf_dims(biomass, vars) - leaf_biomass = biomass - leaf_area = biomass / vars.SLW - leaf_length = sqrt(leaf_area * 4 * vars.LS / pi) - leaf_width = leaf_length / vars.LS - return leaf_length, leaf_width - end - function int_dims(biomass, vars) - int_biomass = biomass - int_volume = biomass / vars.SIW - int_length = cbrt(int_volume * 4 * vars.IS^2 / pi) - int_width = int_length / vars.IS - return int_length, int_width - end +### Growth - #= - Each day, the total biomass of the tree is updated using a simple RGR formula - and the increment of biomass is distributed across the organs proportionally to - their relative sink strength (of leaves or internodes). - - The sink strength of leaves is modelled with a beta distribution scaled to the - `leaf_expansion` argument that determines the duration of leaf growth, whereas - for the internodes it follows a negative exponential distribution. The `pdf` - function computes the probability density of each distribution which is taken as - proportional to the sink strength (the model is actually source-limited since we - imposed a particular growth rate). - =# - sink_strength(leaf, vars) = - leaf.age > vars.leaf_expansion ? 0.0 : - pdf(leaf.sink, leaf.age / vars.leaf_expansion) / 100.0 - # plot( - # 0:1:50, - # x -> sink_strength(TreeTypes.Leaf(age = x), TreeTypes.treeparams()), - # xlabel = "Age", - # ylabel = "Sink strength", - # label = "Leaf", - # ) - - sink_strength(int) = pdf(int.sink, int.age) - #plot!(0:1:50, x -> sink_strength(TreeTypes.Internode(age = x)), label = "Internode") - - #= - Now we need a function that updates the biomass of the tree, allocates it to the - different organs and updates the dimensions of said organs. For simplicity, - we create the functions `leaves()` and `internodes()` that will apply the queries - to the tree required to extract said nodes: - =# - get_leaves(tree) = apply(tree, Query(TreeTypes.Leaf)) - get_internodes(tree) = apply(tree, Query(TreeTypes.Internode)) - - #= - The age of the different organs is updated every time step: - =# - function age!(all_leaves, all_internodes, all_meristems) - for leaf in all_leaves - leaf.age += 1 - end - for int in all_internodes - int.age += 1 - end - for mer in all_meristems - mer.age += 1 - end - return nothing - end +We need some functions to compute the length and width of a leaf or internode +from its biomass - #= - The daily growth is allocated to different organs proportional to their sink - strength. - =# - function grow!(tree, all_leaves, all_internodes) - # Compute total biomass increment - tvars = data(tree) - ΔB = tvars.RGR * tvars.biomass - tvars.biomass += ΔB - # Total sink strength - total_sink = 0.0 - for leaf in all_leaves - total_sink += sink_strength(leaf, tvars) - end - for int in all_internodes - total_sink += sink_strength(int) - end - # Allocate biomass to leaves and internodes - for leaf in all_leaves - leaf.biomass += ΔB * sink_strength(leaf, tvars) / total_sink - end - for int in all_internodes - int.biomass += ΔB * sink_strength(int) / total_sink - end - return nothing - end +=# +function leaf_dims(biomass, vars) + leaf_biomass = biomass + leaf_area = biomass/vars.SLW + leaf_length = sqrt(leaf_area*4*vars.LS/pi) + leaf_width = leaf_length/vars.LS + return leaf_length, leaf_width +end - #= - Finally, we need to update the dimensions of the organs. The leaf dimensions are - =# - function size_leaves!(all_leaves, tvars) - for leaf in all_leaves - leaf.length, leaf.width = leaf_dims(leaf.biomass, tvars) - end - return nothing +function int_dims(biomass, vars) + int_biomass = biomass + int_volume = biomass/vars.SIW + int_length = cbrt(int_volume*4*vars.IS^2/pi) + int_width = int_length/vars.IS + return int_length, int_width +end +#= + +Each day, the total biomass of the tree is updated using a simple RGR formula +and the increment of biomass is distributed across the organs proportionally to +their relative sink strength (of leaves or internodes). + +The sink strength of leaves is modelled with a beta distribution scaled to the +`leaf_expansion` argument that determines the duration of leaf growth, whereas +for the internodes it follows a negative exponential distribution. The `pdf` +function computes the probability density of each distribution which is taken as +proportional to the sink strength (the model is actually source-limited since we +imposed a particular growth rate). + +=# +sink_strength(leaf, vars) = leaf.age > vars.leaf_expansion ? 0.0 : + pdf(leaf.sink, leaf.age/vars.leaf_expansion)/100.0 +plot(0:1:50, x -> sink_strength(TreeTypes.Leaf(age = x), TreeTypes.treeparams()), + xlabel = "Age", ylabel = "Sink strength", label = "Leaf") +#= + +=# +sink_strength(int) = pdf(int.sink, int.age) +plot!(0:1:50, x -> sink_strength(TreeTypes.Internode(age = x)), label = "Internode") +#= + +Now we need a function that updates the biomass of the tree, allocates it to the +different organs and updates the dimensions of said organs. For simplicity, +we create the functions `leaves()` and `internodes()` that will apply the queries +to the tree required to extract said nodes: + +=# +get_leaves(tree) = apply(tree, Query(TreeTypes.Leaf)) +get_internodes(tree) = apply(tree, Query(TreeTypes.Internode)) +#= + +The age of the different organs is updated every time step: + +=# +function age!(all_leaves, all_internodes, all_meristems) + for leaf in all_leaves + leaf.age += 1 end - function size_internodes!(all_internodes, tvars) - for int in all_internodes - int.length, int.width = int_dims(int.biomass, tvars) - end - return nothing + for int in all_internodes + int.age += 1 end - - #= - ### Daily step - - All the growth and developmental functions are combined together into a daily - step function that updates the forest by iterating over the different trees in - parallel. - =# - get_meristems(tree) = apply(tree, Query(TreeTypes.Meristem)) - function daily_step!(forest) - @threads for tree in forest - # Retrieve all the relevant organs - all_leaves = get_leaves(tree) - all_internodes = get_internodes(tree) - all_meristems = get_meristems(tree) - # Update the age of the organs - age!(all_leaves, all_internodes, all_meristems) - # Grow the tree - grow!(tree, all_leaves, all_internodes) - tvars = data(tree) - size_leaves!(all_leaves, tvars) - size_internodes!(all_internodes, tvars) - # Developmental rules - rewrite!(tree) - end + for mer in all_meristems + mer.age += 1 end + return nothing +end +#= - #= - ### Initialization - - The trees are initialized in a regular grid with random values for the initial - orientation and RGR: - =# - RGRs = rand(Normal(0.3, 0.01), 10, 10) - #histogram(vec(RGRs)) - - orientations = [rand() * 360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0] - #histogram(vec(orientations)) - - origins = [Vec(i, j, 0) for i = 1:2.0:20.0, j = 1:2.0:20.0] - - #= - The following initalizes a tree based on the origin, orientation and RGR: - =# - function create_tree(origin, orientation, RGR) - # Initial state and parameters of the tree - vars = TreeTypes.treeparams(RGR = RGR) - # Initial states of the leaves - leaf_length, leaf_width = leaf_dims(vars.LB0, vars) - vleaf = (biomass = vars.LB0, length = leaf_length, width = leaf_width) - # Initial states of the internodes - int_length, int_width = int_dims(vars.LB0, vars) - vint = (biomass = vars.IB0, length = int_length, width = int_width) - # Growth rules - meristem_rule = create_meristem_rule(vleaf, vint) - branch_rule = create_branch_rule(vint) - axiom = - T(origin) + - RH(orientation) + - TreeTypes.Internode( - biomass = vint.biomass, - length = vint.length, - width = vint.width, - ) + - TreeTypes.Meristem() - tree = Graph(axiom = axiom, rules = (meristem_rule, branch_rule), data = vars) - return tree - end +The daily growth is allocated to different organs proportional to their sink +strength. - #= - ## Visualization - - As in the previous example, it makes sense to visualize the forest with a soil - tile beneath it. Unlike in the previous example, we will construct the soil tile - using a dedicated graph and generate a `Scene` object which can later be - merged with the rest of scene generated in daily step: - =# - Base.@kwdef struct Soil <: VirtualPlantLab.Node - length::Float64 - width::Float64 +=# +function grow!(tree, all_leaves, all_internodes) + ## Compute total biomass increment + tvars = data(tree) + ΔB = tvars.RGR*tvars.biomass + tvars.biomass += ΔB + ## Total sink strength + total_sink = 0.0 + for leaf in all_leaves + total_sink += sink_strength(leaf, tvars) + end + for int in all_internodes + total_sink += sink_strength(int) end - function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, vars) - Rectangle!( - turtle, - length = s.length, - width = s.width, - color = RGB(255 / 255, 236 / 255, 179 / 255), - ) + ## Allocate biomass to leaves and internodes + for leaf in all_leaves + leaf.biomass += ΔB*sink_strength(leaf, tvars)/total_sink end - soil_graph = - RA(-90.0) + - T(Vec(0.0, 10.0, 0.0)) + # Moves into position - Soil(length = 20.0, width = 20.0) # Draws the soil tile - soil = Scene(Graph(axiom = soil_graph)) - render(soil, axes = false) - - #= - And the following function renders the entire scene (notice that we need to - use `display()` to force the rendering of the scene when called within a loop - or a function): - =# - function render_forest(forest, soil) - scene = Scene(vec(forest)) # create scene from forest - scene = Scene([scene, soil]) # merges the two scenes - render(scene) + for int in all_internodes + int.biomass += ΔB*sink_strength(int)/total_sink end + return nothing +end +#= - #= - ## Retrieving canopy-level data - - We may want to extract some information at the canopy level such as LAI. This is - best achieved with a query: - =# - function get_LAI(forest) - LAI = 0.0 - @threads for tree in forest - for leaf in get_leaves(tree) - LAI += leaf.length * leaf.width * pi / 4.0 - end - end - return LAI / 400.0 +Finally, we need to update the dimensions of the organs. The leaf dimensions are + +=# +function size_leaves!(all_leaves, tvars) + for leaf in all_leaves + leaf.length, leaf.width = leaf_dims(leaf.biomass, tvars) + end + return nothing +end +function size_internodes!(all_internodes, tvars) + for int in all_internodes + int.length, int.width = int_dims(int.biomass, tvars) + end + return nothing +end +#= + +### Daily step + +All the growth and developmental functions are combined together into a daily +step function that updates the forest by iterating over the different trees in +parallel. + +=# +get_meristems(tree) = apply(tree, Query(TreeTypes.Meristem)) +function daily_step!(forest) + @threads for tree in forest + ## Retrieve all the relevant organs + all_leaves = get_leaves(tree) + all_internodes = get_internodes(tree) + all_meristems = get_meristems(tree) + ## Update the age of the organs + age!(all_leaves, all_internodes, all_meristems) + ## Grow the tree + grow!(tree, all_leaves, all_internodes) + tvars = data(tree) + size_leaves!(all_leaves, tvars) + size_internodes!(all_internodes, tvars) + ## Developmental rules + rewrite!(tree) end +end +#= + +### Initialization + +The trees are initialized in a regular grid with random values for the initial +orientation and RGR: + +=# +RGRs = rand(Normal(0.3,0.01), 10, 10) +histogram(vec(RGRs)) +#= + +=# +orientations = [rand()*360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0] +histogram(vec(orientations)) +#= - #= - ## Simulation +=# +origins = [Vec(i,j,0) for i = 1:2.0:20.0, j = 1:2.0:20.0]; +#= + +The following initalizes a tree based on the origin, orientation and RGR: + +=# +function create_tree(origin, orientation, RGR) + ## Initial state and parameters of the tree + vars = TreeTypes.treeparams(RGR = RGR) + ## Initial states of the leaves + leaf_length, leaf_width = leaf_dims(vars.LB0, vars) + vleaf = (biomass = vars.LB0, length = leaf_length, width = leaf_width) + ## Initial states of the internodes + int_length, int_width = int_dims(vars.LB0, vars) + vint = (biomass = vars.IB0, length = int_length, width = int_width) + ## Growth rules + meristem_rule = create_meristem_rule(vleaf, vint) + branch_rule = create_branch_rule(vint) + axiom = T(origin) + RH(orientation) + + TreeTypes.Internode(biomass = vint.biomass, + length = vint.length, + width = vint.width) + + TreeTypes.Meristem() + tree = Graph(axiom = axiom, rules = (meristem_rule, branch_rule), + data = vars) + return tree +end +#= + + +## Visualization + +As in the previous example, it makes sense to visualize the forest with a soil +tile beneath it. Unlike in the previous example, we will construct the soil tile +using a dedicated graph and generate a `Scene` object which can later be +merged with the rest of scene generated in daily step: + +=# +Base.@kwdef struct Soil <: VirtualPlantLab.Node + length::Float64 + width::Float64 +end +function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, vars) + Rectangle!(turtle, length = s.length, width = s.width, color = RGB(255/255, 236/255, 179/255)) +end +soil_graph = RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + # Moves into position + Soil(length = 20.0, width = 20.0) # Draws the soil tile +soil = Scene(Graph(axiom = soil_graph)); +render(soil, axes = false) +#= + +And the following function renders the entire scene (notice that we need to +use `display()` to force the rendering of the scene when called within a loop +or a function): - We can now create a forest of trees on a regular grid: - =# - forest = create_tree.(origins, orientations, RGRs) - render_forest(forest, soil) - for i = 1:50 - daily_step!(forest) +=# +function render_forest(forest, soil) + scene = Scene(vec(forest)) # create scene from forest + scene = Scene([scene, soil]) # merges the two scenes + render(scene) +end +#= + + +## Retrieving canopy-level data + +We may want to extract some information at the canopy level such as LAI. This is +best achieved with a query: + +=# +function get_LAI(forest) + LAI = 0.0 + @threads for tree in forest + for leaf in get_leaves(tree) + LAI += leaf.length*leaf.width*pi/4.0 + end end - render_forest(forest, soil) + return LAI/400.0 +end +#= - #= - And compute the leaf area index: - =# - get_LAI(forest) +## Simulation +We can now create a forest of trees on a regular grid: + +=# +forest = create_tree.(origins, orientations, RGRs); +render_forest(forest, soil) +for i in 1:50 + daily_step!(forest) end +render_forest(forest, soil) +#= + +And compute the leaf area index: + +=# +get_LAI(forest) diff --git a/test/raytracedforest.jl b/test/raytracedforest.jl index 2c0b3da..0ad6d58 100644 --- a/test/raytracedforest.jl +++ b/test/raytracedforest.jl @@ -1,23 +1,28 @@ #= - Ray-traced forest of binary trees" -Alejandro Morales Sierra +# Ray-traced forest + +Alejandro Morales + Centre for Crop Systems Analysis - Wageningen University + In this example we extend the forest growth model to include PAR interception a radiation use efficiency to compute the daily growth rate. The following packages are needed: + =# using VirtualPlantLab, ColorTypes import GLMakie using Base.Threads: @threads +using Plots import Random using FastGaussQuadrature using Distributions using SkyDomes Random.seed!(123456789) - #= + ## Model definition ### Node types @@ -26,583 +31,558 @@ The data types needed to simulate the trees are given in the following module. The difference with respec to the previous model is that Internodes and Leaves have optical properties needed for ray tracing (they are defined as Lambertian surfaces). + =# -# Data types +## Data types module TreeTypes -using VirtualPlantLab -using Distributions -# Meristem -Base.@kwdef mutable struct Meristem <: VirtualPlantLab.Node - age::Int64 = 0 # Age of the meristem -end -# Bud -struct Bud <: VirtualPlantLab.Node end -# Node -struct Node <: VirtualPlantLab.Node end -# BudNode -struct BudNode <: VirtualPlantLab.Node end -# Internode (needs to be mutable to allow for changes over time) -Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node - age::Int64 = 0 # Age of the internode - biomass::Float64 = 0.0 # Initial biomass - length::Float64 = 0.0 # Internodes - width::Float64 = 0.0 # Internodes - sink::Exponential{Float64} = Exponential(5) - material::Lambertian{1} = Lambertian(τ = 0.1, ρ = 0.05) # Leaf material -end -# Leaf -Base.@kwdef mutable struct Leaf <: VirtualPlantLab.Node - age::Int64 = 0 # Age of the leaf - biomass::Float64 = 0.0 # Initial biomass - length::Float64 = 0.0 # Leaves - width::Float64 = 0.0 # Leaves - sink::Beta{Float64} = Beta(2, 5) - material::Lambertian{1} = Lambertian(τ = 0.1, ρ = 0.05) # Leaf material + using VirtualPlantLab + using Distributions + ## Meristem + Base.@kwdef mutable struct Meristem <: VirtualPlantLab.Node + age::Int64 = 0 ## Age of the meristem + end + ## Bud + struct Bud <: VirtualPlantLab.Node end + ## Node + struct Node <: VirtualPlantLab.Node end + ## BudNode + struct BudNode <: VirtualPlantLab.Node end + ## Internode (needs to be mutable to allow for changes over time) + Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node + age::Int64 = 0 ## Age of the internode + biomass::Float64 = 0.0 ## Initial biomass + length::Float64 = 0.0 ## Internodes + width::Float64 = 0.0 ## Internodes + sink::Exponential{Float64} = Exponential(5) + material::Lambertian{1} = Lambertian(τ = 0.1, ρ = 0.05) ## Leaf material + end + ## Leaf + Base.@kwdef mutable struct Leaf <: VirtualPlantLab.Node + age::Int64 = 0 ## Age of the leaf + biomass::Float64 = 0.0 ## Initial biomass + length::Float64 = 0.0 ## Leaves + width::Float64 = 0.0 ## Leaves + sink::Beta{Float64} = Beta(2,5) + material::Lambertian{1} = Lambertian(τ = 0.1, ρ = 0.05) ## Leaf material + end + ## Graph-level variables -> mutable because we need to modify them during growth + Base.@kwdef mutable struct treeparams + ## Variables + PAR::Float64 = 0.0 ## Total PAR absorbed by the leaves on the tree (MJ) + biomass::Float64 = 2e-3 ## Current total biomass (g) + ## Parameters + RUE::Float64 = 5.0 ## Radiation use efficiency (g/MJ) -> unrealistic to speed up sim + IB0::Float64 = 1e-3 ## Initial biomass of an internode (g) + SIW::Float64 = 0.1e6 ## Specific internode weight (g/m3) + IS::Float64 = 15.0 ## Internode shape parameter (length/width) + LB0::Float64 = 1e-3 ## Initial biomass of a leaf + SLW::Float64 = 100.0 ## Specific leaf weight (g/m2) + LS::Float64 = 3.0 ## Leaf shape parameter (length/width) + budbreak::Float64 = 1/0.5 ## Bud break probability coefficient (in 1/m) + plastochron::Int64 = 5 ## Number of days between phytomer production + leaf_expansion::Float64 = 15.0 ## Number of days that a leaf expands + phyllotaxis::Float64 = 140.0 + leaf_angle::Float64 = 30.0 + branch_angle::Float64 = 45.0 + end end -# Graph-level variables -> mutable because we need to modify them during growth -Base.@kwdef mutable struct treeparams - # Variables - PAR::Float64 = 0.0 # Total PAR absorbed by the leaves on the tree (MJ) - biomass::Float64 = 2e-3 # Current total biomass (g) - # Parameters - RUE::Float64 = 5.0 # Radiation use efficiency (g/MJ) -> unrealistic to speed up sim - IB0::Float64 = 1e-3 # Initial biomass of an internode (g) - SIW::Float64 = 0.1e6 # Specific internode weight (g/m3) - IS::Float64 = 15.0 # Internode shape parameter (length/width) - LB0::Float64 = 1e-3 # Initial biomass of a leaf - SLW::Float64 = 100.0 # Specific leaf weight (g/m2) - LS::Float64 = 3.0 # Leaf shape parameter (length/width) - budbreak::Float64 = 1 / 0.5 # Bud break probability coefficient (in 1/m) - plastochron::Int64 = 5 # Number of days between phytomer production - leaf_expansion::Float64 = 15.0 # Number of days that a leaf expands - phyllotaxis::Float64 = 140.0 - leaf_angle::Float64 = 30.0 - branch_angle::Float64 = 45.0 + +import .TreeTypes +#= + +### Geometry + +The methods for creating the geometry and color of the tree are the same as in +the previous example but include the materials for the ray tracer. + +=# +## Create geometry + color for the internodes +function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, data) + ## Rotate turtle around the head to implement elliptical phyllotaxis + rh!(turtle, data.phyllotaxis) + HollowCylinder!(turtle, length = i.length, height = i.width, width = i.width, + move = true, color = RGB(0.5,0.4,0.0), material = i.material) + return nothing end + +## Create geometry + color for the leaves +function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, data) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -data.leaf_angle) + ## Generate the leaf + Ellipse!(turtle, length = l.length, width = l.width, move = false, + color = RGB(0.2,0.6,0.2), material = l.material) + ## Rotate turtle back to original direction + ra!(turtle, data.leaf_angle) + return nothing end -import .TreeTypes +## Insertion angle for the bud nodes +function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, data) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -data.branch_angle) +end +#= -let - #= - ### Geometry - - The methods for creating the geometry and color of the tree are the same as in - the previous example but include the materials for the ray tracer. - =# - # Create geometry + color for the internodes - function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, data) - # Rotate turtle around the head to implement elliptical phyllotaxis - rh!(turtle, data.phyllotaxis) - HollowCylinder!( - turtle, - length = i.length, - height = i.width, - width = i.width, - move = true, - color = RGB(0.5, 0.4, 0.0), - material = i.material, - ) - return nothing - end +### Development - # Create geometry + color for the leaves - function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, data) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -data.leaf_angle) - # Generate the leaf - Ellipse!( - turtle, - length = l.length, - width = l.width, - move = false, - color = RGB(0.2, 0.6, 0.2), - material = l.material, - ) - # Rotate turtle back to original direction - ra!(turtle, data.leaf_angle) - return nothing - end +The meristem rule is now parameterized by the initial states of the leaves and +internodes and will only be triggered every X days where X is the plastochron. - # Insertion angle for the bud nodes - function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, data) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -data.branch_angle) - end +=# +## Create right side of the growth rule (parameterized by the initial states +## of the leaves and internodes) +function create_meristem_rule(vleaf, vint) + meristem_rule = Rule(TreeTypes.Meristem, + lhs = mer -> mod(data(mer).age, graph_data(mer).plastochron) == 0, + rhs = mer -> TreeTypes.Node() + + (TreeTypes.Bud(), + TreeTypes.Leaf(biomass = vleaf.biomass, + length = vleaf.length, + width = vleaf.width)) + + TreeTypes.Internode(biomass = vint.biomass, + length = vint.length, + width = vint.width) + + TreeTypes.Meristem()) +end +#= - #= - ### Development - - The meristem rule is now parameterized by the initial states of the leaves and - internodes and will only be triggered every X days where X is the plastochron. - =# - # Create right side of the growth rule (parameterized by the initial states - # of the leaves and internodes) - function create_meristem_rule(vleaf, vint) - meristem_rule = Rule( - TreeTypes.Meristem, - lhs = mer -> mod(data(mer).age, graph_data(mer).plastochron) == 0, - rhs = mer -> - TreeTypes.Node() + - ( - TreeTypes.Bud(), - TreeTypes.Leaf( - biomass = vleaf.biomass, - length = vleaf.length, - width = vleaf.width, - ), - ) + - TreeTypes.Internode( - biomass = vint.biomass, - length = vint.length, - width = vint.width, - ) + - TreeTypes.Meristem(), - ) - end +The bud break probability is now a function of distance to the apical meristem +rather than the number of internodes. An adhoc traversal is used to compute this +length of the main branch a bud belongs to (ignoring the lateral branches). - #= - The bud break probability is now a function of distance to the apical meristem - rather than the number of internodes. An adhoc traversal is used to compute this - length of the main branch a bud belongs to (ignoring the lateral branches). - =# - # Compute the probability that a bud breaks as function of distance to the meristem - function prob_break(bud) - # We move to parent node in the branch where the bud was created - node = parent(bud) - # Extract the first internode - child = filter(x -> data(x) isa TreeTypes.Internode, children(node))[1] - data_child = data(child) - # We measure the length of the branch until we find the meristem - distance = 0.0 - while !isa(data_child, TreeTypes.Meristem) - # If we encounter an internode, store the length and move to the next node - if data_child isa TreeTypes.Internode - distance += data_child.length - child = children(child)[1] - data_child = data(child) - # If we encounter a node, extract the next internode - elseif data_child isa TreeTypes.Node +=# +# Compute the probability that a bud breaks as function of distance to the meristem +function prob_break(bud) + ## We move to parent node in the branch where the bud was created + node = parent(bud) + ## Extract the first internode + child = filter(x -> data(x) isa TreeTypes.Internode, children(node))[1] + data_child = data(child) + ## We measure the length of the branch until we find the meristem + distance = 0.0 + while !isa(data_child, TreeTypes.Meristem) + ## If we encounter an internode, store the length and move to the next node + if data_child isa TreeTypes.Internode + distance += data_child.length + child = children(child)[1] + data_child = data(child) + ## If we encounter a node, extract the next internode + elseif data_child isa TreeTypes.Node child = filter(x -> data(x) isa TreeTypes.Internode, children(child))[1] data_child = data(child) - else - error("Should be Internode, Node or Meristem") - end + else + error("Should be Internode, Node or Meristem") end - # Compute the probability of bud break as function of distance and - # make stochastic decision - prob = min(1.0, distance * graph_data(bud).budbreak) - return rand() < prob end + ## Compute the probability of bud break as function of distance and + ## make stochastic decision + prob = min(1.0, distance*graph_data(bud).budbreak) + return rand() < prob +end - # Branch rule parameterized by initial states of internodes - function create_branch_rule(vint) - branch_rule = Rule( - TreeTypes.Bud, +## Branch rule parameterized by initial states of internodes +function create_branch_rule(vint) + branch_rule = Rule(TreeTypes.Bud, lhs = prob_break, - rhs = bud -> - TreeTypes.BudNode() + - TreeTypes.Internode( - biomass = vint.biomass, - length = vint.length, - width = vint.width, - ) + - TreeTypes.Meristem(), - ) - end + rhs = bud -> TreeTypes.BudNode() + + TreeTypes.Internode(biomass = vint.biomass, + length = vint.length, + width = vint.width) + + TreeTypes.Meristem()) +end +#= - #= - ### Light interception - - As growth is now dependent on intercepted PAR via RUE, we now need to simulate - light interception by the trees. We will use a ray-tracing approach to do so. - The first step is to create a scene with the trees and the light sources. As for - rendering, the scene can be created from the `forest` object by simply calling - `Scene(forest)` that will generate the 3D meshes and connect them to their - optical properties. - - However, we also want to add the soil surface as this will affect the light - distribution within the scene due to reflection from the soil surface. This is - similar to the customized scene that we created before for rendering, but now - for the light simulation. - =# - function create_soil() - soil = Rectangle(length = 21.0, width = 21.0) - rotatey!(soil, π / 2) # To put it in the XY plane - VirtualPlantLab.translate!(soil, Vec(0.0, 10.5, 0.0)) # Corner at (0,0,0) - return soil - end - function create_scene(forest) - # These are the trees - scene = Scene(vec(forest)) - # Add a soil surface - soil = create_soil() - soil_material = Lambertian(τ = 0.0, ρ = 0.21) - add!(scene, mesh = soil, material = soil_material) - # Return the scene - return scene - end +### Light interception - #= - Given the scene, we can create the light sources that can approximate the solar - irradiance on a given day, location and time of the day using the functions from - the Sky package (see package documentation for details). Given the latitude, - day of year and fraction of the day (`f = 0` being sunrise and `f = 1` being sunset), - the function `clear_sky()` computes the direct and diffuse solar radiation assuming - a clear sky. These values may be converted to different wavebands and units using - `waveband_conversion()`. Finally, the collection of light sources approximating - the solar irradiance distribution over the sky hemisphere is constructed with the - function `sky()` (this last step requires the 3D scene as input in order to place - the light sources adequately). - =# - function create_sky(; scene, lat = 52.0 * π / 180.0, DOY = 182) - # Fraction of the day and day length - fs = collect(0.1:0.1:0.9) - dec = declination(DOY) - DL = day_length(lat, dec) * 3600 - # Compute solar irradiance - temp = [clear_sky(lat = lat, DOY = DOY, f = f) for f in fs] # W/m2 - Ig = getindex.(temp, 1) - Idir = getindex.(temp, 2) - Idif = getindex.(temp, 3) - # Conversion factors to PAR for direct and diffuse irradiance - f_dir = waveband_conversion(Itype = :direct, waveband = :PAR, mode = :power) - f_dif = waveband_conversion(Itype = :diffuse, waveband = :PAR, mode = :power) - # Actual irradiance per waveband - Idir_PAR = f_dir .* Idir - Idif_PAR = f_dif .* Idif - # Create the dome of diffuse light - dome = sky( - scene, - Idir = 0.0, # No direct solar radiation - Idif = sum(Idir_PAR) / 10 * DL, # Daily Diffuse solar radiation - nrays_dif = 1_000_000, # Total number of rays for diffuse solar radiation - sky_model = StandardSky, # Angular distribution of solar radiation - dome_method = equal_solid_angles, # Discretization of the sky dome - ntheta = 9, # Number of discretization steps in the zenith angle - nphi = 12, - ) # Number of discretization steps in the azimuth angle - # Add direct sources for different times of the day - for I in Idir_PAR - push!(dome, sky(scene, Idir = I / 10 * DL, nrays_dir = 100_000, Idif = 0.0)[1]) - end - return dome - end +As growth is now dependent on intercepted PAR via RUE, we now need to simulate +light interception by the trees. We will use a ray-tracing approach to do so. +The first step is to create a scene with the trees and the light sources. As for +rendering, the scene can be created from the `forest` object by simply calling +`Scene(forest)` that will generate the 3D meshes and connect them to their +optical properties. - #= - The 3D scene and the light sources are then combined into a `RayTracer` object, - together with general settings for the ray tracing simulation chosen via `RTSettings()`. - The most important settings refer to the Russian roulette system and the grid - cloner (see section on Ray Tracing for details). The settings for the Russian - roulette system include the number of times a ray will be traced - deterministically (`maxiter`) and the probability that a ray that exceeds `maxiter` - is terminated (`pkill`). The grid cloner is used to approximate an infinite canopy - by replicating the scene in the different directions (`nx` and `ny` being the - number of replicates in each direction along the x and y axes, respectively). It - is also possible to turn on parallelization of the ray tracing simulation by - setting `parallel = true` (currently this uses Julia's builtin multithreading - capabilities). - - In addition `RTSettings()`, an acceleration structure and a splitting rule can - be defined when creating the `RayTracer` object (see ray tracing documentation - for details). The acceleration structure allows speeding up the ray tracing - by avoiding testing all rays against all objects in the scene. - =# - function create_raytracer(scene, sources) - settings = RTSettings(pkill = 0.9, maxiter = 4, nx = 5, ny = 5, parallel = true) - RayTracer( - scene, - sources, - settings = settings, - acceleration = BVH, - rule = SAH{3}(5, 10), - ) - end +However, we also want to add the soil surface as this will affect the light +distribution within the scene due to reflection from the soil surface. This is +similar to the customized scene that we created before for rendering, but now +for the light simulation. - #= - The actual ray tracing simulation is performed by calling the `trace!()` method - on the ray tracing object. This will trace all rays from all light sources and - update the radiant power absorbed by the different surfaces in the scene inside - the `Material` objects (see `feed!()` above): - =# - function run_raytracer!(forest; DOY = 182) - scene = create_scene(forest) - sources = create_sky(scene = scene, DOY = DOY) - rtobj = create_raytracer(scene, sources) - trace!(rtobj) - return nothing - end +=# +function create_soil() + soil = Rectangle(length = 21.0, width = 21.0) + rotatey!(soil, π/2) ## To put it in the XY plane + VirtualPlantLab.translate!(soil, Vec(0.0, 10.5, 0.0)) ## Corner at (0,0,0) + return soil +end +function create_scene(forest) + ## These are the trees + scene = Scene(vec(forest)) + ## Add a soil surface + soil = create_soil() + soil_material = Lambertian(τ = 0.0, ρ = 0.21) + add!(scene, mesh = soil, material = soil_material) + ## Return the scene + return scene +end +#= - #= - The total PAR absorbed for each tree is calculated from the material objects of - the different internodes (using `power()` on the `Material` object). Note that - the `power()` function returns three different values, one for each waveband, - but they are added together as RUE is defined for total PAR. - =# - # Run the ray tracer, calculate PAR absorbed per tree and add it to the daily - # total using general weighted quadrature formula - function calculate_PAR!(forest; DOY = 182) - # Reset PAR absorbed by the tree (at the start of a new day) - reset_PAR!(forest) - # Run the ray tracer to compute daily PAR absorption - run_raytracer!(forest, DOY = DOY) - # Add up PAR absorbed by each leaf within each tree - @threads for tree in forest - for l in get_leaves(tree) - data(tree).PAR += power(l.material)[1] - end - end - return nothing - end +Given the scene, we can create the light sources that can approximate the solar +irradiance on a given day, location and time of the day using the functions from +the package (see package documentation for details). Given the latitude, +day of year and fraction of the day (`f = 0` being sunrise and `f = 1` being sunset), +the function `clear_sky()` computes the direct and diffuse solar radiation assuming +a clear sky. These values may be converted to different wavebands and units using +`waveband_conversion()`. Finally, the collection of light sources approximating +the solar irradiance distribution over the sky hemisphere is constructed with the +function `sky()` (this last step requires the 3D scene as input in order to place +the light sources adequately). - # Reset PAR absorbed by the tree (at the start of a new day) - function reset_PAR!(forest) - for tree in forest - data(tree).PAR = 0.0 - end - return nothing +=# +function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182) + ## Fraction of the day and day length + fs = collect(0.1:0.1:0.9) + dec = declination(DOY) + DL = day_length(lat, dec)*3600 + ## Compute solar irradiance + temp = [clear_sky(lat = lat, DOY = DOY, f = f) for f in fs] # W/m2 + Ig = getindex.(temp, 1) + Idir = getindex.(temp, 2) + Idif = getindex.(temp, 3) + ## Conversion factors to PAR for direct and diffuse irradiance + f_dir = waveband_conversion(Itype = :direct, waveband = :PAR, mode = :power) + f_dif = waveband_conversion(Itype = :diffuse, waveband = :PAR, mode = :power) + ## Actual irradiance per waveband + Idir_PAR = f_dir.*Idir + Idif_PAR = f_dif.*Idif + ## Create the dome of diffuse light + dome = sky(scene, + Idir = 0.0, ## No direct solar radiation + Idif = sum(Idir_PAR)/10*DL, ## Daily Diffuse solar radiation + nrays_dif = 1_000_000, ## Total number of rays for diffuse solar radiation + sky_model = StandardSky, ## Angular distribution of solar radiation + dome_method = equal_solid_angles, ## Discretization of the sky dome + ntheta = 9, ## Number of discretization steps in the zenith angle + nphi = 12) ## Number of discretization steps in the azimuth angle + ## Add direct sources for different times of the day + for I in Idir_PAR + push!(dome, sky(scene, Idir = I/10*DL, nrays_dir = 100_000, Idif = 0.0)[1]) end + return dome +end +#= - #= - ### Growth - - We need some functions to compute the length and width of a leaf or internode - from its biomass - =# - function leaf_dims(biomass, data) - leaf_biomass = biomass - leaf_area = biomass / data.SLW - leaf_length = sqrt(leaf_area * 4 * data.LS / pi) - leaf_width = leaf_length / data.LS - return leaf_length, leaf_width - end +The 3D scene and the light sources are then combined into a `RayTracer` object, +together with general settings for the ray tracing simulation chosen via `RTSettings()`. +The most important settings refer to the Russian roulette system and the grid +cloner (see section on Ray Tracing for details). The settings for the Russian +roulette system include the number of times a ray will be traced +deterministically (`maxiter`) and the probability that a ray that exceeds `maxiter` +is terminated (`pkill`). The grid cloner is used to approximate an infinite canopy +by replicating the scene in the different directions (`nx` and `ny` being the +number of replicates in each direction along the x and y axes, respectively). It +is also possible to turn on parallelization of the ray tracing simulation by +setting `parallel = true` (currently this uses Julia's builtin multithreading +capabilities). + +In addition `RTSettings()`, an acceleration structure and a splitting rule can +be defined when creating the `RayTracer` object (see ray tracing documentation +for details). The acceleration structure allows speeding up the ray tracing +by avoiding testing all rays against all objects in the scene. - function int_dims(biomass, data) - int_biomass = biomass - int_volume = biomass / data.SIW - int_length = cbrt(int_volume * 4 * data.IS^2 / pi) - int_width = int_length / data.IS - return int_length, int_width - end +=# +function create_raytracer(scene, sources) + settings = RTSettings(pkill = 0.9, maxiter = 4, nx = 5, ny = 5, parallel = true) + RayTracer(scene, sources, settings = settings, acceleration = BVH, + rule = SAH{3}(5, 10)); +end +#= - #= - Each day, the total biomass of the tree is updated using a simple RUE formula - and the increment of biomass is distributed across the organs proportionally to - their relative sink strength (of leaves or internodes). - - The sink strength of leaves is modelled with a beta distribution scaled to the - `leaf_expansion` argument that determines the duration of leaf growth, whereas - for the internodes it follows a negative exponential distribution. The `pdf` - function computes the probability density of each distribution which is taken as - proportional to the sink strength (the model is actually source-limited since we - imposed a particular growth rate). - =# - sink_strength(leaf, data) = - leaf.age > data.leaf_expansion ? 0.0 : - pdf(leaf.sink, leaf.age / data.leaf_expansion) / 100.0 - # plot( - # 0:1:50, - # x -> sink_strength(TreeTypes.Leaf(age = x), TreeTypes.treeparams()), - # xlabel = "Age", - # ylabel = "Sink strength", - # label = "Leaf", - # ) - - sink_strength(int) = pdf(int.sink, int.age) - # plot!(0:1:50, x -> sink_strength(TreeTypes.Internode(age = x)), label = "Internode") - - #= - Now we need a function that updates the biomass of the tree, allocates it to the - different organs and updates the dimensions of said organs. For simplicity, - we create the functions `leaves()` and `internodes()` that will apply the queries - to the tree required to extract said nodes: - =# - get_leaves(tree) = apply(tree, Query(TreeTypes.Leaf)) - get_internodes(tree) = apply(tree, Query(TreeTypes.Internode)) - - #= - The age of the different organs is updated every time step: - =# - function age!(all_leaves, all_internodes, all_meristems) - for leaf in all_leaves - leaf.age += 1 - end - for int in all_internodes - int.age += 1 - end - for mer in all_meristems - mer.age += 1 +The actual ray tracing simulation is performed by calling the `trace!()` method +on the ray tracing object. This will trace all rays from all light sources and +update the radiant power absorbed by the different surfaces in the scene inside +the `Material` objects (see `feed!()` above): + +=# +function run_raytracer!(forest; DOY = 182) + scene = create_scene(forest) + sources = create_sky(scene = scene, DOY = DOY) + rtobj = create_raytracer(scene, sources) + trace!(rtobj) + return nothing +end +#= + +The total PAR absorbed for each tree is calculated from the material objects of +the different internodes (using `power()` on the `Material` object). Note that +the `power()` function returns three different values, one for each waveband, +but they are added together as RUE is defined for total PAR. + + +=# +# Run the ray tracer, calculate PAR absorbed per tree and add it to the daily +# total using general weighted quadrature formula +function calculate_PAR!(forest; DOY = 182) + ## Reset PAR absorbed by the tree (at the start of a new day) + reset_PAR!(forest) + ## Run the ray tracer to compute daily PAR absorption + run_raytracer!(forest, DOY = DOY) + ## Add up PAR absorbed by each leaf within each tree + @threads for tree in forest + for l in get_leaves(tree) + data(tree).PAR += power(l.material)[1] end - return nothing end + return nothing +end - #= - The daily growth is allocated to different organs proportional to their sink - strength. - =# - function grow!(tree, all_leaves, all_internodes) - # Compute total biomass increment - tdata = data(tree) - ΔB = max(0.5, tdata.RUE * tdata.PAR / 1e6) # Trick to emulate reserves in seedling - tdata.biomass += ΔB - # Total sink strength - total_sink = 0.0 - for leaf in all_leaves - total_sink += sink_strength(leaf, tdata) - end - for int in all_internodes - total_sink += sink_strength(int) - end - # Allocate biomass to leaves and internodes - for leaf in all_leaves - leaf.biomass += ΔB * sink_strength(leaf, tdata) / total_sink - end - for int in all_internodes - int.biomass += ΔB * sink_strength(int) / total_sink - end - return nothing +# Reset PAR absorbed by the tree (at the start of a new day) +function reset_PAR!(forest) + for tree in forest + data(tree).PAR = 0.0 end + return nothing +end +#= - #= - Finally, we need to update the dimensions of the organs. The leaf dimensions are - =# - function size_leaves!(all_leaves, tdata) - for leaf in all_leaves - leaf.length, leaf.width = leaf_dims(leaf.biomass, tdata) - end - return nothing +### Growth + +We need some functions to compute the length and width of a leaf or internode +from its biomass + +=# +function leaf_dims(biomass, vars) + leaf_biomass = biomass + leaf_area = biomass/vars.SLW + leaf_length = sqrt(leaf_area*4*vars.LS/pi) + leaf_width = leaf_length/vars.LS + return leaf_length, leaf_width +end + +function int_dims(biomass, vars) + int_biomass = biomass + int_volume = biomass/vars.SIW + int_length = cbrt(int_volume*4*vars.IS^2/pi) + int_width = int_length/vars.IS + return int_length, int_width +end +#= + +Each day, the total biomass of the tree is updated using a simple RUE formula +and the increment of biomass is distributed across the organs proportionally to +their relative sink strength (of leaves or internodes). + +The sink strength of leaves is modelled with a beta distribution scaled to the +`leaf_expansion` argument that determines the duration of leaf growth, whereas +for the internodes it follows a negative exponential distribution. The `pdf` +function computes the probability density of each distribution which is taken as +proportional to the sink strength (the model is actually source-limited since we +imposed a particular growth rate). + +=# +sink_strength(leaf, vars) = leaf.age > vars.leaf_expansion ? 0.0 : + pdf(leaf.sink, leaf.age/vars.leaf_expansion)/100.0 +plot(0:1:50, x -> sink_strength(TreeTypes.Leaf(age = x), TreeTypes.treeparams()), + xlabel = "Age", ylabel = "Sink strength", label = "Leaf") +#= + +=# +sink_strength(int) = pdf(int.sink, int.age) +plot!(0:1:50, x -> sink_strength(TreeTypes.Internode(age = x)), label = "Internode") +#= + +Now we need a function that updates the biomass of the tree, allocates it to the +different organs and updates the dimensions of said organs. For simplicity, +we create the functions `leaves()` and `internodes()` that will apply the queries +to the tree required to extract said nodes: + +=# +get_leaves(tree) = apply(tree, Query(TreeTypes.Leaf)) +get_internodes(tree) = apply(tree, Query(TreeTypes.Internode)) +#= + +The age of the different organs is updated every time step: + +=# +function age!(all_leaves, all_internodes, all_meristems) + for leaf in all_leaves + leaf.age += 1 end - function size_internodes!(all_internodes, tdata) - for int in all_internodes - int.length, int.width = int_dims(int.biomass, tdata) - end - return nothing + for int in all_internodes + int.age += 1 end - - #= - ### Daily step - - All the growth and developmental functions are combined together into a daily - step function that updates the forest by iterating over the different trees in - parallel. - =# - get_meristems(tree) = apply(tree, Query(TreeTypes.Meristem)) - function daily_step!(forest, DOY) - # Compute PAR absorbed by each tree - calculate_PAR!(forest, DOY = DOY) - # Grow the trees - @threads for tree in forest - # Retrieve all the relevant organs - all_leaves = get_leaves(tree) - all_internodes = get_internodes(tree) - all_meristems = get_meristems(tree) - # Update the age of the organs - age!(all_leaves, all_internodes, all_meristems) - # Grow the tree - grow!(tree, all_leaves, all_internodes) - tdata = data(tree) - size_leaves!(all_leaves, tdata) - size_internodes!(all_internodes, tdata) - # Developmental rules - rewrite!(tree) - end + for mer in all_meristems + mer.age += 1 end + return nothing +end +#= - #= - ### Initialization - - The trees are initialized on a regular grid with random values for the initial - orientation and RUE: - =# - RUEs = rand(Normal(1.5, 0.2), 10, 10) - # histogram(vec(RUEs)) - - orientations = [rand() * 360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0] - # histogram(vec(orientations)) - - origins = [Vec(i, j, 0) for i = 1:2.0:20.0, j = 1:2.0:20.0] - - #= - The following initalizes a tree based on the origin, orientation and RUE: - =# - function create_tree(origin, orientation, RUE) - # Initial state and parameters of the tree - data = TreeTypes.treeparams(RUE = RUE) - # Initial states of the leaves - leaf_length, leaf_width = leaf_dims(data.LB0, data) - vleaf = (biomass = data.LB0, length = leaf_length, width = leaf_width) - # Initial states of the internodes - int_length, int_width = int_dims(data.LB0, data) - vint = (biomass = data.IB0, length = int_length, width = int_width) - # Growth rules - meristem_rule = create_meristem_rule(vleaf, vint) - branch_rule = create_branch_rule(vint) - axiom = - T(origin) + - RH(orientation) + - TreeTypes.Internode( - biomass = vint.biomass, - length = vint.length, - width = vint.width, - ) + - TreeTypes.Meristem() - tree = Graph(axiom = axiom, rules = (meristem_rule, branch_rule), data = data) - return tree - end +The daily growth is allocated to different organs proportional to their sink +strength. - #= - ## Visualization - - As in the previous example, it makes sense to visualize the forest with a soil - tile beneath it. Unlike in the previous example, we will construct the soil tile - using a dedicated graph and generate a `Scene` object which can later be - merged with the rest of scene generated in daily step: - =# - Base.@kwdef struct Soil <: VirtualPlantLab.Node - length::Float64 - width::Float64 +=# +function grow!(tree, all_leaves, all_internodes) + ## Compute total biomass increment + tdata = data(tree) + ΔB = max(0.5, tdata.RUE*tdata.PAR/1e6) # Trick to emulate reserves in seedling + tdata.biomass += ΔB + ## Total sink strength + total_sink = 0.0 + for leaf in all_leaves + total_sink += sink_strength(leaf, tdata) + end + for int in all_internodes + total_sink += sink_strength(int) + end + ## Allocate biomass to leaves and internodes + for leaf in all_leaves + leaf.biomass += ΔB*sink_strength(leaf, tdata)/total_sink end - function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, data) - Rectangle!( - turtle, - length = s.length, - width = s.width, - color = RGB(255 / 255, 236 / 255, 179 / 255), - ) + for int in all_internodes + int.biomass += ΔB*sink_strength(int)/total_sink end - soil_graph = - RA(-90.0) + - T(Vec(0.0, 10.0, 0.0)) + # Moves into position - Soil(length = 20.0, width = 20.0) # Draws the soil tile - soil = Scene(Graph(axiom = soil_graph)) - render(soil, axes = false) - - #= - And the following function renders the entire scene (notice that we need to - use `display()` to force the rendering of the scene when called within a loop - or a function): - =# - function render_forest(forest, soil) - scene = Scene(vec(forest)) # create scene from forest - scene = Scene([scene, soil]) # merges the two scenes - display(render(scene)) + return nothing +end +#= + +Finally, we need to update the dimensions of the organs. The leaf dimensions are + +=# +function size_leaves!(all_leaves, tvars) + for leaf in all_leaves + leaf.length, leaf.width = leaf_dims(leaf.biomass, tvars) end + return nothing +end +function size_internodes!(all_internodes, tvars) + for int in all_internodes + int.length, int.width = int_dims(int.biomass, tvars) + end + return nothing +end +#= - #= - ## Simulation - - We can now create a forest of trees on a regular grid: - =# - forest = create_tree.(origins, orientations, RUEs) - render_forest(forest, soil) - start = 180 - for i = 1:5 - println("Day $i") - daily_step!(forest, i + start) - if mod(i, 5) == 0 - soil = Scene(Graph(axiom = soil_graph)) # Otherwise soil is converted into a mesh? - render_forest(forest, soil) - end +### Daily step + +All the growth and developmental functions are combined together into a daily +step function that updates the forest by iterating over the different trees in +parallel. + +=# +get_meristems(tree) = apply(tree, Query(TreeTypes.Meristem)) +function daily_step!(forest, DOY) + ## Compute PAR absorbed by each tree + calculate_PAR!(forest, DOY = DOY) + ## Grow the trees + @threads for tree in forest + ## Retrieve all the relevant organs + all_leaves = get_leaves(tree) + all_internodes = get_internodes(tree) + all_meristems = get_meristems(tree) + ## Update the age of the organs + age!(all_leaves, all_internodes, all_meristems) + ## Grow the tree + grow!(tree, all_leaves, all_internodes) + tdata = data(tree) + size_leaves!(all_leaves, tdata) + size_internodes!(all_internodes, tdata) + ## Developmental rules + rewrite!(tree) end +end +#= + +### Initialization + +The trees are initialized on a regular grid with random values for the initial +orientation and RUE: + +=# +RUEs = rand(Normal(1.5,0.2), 10, 10) +histogram(vec(RUEs)) +#= + +=# +orientations = [rand()*360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0] +histogram(vec(orientations)) +#= + +=# +origins = [Vec(i,j,0) for i = 1:2.0:20.0, j = 1:2.0:20.0]; +#= + +The following initalizes a tree based on the origin, orientation and RUE: + +=# +function create_tree(origin, orientation, RUE) + ## Initial state and parameters of the tree + data = TreeTypes.treeparams(RUE = RUE) + ## Initial states of the leaves + leaf_length, leaf_width = leaf_dims(data.LB0, data) + vleaf = (biomass = data.LB0, length = leaf_length, width = leaf_width) + ## Initial states of the internodes + int_length, int_width = int_dims(data.LB0, data) + vint = (biomass = data.IB0, length = int_length, width = int_width) + ## Growth rules + meristem_rule = create_meristem_rule(vleaf, vint) + branch_rule = create_branch_rule(vint) + axiom = T(origin) + RH(orientation) + + TreeTypes.Internode(biomass = vint.biomass, + length = vint.length, + width = vint.width) + + TreeTypes.Meristem() + tree = Graph(axiom = axiom, rules = (meristem_rule, branch_rule), + data = data) + return tree +end +#= + + +## Visualization + +As in the previous example, it makes sense to visualize the forest with a soil +tile beneath it. Unlike in the previous example, we will construct the soil tile +using a dedicated graph and generate a `Scene` object which can later be +merged with the rest of scene generated in daily step: + +=# +Base.@kwdef struct Soil <: VirtualPlantLab.Node + length::Float64 + width::Float64 +end +function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, data) + Rectangle!(turtle, length = s.length, width = s.width, color = RGB(255/255, 236/255, 179/255)) +end +soil_graph = RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + ## Moves into position + Soil(length = 20.0, width = 20.0) ## Draws the soil tile +soil = Scene(Graph(axiom = soil_graph)); +render(soil, axes = false) +#= + +And the following function renders the entire scene (notice that we need to +use `display()` to force the rendering of the scene when called within a loop +or a function): + +=# +function render_forest(forest, soil) + scene = Scene(vec(forest)) ## create scene from forest + scene = Scene([scene, soil]) ## merges the two scenes + display(render(scene)) +end +#= + +## Simulation + +We can now create a forest of trees on a regular grid: +=# +forest = create_tree.(origins, orientations, RUEs); +render_forest(forest, soil) +start = 180 +for i in 1:20 + println("Day $i") + daily_step!(forest, i + start) + if mod(i, 5) == 0 + render_forest(forest, soil) + end end diff --git a/test/relationalqueries.jl b/test/relationalqueries.jl index c4f87d2..a1c59c7 100644 --- a/test/relationalqueries.jl +++ b/test/relationalqueries.jl @@ -1,325 +1,326 @@ #= -Relational queries -Alejandro Morales Sierra +# Relational queries + +Alejandro Morales + Centre for Crop Systems Analysis - Wageningen University -In this example we illustrate how to test relationships among nodes inside -queries. Relational queries allow to establish relationships between nodes in -the graph, which generally requires a intimiate knowledge of the graph. For this -reason, relational queries are inheretly complex as graphs can become complex -and there may be solutions that do not require relational queries in many -instances. Nevertheless, they are integral part of VirtualPlantLab and can sometimes be -useful. As they can be hard to grasp, this tutorial will illustrate with a -relatively simple graph a series of relational queries with increasing -complexity with the aim that users will get a better understanding of relational -queries. For this purpose, an abstract graph with several branching levels will -be used, so that we can focus on the relations among the nodes without being -distracted by case-specific details. - -The graph will be composed of two types of nodes: the inner nodes (`A` and `C`) -and the leaf nodes (`B`). Each leaf node will be identified uniquely with an -index and the objective is to write queries that can identify a specific subset -of the leaf nodes, without using the data stored in the nodes themselves. That -is, the queries should select the right nodes based on their relationships to -the rest of nodes in the graph. Note that `C` nodes contain a single value that -may be positive or negative, whereas `A` nodes contain no data. +In this example we illustrate how to test relationships among nodes inside queries. +Relational queries allow to establish relationships between nodes in the graph, +which generally requires a intimiate knowledge of the graph. For this reason, +relational queries are inheretly complex as graphs can become complex and there +may be solutions that do not require relational queries in many instances. +Nevertheless, they are integral part of VPL and can sometimes be useful. As they +can be hard to grasp, this tutorial will illustrate with a relatively simple +graph a series of relational queries with increasing complexity with the aim +that users will get a better understanding of relational queries. For this purpose, +an abstract graph with several branching levels will be used, so that we can focus +on the relations among the nodes without being distracted by case-specific details. + +The graph will be composed of two types of nodes: the inner nodes (`A` and `C`) and the +leaf nodes (`B`). Each leaf node will be identified uniquely with an index and +the objective is to write queries that can identify a specific subset of the leaf +nodes, without using the data stored in the nodes themselves. That is, the queries +should select the right nodes based on their relationships to the rest of nodes +in the graph. Note that `C` nodes contain a single value that may be positive or negative, +whereas `A` nodes contain no data. As usual, we start with defining the types of nodes in the graph + =# using VirtualPlantLab import GLMakie module Queries -using VirtualPlantLab -struct A <: Node end + using VirtualPlantLab + struct A <: Node end -struct C <: Node - val::Float64 + struct C <: Node + val::Float64 + end + + struct B <: Node + ID::Int + end end +import .Queries as Q +#= + +We generate the graph directly, rather than with rewriting rules. The graph has +a motif that is repeated three times (with a small variation), so we can create +the graph in a piecewise manner. Note how we can use the function `sum` to add +nodes to the graph (i.e. `sum(A() for i in 1:3)` is equivalent to `A() + A() + A()`) + +=# +motif(n, i = 0) = Q.A() + (Q.C(45.0) + Q.A() + (Q.C(45.0) + Q.A() + Q.B(i + 1), + Q.C(-45.0) + Q.A() + Q.B(i + 2), + Q.A() + Q.B(i + 3)), + Q.C(- 45.0) + sum(Q.A() for i in 1:n) + Q.B(i + 4)) +axiom = motif(3, 0) + motif(2, 4) + motif(1, 8) + Q.A() + Q.A() + Q.B(13) +graph = Graph(axiom = axiom); +draw(graph) +#= + +By default, VPL will use as node label the type of node and the internal ID generated by VPL itself. This ID is useful if we want to +extract a particular node from the graph, but it is not controlled by the user. However, the user can specialized the function `node_label()` +to specify exactly how to label the nodes of a particular type. In this case, we want to just print `A` or `C` for nodes of type `A` and `C`, whereas +for nodes of type `B` we want to use the `ID` field that was stored inside the node during the graph generation. + +=# +VirtualPlantLab.node_label(n::Q.B, id) = "B-$(n.ID)" +VirtualPlantLab.node_label(n::Q.A, id) = "A" +VirtualPlantLab.node_label(n::Q.C, id) = "C" +draw(graph) +#= + +To clarify, the `id` argument of the function `node_label()` refers to the internal id generated by VPL (used by the default method for `node_label()`, whereas the the first argument is the data stored inside a node (in the case of `B` nodes, there is a field called `ID` that will not be modified by VPL as that is user-provided data). + +The goal of this exercise is then to write queries that retrieve specific `B` +nodes without using the data stored in the node in the query. That is, we have +to identify nodes based on their relationships to other nodes. + +## All nodes of type `B` + +First, we create the query object. In this case, there is no special condition as +we want to retrieve all the nodes of type `B` + +=# +Q1 = Query(Q.B) +#= + +Applying the query to the graph returns an array with all the `B` nodes + +=# +A1 = apply(graph, Q1) +#= -struct B <: Node - ID::Int + +For the remainder of this tutorial, the code will be hidden by default to allow users to try on their own. + + +## Node containing value 13 + +Since the `B` node 13 is the leaf node of the main branch of the graph (e.g. this could be the apical meristem of the main stem of a plant), there +are no rotations between the root node of the graph and this node. Therefore, +the only condition require to single out this node is that it has no ancestor +node of type `C`. + +Checking whether a node has an ancestor that meets a certain +condition can be achieved with the function `hasAncestor()`. Similarly to the +condition of the `Query` object, the `hasAncestor()` function also has a condition, +in this case applied to the parent node of the node being tested, and moving +upwards in the graph recursively (until reaching the root node). Note that, in +order to access the object stored inside the node, we need to use the `data()` +function, and then we can test if that object is of type `C`. The `B` node 13 +is the only node for which `hasAncestor()` should return `false`: + +=# +function Q2_fun(n) + check, steps = has_ancestor(n, condition = x -> data(x) isa Q.C) + !check end +#= + +As before, we just need to apply the `Query` object to the graph: + +=# +Q2 = Query(Q.B, condition = Q2_fun) +A2 = apply(graph, Q2) +#= + +## Nodes containing values 1, 2 and 3 + +These three nodes belong to one of the branch motifs repeated through the graph. Thus, +we need to identify the specific motif they belong to and chose all the `B` nodes +inside that motif. The motif is defined by an `A` node that has a `C` child with +a negative `val` and parent node `C` with positive `val`. This `A` node +should then be 2 nodes away from the root node to separate it from upper repetitions +of the motif. +Therefore, we need to test for two conditions, first find those nodes inside a +branch motif, then retrieve the root of the branch motif (i.e., the `A` node +described in the above) and then check the distance of that node from the root: + +=# +function branch_motif(p) + data(p) isa Q.A && + has_descendant(p, condition = x -> data(x) isa Q.C && data(x).val < 0.0)[1] && + has_ancestor(p, condition = x -> data(x) isa Q.C && data(x).val > 0.0)[1] end -import .Queries as Q -let - #= - We generate the graph directly, rather than with rewriting rules. The graph has - a motif that is repeated three times (with a small variation), so we can create - the graph in a piecewise manner. Note how we can use the function `sum` to add - nodes to the graph (i.e. `sum(A() for i in 1:3)` is equivalent to `A() + A() + - A()`) - =# - motif(n, i = 0) = - Q.A() + ( - Q.C(45.0) + - Q.A() + - ( - Q.C(45.0) + Q.A() + Q.B(i + 1), - Q.C(-45.0) + Q.A() + Q.B(i + 2), - Q.A() + Q.B(i + 3), - ), - Q.C(-45.0) + sum(Q.A() for i = 1:n) + Q.B(i + 4), - ) - axiom = motif(3, 0) + motif(2, 4) + motif(1, 8) + Q.A() + Q.A() + Q.B(13) - graph = Graph(axiom = axiom) - draw(graph) - - #= - By default, VirtualPlantLab will use as node label the type of node and the internal ID - generated by VirtualPlantLab itself. This ID is useful if we want to extract a particular - node from the graph, but it is not controlled by the user. However, the user can - specialized the function `node_label()` to specify exactly how to label the - nodes of a particular type. In this case, we want to just print `A` or `C` for - nodes of type `A` and `C`, whereas for nodes of type `B` we want to use the `ID` - field that was stored inside the node during the graph generation. - =# - VirtualPlantLab.node_label(n::Q.B, id) = "B-$(n.ID)" - VirtualPlantLab.node_label(n::Q.A, id) = "A" - VirtualPlantLab.node_label(n::Q.C, id) = "C" - draw(graph) - - #= - To clarify, the `id` argument of the function `node_label()` refers to the - internal id generated by VirtualPlantLab (used by the default method for `node_label()`, - whereas the the first argument is the data stored inside a node (in the case of - `B` nodes, there is a field called `ID` that will not be modified by VirtualPlantLab as that - is user-provided data). The goal of this exercise is then to write queries that - retrieve specific `B` nodes (without using the data stored in the node in the - query, that is, we have to identify nodes based on their topological connections - to other nodes). - - ## All nodes of type `B` - - First, we create the query object. In this case, there is no special condition - as we want to retrieve all the nodes of type `B` - =# - Q1 = Query(Q.B) - - # Applying the query to the graph returns an array with all the `B` nodes - A1 = apply(graph, Q1) - - #= - For the remainder of this tutorial, the code will be hidden by default to allow - users to try on their own. - - ## Node containing value 13 - - Since the `B` node 13 is the leaf node of the main branch of the graph (e.g. - this could be the apical meristem of the main stem of a plant), there are no - rotations between the root node of the graph and this node. Therefore, the only - condition require to single out this node is that it has no ancestor node of - type `C`. - - Checking whether a node has an ancestor that meets a certain condition can be - achieved with the function `has_ancestor()`. Similarly to the condition of the - `Query` object, the `has_ancestor()` function also has a condition, in this case - applied to the parent node of the node being tested, and moving upwards in the - graph recursively (until reaching the root node). Note that, in order to access - the object stored inside the node, we need to use the `data()` function, and - then we can test if that object is of type `C`. The `B` node 13 is the only node - for which `has_ancestor()` should return `false`: - =# - function Q2_fun(n) - check, steps = has_ancestor(n, condition = x -> data(x) isa Q.C) - !check - end +function Q3_fun(n, nsteps) + ## Condition 1 + check, steps = has_ancestor(n, condition = branch_motif) + !check && return false + ## Condition 2 + p = parent(n, nsteps = steps) + check, steps = has_ancestor(p, condition = is_root) + steps != nsteps && return false + return true +end +#= - # As before, we just need to apply the `Query` object to the graph: - Q2 = Query(Q.B, condition = Q2_fun) - A2 = apply(graph, Q2) - - #= - ## Nodes containing values 1, 2 and 3 - - These three nodes belong to one of the branch motifs repeated through the graph. - Thus, we need to identify the specific motif they belong to and chose all the - `B` nodes inside that motif. The motif is defined by an `A` node that has a `C` - child with a negative `val` and parent node `C` with positive `val`. This `A` - node should then be 2 nodes away from the root node to separate it from upper - repetitions of the motif. Therefore, we need to test for two conditions, first - find those nodes inside a branch motif, then retrieve the root of the branch - motif (i.e., the `A` node described in the above) and then check the distance of - that node from the root: - =# - function branch_motif(p) - data(p) isa Q.A && - has_descendant(p, condition = x -> data(x) isa Q.C && data(x).val < 0.0)[1] && - has_ancestor(p, condition = x -> data(x) isa Q.C && data(x).val > 0.0)[1] - end +And applying the query to the object results in the required nodes: - function Q3_fun(n, nsteps) - # Condition 1 - check, steps = has_ancestor(n, condition = branch_motif) - !check && return false - # Condition 2 - p = parent(n, nsteps = steps) - check, steps = has_ancestor(p, condition = is_root) - steps != nsteps && return false - return true - end +=# +Q3 = Query(Q.B, condition = n -> Q3_fun(n, 2)) +A3 = apply(graph, Q3) +#= - # And applying the query to the object results in the required nodes: - Q3 = Query(Q.B, condition = n -> Q3_fun(n, 2)) - A3 = apply(graph, Q3) - - #= - ## Node containing value 4 - - The node `B` with value 4 can be singled-out because there is no branching point - between the root node and this node. This means that no ancestor node should - have more than one children node except the root node. Remember that - `has_ancestor()` returns two values, but we are only interested in the first - value. You do not need to assign the returned object from a Julia function, you - can just index directly the element to be selected from the returned tuple: - =# - function Q4_fun(n) - !has_ancestor(n, condition = x -> is_root(x) && length(children(x)) > 1)[1] - end +## Node containing value 4 + +The node `B` with value 4 can be singled-out because there is no branching point +between the root node and this node. This means that no ancestor node should have +more than one children node except the root node. Remember that `hasAncestor()` +returns two values, but we are only interested in the first value. You do not need to +assign the returned object from a Julia function, you can just index directly the element +to be selected from the returned tuple: + +=# +function Q4_fun(n) + !has_ancestor(n, condition = x -> is_root(x) && length(children(x)) > 1)[1] +end +#= - # And applying the query to the object results in the required node: - Q4 = Query(Q.B, condition = Q4_fun) - A4 = apply(graph, Q4) +And applying the query to the object results in the required node: - #= - ## Node containing value 3 +=# +Q4 = Query(Q.B, condition = Q4_fun) +A4 = apply(graph, Q4) +#= - This node is the only `B` node that is four steps from the root node, which we - can retrieve from the second argument returned by `has_ancestor()`: - =# - function Q5_fun(n) - check, steps = has_ancestor(n, condition = is_root) - steps == 4 - end +## Node containing value 3 - Q5 = Query(Q.B, condition = Q5_fun) - A5 = apply(graph, Q5) - - #= - ## Node containing value 7 - - Node `B` 7 belongs to the second lateral branch motif and the second parent node - is of type `A`. Note that we can reuse the `Q3_fun` from before in the condition - required for this node: - =# - function Q6_fun(n, nA) - check = Q3_fun(n, nA) - !check && return false - p2 = parent(n, nsteps = 2) - data(p2) isa Q.A - end +This node is the only `B` node that is four steps from the root node, which we can +retrieve from the second argument returned by `hasAncestor()`: - Q6 = Query(Q.B, condition = n -> Q6_fun(n, 3)) - A6 = apply(graph, Q6) - - #= - ## Nodes containing values 11 and 13 - - The `B` nodes 11 and 13 actually have different relationships to the rest of the - graph, so we just need to define two different condition functions and combine - them. The condition for the `B` node 11 is similar to the `B` node 7, whereas - the condition for node 13 was already constructed before, so we just need to - combined them with an OR operator: - =# - Q7 = Query(Q.B, condition = n -> Q6_fun(n, 4) || Q2_fun(n)) - A7 = apply(graph, Q7) - - #= - ## Nodes containing values 1, 5 and 9 - - These nodes play the same role in the three lateral branch motifs. They are the - only `B` nodes preceded by the sequence A C+ A. We just need to check the - sequence og types of objects for the the first three parents of each `B` node: - =# - function Q8_fun(n) - p1 = parent(n) - p2 = parent(n, nsteps = 2) - p3 = parent(n, nsteps = 3) - data(p1) isa Q.A && data(p2) isa Q.C && data(p2).val > 0.0 && data(p3) isa Q.A - end +=# +function Q5_fun(n) + check, steps = has_ancestor(n, condition = is_root) + steps == 4 +end - Q8 = Query(Q.B, condition = Q8_fun) - A8 = apply(graph, Q8) - - #= - ## Nodes containing values 2, 6 and 10 - - This exercise is similar to the previous one, but the C node has a negative - `val`. The problem is that node 12 would also match the pattern A C- A. We can - differentiate between this node and the rest by checking for a fourth ancestor - node of class `C`: - =# - function Q9_fun(n) - p1 = parent(n) - p2 = parent(n, nsteps = 2) - p3 = parent(n, nsteps = 3) - p4 = parent(n, nsteps = 4) - data(p1) isa Q.A && - data(p2) isa Q.C && - data(p2).val < 0.0 && - data(p3) isa Q.A && - data(p4) isa Q.C - end +Q5 = Query(Q.B, condition = Q5_fun) +A5 = apply(graph, Q5) +#= - Q9 = Query(Q.B, condition = Q9_fun) - A9 = apply(graph, Q9) - - #= - ## Nodes containing values 6, 7 and 8 - - We already came up with a condition to extract node 7. We can also modify the - previous condition so that it only node 6. Node 8 can be identified by checking - for the third parent node being of type `C` and being 5 nodes from the root of - the graph. - - As always, we can reusing previous conditions since they are just regular Julia - functions: - =# - function Q10_fun(n) - Q6_fun(n, 3) && return true # Check node 7 - Q9_fun(n) && has_ancestor(n, condition = is_root)[2] == 6 && return true # Check node 6 - has_ancestor(n, condition = is_root)[2] == 5 && - data(parent(n, nsteps = 3)) isa Q.C && - return true # Check node 8 (and not 4!) - end +## Node containing value 7 - Q10 = Query(Q.B, condition = Q10_fun) - A10 = apply(graph, Q10) - - #= - ## Nodes containig values 3, 7, 11 and 12 - - We already have conditions to select nodes 3, 7 and 11 so we just need a new - condition for node 12 (similar to the condition for 8). - =# - function Q11_fun(n) - Q5_fun(n) && return true # 3 - Q6_fun(n, 3) && return true # 7 - Q6_fun(n, 4) && return true # 11 - has_ancestor(n, condition = is_root)[2] == 5 && - data(parent(n, nsteps = 2)) isa Q.C && - data(parent(n, nsteps = 4)) isa Q.A && - return true # 12 - end +Node `B` 7 belongs to the second lateral branch motif and the second parent +node is of type `A`. Note that we can reuse the `Q3_fun` from before in the +condition required for this node: - Q11 = Query(Q.B, condition = Q11_fun) - A11 = apply(graph, Q11) +=# +function Q6_fun(n, nA) + check = Q3_fun(n, nA) + !check && return false + p2 = parent(n, nsteps = 2) + data(p2) isa Q.A +end - #= - ## Nodes containing values 7 and 12 +Q6 = Query(Q.B, condition = n -> Q6_fun(n, 3)) +A6 = apply(graph, Q6) +#= + +## Nodes containing values 11 and 13 + +The `B` nodes 11 and 13 actually have different relationships to the rest of the graph, +so we just need to define two different condition functions and combine them. +The condition for the `B` node 11 is similar to the `B` node 7, whereas the condition +for node 13 was already constructed before, so we just need to combined them with an +OR operator: + +=# +Q7 = Query(Q.B, condition = n -> Q6_fun(n, 4) || Q2_fun(n)) +A7 = apply(graph, Q7) +#= - We just need to combine the conditions for the nodes 7 and 12 - =# - function Q12_fun(n) - Q6_fun(n, 3) && return true # 7 - has_ancestor(n, condition = is_root)[2] == 5 && - data(parent(n, nsteps = 2)) isa Q.C && - data(parent(n, nsteps = 4)) isa Q.A && - return true # 12 - end - Q12 = Query(Q.B, condition = Q12_fun) - A12 = apply(graph, Q12) +## Nodes containing values 1, 5 and 9 +These nodes play the same role in the three lateral branch motifs. They are the +only `B` nodes preceded by the sequence A C+ A. We just need to check the +sequence og types of objects for the the first three parents of each `B` node: + +=# +function Q8_fun(n) + p1 = parent(n) + p2 = parent(n, nsteps = 2) + p3 = parent(n, nsteps = 3) + data(p1) isa Q.A && data(p2) isa Q.C && data(p2).val > 0.0 && data(p3) isa Q.A end + +Q8 = Query(Q.B, condition = Q8_fun) +A8 = apply(graph, Q8) +#= + +## Nodes contaning values 2, 6 and 10 + +This exercise is similar to the previous one, but the C node has a negative +`val`. The problem is that node 12 would also match the pattern A C- A. We +can differentiate between this node and the rest by checking for a fourth +ancestor node of class `C`: + +=# +function Q9_fun(n) + p1 = parent(n) + p2 = parent(n, nsteps = 2) + p3 = parent(n, nsteps = 3) + p4 = parent(n, nsteps = 4) + data(p1) isa Q.A && data(p2) isa Q.C && data(p2).val < 0.0 && + data(p3) isa Q.A && data(p4) isa Q.C +end + +Q9 = Query(Q.B, condition = Q9_fun) +A9 = apply(graph, Q9) +#= + +## Nodes containg values 6, 7 and 8 + +We already came up with a condition to extract node 7. We can also modify the previous +condition so that it only node 6. Node 8 can be identified by checking for the third +parent node being of type `C` and being 5 nodes from the root of the graph. + +As always, we can reusing previous conditions since they are just regular Julia functions: + +=# +function Q10_fun(n) + Q6_fun(n, 3) && return true ## Check node 7 + Q9_fun(n) && has_ancestor(n, condition = is_root)[2] == 6 && return true ## Check node 6 + has_ancestor(n, condition = is_root)[2] == 5 && data(parent(n, nsteps = 3)) isa Q.C && return true ## Check node 8 (and not 4!) +end + +Q10 = Query(Q.B, condition = Q10_fun) +A10 = apply(graph, Q10) +#= + +## Nodes containig values 3, 7, 11 and 12 + +We already have conditions to select nodes 3, 7 and 11 so we just need a new condition +for node 12 (similar to the condition for 8). + +=# +function Q11_fun(n) + Q5_fun(n) && return true ## 3 + Q6_fun(n, 3) && return true ## 7 + Q6_fun(n, 4) && return true ## 11 + has_ancestor(n, condition = is_root)[2] == 5 && data(parent(n, nsteps = 2)) isa Q.C && + data(parent(n, nsteps = 4)) isa Q.A && return true ## 12 +end + +Q11 = Query(Q.B, condition = Q11_fun) +A11 = apply(graph, Q11) +#= + +## Nodes containing values 7 and 12 + +We just need to combine the conditions for the nodes 7 and 12 + +=# +function Q12_fun(n) + Q6_fun(n, 3) && return true # 7 + has_ancestor(n, condition = is_root)[2] == 5 && data(parent(n, nsteps = 2)) isa Q.C && + data(parent(n, nsteps = 4)) isa Q.A && return true ## 12 +end + +Q12 = Query(Q.B, condition = Q12_fun) +A12 = apply(graph, Q12) diff --git a/test/runtests.jl b/test/runtests.jl index 81bfe50..3d82f9b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,27 +11,43 @@ end # Tutorials @testset "Tutorials" begin @testset "Algae" begin - include("algae.jl") + let + include("algae.jl") + end end @testset "Snowflakes" begin - include("snowflakes.jl") + let + include("snowflakes.jl") + end end @testset "Tree" begin - include("tree.jl") + let + include("tree.jl") + end end @testset "Forest" begin - include("forest.jl") + let + include("forest.jl") + end end @testset "Growth forest" begin - include("growthforest.jl") + let + include("growthforest.jl") + end end @testset "Raytraced forest" begin - include("raytracedforest.jl") + let + include("raytracedforest.jl") + end end @testset "Context" begin - include("context.jl") + let + include("context.jl") + end end @testset "Relational queries" begin - include("relationalqueries.jl") + let + include("relationalqueries.jl") + end end end diff --git a/test/snowflakes.jl b/test/snowflakes.jl index 4d6bfb3..fa10ea7 100644 --- a/test/snowflakes.jl +++ b/test/snowflakes.jl @@ -1,8 +1,11 @@ #= -The Koch snowflake -Alejandro Morales Sierra +# The Koch snowflake + +Alejandro Morales + Centre for Crop Systems Analysis - Wageningen University + In this example, we create a Koch snowflake, which is one of the earliest fractals to be described. The Koch snowflake is a closed curve composed on multiple of segments of different lengths. Starting with an equilateral @@ -11,11 +14,11 @@ length arrange in a specific manner. Graphically, the first four iterations of the Koch snowflake construction process result in the following figures (the green segments are shown as guides but they are not part of the snowflake): -![First four iterations fo Koch snowflake fractal](./KochWikipedia.png) +![First four iterations fo Koch snowflake fractal](https://upload.wikimedia.org/wikipedia/commons/8/8e/KochFlake.png) -In order to implement the construction process of a Koch snowflake in VirtualPlantLab we +In order to implement the construction process of a Koch snowflake in VPL we need to understand how a 3D structure can be generated from a graph of nodes. -VirtualPlantLab uses a procedural approach to generate of structure based on the concept of +VPL uses a procedural approach to generate of structure based on the concept of turtle graphics. The idea behind this approach is to imagine a turtle located in space with a @@ -26,7 +29,7 @@ node may also include instructions to move and/or rotate the turtle, which allows to alter the relative position of the different 3D structures described by a graph. -The construction process of the Koch snowflake in VirtualPlantLab could then be represented +The construction process of the Koch snowflake in VPL could then be represented by the following axiom and rewriting rule: axiom: E(L) + RU(120) + E(L) + RU(120) + E(L) @@ -37,165 +40,168 @@ represents a rotation of the turtle around the upward axis, with angle of rotation given in parenthesis in hexadecimal degrees. The rule can be visualized as follows: -![Koch construction rule](./Koch_order_1.png) +![Koch construction rule](https://python-with-science.readthedocs.io/en/latest/_images/koch_order_1.png) -Note that VirtualPlantLab already provides several classes for common turtle movements and +Note that VPL already provides several classes for common turtle movements and rotations, so our implementation of the Koch snowflake only needs to define a class to implement the edges of the snowflake. This can be achieved as follows: + =# using VirtualPlantLab -import GLMakie # Import rather than "using" to avoid masking Scene -using ColorTypes # To define colors for the rendering +import GLMakie ## Import rather than "using" to avoid masking Scene +using ColorTypes ## To define colors for the rendering module sn -import VirtualPlantLab -struct E <: VirtualPlantLab.Node - length::Float64 -end + import VirtualPlantLab + struct E <: VirtualPlantLab.Node + length::Float64 + end end import .sn +#= -let - - #Note that nodes of type E need to keep track of the length as illustrated in - #the above. The axiom is straightforward: - L = 1.0 - axiom = sn.E(L) + VirtualPlantLab.RU(120.0) + sn.E(L) + VirtualPlantLab.RU(120.0) + sn.E(L) - - #The rule is also straightforward to implement as all the nodes of type E will - #be replaced in each iteration. However, we need to ensure that the length of - #the new edges is a calculated from the length of the edge being replaced. In - #order to extract the data stored in the node being replaced we can simply use - #the function data. In this case, the replacement function is defined and then - #added to the rule. This can make the code more readable but helps debugging and - #testing the replacement function. - function Kochsnowflake(x) - L = data(x).length - sn.E(L / 3) + - RU(-60.0) + - sn.E(L / 3) + - RU(120.0) + - sn.E(L / 3) + - RU(-60.0) + - sn.E(L / 3) - end - rule = Rule(sn.E, rhs = Kochsnowflake) - - # The model is then created by constructing the graph - Koch = Graph(axiom = axiom, rules = Tuple(rule)) - - #= - In order to be able to generate a 3D structure we need to define a method for - the function `VirtualPlantLab.feed!` (notice the need to prefix it with `VirtualPlantLab.` as we are - going to define a method for this function). The method needs to two take two - arguments, the first one is always an object of type Turtle and the second is an - object of the type for which the method is defined (in this case, E). - - The body of the method should generate the 3D structures using the geometry - primitives provided by VirtualPlantLab and feed them to the turtle that is being passed to - the method as first argument. In this case, we are going to represent the edges - of the Koch snowflakes with cylinders, which can be generated with the - `HollowCylinder!` function from VirtualPlantLab. Note that the `feed!` should return - `nothing`, the turtle will be modified in place (hence the use of `!` at the end - of the function as customary in the VirtualPlantLab community). - - In order to render the geometry, we need assign a `color` (i.e., any type of - color support by the package ColorTypes.jl). In this case, we just feed a basic - `RGB` color defined by the proportion of red, green and blue. To make the - figures more appealing, we can assign random values to each channel of the color - to generate random colors. - =# - function VirtualPlantLab.feed!(turtle::Turtle, e::sn.E, vars) - HollowCylinder!( - turtle, - length = e.length, - width = e.length / 10, - height = e.length / 10, - move = true, - color = RGB(rand(), rand(), rand()), - ) - return nothing - end +Note that nodes of type E need to keep track of the length as illustrated in the +above. The axiom is straightforward: + +=# +const L = 1.0 +axiom = sn.E(L) + VirtualPlantLab.RU(120.0) + sn.E(L) + VirtualPlantLab.RU(120.0) + sn.E(L) +#= + +The rule is also straightforward to implement as all the nodes of type E will be +replaced in each iteration. However, we need to ensure that the length of the +new edges is a calculated from the length of the edge being replaced. In order +to extract the data stored in the node being replaced we can simply use the +function data. In this case, the replacement function is defined and then added +to the rule. This can make the code more readable but helps debugging and +testing the replacement function. + +=# +function Kochsnowflake(x) + L = data(x).length + sn.E(L/3) + RU(-60.0) + sn.E(L/3) + RU(120.0) + sn.E(L/3) + RU(-60.0) + sn.E(L/3) +end +rule = Rule(sn.E, rhs = Kochsnowflake) +#= + +The model is then created by constructing the graph + +=# +Koch = Graph(axiom = axiom, rules = Tuple(rule)) +#= + +In order to be able to generate a 3D structure we need to define a method for +the function `VirtualPlantLab.feed!` (notice the need to prefix it with `VirtualPlantLab.` as we are +going to define a method for this function). The method needs to two take two +arguments, the first one is always an object of type Turtle and the second is an +object of the type for which the method is defined (in this case, E). + +The body of the method should generate the 3D structures using the geometry +primitives provided by VPL and feed them to the turtle that is being passed to +the method as first argument. In this case, we are going to represent the edges +of the Koch snowflakes with cylinders, which can be generated with the +`HollowCylinder!` function from VirtualPlantLab. Note that the `feed!` should return +`nothing`, the turtle will be modified in place (hence the use of `!` at the end +of the function as customary in the VPL community). + +In order to render the geometry, we need assign a `color` (i.e., any type of +color support by the package ColorTypes.jl). In this case, we just feed a basic +`RGB` color defined by the proportion of red, green and blue. To make the +figures more appealing, we can assign random values to each channel of the color +to generate random colors. + +=# +function VirtualPlantLab.feed!(turtle::Turtle, e::sn.E, vars) + HollowCylinder!(turtle, length = e.length, width = e.length/10, + height = e.length/10, move = true, + color = RGB(rand(), rand(), rand())) + return nothing +end +#= + +Note that the argument `move = true` indicates that the turtle should move +forward as the cylinder is generated a distance equal to the length of the +cylinder. Also, the `feed!` method has a third argument called `vars`. This +gives acess to the shared variables stored within the graph (such that they can +be accessed by any node). In this case, we are not using this argument. + +After defining the method, we can now call the function render on the graph to +generate a 3D interactive image of the Koch snowflake in the current state - #= - Note that the argument `move = true` indicates that the turtle should move - forward as the cylinder is generated a distance equal to the length of the - cylinder. Also, the `feed!` method has a third argument called `vars`. This - gives acess to the shared variables stored within the graph (such that they can - be accessed by any node). In this case, we are not using this argument. - - After defining the method, we can now call the function render on the graph to - generate a 3D interactive image of the Koch snowflake in the current state - =# - sc = Scene(Koch) - render(sc, axes = false) - - #This renders the initial triangle of the construction procedure of the Koch - #snowflake. Let's execute the rules once to verify that we get the 2nd iteration - #(check the figure at the beginning of this document): +=# +sc = Scene(Koch) +render(sc, axes = false) +#= + +This renders the initial triangle of the construction procedure of the Koch +snowflake. Let's execute the rules once to verify that we get the 2nd iteration +(check the figure at the beginning of this document): + +=# +rewrite!(Koch) +render(Scene(Koch), axes = false) +#= + +And two more times + +=# +for i in 1:3 rewrite!(Koch) - render(Scene(Koch), axes = false) +end +render(Scene(Koch), axes = false) +#= - # And two more times - for i = 1:3 - rewrite!(Koch) - end - render(Scene(Koch), axes = false) - - #= - # Other snowflake fractals - - To demonstrate the power of this approach, let's create an alternative - snowflake. We will simply invert the rotations of the turtle in the rewriting - rule - =# - function Kochsnowflake2(x) - L = data(x).length - sn.E(L / 3) + - RU(60.0) + - sn.E(L / 3) + - RU(-120.0) + - sn.E(L / 3) + - RU(60.0) + - sn.E(L / 3) - end - rule2 = Rule(sn.E, rhs = Kochsnowflake2) - Koch2 = Graph(axiom = axiom, rules = Tuple(rule2)) - - # The axiom is the same, but now the edges added by the rule will generate the - # edges towards the inside of the initial triangle. Let's execute the first - # three iterations and render the results - - # First iteration - rewrite!(Koch2) - render(Scene(Koch2), axes = false) - # Second iteration - rewrite!(Koch2) - render(Scene(Koch2), axes = false) - # Third iteration - rewrite!(Koch2) - render(Scene(Koch2), axes = false) - - # This is know as [Koch - # antisnowflake](https://mathworld.wolfram.com/KochAntisnowflake.html). We could - # also easily generate a [Cesàro - # fractal](https://mathworld.wolfram.com/CesaroFractal.html) by also changing - # the axiom: - axiomCesaro = sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L) - Cesaro = Graph(axiom = axiomCesaro, rules = (rule2,)) - render(Scene(Cesaro), axes = false) - - # And, as before, let's go through the first three iterations - - # First iteration - rewrite!(Cesaro) - render(Scene(Cesaro), axes = false) - - # Second iteration - rewrite!(Cesaro) - render(Scene(Cesaro), axes = false) - - # Third iteration - rewrite!(Cesaro) - render(Scene(Cesaro), axes = false) +# Other snowflake fractals + +To demonstrate the power of this approach, let's create an alternative +snowflake. We will simply invert the rotations of the turtle in the rewriting +rule +=# +function Kochsnowflake2(x) + L = data(x).length + sn.E(L/3) + RU(60.0) + sn.E(L/3) + RU(-120.0) + sn.E(L/3) + RU(60.0) + sn.E(L/3) end +rule2 = Rule(sn.E, rhs = Kochsnowflake2) +Koch2 = Graph(axiom = axiom, rules = Tuple(rule2)) +#= + +The axiom is the same, but now the edges added by the rule will generate the +edges towards the inside of the initial triangle. Let's execute the first three +iterations and render the results + +=# +## First iteration +rewrite!(Koch2) +render(Scene(Koch2), axes = false) +## Second iteration +rewrite!(Koch2) +render(Scene(Koch2), axes = false) +## Third iteration +rewrite!(Koch2) +render(Scene(Koch2), axes = false) +#= + +This is know as [Koch +antisnowflake](https://mathworld.wolfram.com/KochAntisnowflake.html). We could +also easily generate a [Cesàro +fractal](https://mathworld.wolfram.com/CesaroFractal.html) by also changing the +axiom: + +=# +axiomCesaro = sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L) +Cesaro = Graph(axiom = axiomCesaro, rules = (rule2,)) +render(Scene(Cesaro), axes = false) +#= + +And, as before, let's go through the first three iterations + +=# +## First iteration +rewrite!(Cesaro) +render(Scene(Cesaro), axes = false) +## Second iteration +rewrite!(Cesaro) +render(Scene(Cesaro), axes = false) +## Third iteration +rewrite!(Cesaro) +render(Scene(Cesaro), axes = false) diff --git a/test/tree.jl b/test/tree.jl index 427b47b..b3cd0bd 100644 --- a/test/tree.jl +++ b/test/tree.jl @@ -1,276 +1,196 @@ #= -Tree Alejandro Morales Sierra Centre for Crop Systems Analysis - Wageningen -University +# Tree -In this example we build a 3D representation of a binary TreeTypes. Although -this will not look like a real plant, this example will help introduce -additional features of VirtualPlantLab. +Alejandro Morales + +Centre for Crop Systems Analysis - Wageningen University + + +In this example we build a 3D representation of a binary TreeTypes. Although this will not look like a real plant, this example will help introduce additional features of VPL. The model requires five types of nodes: -*Meristem*: These are the nodes responsible for growth of new organs in our -binary TreeTypes. They contain no data or geometry (i.e. they are a point in the -3D structure). - -*Internode*: The result of growth of a branch, between two nodes. Internodes are -represented by cylinders with a fixed width but variable length. - -*Node*: What is left after a meristem produces a new organ (it separates -internodes). They contain no data or geometry (so also a point) but are required -to keep the branching structure of the tree as well as connecting leaves. - -*Bud*: These are dormant meristems associated to tree nodes. When they are -activated, they become an active meristem that produces a branch. They contain -no data or geometry but they change the orientation of the turtle. - -*BudNode*: The node left by a bud after it has been activated. They contain no -data or geometry but they change the orientation of the turtle. - -*Leaf*: These are the nodes associated to leaves in the TreeTypes. They are -represented by ellipses with a particular orientation and insertion angle. The -insertion angle is assumed constant, but the orientation angle varies according -to an elliptical phyllotaxis rule. - -In this first simple model, only internodes grow over time according to a -relative growth rate, whereas leaves are assumed to be of fixed sized determined -at their creation. For simplicity, all active meristems will produce an phytomer -(combination of node, internode, leaves and buds) per time step. Bud break is -assumed stochastic, with a probability that increases proportional to the number -of phytomers from the apical meristem (up to 1). In the following tutorials, -these assumptions are replaced by more realistic models of light interception, -photosynthesis, etc. - -In order to simulate growth of the 3D binary tree, we need to define a parameter -describing the relative rate at which each internode elongates in each iteration -of the simulation, a coefficient to compute the probability of bud break as well -as the insertion and orientation angles of the leaves. We could stored these -values as global constants, but VirtualPlantLab offers to opportunity to store them per -plant. This makes it easier to manage multiple plants in the same simulation -that may belong to different species, cultivars, ecotypes or simply to simulate -plant-to-plant variation. Graphs in VirtualPlantLab can store an object of any user-defined -type that will me made accessible to graph rewriting rules and queries. For this -example, we define a data type `treeparams` that holds the relevant parameters. -We use `Base.@kwdef` to assign default values to all parameters and allow to -assign them by name. +*Meristem*: These are the nodes responsible for growth of new organs in our binary TreeTypes. They contain no data or geometry (i.e. they are a point in the 3D structure). + +*Internode*: The result of growth of a branch, between two nodes. Internodes are represented by cylinders with a fixed width but variable length. + +*Node*: What is left after a meristem produces a new organ (it separates internodes). They contain no data or geometry (so also a point) but are required to keep the branching structure of the tree as well as connecting leaves. + +*Bud*: These are dormant meristems associated to tree nodes. When they are activated, they become an active meristem that produces a branch. They contain no data or geometry but they change the orientation of the turtle. + +*BudNode*: The node left by a bud after it has been activated. They contain no data or geometry but they change the orientation of the turtle. + +*Leaf*: These are the nodes associated to leaves in the TreeTypes. They are represented by ellipses with a particular orientation and insertion angle. The insertion angle is assumed constant, but the orientation angle varies according to an elliptical phyllotaxis rule. + +In this first simple model, only internodes grow over time according to a relative growth rate, whereas leaves are assumed to be of fixed sized determined at their creation. For simplicity, all active meristems will produce an phytomer (combination of node, internode, leaves and buds) per time step. Bud break is assumed stochastic, with a probability that increases proportional to the number of phytomers from the apical meristem (up to 1). In the following tutorials, these assumptions are replaced by more realistic models of light interception, photosynthesis, etc. + +In order to simulate growth of the 3D binary tree, we need to define a parameter describing the relative rate at which each internode elongates in each iteration of the simulation, a coefficient to compute the probability of bud break as well as the insertion and orientation angles of the leaves. We could stored these values as global constants, but VPL offers to opportunity to store them per plant. This makes it easier to manage multiple plants in the same simulation that may belong to different species, cultivars, ecotypes or simply to simulate plant-to-plant variation. Graphs in VPL can store an object of any user-defined type that will me made accessible to graph rewriting rules and queries. For this example, we define a data type `treeparams` that holds the relevant parameters. We use `Base.@kwdef` to assign default values to all parameters and allow to assign them by name. + =# using VirtualPlantLab using ColorTypes import GLMakie module TreeTypes -import VirtualPlantLab -# Meristem -struct Meristem <: VirtualPlantLab.Node end -# Bud -struct Bud <: VirtualPlantLab.Node end -# Node -struct Node <: VirtualPlantLab.Node end -# BudNode -struct BudNode <: VirtualPlantLab.Node end -# Internode (needs to be mutable to allow for changes over time) -Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node - length::Float64 = 0.10 # Internodes start at 10 cm + import VirtualPlantLab + ## Meristem + struct Meristem <: VirtualPlantLab.Node end + ## Bud + struct Bud <: VirtualPlantLab.Node end + ## Node + struct Node <: VirtualPlantLab.Node end + ## BudNode + struct BudNode <: VirtualPlantLab.Node end + ## Internode (needs to be mutable to allow for changes over time) + Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node + length::Float64 = 0.10 ## Internodes start at 10 cm + end + ## Leaf + Base.@kwdef struct Leaf <: VirtualPlantLab.Node + length::Float64 = 0.20 ## Leaves are 20 cm long + width::Float64 = 0.1 ## Leaves are 10 cm wide + end + ## Graph-level variables + Base.@kwdef struct treeparams + growth::Float64 = 0.1 + budbreak::Float64 = 0.25 + phyllotaxis::Float64 = 140.0 + leaf_angle::Float64 = 30.0 + branch_angle::Float64 = 45.0 + end end -# Leaf -Base.@kwdef struct Leaf <: VirtualPlantLab.Node - length::Float64 = 0.20 # Leaves are 20 cm long - width::Float64 = 0.1 # Leaves are 10 cm wide +#= + +As always, the 3D structure and the color of each type of node are implemented with the `feed!` method. In this case, the internodes and leaves have a 3D representation, whereas bud nodes rotate the turtle. The rest of the elements of the trees are just points in the 3D structure, and hence do not have an explicit geometry: + +=# +## Create geometry + color for the internodes +function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, vars) + ## Rotate turtle around the head to implement elliptical phyllotaxis + rh!(turtle, vars.phyllotaxis) + HollowCylinder!(turtle, length = i.length, height = i.length/15, width = i.length/15, + move = true, color = RGB(0.5,0.4,0.0)) + return nothing end -# Graph-level variables -Base.@kwdef struct treeparams - growth::Float64 = 0.1 - budbreak::Float64 = 0.25 - phyllotaxis::Float64 = 140.0 - leaf_angle::Float64 = 30.0 - branch_angle::Float64 = 45.0 + +## Create geometry + color for the leaves +function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -vars.leaf_angle) + ## Generate the leaf + Ellipse!(turtle, length = l.length, width = l.width, move = false, + color = RGB(0.2,0.6,0.2)) + ## Rotate turtle back to original direction + ra!(turtle, vars.leaf_angle) + return nothing end + +## Insertion angle for the bud nodes +function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars) + ## Rotate turtle around the arm for insertion angle + ra!(turtle, -vars.branch_angle) end +#= -let - #= - As always, the 3D structure and the color of each type of node are implemented - with the `feed!` method. In this case, the internodes and leaves have a 3D - representation, whereas bud nodes rotate the turtle. The rest of the elements of - the trees are just points in the 3D structure, and hence do not have an explicit - geometry: - =# - # Create geometry + color for the internodes - function VirtualPlantLab.feed!(turtle::Turtle, i::TreeTypes.Internode, vars) - # Rotate turtle around the head to implement elliptical phyllotaxis - rh!(turtle, vars.phyllotaxis) - HollowCylinder!( - turtle, - length = i.length, - height = i.length / 15, - width = i.length / 15, - move = true, - color = RGB(0.5, 0.4, 0.0), - ) - return nothing - end +The growth rule for a branch within a tree is simple: a phytomer (or basic unit of morphology) is composed of a node, a leaf, a bud node, an internode and an active meristem at the end. Each time step, the meristem is replaced by a new phytomer, allowing for developmemnt within a branch. - # Create geometry + color for the leaves - function VirtualPlantLab.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -vars.leaf_angle) - # Generate the leaf - Ellipse!( - turtle, - length = l.length, - width = l.width, - move = false, - color = RGB(0.2, 0.6, 0.2), - ) - # Rotate turtle back to original direction - ra!(turtle, vars.leaf_angle) - return nothing - end +=# +meristem_rule = Rule(TreeTypes.Meristem, rhs = mer -> TreeTypes.Node() + + (TreeTypes.Bud(), TreeTypes.Leaf()) + + TreeTypes.Internode() + TreeTypes.Meristem()) +#= - # Insertion angle for the bud nodes - function VirtualPlantLab.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars) - # Rotate turtle around the arm for insertion angle - ra!(turtle, -vars.branch_angle) - end +In addition, every step of the simulation, each bud may break, creating a new branch. The probability of bud break is proportional to the number of phytomers from the apical meristem (up to 1), which requires a relational rule to count the number of internodes in the graph up to reaching a meristem. When a bud breaks, it is replaced by a bud node, an internode and a new meristem. This new meristem becomes the apical meristem of the new branch, such that `meristem_rule` would apply. Note how we create an external function to compute whether a bud breaks or not. This is useful to keep the `branch_rule` rule simple and readable, while allow for a relatively complex bud break model. It also makes it easier to debug the bud break model, since it can be tested independently of the graph rewriting. - #= - The growth rule for a branch within a tree is simple: a phytomer (or basic unit - of morphology) is composed of a node, a leaf, a bud node, an internode and an - active meristem at the end. Each time step, the meristem is replaced by a new - phytomer, allowing for developmemnt within a branch. - =# - meristem_rule = Rule( - TreeTypes.Meristem, - rhs = mer -> - TreeTypes.Node() + - (TreeTypes.Bud(), TreeTypes.Leaf()) + - TreeTypes.Internode() + - TreeTypes.Meristem(), - ) - - #= - In addition, every step of the simulation, each bud may break, creating a new - branch. The probability of bud break is proportional to the number of phytomers - from the apical meristem (up to 1), which requires a relational rule to count - the number of internodes in the graph up to reaching a meristem. When a bud - breaks, it is replaced by a bud node, an internode and a new meristem. This new - meristem becomes the apical meristem of the new branch, such that - `meristem_rule` would apply. Note how we create an external function to compute - whether a bud breaks or not. This is useful to keep the `branch_rule` rule - simple and readable, while allow for a relatively complex bud break model. It - also makes it easier to debug the bud break model, since it can be tested - independently of the graph rewriting. - =# - function prob_break(bud) - # We move to parent node in the branch where the bud was created - node = parent(bud) - # We count the number of internodes between node and the first Meristem - # moving down the graph - check, steps = has_descendant(node, condition = n -> data(n) isa TreeTypes.Meristem) - steps = Int(ceil(steps / 2)) # Because it will count both the nodes and the internodes - # Compute probability of bud break and determine whether it happens - if check - prob = min(1.0, steps * graph_data(bud).budbreak) - return rand() < prob - # If there is no meristem, an error happened since the model does not allow - # for this - else - error("No meristem found in branch") - end - end - branch_rule = Rule( - TreeTypes.Bud, - lhs = prob_break, - rhs = bud -> TreeTypes.BudNode() + TreeTypes.Internode() + TreeTypes.Meristem(), - ) - - #= - A binary tree initializes as an internode followed by a meristem, so the axiom - can be constructed simply as: - =# - axiom = TreeTypes.Internode() + TreeTypes.Meristem() - - #= - And the object for the tree can be constructed as in previous examples, by - passing the axiom and the graph rewriting rules, but in this case also with the - object with growth-related parameters. - =# - tree = Graph( - axiom = axiom, - rules = (meristem_rule, branch_rule), - data = TreeTypes.treeparams(), - ) - - #= - Note that so far we have not included any code to simulate growth of the - internodes. The reason is that, as elongation of internotes does not change the - topology of the graph (it simply changes the data stored in certain nodes), this - process does not need to be implemented with graph rewriting rules. Instead, we - will use a combination of a query (to identify which nodes need to be altered) - and direct modification of these nodes: - =# - getInternode = Query(TreeTypes.Internode) - - #= - If we apply the query to a graph using the `apply` function, we will get an - array of all the nodes that match the query, allow for direct manipulation of - their contents. To help organize the code, we will create a function that - simulates growth by multiplying the `length` argument of all internodes in a - tree by the `growth` parameter defined in the above: - =# - function elongate!(tree, query) - for x in apply(tree, query) - x.length = x.length * (1.0 + data(tree).growth) - end +=# +function prob_break(bud) + ## We move to parent node in the branch where the bud was created + node = parent(bud) + ## We count the number of internodes between node and the first Meristem + ## moving down the graph + check, steps = has_descendant(node, condition = n -> data(n) isa TreeTypes.Meristem) + steps = Int(ceil(steps/2)) ## Because it will count both the nodes and the internodes + ## Compute probability of bud break and determine whether it happens + if check + prob = min(1.0, steps*graph_data(bud).budbreak) + return rand() < prob + ## If there is no meristem, an error happened since the model does not allow + ## for this + else + error("No meristem found in branch") end +end +branch_rule = Rule(TreeTypes.Bud, + lhs = prob_break, + rhs = bud -> TreeTypes.BudNode() + TreeTypes.Internode() + TreeTypes.Meristem()) +#= - #= - Note that we use `vars` on the `Graph` object to extract the object that was - stored inside of it. Also, as this function will modify the graph which is - passed as input, we append an `!` to the name (this not a special syntax of the - language, its just a convention in the Julia community). Also, in this case, the - query object is kept separate from the graph. We could have also stored it - inside the graph like we did for the parameter `growth`. We could also have - packaged the graph and the query into another type representing an individual - TreeTypes. This is entirely up to the user and indicates that a model can be - implemented in many differences ways with VirtualPlantLab. - - Simulating the growth a tree is a matter of elongating the internodes and - applying the rules to create new internodes: - =# - function growth!(tree, query) - elongate!(tree, query) - rewrite!(tree) - end - #= - and a simulation for n steps is achieved with a simple loop: - =# - function simulate(tree, query, nsteps) - new_tree = deepcopy(tree) - for i = 1:nsteps - growth!(new_tree, query) - end - return new_tree +A binary tree initializes as an internode followed by a meristem, so the axiom can be constructed simply as: + +=# +axiom = TreeTypes.Internode() + TreeTypes.Meristem() +#= + +And the object for the tree can be constructed as in previous examples, by passing the axiom and the graph rewriting rules, but in this case also with the object with growth-related parameters. + +=# +tree = Graph(axiom = axiom, rules = (meristem_rule, branch_rule), data = TreeTypes.treeparams()) +#= + +Note that so far we have not included any code to simulate growth of the internodes. The reason is that, as elongation of internotes does not change the topology of the graph (it simply changes the data stored in certain nodes), this process does not need to be implemented with graph rewriting rules. Instead, we will use a combination of a query (to identify which nodes need to be altered) and direct modification of these nodes: + +=# +getInternode = Query(TreeTypes.Internode) +#= + +If we apply the query to a graph using the `apply` function, we will get an array of all the nodes that match the query, allow for direct manipulation of their contents. To help organize the code, we will create a function that simulates growth by multiplying the `length` argument of all internodes in a tree by the `growth` parameter defined in the above: + +=# +function elongate!(tree, query) + for x in apply(tree, query) + x.length = x.length*(1.0 + data(tree).growth) end +end +#= - #= - Notice that the `simulate` function creates a copy of the object to avoid - overwriting it. If we run the simulation for a couple of steps - =# - newtree = simulate(tree, getInternode, 2) - - #= - The binary tree after two iterations has two branches, as expected: - =# - render(Scene(newtree)) - - #= - Notice how the lengths of the prisms representing internodes decreases as the - branching order increases, as the internodes are younger (i.e. were generated - fewer generations ago). Further steps will generate a structure that is more - tree-like. - =# - newtree = simulate(newtree, getInternode, 15) - render(Scene(newtree)) +Note that we use `vars` on the `Graph` object to extract the object that was stored inside of it. Also, as this function will modify the graph which is passed as input, we append an `!` to the name (this not a special syntax of the language, its just a convention in the Julia community). Also, in this case, the query object is kept separate from the graph. We could have also stored it inside the graph like we did for the parameter `growth`. We could also have packaged the graph and the query into another type representing an individual TreeTypes. This is entirely up to the user and indicates that a model can be implemented in many differences ways with VPL. +Simulating the growth a tree is a matter of elongating the internodes and applying the rules to create new internodes: + +=# +function growth!(tree, query) + elongate!(tree, query) + rewrite!(tree) end +#= + +and a simulation for n steps is achieved with a simple loop: + +=# +function simulate(tree, query, nsteps) + new_tree = deepcopy(tree) + for i in 1:nsteps + growth!(new_tree, query) + end + return new_tree +end +#= + +Notice that the `simulate` function creates a copy of the object to avoid overwriting it. If we run the simulation for a couple of steps + +=# +newtree = simulate(tree, getInternode, 2) +#= + +The binary tree after two iterations has two branches, as expected: + +=# +render(Scene(newtree)) +#= + +Notice how the lengths of the prisms representing internodes decreases as the branching order increases, as the internodes are younger (i.e. were generated fewer generations ago). Further steps will generate a structure that is more tree-like. + +=# +newtree = simulate(newtree, getInternode, 15) +render(Scene(newtree))