FF16 Rainfall API

[devmitch]

This PR begins to implement the rainfall model in the FF16 environment. Currently only a few necessary getters, setters and spline computation functions are exposed, we can add/remove depending on what will be required for the model.

A basic test is included, and some slight refactoring was done to the Interpolator class.

This may need to wait for #304 to be merged.

[codecov-commenter]

Codecov Report

❗ No coverage uploaded for pull request base (develop@a8f8a7b). Click here to learn what that means.
The diff coverage is n/a.

@@            Coverage Diff             @@
##             develop     #324   +/-   ##
==========================================
  Coverage           ?   79.54%           
==========================================
  Files              ?       95           
  Lines              ?     8600           
  Branches           ?        0           
==========================================
  Hits               ?     6841           
  Misses             ?     1759           
  Partials           ?        0           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a8f8a7b...846c0d8. Read the comment docs.

[aornugent]

Nice one @devmitch! I think the next step is to show that it accumulates correctly - i.e. accessing the spline from environment.compute_rates() and checking that the total input rainfall for the duration of a simulation is equivalent to integrating the area under the spline. The code on L57 is a bit nonsense but is testing the same thing.

We now want to replace soil_infiltration_rate with rainfall_spline and update the FF16w tests. I think FF16_Environment holds a time variable that can be accessed from within compute_rates (as the x value/eval point) but you might like to check the call stack and verify how compute_rates is called.

[dfalster]

Hi @devmitch

great work. I'm really happy to see this coming together. Well done setting it up and extending in a way that keeps all tests passing (and adds new tests)

I agree with Andrew's suggested next steps.

[devmitch]

I've set a default spline of y = 0 to stay in line with the default soil infiltration rate of 0.0 in FF16w_make_environment . I used domain [0, 10] to construct it, I hope the spline isn't volatile when you try extrapolate points outside the domain like x = 1000, might investigate that later.

Now the rainfall spline will be set in either of 3 ways:

• Through the default argument (y=0) of FF16w_make_environment
• Through passing in customs arguments of x, y points to FF16w_make_environment
• Doing either of the above and after construction, setting the spline through rainfall_init(x, y)

[dfalster]

Really great progress @devmitch .

Some minor suggestions included in this review. @aornugent may have more to mention

[dfalster]

You're right that need x & y. I would rather pass in a single data frame than an array, so perhaps change this accordingly?

[devmitch]

So single data frame for FF15w_make_environment() argument which is split into two arrays to pass into rainfall_add_points()?

[dfalster]

probably. We will occasionally want to specify an exact environment for a non-dynamic use

[devmitch]

After discussions with @dfalster and @itowers1, the rainfall spline (which now exists in the Environment class) has been replaced with a set of splines called extrinsic drivers.

Construction of the rainfall spline is now done in either of 3 ways:

# default rainfall spline of y = 1
env <- make_environment("FF16w")

# rainfall spline of y = 5
env <- make_environment("FF16w", rainfall=5)

# custom rainfall spline based on set of control points
env <- make_environment("FF16w")
x <- seq(0, 10, 1)
y <- x^2
env$set_extrinsic_driver("rainfall", x, y) # overwrites previous spline of y = 1

Querying the spline now requires the user to pass in the name of the extrinsic driver they wish to query:

# determine the rainfall value at time 84
res <- env$extrinsic_driver_evaluate("rainfall", 84)
# determine the rainfall values over the range of integer times from 0 to 100
res <- env$extrinsic_driver_evaluate_range("rainfall", seq(0, 100, 1))

The list of potentially active extrinsic_drivers is now also available from any environment:

env <- make_environment("FF16w")
expect_equal(env$get_extrinsic_driver_names(), c("rainfall"))

Two more things that need to be looked at:

1. One test is still failing in test-patch.R:

# context: env is an FF16w environment  
env_size <- env$ode_size  
env_state <- patch$ode_state  
env_rates <- patch$ode_rates  
expect_equal(patch$ode_size, env_size)  
  
# either 0 or numeric(0)  
expect_identical(patch$ode_state, numeric(env_size))  
expect_identical(patch$ode_rates, numeric(env_size)) # FAILS, is actually 1  

I'm not exactly sure what the relationship between ode_rates and env_size is, and I'm not sure how to resolve this. This was previously passing as ode_rates was 0 because initially there was no rainfall / no soil infiltration.

2. The rainfall spline is being constructed in the FF16_Environment constructor, though I'm not sure it should be however since rainfall is an FF16w specific feature. However as compute_rates needs the rainfall spline, FF16_Environment needs to have it as well. Bit confused here... this means all FF16 environments will have a possibly uninitialised rainfall spline.

Thoughts @dfalster @aornugent ?

[aornugent]

Tackling the failing test first!

# context: env is an FF16w environment  
env_size <- env$ode_size  
env_state <- patch$ode_state  
env_rates <- patch$ode_rates  

