Making a Model

This tutorial will cover the topic of creating a Cropbox model.

• Making a Model
• • Growing Degree-Day
  • • System
    • • Variables
      • Mixins
    • Configuration
    • Simulation
    • Visualization
  • Lotka-Volterra Equations
  • • System
    • Configuration
    • Visualization
    • Density-Dependent Lotka-Volterra Equations
    • • System
      • Configuration
      • Visualization
      • Calibration
      • Evaluation

Growing Degree-Day

You might have heard the terms like growing degree days (GDD), thermal units, heat units, heat sums, temperature sums, and thermal-time that are used to relate the rate of plant or insect development to temperature. They are all synonymous. The concept of thermal-time or thermal-units derives from the long-standing observation and assumption that timing of development is primarily driven by temperature in plants and the relationship is largely linear. The linear relationship is generally held true over normal growing temperatures that are bracketed by the base temperature (Tb) and optimal temperature (Topt). Many existing crop models and tree growth models use thermal-unit approaches (e.g., GDD) for modeling phenology with some modifications to account for other factors like photoperiod, vernalization, dormancy, and stress. The growing degree days (GDD) is defined as the difference between the average daily air temperature (T) and the base temperature below which the developmental process stops. The bigger the difference in a day, the faster the development takes place up to a certain optimal temperature (Topt). The Cumulative GDD (cGDD) since the growth initiation (e.g., sowing, imbibition for germination) is then calculated by:

\[\begin{align} \mathrm{GDD}(T) &= \max \{ 0, \min \{ T, T_{opt} \} - T_b \} \\ \mathrm{cGDD} &= \sum_i^n \mathrm{GDD}(T_i) \\ \end{align}\]

In this section, we will create a model that simulates GDD and cGDD.

System

Let us start by making a system called GrowingDegreeDay. This can be done using a simple Cropbox macro @system.

@system GrowingDegreeDay

Variables

From the equation, let's identify the variables we need to declare in our system. In the equation for GDD, we have two parameters Topt and Tb. Since they are fixed values, we will declare them as preserve variables, which are variables that remain constant throughout a simulation.

@system GrowingDegreeDay begin
    Tb ~ preserve
    To ~ preserve
end

Tb and To are parameters that we may want to change depending on the simulation. To make this possible, we will assign them the parameter tag, which allows the tagged variables to be altered through a configuration for each simulation. Note that we will not assign values at declaration because we will configure them when we run the simulation.

@system GrowingDegreeDay begin
    Tb ~ preserve(parameter)
    To ~ preserve(parameter)
end

Lastly, we will tag the variables with units. Tagging units is the recommended practice for many reasons, one of which is to catch mismatching units during calculations.

@system GrowingDegreeDay begin
    Tb ~ preserve(parameter, u"°C")
    To ~ preserve(parameter, u"°C")
end

In the GDD equation, T represents the average daily temperature value necessary to calculate the GDD. Likewise, the variable in our system will represent a series of daily average temperatures. The series of temperature values will be driven from an external data source, for which we will create a separate system later on for data extraction. For the GrowingDegreeDay system, we will declare T as a hold variable, which represents a placeholder that will be replaced by a T from another system.

@system GrowingDegreeDay begin
    T ~ hold
    Tb ~ preserve(parameter, u"°C")
    To ~ preserve(parameter, u"°C")
end

We declared all the necessary variables required to calculate GDD. Now it is time to declare GDD as a variable in the system. Because GDD is a variable that we want to evaluate and store in each update, we will declare it as a track variable with T, Tb, and To as its depending variables.

@system GrowingDegreeDay begin
    T ~ hold
    Tb ~ preserve(parameter, u"°C")
    To ~ preserve(parameter, u"°C")

    GDD(T, Tb, To) => begin
        min(T, To) - Tb
    end ~ track(min = 0, u"K")
end

Note that we have tagged the unit for GDD as u"K". This is to avoid incompatibilities that u"°C" has with certain operations.

Now that GDD is declared in the system, we will declare cGDD as an accumulate variable with GDD as its depending variable. Recall that accumulate variables perform the Euler method of integration.

@system GrowingDegreeDay begin
    T ~ hold
    Tb ~ preserve(parameter, u"°C")
    To ~ preserve(parameter, u"°C")

    GDD(T, Tb, To) => begin
        min(T, To) - Tb
    end ~ track(min = 0, u"K")

    cGDD(GDD) ~ accumulate(u"K*d")
end

We have declared all the necessary variables for GrowingDegreeDay.

Mixins

