devmitch commented 2 years ago

This PR is an attempt to apply the flexibility of the extrinsic drivers API that was recently applied in the Environment class to the Species class through the birth_rate variable. Now, the birth_rate for each species can be either a constant or a custom function over time defined by a spline. The control points are set and stored in the Parameters object which then passes them down to each Species object through the Patch constructor when the runner starts.

The Parameters object is now constructed as:

p0 <- scm_base_parameters("FF16")
p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16w_hyperpar,FALSE)
p2 <- set_constant_birth_rate(p1, 1, 20) # set 1st species to have birth_rate of 20

To set the birth rate as either constant or interpolated with a spline, two functions are defined in util_model.R:

set_constant_birth_rate <- function(p, i, k) {
  #p <- validate(p)
  p$birth_rate_x[[i]] <- seq(1, 10)
  p$birth_rate_y[[i]] <- rep(k, 10)
  p$is_constant_birth_rate[[i]] = TRUE
  p
}

set_interpolated_birth_rate <- function(p, i, x, y) {
  #p <- validate(p)
  p$birth_rate_x[i] <- x
  p$birth_rate_y[i] <- y
  p$is_constant_birth_rate[i] = FALSE
  p
}

Then, the runner can be used as normal and we can check the offspring_production is as expected:

env <- make_environment("FF16")
ctrl <- scm_base_control()

out <- run_scm(p2, env, ctrl)

expect_equal(out$patch$species[[1]]$extrinsic_driver_evaluate("birth_rate", 0), 20)
expect_equal(out$patch$species[[1]]$extrinsic_driver_evaluate("birth_rate", 10), 20)
expect_equal(out$patch$species[[1]]$extrinsic_driver_evaluate("birth_rate", 100), 20)
expect_equal(out$offspring_production, 16.88946, tolerance=1e-5)

Some things to consider and work on before merging:

[x] Test the results of net_reproduction_ratios and compute_rates as these both depend on the birth rate
[ ] It might be possible to be consistent with the current way of setting birth rates by being smarter with how things are set up in Parameters
[x] Currently no defaults are set with the standard Parameters constructor
[x] Since the exact same extrinsic drivers API is used in Species and Environment, its probably best to abstract these functions into their own classes, but I'm not sure exactly how. One possibility that comes to mind is to have both classes hold a public ExtrinsicDrivers object, and then call for example env.extrinsic_drivers.evaluate("blah", 42)
[ ] Perhaps the control point setting functions should be on the C++ side (through Parameters) rather than the R side

aornugent commented 2 years ago

Nice one @devmitch! I'll take a look at this shortly.

devmitch commented 2 years ago

Currently no defaults are set with the standard Parameters constructor

With the newest commit, birth rates for each species are now set at a default constant of 1, which was the previous behaviour.

Perhaps the control point setting functions should be on the C++ side (through Parameters) rather than the R side

It seems that RcppR6 doesn't allow methods in structs:

Parameters:
  name_cpp: "plant::Parameters<T,E>"
  templates:
    parameters: [T, E]
    concrete:
      - ["FF16": "plant::FF16_Strategy", "FF16_Env": "plant::FF16_Environment"]
      - ["FF16w": "plant::FF16w_Strategy", "FF16_Env": "plant::FF16_Environment"]
      - ["FF16r": "plant::FF16r_Strategy", "FF16_Env": "plant::FF16_Environment"]
      - ["K93": "plant::K93_Strategy", "K93_Env": "plant::K93_Environment"]
  validator:
    name_cpp: validate
  list:
    - patch_area: double
    - n_patches: size_t
    - patch_type: "std::string"
    - max_patch_lifetime: double
    - strategies: "std::vector<T>"
    - birth_rate_x: "std::vector<std::vector<double> >"
    - birth_rate_y: "std::vector<std::vector<double> >"
    - is_interpolated_birth_rate: "std::vector<bool>"
    - is_resident: "std::vector<bool>"
    - strategy_default: "T"
    - cohort_schedule_times_default: "std::vector<double>"
    - cohort_schedule_times: "std::vector<std::vector<double> >"
    - cohort_schedule_ode_times: "std::vector<double>"
  methods:
    set_constant_birth_rate:
      args: [ species_index: int, birth_rate: double ]
      return_type: void
    set_interpolated_birth_rate:
      args: [ species_index: int, x: std::vector<double>, y: std::vector<double> ]
      return_type: void