expect_equal(patch$ode_size, env_size)  

# these two are typos that previously passed because we initialised at 0
-expect_identical(patch$ode_state, numeric(env_size))  
+expect_identical(patch$ode_state, env_state)
-expect_identical(patch$ode_rates, numeric(env_size)) 
+expect_identical(patch$ode_rates, env_rates)

[aornugent]

This is awesome work @devmitch. I really like the approach of storing multiple drivers in a map- it makes it much more general and hopefully we can apply the same pattern to a variety of new model processes.

To your second point:

The rainfall spline is being constructed in the FF16_Environment constructor, though I'm not sure it should be

I wasn't able to find a way to inherit from FF16_Environment and extend it with custom methods, like rainfall, while also inheriting from FF16_Strategy as we do in FF16w_Strategy due to templating. In the end, we decided that it was ok to extend FF16_Environment as long as it maintained backwards compatibility with the FF16_Strategy - that is, unused methods were ok as long as we were still able to reproduce the original results published by Dan.

In this case, the uninitialised driver fails loudly, which I think is a good thing:

Error in FF16_Environment__extrinsic_driver_evaluate(self, driver_name,  : 
  Interpolator not initialised -- cannot evaluate 

My only outstanding question is what happens when a spline has been initialised for a shorter duration than the model simulation? I've done a little tinkering and the spline will continue to interpolate outside the bounds of which it has been initialised, but often with strange behaviour:

env$set_extrinsic_driver("rainfall", x = c(1, 2, 3), y = c(10, 1, 0))
env$extrinsic_driver_evaluate_range('rainfall', c(1, 2, 3, 10, 50, 100, 1000))
> 10   1   0   7  47  97 997

I suggest then that we stick with the philosophy of 'fail loudly' and have a simulation stop if it goes beyond the length of the spline.

This might require modifying interpolator.h to store min_x and max_x, then adding a check in environment::extrinsic_driver_evaluate_*. I think (but should confirm) that we don't want to do this in interpolator::evaluate_* because we sometimes do extrapolate outside of the range of control points in canopy.h.

I'd also suggest running this check in the Patch constructor which has access to parameters.max_patch_lifetime and the environment, which provides an opportunity to add additional context to the error message.

Let me know if you'd like to discuss!

[devmitch]

I wasn't able to find a way to inherit from FF16_Environment and extend it with custom methods, like rainfall, while also inheriting from FF16_Strategy as we do in FF16w_Strategy due to templating. In the end, we decided that it was ok to extend FF16_Environment as long as it maintained backwards compatibility with the FF16_Strategy - that is, unused methods were ok as long as we were still able to reproduce the original results published by Dan.

Ah I see, makes sense. I think I was confusing myself a little as well when I was considering this.

My only outstanding question is what happens when a spline has been initialised for a shorter duration than the model simulation? I've done a little tinkering and the spline will continue to interpolate outside the bounds of which it has been initialised, but often with strange behaviour

Yeah it seems to be smart enough to extrapolate correctly for simple curves like the linear y = k spline constructed for FF16w_make_environment() which is constructed with points (0, k), (1, k), ... , (10, k) yet correctly places (1000, k) on the spline. However for curves like your example with low control points and no real pattern it seems to be volatile.

This might require modifying interpolator.h to store min_x and max_x, then adding a check in environment::extrinsic_driver_evaluate_*. I think (but should confirm) that we don't want to do this in interpolator::evaluate_* because we sometimes do extrapolate outside of the range of control points in canopy.h.

Another solution that comes to me is adding a boolean flag in Interpolator called extrapolate, which is false by default and can be changed with setExtrapolation() for example. The logic for interpolator::eval() would then look like:

double Interpolator::eval(double u) const {
  check_active();
  if (not extrapolate and (u < min() or u > max())) { // min() & max() already implemented
    util::stop("Extrapolation disabled and evaluation point of " + u + " outside of interpolated domain.");
  }
  return tk_spline(u);
}

It then seems ideal to also provide a faster, less safe method for spline evaluation as well (this is how it is done in the STL with maps for example - at() checks the key exists before getting the value, operator[] does not)

double interpolator::operator()(double u) const {
  return tk_spline(u);
}

However, if we were to disable extrapolation for rainfall splines, how do we know what domain to construct the default spline of y = k in FF16w_make_environment with if we don't have access to parameters.max_patch_lifetime? Should max_lifetime also be an argument in FF16w_make_environment? @aornugent

[aornugent]

Another solution that comes to me is adding a boolean flag in Interpolator called extrapolate, which is false by default and can be changed with setExtrapolation() for example.

Sounds very tidy! Could we then enable extrapolation in the default constructor, but reset the flag to FALSE when a custom spline is created using set_extrinsic_driver?

Should max_lifetime also be an argument in FF16w_make_environment?

Not if we can avoid it. It'd be nice to re-use the same environment for a variety of simulations. My suggestion to run the extrapolation check in Patch was only really to make the code clear that it should test for compatibility between environment and parameters at the start of a run, if extrapolate = FALSE then it should fail loudly regardless.

I see the tests are failing!

...
40: test_check("plant")
An irrecoverable exception occurred. R is aborting now ...
Segmentation fault (core dumped)

My intuition is it doesn't like calling things that haven't been initialised on Windows and Ubuntu (strange that it passes on OSX). All tests passed locally on my Ubuntu machine the other day, which is a bit curly. Maybe try make rebuild and push the Rcpp and NAMESPACE files? 🤷‍♂️

[dfalster]

Good work @devmitch . I agree with @aornugent 's suggestions above.

[devmitch]

Could we then enable extrapolation in the default constructor, but reset the flag to FALSE when a custom spline is created using set_extrinsic_driver?

Yeah I originally meant this but accidentally said false (disabled) by default.

Not if we can avoid it. It'd be nice to re-use the same environment for a variety of simulations. My suggestion to run the extrapolation check in Patch was only really to make the code clear that it should test for compatibility between environment and parameters at the start of a run, if extrapolate = FALSE then it should fail loudly regardless.

The following is currently how the default spline of y = k is generated in FF16w_make_environment:

FF16w_make_environment <- function(canopy_light_tol = 1e-4, 
                                   canopy_light_nbase = 17,
                                   canopy_light_max_depth = 16, 
                                   canopy_rescale_usually = TRUE,
                                   soil_number_of_depths = 1,
                                   soil_initial_state = 0.0,
                                   rainfall = 1) {
  # ...
  x <- seq(0, 10, 1)
  y <- rep(rainfall, 11)
  e$set_extrinsic_driver("rainfall", x, y)
}

If we disable extrapolation for the rainfall spline, then evaluating it at any point > 10 will throw an error as extrapolation is disabled and we are attempting to evaluate the spline above the maximum control point. Since the normal patch lifetime is up to 105 years, this obviously won't work well when we try to evaluate these larger numbers. Since I'm guessing patch lifetime is arbitrary, how can we construct the spline with a large enough domain of control points to account for any arbitrary patch lifetime?

Maybe try make rebuild and push the Rcpp and NAMESPACE files?

Git tells me there are no changes to these files after running this command...

[aornugent]

Ah ha- if we're handling construction on the R side:

FF16w_make_environment <- function(canopy_light_tol = 1e-4, 
                                   canopy_light_nbase = 17,
                                   canopy_light_max_depth = 16, 
                                   canopy_rescale_usually = TRUE,
                                   soil_number_of_depths = 1,
                                   soil_initial_state = 0.0,
                                   rainfall = 1) {
  # ...
  x <- seq(0, 10, 1)
  y <- rep(rainfall, 11)
- e$set_extrinsic_driver("rainfall", x, y)
+ e$set_extrinsic_driver("rainfall", x, y, extrapolate = TRUE)
}

Users can then update the environment with more complex rainfall scenarios:

env <- FF16w_make_environment(rainfall = 1)
x <- seq(0, 105, len = 100)
y <- x^2

# extrapolate = FALSE by default
env$set_extrinsic_driver("rainfall, x, y)

Then, if the patch lifetime exceeds the length of the control points, the simulation will fail (hopefully with a useful error). This is particularly useful where the extrinisic driver comes from real data- we want to prevent the case where people start analysing results that go beyond their domain of interest. If they want to forecast (i.e. extrapolate beyond the available data) they can make an informed decision to set extrapolate = TRUE or create synthetic data that represents a rainfall scenario of their choice.

---

Maybe try make rebuild and push the Rcpp and NAMESPACE files?
Git tells me there are no changes to these files after running this command...

Darn, that's more concerning then. I'll try to check out these changes tonight and see what might be happening, please let me know if you find anything that might address the problem. My hunch is still that the uninitialised splines will be causing trouble in FF16 and K93, but my hunches are nearly always wrong so more methodical investigation is needed.

[devmitch]

Darn, that's more concerning then. I'll try to check out these changes tonight and see what might be happening, please let me know if you find anything that might address the problem. My hunch is still that the uninitialised splines will be causing trouble in FF16 and K93, but my hunches are nearly always wrong so more methodical investigation is needed.

I wasn't able to figure out much, google tells me that it might have something to do with the R packages installed, not sure.

Anyway, it seems like Rcpp can't handle default arguments, so I've gone with the following for the default linear spline in FF16w_make_environment:

x <- seq(0, 10, 1)
y <- rep(rainfall, 11)
e$set_extrinsic_driver("rainfall", x, y)
e$extrinsic_driver_extrapolation("rainfall", TRUE);

(unless there is another way you know of to get them working...)