Now let's address the issue of the missing temperature values. We will make a new system that will provide the missing temperature data we need for simulating GrowingDegreeDay. We will call this system Temperature. The purpose of Temperature will be to obtain a time series of daily average temperature values from an external data source.

@system Temperature

For this tutorial, we will be using weather data from Beltsville, Maryland in 2002. The data is available here.

If you are unsure how to read CSV files, check the following documentation. You can also follow the template below (make sure path corresponds to the path to your data file).

using CSV, DataFrames

weather = CSV.read(path, DataFrame)

The data should look something like the following:

first(weather, 3)

Notice that the column names have units in parentheses. The unitfy() function in Cropbox automatically assigns units to values based on names of the columns (if the unit is specified).

weather = unitfy(weather)
first(weather, 3)

In the Temperature system, there is one variable that we will declare before declaring any other variable. We will name this variable calendar.

@system Temperature begin
    calendar(context) ~ ::Calendar
end

The purpose of calendar is to have access to variables inside the Calendar system such as init, last, and date, which represent initial, last, and current date, respectively.

The next variable we will add is a variable storing the weather data as a DataFrame. This variable will be a provide variable named data.

@system Temperature begin
    calendar(context) ~ ::Calendar
    data ~ provide(parameter, index=:date, init=calendar.date)
end

Note that we have tagged the variable with a parameter tag so that we can assign a DataFrame during the configuration. We will set the index of the extracted DataFrame as the "date" column of the data source. The init tag is used to specify the starting row of the data that we want to store. calendar.date refers to the date variable in the Calendar system, and is a track variable that keeps track of the dates of simulation. The initial value of date is dependent on calendar.init which we will assign during configuration. By setting init to calendar.date, we are making sure that the provide variable extracts data from the correct starting row corresponding to the desired initial date of simulation.

Now we can finally declare the temperature variable using one of the columns of the DataFrame represented by data. Because this variable is driven from a source, we will be declaring a drive variable named T. The from tag specifies the DataFrame source and the by tag specifies which column to take the values from.

@system Temperature begin
    calendar(context) ~ ::Calendar
    data ~ provide(parameter, index=:date, init=calendar.date)
    T ~ drive(from=data, by=:Tavg, u"°C")
end

Main.Temperature

We finally have all the components to define our model. Because GrowingDegreeDay requires values for T from Temperature, let's redeclare GrowingDegreeDay with Temperature as a mixin. Because we want to run a simulation of GrowingDegreeDay, we also want to include Controller as a mixin. Recall that Controller must be included as a mixin for any system that you want to simulate.

@system GrowingDegreeDay(Temperature, Controller) begin
    T ~ hold
    Tb ~ preserve(parameter, u"°C")
    To ~ preserve(parameter, u"°C")

    GDD(T, Tb, To) => begin
        min(T, To) - Tb
    end ~ track(min = 0, u"K")

    cGDD(GDD) ~ accumulate(u"K*d")
end

Main.GrowingDegreeDay

Configuration

The next step is to create a configuration object to assign the values of parameters. Recall that data, T, Tb, and To are empty variables at the moment.

As covered in the Configuration section, we can make a single Config object with all the configurations we need for our systems.

Given the nature of GDD, this model is a daily model. To run a daily simulation, we need to configure the step variable in the Clock system from 1u"hr" to 1u"d". This will change the time interval of the simulation from hourly (default) to daily.

c = @config :Clock => :step => 1u"d"

Next we will add the configurations for GrowingDegreeDay. The only parameters we have to configure are Tb and To.

c = @config (
    :Clock => (
        :step => 1u"d"
    ),
    :GrowingDegreeDay => (
        :Tb => 8.0u"°C",
        :To => 32.0u"°C"
    )
)

Next we will pair the aforementioned DataFrame weather to data in Temperature

c = @config(
    :Clock => (
        :step => 1u"d"
    ),
    :GrowingDegreeDay => (
        :Tb => 8.0u"°C",
        :To => 32.0u"°C"
    ),
    :Temperature => (
        :data => weather
    )
)

Lastly, we will configure the init and last parameters of the Calendar system, which will define the time range of our simulation.

c = @config(
    :Clock => (
        :step => 1u"d"
    ),
    :GrowingDegreeDay  => (
        :Tb => 8.0u"°C",
        :To => 32.0u"°C"
    ),
    :Temperature => (
        :data => weather
    ),
    :Calendar => (
        :init => ZonedDateTime(2002, 5, 15, tz"America/New_York"),
        :last => ZonedDateTime(2002, 9, 30, tz"America/New_York")
    )
)

Config for 4 systems:

Simulation