$ make RcppR6
Rscript -e "library(methods); RcppR6::RcppR6()"
Reading classes from inst/RcppR6_classes.yml
Warning in warn_unknown(ret$name_r, ret$defn, valid) :
  Unknown fields in Parameters: methods

So I'm not sure if this is possible or not. It would be more convenient as you wouldn't need to reassign the parameters object every time you wanted to set the birth rate of one species:

# suppose there are 3 species
p <- set_constant_birth_rate(p, 1, 20)
p <- set_constant_birth_rate(p, 2, 20)
p <- set_constant_birth_rate(p, 3, 20)

compared to

p$set_constant_birth_rate(p, 1, 20)
p$set_constant_birth_rate(p, 2, 20)
p$set_constant_birth_rate(p, 3, 20)

But otherwise, once the design is agreed upon I can update the tests and merge these changes.

dfalster commented 2 years ago

Hi @devmitch ,

Great progress here! It's moving in the right direction but still WIP.

Since the exact same extrinsic drivers API is used in Species and Environment, its probably best to abstract these functions into their own classes, but I'm not sure exactly how.

This seems sensible. A stand-alone class would be easy to test. When you say you don't know how, what do you mean?

Further, I had an idea. The driver class could have a function set_constant(y) and set_variable(t,y) and a boolean flag is_constant. When setting as a constant, you would just store a constant (currently you fit a spline with constant y). Then when it's evaluated:

get_driver(t){
   if(is_constant)
      return(constant)
   else
     return(evaulate_spline(t))
}

It seems that RcppR6 doesn't allow methods in structs Correct. But we could turn this into another type of object, which can have methods. @aornugent is better at this than me. Advice @aornugent ?

I see the driver is currently stored in species. Ideally we'd add it to the species strategy, where it currently resides. But now stored as an instance of driver class than species. Again, this would require changing the type of object the strategy is in RcppR6. But this is something I want to do anyway, to enable access to the functions in strategy directly from R.

So perhaps we should first

[ ] make stand alone driver class
[ ] modify the object types for strategy in RcppR6 to make it a fully functional class?

devmitch commented 2 years ago

When you say you don't know how, what do you mean?

Well, aside from separating the functions into their own class and storing an object of this class in both Species and Environment, another option could be to make both these classes inherit from the drivers class, as C++ allows multiple inheritance. Though I've seen the idea that software design should favour "composition over inheritance" and it probably seems more difficult to go the inheritance route, so I suppose not.

Further, I had an idea. The driver class could have a function set_constant(y) and set_variable(t,y) and a boolean flag is_constant. When setting as a constant, you would just store a constant (currently you fit a spline with constant y). Then when it's evaluated:

I had a similar idea for the extrinsic driver functions in the species class, but it makes even more sense to put this in the separated drivers class. Doing this will also allow us to simplify some of the extrinsic driver code for the FF16w helper functions.

I'll start to work on the standalone drivers class now.

dfalster commented 2 years ago

Great. Just to clarify, my suggestions that we

create a generic driver class
use an instance of the driver class to store climate drivers, i.e. updating the work for rainfall api
use another instance of the driver class to store birth_rates in the parameters object
modify the object types for strategies in Rcpp to enable better interaction and access to internal functions

aornugent commented 2 years ago

Attempting to summarise a conversation with @devmitch:

A general implementation of the drivers class is nearly done :rocket:. Both Species and Environments will have an ExtrinsicDrivers object, parameterised as required by a Strategy.
Static fields for control point inputs will be added to each Strategy as needed. This means that the existing birth_rate list can be removed from Parameters.
On the R side, the input of control point values (or constants) will be moved into expand_parameters by adding an extra birth_rate argument that holds x and y values for each species specified in the trait_matrix input. See code excerpt below.
The control point inputs will be accessed during the initialisation of a Species extrinsic drivers object as defined in *_Strategy::prepare_strategy, keeping the different drivers closely tied to a given model.
The Drivers class will be responsible for validating the control point inputs and decide whether to use a constant or interpolate function.

@dfalster - this means that we don't need to change the Strategy Rcpp object type yet. Ok with you if we punt this into a later PR?

Also curious if we can assume that birth_rate is general to all size structured population models? If so, it could be pushed down into the base Strategy class.

@devmitch - did I get this right? Holla if there's any questions.

util_model.R

-expand_parameters <- function(trait_matrix, p, hyperpar=param_hyperpar(p), mutant=TRUE ) {
+expand_parameters <- function(trait_matrix, p, hyperpar=param_hyperpar(p), mutant=TRUE, birth_rate_list) {
  ...
+ if(nrow(trait_matrix) != length(birth_rate_list)
+   stop("Must provide a birth rate input for each species")
-  extra <- strategy_list(trait_matrix, p, hyperpar)
+  extra <- strategy_list(trait_matrix, p, hyperpar, birth_rate_list)
  n_extra <- length(extra)

  ret <- p <- validate(p) # Ensure times are set up correctly.
  ret$strategies <- c(p$strategies, extra)
  ret$is_resident <- c(p$is_resident, rep(!mutant, n_extra))
-  ret$birth_rate <- c(p$birth_rate, rep(1.0, n_extra))
   ...

 ret
}

-strategy_list <- function(x, parameters, hyperpar=param_hyperpar(parameters)) {
+strategy_list <- function(x, parameters, hyperpar=param_hyperpar(parameters), birth_rate_list) {
  if (!is.matrix(x)) {
    stop("Invalid type x -- expected a matrix")
  }

  strategy <- parameters$strategy_default
  x <- hyperpar(x, strategy)

  trait_names <- colnames(x)

  # we could be smarter and generalise birth_rate_list -> species_drivers and map to strategies 
  # using the same approach as `trait_names`, but let's start with one for now
-  f <- function(xi) {
+  f <- function(xi, br) {
    strategy[trait_names] <- xi
+   strategy$birth_rate_x <- br$x
+   strategy$birth_rate_y < br$y
    strategy
  }

-  lapply(matrix_to_list(x), f)
+  mapply(f, matrix_to_list(x), birth_rate_list)
}

devmitch commented 2 years ago

I've made changes to the strategy_list function such that it is slightly easier and less redundant to set up constant birth_rate splines:

x <- seq(0, 200)
birth_rates <- list(
  species1 = list(x = x, y = 1 + sin(x)), # variable
  species2 = 2 # constant
)
p1 <- expand_parameters(lmas, p0, FF16_hyperpar, FALSE, birth_rates)

I have also updated the test-birth-rate-splines.R tests with actual expect statements, however I'm not sure what to test after the SCM run has finished (probably anything that uses birth rates - net_reproduction_ratios for example). There are still a few other tests not passing that I am still working on.

codecov-commenter commented 2 years ago

Codecov Report

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

@@            Coverage Diff             @@
##             develop     #334   +/-   ##
==========================================
  Coverage           ?   79.65%           
==========================================
  Files              ?       97           
  Lines              ?     8818           
  Branches           ?        0           
==========================================
  Hits               ?     7024           
  Misses             ?     1794           
  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 425004d...ecc4ef2. Read the comment docs.

traitecoevo / plant

Further extrinsic drivers #334

Codecov Report