FF16 Rainfall API (#324)

* Implemented functions to explose auxillary variables in R API * Resized auxillary vector to appropriate size * Exposed auxillary variables in R API + basic test * Fixed all typos of 'auxillary' to 'auxiliary' * Spline added to FF16 environment with basic API exposed + tests * basic test for rainfall spline over run_scm * Default rainfall spline of y = 0 * Removed depricated infiltration rate settings, updated documentation * some changes to testing with regards to tolerance * fixed incorrect tolerance * fixed comp function for std::is_sorted + more specific testsasdasd * less lenient tolerance to pass windows tests * Refactored FF16w rainfall splines in favour of more generic 'extrinsic drivers' * Ability to get extrinsic driver names in R + small refactoring * fixed old typo in patch tests * extrapolation setting in interpolator * trying with .at() rather than operator[] * Remove rainfall_spline; update Rcpp * removed stringstream in favour of classic std::string for error handling * more useful error messages for extrapolation * adding return statement in evaluate * Removed comments + minor refactor Co-authored-by: aornugent <[email protected]> Co-authored-by: Daniel Falster <[email protected]>
traitecoevo · Dec 13, 2021 · 930b9c9 · 930b9c9
1 parent a8f8a7b
commit 930b9c9
Show file tree

Hide file tree

Showing 14 changed files with 299 additions and 72 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -3293,6 +3293,26 @@ FF16_Environment__ctor <- function(canopy_rescale_usually, soil_number_of_depths
     .Call('_plant_FF16_Environment__ctor', PACKAGE = 'plant', canopy_rescale_usually, soil_number_of_depths)
 }
 
+FF16_Environment__set_extrinsic_driver <- function(obj_, driver_name, x, y) {
+    invisible(.Call('_plant_FF16_Environment__set_extrinsic_driver', PACKAGE = 'plant', obj_, driver_name, x, y))
+}
+
+FF16_Environment__extrinsic_driver_extrapolation <- function(obj_, driver_name, extrapolate) {
+    invisible(.Call('_plant_FF16_Environment__extrinsic_driver_extrapolation', PACKAGE = 'plant', obj_, driver_name, extrapolate))
+}
+
+FF16_Environment__extrinsic_driver_evaluate <- function(obj_, driver_name, u) {
+    .Call('_plant_FF16_Environment__extrinsic_driver_evaluate', PACKAGE = 'plant', obj_, driver_name, u)
+}
+
+FF16_Environment__extrinsic_driver_evaluate_range <- function(obj_, driver_name, u) {
+    .Call('_plant_FF16_Environment__extrinsic_driver_evaluate_range', PACKAGE = 'plant', obj_, driver_name, u)
+}
+
+FF16_Environment__get_extrinsic_driver_names <- function(obj_) {
+    .Call('_plant_FF16_Environment__get_extrinsic_driver_names', PACKAGE = 'plant', obj_)
+}
+
 FF16_Environment__canopy_openness <- function(obj_, height) {
     .Call('_plant_FF16_Environment__canopy_openness', PACKAGE = 'plant', obj_, height)
 }
@@ -3305,10 +3325,6 @@ FF16_Environment__set_fixed_environment <- function(obj_, value, height_max) {
     invisible(.Call('_plant_FF16_Environment__set_fixed_environment', PACKAGE = 'plant', obj_, value, height_max))
 }
 
-FF16_Environment__set_soil_infiltration_rate <- function(obj_, rate) {
-    invisible(.Call('_plant_FF16_Environment__set_soil_infiltration_rate', PACKAGE = 'plant', obj_, rate))
-}
-
 FF16_Environment__set_soil_water_state <- function(obj_, state) {
     invisible(.Call('_plant_FF16_Environment__set_soil_water_state', PACKAGE = 'plant', obj_, state))
 }

diff --git a/R/RcppR6.R b/R/RcppR6.R
@@ -1,6 +1,6 @@
 ## Generated by RcppR6: do not edit by hand
 ## Version: 0.2.4
-## Hash: f6b511a10ca3bb03a599e71a6fa83092
+## Hash: 2cb4fee1111abff105bb05242a8b85a0
 
 ##' @importFrom Rcpp evalCpp
 ##' @importFrom R6 R6Class
@@ -4814,6 +4814,21 @@ StochasticPatchRunner <- function(T, E) {
       initialize = function(ptr) {
         self$.ptr <- ptr
       },
+      set_extrinsic_driver = function(driver_name, x, y) {
+        FF16_Environment__set_extrinsic_driver(self, driver_name, x, y)
+      },
+      extrinsic_driver_extrapolation = function(driver_name, extrapolate) {
+        FF16_Environment__extrinsic_driver_extrapolation(self, driver_name, extrapolate)
+      },
+      extrinsic_driver_evaluate = function(driver_name, u) {
+        FF16_Environment__extrinsic_driver_evaluate(self, driver_name, u)
+      },
+      extrinsic_driver_evaluate_range = function(driver_name, u) {
+        FF16_Environment__extrinsic_driver_evaluate_range(self, driver_name, u)
+      },
+      get_extrinsic_driver_names = function() {
+        FF16_Environment__get_extrinsic_driver_names(self)
+      },
       canopy_openness = function(height) {
         FF16_Environment__canopy_openness(self, height)
       },
@@ -4823,9 +4838,6 @@ StochasticPatchRunner <- function(T, E) {
       set_fixed_environment = function(value, height_max) {
         FF16_Environment__set_fixed_environment(self, value, height_max)
       },
-      set_soil_infiltration_rate = function(rate) {
-        FF16_Environment__set_soil_infiltration_rate(self, rate)
-      },
       set_soil_water_state = function(state) {
         FF16_Environment__set_soil_water_state(self, state)
       },

diff --git a/R/ff16w.R b/R/ff16w.R
@@ -64,16 +64,15 @@ FF16w_StochasticPatchRunner <- function(p) {
 ## Helper:
 ##' @export
 ##' @rdname FF16_Environment
-##' @param infil_rate rate of water entering the first layer
-##' @param n_layers the number of layers
-##' @param init starting conditions
+##' @param soil_number_of_depths the number of soil layers
+##' @param rainfall constant function value for rainfall driver, y = rainfall
 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,
-                                   soil_infiltration_rate = 0.0) {
+                                   rainfall = 1) {
 
   if(soil_number_of_depths < 1)
     stop("FF16w Environment must have at least one soil layer")
@@ -85,15 +84,18 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4,
                      canopy_light_nbase, 
                      canopy_light_max_depth)
 
-  # there might be a better way to skip this if using defaultss
+  # there might be a better way to skip this if using defaults
   if(sum(soil_initial_state) > 0.0) {
     if(soil_number_of_depths != length(soil_initial_state))
       stop("Not enough starting points for all layers")
 
     e$set_soil_water_state(soil_initial_state)
   }
 
-  e$set_soil_infiltration_rate(soil_infiltration_rate)
+  x <- seq(0, 10, 1)
+  y <- rep(rainfall, 11)
+  e$set_extrinsic_driver("rainfall", x, y)
+  e$extrinsic_driver_extrapolation("rainfall", TRUE);
 
   return(e)
 }

diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml
@@ -871,6 +871,22 @@ FF16_Environment:
       FF16_Environment object
       @export
   methods:
+    # unforunately need to expose these base Environment functions for every derived class...
+    set_extrinsic_driver:
+      args: [driver_name: "std::string", x: "std::vector<double>", y: "std::vector<double>"]
+      return_type: void
+    extrinsic_driver_extrapolation:
+      args: [driver_name: "std::string", extrapolate: bool]
+      return_type: void
+    extrinsic_driver_evaluate:
+      args: [driver_name: "std::string", u: double]
+      return_type: double
+    extrinsic_driver_evaluate_range:
+      args: [driver_name: "std::string", u: "std::vector<double>"]
+      return_type: "std::vector<double>"
+    get_extrinsic_driver_names:
+      args: []
+      return_type: "std::vector<std::string>"
     canopy_openness:
       args: [height: double]
       return_type: double
@@ -879,9 +895,6 @@ FF16_Environment:
     set_fixed_environment:
       args: [value: double, height_max: double]
       return_type: void
-    set_soil_infiltration_rate:
-      args: [rate: double]
-      return_type: void
     set_soil_water_state:
       args: [state: "std::vector<double>"]
       return_type: void

diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h
@@ -8,6 +8,7 @@
 #include <plant/ode_interface.h>
 #include <plant/internals.h>
 #include <plant/util.h>
+#include <unordered_map>
 #include <Rcpp.h>
 
 using namespace Rcpp;
@@ -69,6 +70,44 @@ class Environment {
   size_t species_arriving_index;
 
   Internals vars;
+
+  /* EXSTRINSIC DRIVERS API*/
+
+  // initialise spline of driver with x, y control points
+  void set_extrinsic_driver(std::string driver_name, std::vector<double> const& x, std::vector<double> const& y) {
+    // if we wanted to be faster we could skip this check (but less safe)
+    if (extrinsic_drivers.count(driver_name) == 0) {
+      util::stop(driver_name + " doesn't exist in the list of extrinsic_drivers.");
+    } else {
+      extrinsic_drivers.at(driver_name).init(x, y);
+      extrinsic_driver_extrapolation(driver_name, false); // default no extrapolation
+    }
+  }
+
+  void extrinsic_driver_extrapolation(std::string driver_name, bool extrapolate) {
+    extrinsic_drivers.at(driver_name).setExtrapolate(extrapolate);
+  }
+
+  // evaluate/query interpolated spline for driver at point u, return s(u), where s is interpolated function
+  double extrinsic_driver_evaluate(std::string driver_name, double u) const {
+    return extrinsic_drivers.at(driver_name).eval(u);
+  }
+
+  // // evaluate/query interpolated spline for driver at vector of points, return vector of values
+  std::vector<double> extrinsic_driver_evaluate_range(std::string driver_name, std::vector<double> u) const {
+    return extrinsic_drivers.at(driver_name).r_eval(u);
+  }
+
+  std::vector<std::string> get_extrinsic_driver_names() {
+    auto ret = std::vector<std::string>();
+    for (auto const& driver : extrinsic_drivers) {
+      ret.push_back(driver.first);
+    }
+    return ret;
+  }
+
+protected:
+  std::unordered_map<std::string, interpolator::Interpolator> extrinsic_drivers;
 };
 
 }

diff --git a/inst/include/plant/interpolator.h b/inst/include/plant/interpolator.h
@@ -16,16 +16,20 @@ class Interpolator {
   void initialise();
 
   void add_point(double xi, double yi);
+  void add_point_sorted(double xi, double yi);
   void clear();
 
   double eval(double u) const;
+  double operator()(double u) const;
   size_t size() const;
 
   double min() const;
   double max() const;
+  void setExtrapolate(bool e);
 
   std::vector<double> get_x() const;
   std::vector<double> get_y() const;
+
 
   // * R interface
   SEXP r_get_xy() const;
@@ -35,7 +39,8 @@ class Interpolator {
   void check_active() const;
   std::vector<double> x, y;
   tk::spline tk_spline;
-  bool active;
+  bool active = false;
+  bool extrapolate = true;
 };
 
 }

diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h
@@ -4,6 +4,8 @@
 
 #include <plant/environment.h>
 #include <plant/canopy.h>
+#include <plant/interpolator.h>
+#include <iostream>
 
 using namespace Rcpp;
 
@@ -19,14 +21,16 @@ class FF16_Environment : public Environment {
       : canopy_rescale_usually(canopy_rescale_usually) {
     time = 0.0;
     canopy = Canopy();
-    soil_infiltration_rate = 0.0;
     vars = Internals(soil_number_of_depths);
     set_soil_water_state(std::vector<double>(soil_number_of_depths, 0.0));
+    extrinsic_drivers["rainfall"] = interpolator::Interpolator();
   };
 
+
   // Light interface
   bool canopy_rescale_usually;
 
+  // private?
   Canopy canopy;
 
   // Should this be here or in canopy?
@@ -51,12 +55,9 @@ class FF16_Environment : public Environment {
     canopy.r_init_interpolators(state);
   }
 
-  // Soil interface
-  double soil_infiltration_rate;
-
   virtual void compute_rates() {
     for (size_t i = 0; i < vars.state_size; i++) {
-      vars.set_rate(i, soil_infiltration_rate / (i+1));
+      vars.set_rate(i, extrinsic_driver_evaluate("rainfall", time) / (i+1));
     }
   }
 
@@ -74,10 +75,6 @@ class FF16_Environment : public Environment {
     }
   }
 
-  void set_soil_infiltration_rate(double rate) {
-    soil_infiltration_rate = rate;
-  }
-
   // Core functions
   template <typename Function>
   void compute_environment(Function f_compute_competition, double height_max) {
@@ -92,7 +89,6 @@ class FF16_Environment : public Environment {
   void clear_environment() {
     canopy.clear();
   }
-
 };
 
 inline Rcpp::NumericMatrix get_state(const FF16_Environment environment) {