Now that we have fully defined GrowingDegreeDay and created a configuration for it, we can finally simulate the model.

s = simulate(GrowingDegreeDay;
    config = c,
    stop = "calendar.stop",
    index = "calendar.date",
    target = [:GDD, :cGDD]
)

first(s, 10)

Visualization

To end the tutorial, let's visualize the simulation using the plot() and visualize() functions.

We can input the DataFrame from our simulation in the plot() function to create a plot.

Here is a plot of GDD over time.

plot(s, "calendar.date", :GDD; kind=:line)

We can also simultaneously run a new simulation and plot its result using the visualize() function.

Here is a plot of cGDD over time.

visualize(GrowingDegreeDay, "calendar.date", :cGDD; config=c, stop="calendar.stop", kind=:line)

Lotka-Volterra Equations

In this tutorial, we will create a model that simulates population dynamics between prey and predator using the Lotka-Volterra equations. The Lotka-Volterra equations are as follows:

\[\begin{align} \frac{dN}{dt} &= bN - aNP \\ \frac{dP}{dt} &= caNP - mP \\ \end{align}\]

Here is a list of variables used in the system:

System

Let's begin by creating a system called LotkaVolterra. Since this is a system that we want to simulate later on, we must include Controller as a mixin.

@system LotkaVolterra(Controller)

We will first declare a time variable with a yearly unit, which we will use for plotting the model simulations later on. Recall that context.clock.time is a variable that keeps track of the progression of time. We are simply declaring a variable to keep track of the time in years.

@system LotkaVolterra(Controller) begin
    t(context.clock.time) ~ track(u"yr")
end

Next, we will declare the parameters in the equations as preserve variables. preserve variables are variables that remain constant throughout a simulation.

@system LotkaVolterra(Controller) begin
    t(context.clock.time) ~ track(u"yr")

    b: prey_birth_rate            ~ preserve(parameter, u"yr^-1")
    a: predation_rate             ~ preserve(parameter, u"yr^-1")
    c: predator_reproduction_rate ~ preserve(parameter)
    m: predator_mortality_rate    ~ preserve(parameter, u"yr^-1")
end

Now let's declare the prey and predator populations as variables. The Lotka-Volterra equations describe the rates of change for the two populations. As we want to track the actual number of the two populations, we will declare the two populations as accumulate variables, which are simply Euler integrations of the two population rates. Note that a variable can be used as its own depending variable.

@system LotkaVolterra(Controller) begin
    t(context.clock.time) ~ track(u"yr")

    b: prey_birth_rate            ~ preserve(parameter, u"yr^-1")
    a: predation_rate             ~ preserve(parameter, u"yr^-1")
    c: predator_reproduction_rate ~ preserve(parameter)
    m: predator_mortality_rate    ~ preserve(parameter, u"yr^-1")

    N(N, P, b, a):    prey_population     =>     b*N - a*N*P ~ accumulate
    P(N, P, c, a, m): predator_population => c*a*N*P -   m*P ~ accumulate
end

By default, accumulate variables initialize at a value of zero. In our current model, that would result in two populations remaining at zero indefinitely. To address this, we will define the initial values for the two accumulate variables using the init tag. We can specify a particular value, or we can also create and reference new parameters representing the two initial populations. We will go with the latter option as it allows us to flexibly change the initial populations with a configuration.

@system LotkaVolterra(Controller) begin
    t(context.clock.time) ~ track(u"yr")

    b: prey_birth_rate            ~ preserve(parameter, u"yr^-1")
    a: predation_rate             ~ preserve(parameter, u"yr^-1")
    c: predator_reproduction_rate ~ preserve(parameter)
    m: predator_mortality_rate    ~ preserve(parameter, u"yr^-1")

    N0: prey_initial_population     ~ preserve(parameter)
    P0: predator_initial_population ~ preserve(parameter)

    N(N, P, b, a):    prey_population     =>     b*N - a*N*P ~ accumulate(init=N0)
    P(N, P, c, a, m): predator_population => c*a*N*P -   m*P ~ accumulate(init=P0)
end

Main.LotkaVolterra

Configuration

With the system now defined, we will create a Config object to fill or adjust the parameters.

First, we will change the step variable in the Clock system to 1u"d", which will make the system update at a daily interval. Recall that Clock is a system that is referenced in all systems by default. You can technically run the model with any timestep.

lvc = @config (:Clock => :step => 1u"d")

Config for 1 system:

Next, we will configure the parameters in the LotkaVolterra system that we defined. Note that we can easily combine configurations by providing multiple elements.

lvc = @config (lvc,
    :LotkaVolterra => (
        b = 0.6,
        a = 0.02,
        c = 0.5,
        m = 0.5,
        N0 = 20,
        P0 = 30
    )
)

Config for 2 systems:

Visualization

Let's visualize the LotkaVolterra system with the configuration that we just created, using the visualize() function. The visualize() function both runs a simulation and plots the resulting DataFrame.

visualize(LotkaVolterra, :t, [:N, :P]; config = lvc, stop = 100u"yr", kind = :line)

Density-Dependent Lotka-Volterra Equations

Now let's try to make a density-dependent version of the original Lotka-Volterra model which incorporates a new term in the prey population rate. The new variable K represents the carrying capacity of the prey population.

\[\begin{align} \frac{dN}{dt} &= bN-\frac{b}{K}N^2-aNP \\ \frac{dP}{dt} &= caNP-mP \\ \end{align}\]

System

We will call this new system LotkaVolterraDD.

@system LotkaVolterraDD(Controller)

Since we already defined the LotkaVolterra system, which already has most of the variables we require, we can use LotkaVolterra as a mixin for LotkaVolterraDD. This makes our task a lot simpler, as all that remains is to declare the variable K for carrying capacity and redeclare the variable N for prey population. The variable N in the new system will automatically overwrite the N from LotkaVolterra.

@system LotkaVolterraDD(LotkaVolterra, Controller) begin
    N(N, P, K, b, a): prey_population => begin
        b*N - b/K*N^2 - a*N*P
    end ~ accumulate(init = N0)

    K: carrying_capacity ~ preserve(parameter)
end

Main.LotkaVolterraDD

Configuration

Much like the new system, the new configuration can be created by reusing the old configuration. All we need to configure is the new variable K.

lvddc = @config(lvc, (:LotkaVolterraDD => :K => 1000))

Config for 3 systems:

Visualization

Once again, let's visualize the system using the visualize() function.

visualize(LotkaVolterraDD, :t, [:N, :P]; config = lvddc, stop = 100u"yr", kind = :line)

Calibration

If you want to calibrate the parameters according to a particular dataset, Cropbox provides the calibrate() function, which relies on BlackBoxOptim.jl for global optimization methods. If you are interested in local optimization methods, refer to Optim.jl package for more information.

For this tutorial, we will use a dataset containing the number of pelts (in thousands) of Canadian lynx and snowshoe hare traded by the Hudson Bay Trading Company in Canada from 1845 to 1935. The data is available here.

If you are unsure how to read CSV files, check the following documentation. You can also follow the template below (make sure path corresponds to the path to your data file).

using CSV, DataFrames

pelts = CSV.read(path, DataFrame)

The data should look something like the this:

first(pelts, 3)

Recall that we can use the unitfy() function in Cropbox to automatically assign units when they are specified in the column headers.

pelts = unitfy(pelts)
first(pelts, 3)

Let's plot the data and see what it looks like.

visualize(pelts, :Year, [:Hare, :Lynx], kind = :scatterline)

For our calibration, we will use a subset of the data covering years 1900 to 1920.

pelts_subset = @subset(pelts, 1900u"yr" .<= :Year .<= 1920u"yr")

Before we calibrate the parameters for LotkaVolterra, let's add one new variable to the system. We will name this variable y for year. The purpose of y is to keep track of the year in the same manner as the dataset.

@system LotkaVolterra(Controller) begin
    t(context.clock.time) ~ track(u"yr")
    y(t): year            ~ track::Int(u"yr", round)

    b: prey_birth_rate            ~ preserve(parameter, u"yr^-1")
    a: predation_rate             ~ preserve(parameter, u"yr^-1")
    c: predator_reproduction_rate ~ preserve(parameter)
    m: predator_mortality_rate    ~ preserve(parameter, u"yr^-1")

    N0: prey_initial_population     ~ preserve(parameter)
    P0: predator_initial_population ~ preserve(parameter)

    N(N, P, b, a):    prey_population     =>     b*N - a*N*P ~ accumulate(init=N0)
    P(N, P, c, a, m): predator_population => c*a*N*P -   m*P ~ accumulate(init=P0)
end

Main.LotkaVolterra

We will now use the calibrate() function to find parameters that fit the data. Keep in mind that the search range for each parameter will be determined by you. We will use the snap option to explicitly indicate that the output should be recorded by 365-day intervals to avoid excessive rows in the DataFrame causing unnecessary slowdown. Note that we will use 365u"d" instead of 1u"yr" which is technically equivalent to 365.25u"d" following the convention in astronomy. For information regarding syntax, please check the reference.

lvcc = calibrate(LotkaVolterra, pelts_subset;
    index = :Year => :y,
    target = [:Hare => :N, :Lynx => :P],
    config = :Clock => (:init => 1900u"yr", :step => 1u"d"),
    parameters = LotkaVolterra => (;
        b = (0, 2),
        a = (0, 2),
        c = (0, 2),
        m = (0, 2),
        N0 = (0, 100),
        P0 = (0, 100),
    ),
    stop = 20u"yr",
    snap = 365u"d"
)

Config for 1 system:

As you can see above, the calibrate() function will return a Config object for the system.

Using the new configuration, let's make a comparison plot to visualize how well the simualation with the new parameters fits the data.

p1 = visualize(pelts_subset, :Year, [:Hare, :Lynx]; kind = :scatterline)
visualize!(p1, LotkaVolterra, :t, [:N, :P];
    config = (lvcc, :Clock => (:init => 1900u"yr", :step => 1u"d")),
    stop = 20u"yr",
    kind = :line,
    colors = [1, 2],
    names = [],
)

Now let's try calibrating the density-dependent version of the model. Since we made a slight change to LotkaVolterra, let's make sure to define LotkaVolterraDD again.

@system LotkaVolterraDD(LotkaVolterra, Controller) begin
    N(N, P, K, b, a): prey_population => begin
        b*N - b/K*N^2 - a*N*P
    end ~ accumulate(init = N0)

    K: carrying_capacity ~ preserve(parameter)
end

Main.LotkaVolterraDD

Don't forget to add K among the parameters that we want to calibrate.

lvddcc = calibrate(LotkaVolterraDD, pelts_subset;
    index = :Year => :y,
    target = [:Hare => :N, :Lynx => :P],
    config = :Clock => (:init => 1900u"yr", :step => 1u"d"),
    parameters = LotkaVolterraDD => (;
        b = (0, 2),
        a = (0, 2),
        c = (0, 2),
        m = (0, 2),
        N0 = (0, 100),
        P0 = (0, 100),
        K = (0, 1000)
    ),
    stop = 20u"yr",
    snap = 365u"d"
)

Config for 1 system:

Once again, let us make a comparison plot to see how the density-dependent version of the model fares against the original dataset.

p2 = visualize(pelts_subset, :Year, [:Hare, :Lynx]; kind = :scatterline)
visualize!(p2, LotkaVolterraDD, :t, [:N, :P];
    config = (lvddcc, :Clock => (:init => 1900u"yr", :step => 1u"d")),
    stop = 20u"yr",
    kind = :line,
    colors = [1, 2],
    names = [],
)

Evaluation

We have visualized how the simulated LotkaVolterra and LotkaVolterraDD systems compare to the the original dataset. Let us obtain a metric for how well the simulations fit the original dataset using the evaluate() function in Cropbox. The evaluate() function supports numerous different metrics for evaluation. Here, we will calculate the root-mean-square error (RMSE) and modeling efficiency (EF).

Here are the evaluation metrics for LotkaVolterra. The numbers in the tuples correspond to hare and lynx, respectively.

evaluate(LotkaVolterra, pelts_subset;
    index = :Year => :y,
    target = [:Hare => :N, :Lynx => :P],
    config = (lvcc, :Clock => (:init => 1900u"yr", :step => 1u"d")),
    stop = 20u"yr",
    snap = 365u"d",
    metric = :rmse,
)

(17.37148524278938, 9.760454951690196)

evaluate(LotkaVolterra, pelts_subset;
    index = :Year => :y,
    target = [:Hare => :N, :Lynx => :P],
    config = (lvcc, :Clock => (:init => 1900u"yr", :step => 1u"d")),
    stop = 20u"yr",
    snap = 365u"d",
    metric = :ef
)

(0.45103924288219865, 0.7479251686137505)

Here are the evaluation metrics for LotkaVolterraDD:

evaluate(LotkaVolterraDD, pelts_subset;
    index = :Year => :y,
    target = [:Hare => :N, :Lynx => :P],
    config = (lvddcc, :Clock => (:init => 1900u"yr", :step => 1u"d")),
    stop = 20u"yr",
    snap = 365u"d",
    metric = :rmse
)

(16.17943364709707, 11.70354851591861)

evaluate(LotkaVolterraDD, pelts_subset;
    index = :Year => :y,
    target = [:Hare => :N, :Lynx => :P],
    config = (lvddcc, :Clock => (:init => 1900u"yr", :step => 1u"d")),
    stop = 20u"yr",
    snap = 365u"d",
    metric = :ef
)

(0.5237949169511267, 0.6375697135890357)