From d4472156f94b1075e9e9a94eb7a4b5370e4e5fff Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Mon, 8 Nov 2021 12:44:08 +1100 Subject: [PATCH 01/22] Implemented functions to explose auxillary variables in R API --- DESCRIPTION | 2 +- R/RcppExports.R | 16 +++++++++ R/RcppR6.R | 30 ++++++++++++++++- inst/RcppR6_classes.yml | 1 + inst/include/plant/cohort.h | 9 +++++ inst/include/plant/individual.h | 10 +++++- inst/include/plant/internals.h | 2 ++ inst/include/plant/ode_interface.h | 17 ++++++++++ inst/include/plant/patch.h | 8 +++++ inst/include/plant/species.h | 6 ++++ src/RcppExports.cpp | 53 ++++++++++++++++++++++++++++++ src/RcppR6.cpp | 20 +++++++++++ 12 files changed, 171 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 8e97382ca..9f1ea81ac 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -43,4 +43,4 @@ Remotes: VignetteBuilder: knitr URL: https://github.com/traitecoevo/plant BugReports: https://github.com/traitecoevo/plant/issues -RoxygenNote: 7.1.1 +RoxygenNote: 7.1.2 diff --git a/R/RcppExports.R b/R/RcppExports.R index 1d175ab7f..d506b111c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1685,6 +1685,10 @@ Patch___FF16__FF16_Env__ode_rates__get <- function(obj_) { .Call('_plant_Patch___FF16__FF16_Env__ode_rates__get', PACKAGE = 'plant', obj_) } +Patch___FF16__FF16_Env__ode_aux__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__ode_aux__get', PACKAGE = 'plant', obj_) +} + Patch___FF16__FF16_Env__cohort_ode_size__get <- function(obj_) { .Call('_plant_Patch___FF16__FF16_Env__cohort_ode_size__get', PACKAGE = 'plant', obj_) } @@ -1789,6 +1793,10 @@ Patch___FF16w__FF16_Env__ode_rates__get <- function(obj_) { .Call('_plant_Patch___FF16w__FF16_Env__ode_rates__get', PACKAGE = 'plant', obj_) } +Patch___FF16w__FF16_Env__ode_aux__get <- function(obj_) { + .Call('_plant_Patch___FF16w__FF16_Env__ode_aux__get', PACKAGE = 'plant', obj_) +} + Patch___FF16w__FF16_Env__cohort_ode_size__get <- function(obj_) { .Call('_plant_Patch___FF16w__FF16_Env__cohort_ode_size__get', PACKAGE = 'plant', obj_) } @@ -1893,6 +1901,10 @@ Patch___FF16r__FF16_Env__ode_rates__get <- function(obj_) { .Call('_plant_Patch___FF16r__FF16_Env__ode_rates__get', PACKAGE = 'plant', obj_) } +Patch___FF16r__FF16_Env__ode_aux__get <- function(obj_) { + .Call('_plant_Patch___FF16r__FF16_Env__ode_aux__get', PACKAGE = 'plant', obj_) +} + Patch___FF16r__FF16_Env__cohort_ode_size__get <- function(obj_) { .Call('_plant_Patch___FF16r__FF16_Env__cohort_ode_size__get', PACKAGE = 'plant', obj_) } @@ -1997,6 +2009,10 @@ Patch___K93__K93_Env__ode_rates__get <- function(obj_) { .Call('_plant_Patch___K93__K93_Env__ode_rates__get', PACKAGE = 'plant', obj_) } +Patch___K93__K93_Env__ode_aux__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__ode_aux__get', PACKAGE = 'plant', obj_) +} + Patch___K93__K93_Env__cohort_ode_size__get <- function(obj_) { .Call('_plant_Patch___K93__K93_Env__cohort_ode_size__get', PACKAGE = 'plant', obj_) } diff --git a/R/RcppR6.R b/R/RcppR6.R index 6d861d52b..19c9c346e 100644 --- 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: 4458bbba521dd4dfd30c08f98140aed9 +## Hash: ae9874a37766c323fc0f943ef595d917 ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -2518,6 +2518,13 @@ Patch <- function(T, E) { stop("Patch$ode_rates is read-only") } }, + ode_aux = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__ode_aux__get(self) + } else { + stop("Patch$ode_aux is read-only") + } + }, cohort_ode_size = function(value) { if (missing(value)) { Patch___FF16__FF16_Env__cohort_ode_size__get(self) @@ -2653,6 +2660,13 @@ Patch <- function(T, E) { stop("Patch$ode_rates is read-only") } }, + ode_aux = function(value) { + if (missing(value)) { + Patch___FF16w__FF16_Env__ode_aux__get(self) + } else { + stop("Patch$ode_aux is read-only") + } + }, cohort_ode_size = function(value) { if (missing(value)) { Patch___FF16w__FF16_Env__cohort_ode_size__get(self) @@ -2788,6 +2802,13 @@ Patch <- function(T, E) { stop("Patch$ode_rates is read-only") } }, + ode_aux = function(value) { + if (missing(value)) { + Patch___FF16r__FF16_Env__ode_aux__get(self) + } else { + stop("Patch$ode_aux is read-only") + } + }, cohort_ode_size = function(value) { if (missing(value)) { Patch___FF16r__FF16_Env__cohort_ode_size__get(self) @@ -2923,6 +2944,13 @@ Patch <- function(T, E) { stop("Patch$ode_rates is read-only") } }, + ode_aux = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__ode_aux__get(self) + } else { + stop("Patch$ode_aux is read-only") + } + }, cohort_ode_size = function(value) { if (missing(value)) { Patch___K93__K93_Env__cohort_ode_size__get(self) diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index 62236341b..bf63fbe5e 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -470,6 +470,7 @@ Patch: ode_time: {type: double, access: function, name_cpp: "plant::ode::r_ode_time"} ode_state: {type: "plant::ode::state_type", access: function, name_cpp: "plant::ode::r_ode_state"} ode_rates: {type: "plant::ode::state_type", access: function, name_cpp: "plant::ode::r_ode_rates"} + ode_aux: {type: "plant::ode::state_type", access: function, name_cpp: "plant::ode::r_ode_aux"} cohort_ode_size: {type: size_t, access: member} methods: introduce_new_cohort: diff --git a/inst/include/plant/cohort.h b/inst/include/plant/cohort.h index 27ecfbb69..700da37c5 100644 --- a/inst/include/plant/cohort.h +++ b/inst/include/plant/cohort.h @@ -44,6 +44,7 @@ class Cohort { ode::const_iterator set_ode_state(ode::const_iterator it); ode::iterator ode_state(ode::iterator it) const; ode::iterator ode_rates(ode::iterator it) const; + ode::iterator ode_aux(ode::iterator it) const; static std::vector ode_names() { std::vector plant_names = strategy_type::state_names(); @@ -201,6 +202,14 @@ ode::iterator Cohort::ode_rates(ode::iterator it) const { return it; } +template +ode::iterator Cohort::ode_aux(ode::iterator it) const { + for (size_t i = 0; i < plant.aux_size(); i++) { + *it++ = plant.aux(i); + } + return it; +} + template Cohort make_cohort(typename Cohort::strategy_type s) { return Cohort(make_strategy_ptr(s)); diff --git a/inst/include/plant/individual.h b/inst/include/plant/individual.h index f1454ce8e..dcb10f7bc 100644 --- a/inst/include/plant/individual.h +++ b/inst/include/plant/individual.h @@ -76,7 +76,8 @@ template class Individual { static size_t ode_size() { return strategy_type::state_size(); } static std::vector ode_names() { return strategy_type::state_names(); } - size_t aux_size() { return strategy->aux_size(); } + // why doesn't strategy->aux_size() also need to be const-qualified? is it because strategy is a pointer? + size_t aux_size() const { return strategy->aux_size(); } std::vector aux_names() { return strategy->aux_names(); } ode::const_iterator set_ode_state(ode::const_iterator it) { @@ -99,6 +100,13 @@ template class Individual { return it; } + ode::iterator ode_aux(ode::iterator it) const { + for (size_t i = 0; i < vars.aux_size; i++) { + *it++ = vars.auxs[i]; + } + return it; + } + // Single individual methods // Used in the stochastic model: diff --git a/inst/include/plant/internals.h b/inst/include/plant/internals.h index ef26f9ab3..828fb7117 100644 --- a/inst/include/plant/internals.h +++ b/inst/include/plant/internals.h @@ -27,6 +27,8 @@ class Internals { size_t state_size; size_t aux_size; + + // Perhaps make these private so the () overloads below have some use std::vector states; std::vector rates; std::vector auxs; diff --git a/inst/include/plant/ode_interface.h b/inst/include/plant/ode_interface.h index 3ff2deab6..0921c3adb 100644 --- a/inst/include/plant/ode_interface.h +++ b/inst/include/plant/ode_interface.h @@ -65,6 +65,16 @@ iterator ode_rates(ForwardIterator first, ForwardIterator last, return it; } +template +iterator ode_aux(ForwardIterator first, ForwardIterator last, + iterator it) { + while (first != last) { + it = first->ode_aux(it); + ++first; + } + return it; +} + template typename std::enable_if::value, double>::type ode_time(const T& obj) { @@ -146,6 +156,13 @@ state_type r_ode_rates(const T& obj) { return dydt; } +template +state_type r_ode_aux(const T& obj) { + state_type dydt(obj.ode_size()); + obj.ode_aux(dydt.begin()); + return dydt; +} + } } diff --git a/inst/include/plant/patch.h b/inst/include/plant/patch.h index c7ac34fac..6bbf66b5f 100644 --- a/inst/include/plant/patch.h +++ b/inst/include/plant/patch.h @@ -56,6 +56,7 @@ class Patch { ode::const_iterator set_ode_state(ode::const_iterator it, double time); ode::iterator ode_state(ode::iterator it) const; ode::iterator ode_rates(ode::iterator it) const; + ode::iterator ode_aux(ode::iterator it) const; // * R interface // Data accessors: @@ -275,6 +276,13 @@ ode::iterator Patch::ode_rates(ode::iterator it) const { return it; } +template +ode::iterator Patch::ode_aux(ode::iterator it) const { + it = ode::ode_aux(species.begin(), species.end(), it); + //it = environment.ode_rates(it); + return it; +} + } #endif diff --git a/inst/include/plant/species.h b/inst/include/plant/species.h index 2acdc7f54..59b9f0edb 100644 --- a/inst/include/plant/species.h +++ b/inst/include/plant/species.h @@ -39,6 +39,7 @@ class Species { ode::const_iterator set_ode_state(ode::const_iterator it); ode::iterator ode_state(ode::iterator it) const; ode::iterator ode_rates(ode::iterator it) const; + ode::iterator ode_aux(ode::iterator it) const; // * R interface std::vector r_heights() const; @@ -208,6 +209,11 @@ ode::iterator Species::ode_rates(ode::iterator it) const { } //double sum_aux(int index) {} +template +ode::iterator Species::ode_aux(ode::iterator it) const { + return ode::ode_aux(cohorts.begin(), cohorts.end(), it); +} + template std::vector Species::r_heights() const { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 435037368..60ed79d96 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -6,6 +6,11 @@ using namespace Rcpp; +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + // Lorenz__ctor plant::ode::test::Lorenz Lorenz__ctor(double sigma, double R, double b); RcppExport SEXP _plant_Lorenz__ctor(SEXP sigmaSEXP, SEXP RSEXP, SEXP bSEXP) { @@ -4734,6 +4739,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// Patch___FF16__FF16_Env__ode_aux__get +plant::ode::state_type Patch___FF16__FF16_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_Patch___FF16__FF16_Env__ode_aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Patch___FF16__FF16_Env__ode_aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // Patch___FF16__FF16_Env__cohort_ode_size__get size_t Patch___FF16__FF16_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_Patch___FF16__FF16_Env__cohort_ode_size__get(SEXP obj_SEXP) { @@ -5030,6 +5046,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// Patch___FF16w__FF16_Env__ode_aux__get +plant::ode::state_type Patch___FF16w__FF16_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_Patch___FF16w__FF16_Env__ode_aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Patch___FF16w__FF16_Env__ode_aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // Patch___FF16w__FF16_Env__cohort_ode_size__get size_t Patch___FF16w__FF16_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_Patch___FF16w__FF16_Env__cohort_ode_size__get(SEXP obj_SEXP) { @@ -5326,6 +5353,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// Patch___FF16r__FF16_Env__ode_aux__get +plant::ode::state_type Patch___FF16r__FF16_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_Patch___FF16r__FF16_Env__ode_aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Patch___FF16r__FF16_Env__ode_aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // Patch___FF16r__FF16_Env__cohort_ode_size__get size_t Patch___FF16r__FF16_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_Patch___FF16r__FF16_Env__cohort_ode_size__get(SEXP obj_SEXP) { @@ -5622,6 +5660,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// Patch___K93__K93_Env__ode_aux__get +plant::ode::state_type Patch___K93__K93_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_Patch___K93__K93_Env__ode_aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Patch___K93__K93_Env__ode_aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // Patch___K93__K93_Env__cohort_ode_size__get size_t Patch___K93__K93_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_Patch___K93__K93_Env__cohort_ode_size__get(SEXP obj_SEXP) { @@ -10083,6 +10132,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_Patch___FF16__FF16_Env__ode_time__get", (DL_FUNC) &_plant_Patch___FF16__FF16_Env__ode_time__get, 1}, {"_plant_Patch___FF16__FF16_Env__ode_state__get", (DL_FUNC) &_plant_Patch___FF16__FF16_Env__ode_state__get, 1}, {"_plant_Patch___FF16__FF16_Env__ode_rates__get", (DL_FUNC) &_plant_Patch___FF16__FF16_Env__ode_rates__get, 1}, + {"_plant_Patch___FF16__FF16_Env__ode_aux__get", (DL_FUNC) &_plant_Patch___FF16__FF16_Env__ode_aux__get, 1}, {"_plant_Patch___FF16__FF16_Env__cohort_ode_size__get", (DL_FUNC) &_plant_Patch___FF16__FF16_Env__cohort_ode_size__get, 1}, {"_plant_Patch___FF16w__FF16_Env__ctor", (DL_FUNC) &_plant_Patch___FF16w__FF16_Env__ctor, 3}, {"_plant_Patch___FF16w__FF16_Env__introduce_new_cohort", (DL_FUNC) &_plant_Patch___FF16w__FF16_Env__introduce_new_cohort, 2}, @@ -10109,6 +10159,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_Patch___FF16w__FF16_Env__ode_time__get", (DL_FUNC) &_plant_Patch___FF16w__FF16_Env__ode_time__get, 1}, {"_plant_Patch___FF16w__FF16_Env__ode_state__get", (DL_FUNC) &_plant_Patch___FF16w__FF16_Env__ode_state__get, 1}, {"_plant_Patch___FF16w__FF16_Env__ode_rates__get", (DL_FUNC) &_plant_Patch___FF16w__FF16_Env__ode_rates__get, 1}, + {"_plant_Patch___FF16w__FF16_Env__ode_aux__get", (DL_FUNC) &_plant_Patch___FF16w__FF16_Env__ode_aux__get, 1}, {"_plant_Patch___FF16w__FF16_Env__cohort_ode_size__get", (DL_FUNC) &_plant_Patch___FF16w__FF16_Env__cohort_ode_size__get, 1}, {"_plant_Patch___FF16r__FF16_Env__ctor", (DL_FUNC) &_plant_Patch___FF16r__FF16_Env__ctor, 3}, {"_plant_Patch___FF16r__FF16_Env__introduce_new_cohort", (DL_FUNC) &_plant_Patch___FF16r__FF16_Env__introduce_new_cohort, 2}, @@ -10135,6 +10186,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_Patch___FF16r__FF16_Env__ode_time__get", (DL_FUNC) &_plant_Patch___FF16r__FF16_Env__ode_time__get, 1}, {"_plant_Patch___FF16r__FF16_Env__ode_state__get", (DL_FUNC) &_plant_Patch___FF16r__FF16_Env__ode_state__get, 1}, {"_plant_Patch___FF16r__FF16_Env__ode_rates__get", (DL_FUNC) &_plant_Patch___FF16r__FF16_Env__ode_rates__get, 1}, + {"_plant_Patch___FF16r__FF16_Env__ode_aux__get", (DL_FUNC) &_plant_Patch___FF16r__FF16_Env__ode_aux__get, 1}, {"_plant_Patch___FF16r__FF16_Env__cohort_ode_size__get", (DL_FUNC) &_plant_Patch___FF16r__FF16_Env__cohort_ode_size__get, 1}, {"_plant_Patch___K93__K93_Env__ctor", (DL_FUNC) &_plant_Patch___K93__K93_Env__ctor, 3}, {"_plant_Patch___K93__K93_Env__introduce_new_cohort", (DL_FUNC) &_plant_Patch___K93__K93_Env__introduce_new_cohort, 2}, @@ -10161,6 +10213,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_Patch___K93__K93_Env__ode_time__get", (DL_FUNC) &_plant_Patch___K93__K93_Env__ode_time__get, 1}, {"_plant_Patch___K93__K93_Env__ode_state__get", (DL_FUNC) &_plant_Patch___K93__K93_Env__ode_state__get, 1}, {"_plant_Patch___K93__K93_Env__ode_rates__get", (DL_FUNC) &_plant_Patch___K93__K93_Env__ode_rates__get, 1}, + {"_plant_Patch___K93__K93_Env__ode_aux__get", (DL_FUNC) &_plant_Patch___K93__K93_Env__ode_aux__get, 1}, {"_plant_Patch___K93__K93_Env__cohort_ode_size__get", (DL_FUNC) &_plant_Patch___K93__K93_Env__cohort_ode_size__get, 1}, {"_plant_SCM___FF16__FF16_Env__ctor", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__ctor, 3}, {"_plant_SCM___FF16__FF16_Env__run", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__run, 1}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index 28a382766..48f4a802f 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -1922,6 +1922,11 @@ plant::ode::state_type Patch___FF16__FF16_Env__ode_rates__get(plant::RcppR6::Rcp return plant::ode::r_ode_rates(*obj_); } +// [[Rcpp::export]] +plant::ode::state_type Patch___FF16__FF16_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_) { + return plant::ode::r_ode_aux(*obj_); +} + // [[Rcpp::export]] size_t Patch___FF16__FF16_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_) { return obj_->cohort_ode_size(); @@ -2038,6 +2043,11 @@ plant::ode::state_type Patch___FF16w__FF16_Env__ode_rates__get(plant::RcppR6::Rc return plant::ode::r_ode_rates(*obj_); } +// [[Rcpp::export]] +plant::ode::state_type Patch___FF16w__FF16_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_) { + return plant::ode::r_ode_aux(*obj_); +} + // [[Rcpp::export]] size_t Patch___FF16w__FF16_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_) { return obj_->cohort_ode_size(); @@ -2154,6 +2164,11 @@ plant::ode::state_type Patch___FF16r__FF16_Env__ode_rates__get(plant::RcppR6::Rc return plant::ode::r_ode_rates(*obj_); } +// [[Rcpp::export]] +plant::ode::state_type Patch___FF16r__FF16_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_) { + return plant::ode::r_ode_aux(*obj_); +} + // [[Rcpp::export]] size_t Patch___FF16r__FF16_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_) { return obj_->cohort_ode_size(); @@ -2270,6 +2285,11 @@ plant::ode::state_type Patch___K93__K93_Env__ode_rates__get(plant::RcppR6::RcppR return plant::ode::r_ode_rates(*obj_); } +// [[Rcpp::export]] +plant::ode::state_type Patch___K93__K93_Env__ode_aux__get(plant::RcppR6::RcppR6 > obj_) { + return plant::ode::r_ode_aux(*obj_); +} + // [[Rcpp::export]] size_t Patch___K93__K93_Env__cohort_ode_size__get(plant::RcppR6::RcppR6 > obj_) { return obj_->cohort_ode_size(); From 157bb29881fd568e9a485a46c1d51f8c37d60670 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Mon, 8 Nov 2021 15:01:59 +1100 Subject: [PATCH 02/22] Resized auxillary vector to appropriate size --- inst/include/plant/cohort.h | 2 ++ inst/include/plant/ode_interface.h | 12 +++++++++++- inst/include/plant/patch.h | 7 +++++++ inst/include/plant/species.h | 6 ++++++ 4 files changed, 26 insertions(+), 1 deletion(-) diff --git a/inst/include/plant/cohort.h b/inst/include/plant/cohort.h index 700da37c5..7617d9e4d 100644 --- a/inst/include/plant/cohort.h +++ b/inst/include/plant/cohort.h @@ -40,7 +40,9 @@ class Cohort { // NOTE: We are a time-independent model here so no need to pass // time in as an argument. All the bits involving time are taken // care of by Environment for us. + // +2 for log_density and offspring_production_dt static size_t ode_size() { return strategy_type::state_size() + 2; } + size_t aux_size() { return plant.aux_size(); } ode::const_iterator set_ode_state(ode::const_iterator it); ode::iterator ode_state(ode::iterator it) const; ode::iterator ode_rates(ode::iterator it) const; diff --git a/inst/include/plant/ode_interface.h b/inst/include/plant/ode_interface.h index 0921c3adb..6bf5b2aa1 100644 --- a/inst/include/plant/ode_interface.h +++ b/inst/include/plant/ode_interface.h @@ -35,6 +35,16 @@ size_t ode_size(ForwardIterator first, ForwardIterator last) { return ret; } +template +size_t aux_size(ForwardIterator first, ForwardIterator last) { + size_t ret = 0; + while (first != last) { + ret += first->aux_size(); + ++first; + } + return ret; +} + template const_iterator set_ode_state(ForwardIterator first, ForwardIterator last, const_iterator it) { @@ -158,7 +168,7 @@ state_type r_ode_rates(const T& obj) { template state_type r_ode_aux(const T& obj) { - state_type dydt(obj.ode_size()); + state_type dydt(obj.aux_size()); obj.ode_aux(dydt.begin()); return dydt; } diff --git a/inst/include/plant/patch.h b/inst/include/plant/patch.h index 6bbf66b5f..b8c728a19 100644 --- a/inst/include/plant/patch.h +++ b/inst/include/plant/patch.h @@ -52,6 +52,7 @@ class Patch { // * ODE interface size_t ode_size() const; + size_t aux_size() const; double ode_time() const; ode::const_iterator set_ode_state(ode::const_iterator it, double time); ode::iterator ode_state(ode::iterator it) const; @@ -241,6 +242,12 @@ size_t Patch::ode_size() const { return ode::ode_size(species.begin(), species.end()) + environment.ode_size(); } +template +size_t Patch::aux_size() const { + // no use for auxillary environment variables (yet) + return ode::aux_size(species.begin(), species.end());// + environment.ode_size(); +} + template double Patch::ode_time() const { return time(); diff --git a/inst/include/plant/species.h b/inst/include/plant/species.h index 59b9f0edb..79ac8ba7e 100644 --- a/inst/include/plant/species.h +++ b/inst/include/plant/species.h @@ -36,6 +36,7 @@ class Species { // time in as an argument. All the bits involving time are taken // care of by Environment for us. size_t ode_size() const; + size_t aux_size() const; ode::const_iterator set_ode_state(ode::const_iterator it); ode::iterator ode_state(ode::iterator it) const; ode::iterator ode_rates(ode::iterator it) const; @@ -193,6 +194,11 @@ size_t Species::ode_size() const { return size() * cohort_type::ode_size(); } +template +size_t Species::aux_size() const { + return size() * strategy->aux_size(); +} + template ode::const_iterator Species::set_ode_state(ode::const_iterator it) { return ode::set_ode_state(cohorts.begin(), cohorts.end(), it); From d64bc0dafa25464137ded817ef85bfa0a752b5f5 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Mon, 8 Nov 2021 17:13:22 +1100 Subject: [PATCH 03/22] Exposed auxillary variables in R API + basic test --- R/RcppExports.R | 16 ++++++++++ R/RcppR6.R | 30 +++++++++++++++++- R/scm_support.R | 28 +++++------------ inst/RcppR6_classes.yml | 4 +++ inst/include/plant.h | 1 + inst/include/plant/cohort.h | 2 +- inst/include/plant/get_aux.h | 52 +++++++++++++++++++++++++++++++ inst/include/plant/species.h | 19 +++++++++-- src/RcppExports.cpp | 48 ++++++++++++++++++++++++++++ src/RcppR6.cpp | 20 ++++++++++++ tests/testthat/test-scm-support.R | 15 +++++++++ 11 files changed, 211 insertions(+), 24 deletions(-) create mode 100644 inst/include/plant/get_aux.h diff --git a/R/RcppExports.R b/R/RcppExports.R index d506b111c..93520faaa 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2085,6 +2085,10 @@ SCM___FF16__FF16_Env__state__get <- function(obj_) { .Call('_plant_SCM___FF16__FF16_Env__state__get', PACKAGE = 'plant', obj_) } +SCM___FF16__FF16_Env__aux__get <- function(obj_) { + .Call('_plant_SCM___FF16__FF16_Env__aux__get', PACKAGE = 'plant', obj_) +} + SCM___FF16__FF16_Env__use_ode_times__get <- function(obj_) { .Call('_plant_SCM___FF16__FF16_Env__use_ode_times__get', PACKAGE = 'plant', obj_) } @@ -2165,6 +2169,10 @@ SCM___FF16w__FF16_Env__state__get <- function(obj_) { .Call('_plant_SCM___FF16w__FF16_Env__state__get', PACKAGE = 'plant', obj_) } +SCM___FF16w__FF16_Env__aux__get <- function(obj_) { + .Call('_plant_SCM___FF16w__FF16_Env__aux__get', PACKAGE = 'plant', obj_) +} + SCM___FF16w__FF16_Env__use_ode_times__get <- function(obj_) { .Call('_plant_SCM___FF16w__FF16_Env__use_ode_times__get', PACKAGE = 'plant', obj_) } @@ -2245,6 +2253,10 @@ SCM___FF16r__FF16_Env__state__get <- function(obj_) { .Call('_plant_SCM___FF16r__FF16_Env__state__get', PACKAGE = 'plant', obj_) } +SCM___FF16r__FF16_Env__aux__get <- function(obj_) { + .Call('_plant_SCM___FF16r__FF16_Env__aux__get', PACKAGE = 'plant', obj_) +} + SCM___FF16r__FF16_Env__use_ode_times__get <- function(obj_) { .Call('_plant_SCM___FF16r__FF16_Env__use_ode_times__get', PACKAGE = 'plant', obj_) } @@ -2325,6 +2337,10 @@ SCM___K93__K93_Env__state__get <- function(obj_) { .Call('_plant_SCM___K93__K93_Env__state__get', PACKAGE = 'plant', obj_) } +SCM___K93__K93_Env__aux__get <- function(obj_) { + .Call('_plant_SCM___K93__K93_Env__aux__get', PACKAGE = 'plant', obj_) +} + SCM___K93__K93_Env__use_ode_times__get <- function(obj_) { .Call('_plant_SCM___K93__K93_Env__use_ode_times__get', PACKAGE = 'plant', obj_) } diff --git a/R/RcppR6.R b/R/RcppR6.R index 19c9c346e..60bf49c77 100644 --- 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: ae9874a37766c323fc0f943ef595d917 +## Hash: edd85bd671370b151c2bba310b3f12e3 ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -3064,6 +3064,13 @@ SCM <- function(T, E) { stop("SCM$state is read-only") } }, + aux = function(value) { + if (missing(value)) { + SCM___FF16__FF16_Env__aux__get(self) + } else { + stop("SCM$aux is read-only") + } + }, use_ode_times = function(value) { if (missing(value)) { SCM___FF16__FF16_Env__use_ode_times__get(self) @@ -3175,6 +3182,13 @@ SCM <- function(T, E) { stop("SCM$state is read-only") } }, + aux = function(value) { + if (missing(value)) { + SCM___FF16w__FF16_Env__aux__get(self) + } else { + stop("SCM$aux is read-only") + } + }, use_ode_times = function(value) { if (missing(value)) { SCM___FF16w__FF16_Env__use_ode_times__get(self) @@ -3286,6 +3300,13 @@ SCM <- function(T, E) { stop("SCM$state is read-only") } }, + aux = function(value) { + if (missing(value)) { + SCM___FF16r__FF16_Env__aux__get(self) + } else { + stop("SCM$aux is read-only") + } + }, use_ode_times = function(value) { if (missing(value)) { SCM___FF16r__FF16_Env__use_ode_times__get(self) @@ -3397,6 +3418,13 @@ SCM <- function(T, E) { stop("SCM$state is read-only") } }, + aux = function(value) { + if (missing(value)) { + SCM___K93__K93_Env__aux__get(self) + } else { + stop("SCM$aux is read-only") + } + }, use_ode_times = function(value) { if (missing(value)) { SCM___K93__K93_Env__use_ode_times__get(self) diff --git a/R/scm_support.R b/R/scm_support.R index c2d4e00e0..e56ffe2a2 100644 --- a/R/scm_support.R +++ b/R/scm_support.R @@ -100,31 +100,25 @@ run_scm <- function(p, env = make_environment(parameters = p), ##' @param p A \code{Parameters} object ##' @param env Environment object (defaults to FF16_Environment) ##' @param ctrl Control object -##' @param include_competition_effect Include total leaf area (will change; see -##' issue #138) +##' @param collect_auxillary_variables Return additional strategy variables (eg +##' competition_effect) ##' @author Rich FitzJohn ##' @export run_scm_collect <- function(p, env = make_environment(parameters = p), ctrl = scm_base_control(), - include_competition_effect=FALSE) { + collect_auxillary_variables=FALSE) { collect_default <- function(scm) { scm$state } - collect_competition_effect <- function(scm) { + collect_aux <- function(scm) { ret <- scm$state - competition_effect <- numeric(length(ret$species)) - for (i in seq_along(ret$species)) { - ## ret$species[[i]] <- rbind( - ## ret$species[[i]] - ## competition_effect=c(scm$patch$species[[i]]$competition_effects, 0.0)) - competition_effect[i] <- scm$patch$species[[i]]$compute_competition(0.0) - } - ret$competition_effect <- competition_effect + aux <- scm$aux + ret$species <- mapply(rbind, ret$species, aux$species, SIMPLIFY=FALSE) ret } - collect <- if (include_competition_effect) collect_competition_effect else collect_default + collect <- if (collect_auxillary_variables) collect_aux else collect_default types <- extract_RcppR6_template_types(p, "Parameters") - + scm <- do.call('SCM', types)(p, env, ctrl) res <- list(collect(scm)) @@ -157,12 +151,6 @@ run_scm_collect <- function(p, env = make_environment(parameters = p), net_reproduction_ratios=scm$net_reproduction_ratios, patch_density=patch_density, p=p) - - if (include_competition_effect) { - ret$competition_effect <- do.call("rbind", lapply(res, "[[", "competition_effect")) - } - - ret } ##' Functions for reconstructing a Patch from an SCM diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index bf63fbe5e..cc7652e7d 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -574,6 +574,10 @@ SCM: type: "Rcpp::List" access: function name_cpp: "plant::get_state" + aux: + type: "Rcpp::List" + access: function + name_cpp: "plant::get_aux" use_ode_times: type: bool access: member diff --git a/inst/include/plant.h b/inst/include/plant.h index b240d095c..0fd634a05 100644 --- a/inst/include/plant.h +++ b/inst/include/plant.h @@ -71,6 +71,7 @@ #include #include #include +#include #include #endif diff --git a/inst/include/plant/cohort.h b/inst/include/plant/cohort.h index 7617d9e4d..e2eb73a6f 100644 --- a/inst/include/plant/cohort.h +++ b/inst/include/plant/cohort.h @@ -42,7 +42,7 @@ class Cohort { // care of by Environment for us. // +2 for log_density and offspring_production_dt static size_t ode_size() { return strategy_type::state_size() + 2; } - size_t aux_size() { return plant.aux_size(); } + size_t aux_size() const { return plant.aux_size(); } ode::const_iterator set_ode_state(ode::const_iterator it); ode::iterator ode_state(ode::iterator it) const; ode::iterator ode_rates(ode::iterator it) const; diff --git a/inst/include/plant/get_aux.h b/inst/include/plant/get_aux.h new file mode 100644 index 000000000..5d9410881 --- /dev/null +++ b/inst/include/plant/get_aux.h @@ -0,0 +1,52 @@ +// -*-c++-*- +#ifndef PLANT_GET_AUX_H_ +#define PLANT_GET_AUX_H_ + +#include +#include +#include +#include + +namespace plant { + +template +Rcpp::NumericMatrix::iterator get_aux(const Cohort& cohort, + Rcpp::NumericMatrix::iterator it) { + std::vector tmp = ode::r_ode_aux(cohort); + return std::copy(tmp.begin(), tmp.end(), it); +} + +template +Rcpp::NumericMatrix get_aux(const Species& species) { + typedef Cohort cohort_type; + size_t aux_size = species.strategy_aux_size(), np = species.size(); + Rcpp::NumericMatrix ret(static_cast(aux_size), np + 1); // +1 is seed + Rcpp::NumericMatrix::iterator it = ret.begin(); + for (size_t i = 0; i < np; ++i) { + it = get_aux(species.r_cohort_at(i), it); + } + it = get_aux(species.r_new_cohort(), it); + ret.attr("dimnames") = + Rcpp::List::create(species.aux_names(), R_NilValue); + return ret; +} + +template +Rcpp::List get_aux(const Patch& patch) { + Rcpp::List ret; + for (size_t i = 0; i < patch.size(); ++i) { + ret.push_back(get_aux(patch.at(i))); + } + return ret; +} + +template +Rcpp::List get_aux(const SCM& scm) { + using namespace Rcpp; + const Patch& patch = scm.r_patch(); + return List::create(_["species"] = get_aux(patch)); +} + +} + +#endif \ No newline at end of file diff --git a/inst/include/plant/species.h b/inst/include/plant/species.h index 79ac8ba7e..07d04d7cd 100644 --- a/inst/include/plant/species.h +++ b/inst/include/plant/species.h @@ -37,6 +37,9 @@ class Species { // care of by Environment for us. size_t ode_size() const; size_t aux_size() const; + size_t strategy_aux_size() const; + std::vector aux_names() const; + ode::const_iterator set_ode_state(ode::const_iterator it); ode::iterator ode_state(ode::iterator it) const; ode::iterator ode_rates(ode::iterator it) const; @@ -116,7 +119,7 @@ double Species::height_max() const { // height is zero, because it will be zero for all cohorts further down // the list. // -// NOTE: This is simply performing numerical integration, via the +// NOTE: This is simply performing numerical integration, via the // trapezium rule, of the compute_competition with respect to plant // height. You'd think that this would be nicer to do in terms of a // call to an external trapezium integration function, but building @@ -194,9 +197,21 @@ size_t Species::ode_size() const { return size() * cohort_type::ode_size(); } +// bit clunky... template size_t Species::aux_size() const { - return size() * strategy->aux_size(); + return size() * strategy_aux_size(); +} + +// these 2 only really used in get_aux.h +template +size_t Species::strategy_aux_size() const { + return strategy->aux_size(); +} + +template +std::vector Species::aux_names() const { + return strategy->aux_names(); } template diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 60ed79d96..2ce5f7266 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -5871,6 +5871,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// SCM___FF16__FF16_Env__aux__get +Rcpp::List SCM___FF16__FF16_Env__aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_SCM___FF16__FF16_Env__aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(SCM___FF16__FF16_Env__aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // SCM___FF16__FF16_Env__use_ode_times__get bool SCM___FF16__FF16_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_SCM___FF16__FF16_Env__use_ode_times__get(SEXP obj_SEXP) { @@ -6093,6 +6104,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// SCM___FF16w__FF16_Env__aux__get +Rcpp::List SCM___FF16w__FF16_Env__aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_SCM___FF16w__FF16_Env__aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(SCM___FF16w__FF16_Env__aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // SCM___FF16w__FF16_Env__use_ode_times__get bool SCM___FF16w__FF16_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_SCM___FF16w__FF16_Env__use_ode_times__get(SEXP obj_SEXP) { @@ -6315,6 +6337,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// SCM___FF16r__FF16_Env__aux__get +Rcpp::List SCM___FF16r__FF16_Env__aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_SCM___FF16r__FF16_Env__aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(SCM___FF16r__FF16_Env__aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // SCM___FF16r__FF16_Env__use_ode_times__get bool SCM___FF16r__FF16_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_SCM___FF16r__FF16_Env__use_ode_times__get(SEXP obj_SEXP) { @@ -6537,6 +6570,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// SCM___K93__K93_Env__aux__get +Rcpp::List SCM___K93__K93_Env__aux__get(plant::RcppR6::RcppR6 > obj_); +RcppExport SEXP _plant_SCM___K93__K93_Env__aux__get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 > >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(SCM___K93__K93_Env__aux__get(obj_)); + return rcpp_result_gen; +END_RCPP +} // SCM___K93__K93_Env__use_ode_times__get bool SCM___K93__K93_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_); RcppExport SEXP _plant_SCM___K93__K93_Env__use_ode_times__get(SEXP obj_SEXP) { @@ -10232,6 +10276,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_SCM___FF16__FF16_Env__cohort_schedule__set", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__cohort_schedule__set, 2}, {"_plant_SCM___FF16__FF16_Env__ode_times__get", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__ode_times__get, 1}, {"_plant_SCM___FF16__FF16_Env__state__get", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__state__get, 1}, + {"_plant_SCM___FF16__FF16_Env__aux__get", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__aux__get, 1}, {"_plant_SCM___FF16__FF16_Env__use_ode_times__get", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__use_ode_times__get, 1}, {"_plant_SCM___FF16__FF16_Env__use_ode_times__set", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__use_ode_times__set, 2}, {"_plant_SCM___FF16__FF16_Env__net_reproduction_ratio_errors__get", (DL_FUNC) &_plant_SCM___FF16__FF16_Env__net_reproduction_ratio_errors__get, 1}, @@ -10252,6 +10297,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_SCM___FF16w__FF16_Env__cohort_schedule__set", (DL_FUNC) &_plant_SCM___FF16w__FF16_Env__cohort_schedule__set, 2}, {"_plant_SCM___FF16w__FF16_Env__ode_times__get", (DL_FUNC) &_plant_SCM___FF16w__FF16_Env__ode_times__get, 1}, {"_plant_SCM___FF16w__FF16_Env__state__get", (DL_FUNC) &_plant_SCM___FF16w__FF16_Env__state__get, 1}, + {"_plant_SCM___FF16w__FF16_Env__aux__get", (DL_FUNC) &_plant_SCM___FF16w__FF16_Env__aux__get, 1}, {"_plant_SCM___FF16w__FF16_Env__use_ode_times__get", (DL_FUNC) &_plant_SCM___FF16w__FF16_Env__use_ode_times__get, 1}, {"_plant_SCM___FF16w__FF16_Env__use_ode_times__set", (DL_FUNC) &_plant_SCM___FF16w__FF16_Env__use_ode_times__set, 2}, {"_plant_SCM___FF16w__FF16_Env__net_reproduction_ratio_errors__get", (DL_FUNC) &_plant_SCM___FF16w__FF16_Env__net_reproduction_ratio_errors__get, 1}, @@ -10272,6 +10318,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_SCM___FF16r__FF16_Env__cohort_schedule__set", (DL_FUNC) &_plant_SCM___FF16r__FF16_Env__cohort_schedule__set, 2}, {"_plant_SCM___FF16r__FF16_Env__ode_times__get", (DL_FUNC) &_plant_SCM___FF16r__FF16_Env__ode_times__get, 1}, {"_plant_SCM___FF16r__FF16_Env__state__get", (DL_FUNC) &_plant_SCM___FF16r__FF16_Env__state__get, 1}, + {"_plant_SCM___FF16r__FF16_Env__aux__get", (DL_FUNC) &_plant_SCM___FF16r__FF16_Env__aux__get, 1}, {"_plant_SCM___FF16r__FF16_Env__use_ode_times__get", (DL_FUNC) &_plant_SCM___FF16r__FF16_Env__use_ode_times__get, 1}, {"_plant_SCM___FF16r__FF16_Env__use_ode_times__set", (DL_FUNC) &_plant_SCM___FF16r__FF16_Env__use_ode_times__set, 2}, {"_plant_SCM___FF16r__FF16_Env__net_reproduction_ratio_errors__get", (DL_FUNC) &_plant_SCM___FF16r__FF16_Env__net_reproduction_ratio_errors__get, 1}, @@ -10292,6 +10339,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_SCM___K93__K93_Env__cohort_schedule__set", (DL_FUNC) &_plant_SCM___K93__K93_Env__cohort_schedule__set, 2}, {"_plant_SCM___K93__K93_Env__ode_times__get", (DL_FUNC) &_plant_SCM___K93__K93_Env__ode_times__get, 1}, {"_plant_SCM___K93__K93_Env__state__get", (DL_FUNC) &_plant_SCM___K93__K93_Env__state__get, 1}, + {"_plant_SCM___K93__K93_Env__aux__get", (DL_FUNC) &_plant_SCM___K93__K93_Env__aux__get, 1}, {"_plant_SCM___K93__K93_Env__use_ode_times__get", (DL_FUNC) &_plant_SCM___K93__K93_Env__use_ode_times__get, 1}, {"_plant_SCM___K93__K93_Env__use_ode_times__set", (DL_FUNC) &_plant_SCM___K93__K93_Env__use_ode_times__set, 2}, {"_plant_SCM___K93__K93_Env__net_reproduction_ratio_errors__get", (DL_FUNC) &_plant_SCM___K93__K93_Env__net_reproduction_ratio_errors__get, 1}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index 48f4a802f..41566f8ec 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -2373,6 +2373,11 @@ Rcpp::List SCM___FF16__FF16_Env__state__get(plant::RcppR6::RcppR6 > obj_) { + return plant::get_aux(*obj_); +} + // [[Rcpp::export]] bool SCM___FF16__FF16_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_) { return obj_->r_use_ode_times(); @@ -2465,6 +2470,11 @@ Rcpp::List SCM___FF16w__FF16_Env__state__get(plant::RcppR6::RcppR6 > obj_) { + return plant::get_aux(*obj_); +} + // [[Rcpp::export]] bool SCM___FF16w__FF16_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_) { return obj_->r_use_ode_times(); @@ -2557,6 +2567,11 @@ Rcpp::List SCM___FF16r__FF16_Env__state__get(plant::RcppR6::RcppR6 > obj_) { + return plant::get_aux(*obj_); +} + // [[Rcpp::export]] bool SCM___FF16r__FF16_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_) { return obj_->r_use_ode_times(); @@ -2649,6 +2664,11 @@ Rcpp::List SCM___K93__K93_Env__state__get(plant::RcppR6::RcppR6 > obj_) { + return plant::get_aux(*obj_); +} + // [[Rcpp::export]] bool SCM___K93__K93_Env__use_ode_times__get(plant::RcppR6::RcppR6 > obj_) { return obj_->r_use_ode_times(); diff --git a/tests/testthat/test-scm-support.R b/tests/testthat/test-scm-support.R index b490b10c7..6611f646e 100644 --- a/tests/testthat/test-scm-support.R +++ b/tests/testthat/test-scm-support.R @@ -64,3 +64,18 @@ test_that("expand_parameters", { p1$max_patch_lifetime <- 100 expect_silent(p2 <- expand_parameters(trait_matrix(0.2, "lma"), p1, mutant=FALSE)) }) + +test_that("collect_auxillary_variables", { + env <- make_environment("FF16") + ctrl <- scm_base_control() + p0 <- scm_base_parameters("FF16") + hyperpar <- make_FF16_hyperpar() + p0$disturbance_mean_interval <- 30.0 + p1 <- expand_parameters(trait_matrix(0.08, "lma"), p0, hyperpar, FALSE) + + res <- run_scm_collect(p1, env, ctrl, collect_auxillary_variables=TRUE) + state <- res$species[[1]] + expect_equal(nrow(state), 9) + expect_equal(rownames(state)[8:9], c("competition_effect", "net_mass_production_dt")) + expect_equal(as.numeric(state[8:9, 142, 141]), c(0.0007211209, 0.0001707292)) +}) From 98433deedf5d2a1237a7305638dd46f2cb407c96 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Mon, 8 Nov 2021 17:35:27 +1100 Subject: [PATCH 04/22] Fixed all typos of 'auxillary' to 'auxiliary' --- R/RcppR6.R | 2 +- R/scm_support.R | 6 +++--- inst/RcppR6_classes.yml | 6 +++--- inst/include/plant/RcppR6_post.hpp | 18 +++++++++--------- inst/include/plant/models/ff16_strategy.h | 2 +- inst/include/plant/patch.h | 2 +- inst/include/plant/strategy.h | 2 +- src/ff16_strategy.cpp | 6 +++--- src/ff16r_strategy.cpp | 2 +- src/ff16w_strategy.cpp | 2 +- tests/testthat/test-scm-support.R | 4 ++-- .../test-strategy-ff16-reference-comparison.R | 4 ++-- tests/testthat/test-strategy-ff16.R | 8 ++++---- tests/testthat/test-strategy-ff16r.R | 8 ++++---- tests/testthat/test-strategy-k93.R | 2 +- 15 files changed, 37 insertions(+), 37 deletions(-) diff --git a/R/RcppR6.R b/R/RcppR6.R index 60bf49c77..67d57ccdc 100644 --- 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: edd85bd671370b151c2bba310b3f12e3 +## Hash: f6b511a10ca3bb03a599e71a6fa83092 ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class diff --git a/R/scm_support.R b/R/scm_support.R index e56ffe2a2..a9efcf951 100644 --- a/R/scm_support.R +++ b/R/scm_support.R @@ -100,13 +100,13 @@ run_scm <- function(p, env = make_environment(parameters = p), ##' @param p A \code{Parameters} object ##' @param env Environment object (defaults to FF16_Environment) ##' @param ctrl Control object -##' @param collect_auxillary_variables Return additional strategy variables (eg +##' @param collect_auxiliary_variables Return additional strategy variables (eg ##' competition_effect) ##' @author Rich FitzJohn ##' @export run_scm_collect <- function(p, env = make_environment(parameters = p), ctrl = scm_base_control(), - collect_auxillary_variables=FALSE) { + collect_auxiliary_variables=FALSE) { collect_default <- function(scm) { scm$state } @@ -116,7 +116,7 @@ run_scm_collect <- function(p, env = make_environment(parameters = p), ret$species <- mapply(rbind, ret$species, aux$species, SIMPLIFY=FALSE) ret } - collect <- if (collect_auxillary_variables) collect_aux else collect_default + collect <- if (collect_auxiliary_variables) collect_aux else collect_default types <- extract_RcppR6_template_types(p, "Parameters") scm <- do.call('SCM', types)(p, env, ctrl) diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index cc7652e7d..815517dbe 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -859,7 +859,7 @@ FF16_Strategy: - a_dG2: double - k_I: double - control: "plant::Control" - - collect_all_auxillary: bool + - collect_all_auxiliary: bool FF16_Environment: name_cpp: "plant::FF16_Environment" @@ -941,7 +941,7 @@ FF16r_Strategy: - a_dG2: double - k_I: double - control: "plant::Control" - - collect_all_auxillary: bool + - collect_all_auxiliary: bool # The following strategy was built from FF16r on Wed Aug 12 15:33:08 2020 K93_Strategy: @@ -1030,4 +1030,4 @@ FF16w_Strategy: - a_dG2: double - k_I: double - control: "plant::Control" - - collect_all_auxillary: bool + - collect_all_auxiliary: bool diff --git a/inst/include/plant/RcppR6_post.hpp b/inst/include/plant/RcppR6_post.hpp index 170d7946b..b6bfd99bb 100644 --- a/inst/include/plant/RcppR6_post.hpp +++ b/inst/include/plant/RcppR6_post.hpp @@ -995,7 +995,7 @@ template <> inline SEXP wrap(const plant::FF16_Strategy& x) { ret["a_dG2"] = Rcpp::wrap(x.a_dG2); ret["k_I"] = Rcpp::wrap(x.k_I); ret["control"] = Rcpp::wrap(x.control); - ret["collect_all_auxillary"] = Rcpp::wrap(x.collect_all_auxillary); + ret["collect_all_auxiliary"] = Rcpp::wrap(x.collect_all_auxiliary); ret.attr("class") = "FF16_Strategy"; return ret; } @@ -1072,8 +1072,8 @@ template <> inline plant::FF16_Strategy as(SEXP x) { ret.k_I = Rcpp::as(xl["k_I"]); // ret.control = Rcpp::as(xl["control"]); ret.control = Rcpp::as(xl["control"]); - // ret.collect_all_auxillary = Rcpp::as(xl["collect_all_auxillary"]); - ret.collect_all_auxillary = Rcpp::as(xl["collect_all_auxillary"]); + // ret.collect_all_auxiliary = Rcpp::as(xl["collect_all_auxiliary"]); + ret.collect_all_auxiliary = Rcpp::as(xl["collect_all_auxiliary"]); return ret; } template <> inline SEXP wrap(const plant::FF16_Environment& x) { @@ -1116,7 +1116,7 @@ template <> inline SEXP wrap(const plant::FF16r_Strategy& x) { ret["a_dG2"] = Rcpp::wrap(x.a_dG2); ret["k_I"] = Rcpp::wrap(x.k_I); ret["control"] = Rcpp::wrap(x.control); - ret["collect_all_auxillary"] = Rcpp::wrap(x.collect_all_auxillary); + ret["collect_all_auxiliary"] = Rcpp::wrap(x.collect_all_auxiliary); ret.attr("class") = "FF16r_Strategy"; return ret; } @@ -1193,8 +1193,8 @@ template <> inline plant::FF16r_Strategy as(SEXP x) { ret.k_I = Rcpp::as(xl["k_I"]); // ret.control = Rcpp::as(xl["control"]); ret.control = Rcpp::as(xl["control"]); - // ret.collect_all_auxillary = Rcpp::as(xl["collect_all_auxillary"]); - ret.collect_all_auxillary = Rcpp::as(xl["collect_all_auxillary"]); + // ret.collect_all_auxiliary = Rcpp::as(xl["collect_all_auxiliary"]); + ret.collect_all_auxiliary = Rcpp::as(xl["collect_all_auxiliary"]); return ret; } template <> inline SEXP wrap(const plant::K93_Strategy& x) { @@ -1289,7 +1289,7 @@ template <> inline SEXP wrap(const plant::FF16w_Strategy& x) { ret["a_dG2"] = Rcpp::wrap(x.a_dG2); ret["k_I"] = Rcpp::wrap(x.k_I); ret["control"] = Rcpp::wrap(x.control); - ret["collect_all_auxillary"] = Rcpp::wrap(x.collect_all_auxillary); + ret["collect_all_auxiliary"] = Rcpp::wrap(x.collect_all_auxiliary); ret.attr("class") = "FF16w_Strategy"; return ret; } @@ -1366,8 +1366,8 @@ template <> inline plant::FF16w_Strategy as(SEXP x) { ret.k_I = Rcpp::as(xl["k_I"]); // ret.control = Rcpp::as(xl["control"]); ret.control = Rcpp::as(xl["control"]); - // ret.collect_all_auxillary = Rcpp::as(xl["collect_all_auxillary"]); - ret.collect_all_auxillary = Rcpp::as(xl["collect_all_auxillary"]); + // ret.collect_all_auxiliary = Rcpp::as(xl["collect_all_auxiliary"]); + ret.collect_all_auxiliary = Rcpp::as(xl["collect_all_auxiliary"]); return ret; } } diff --git a/inst/include/plant/models/ff16_strategy.h b/inst/include/plant/models/ff16_strategy.h index 1b957fc96..6b61080b3 100644 --- a/inst/include/plant/models/ff16_strategy.h +++ b/inst/include/plant/models/ff16_strategy.h @@ -36,7 +36,7 @@ class FF16_Strategy: public Strategy { "net_mass_production_dt" }); // add the associated computation to compute_rates and compute there - if (collect_all_auxillary) { + if (collect_all_auxiliary) { ret.push_back("area_sapwood"); } return ret; diff --git a/inst/include/plant/patch.h b/inst/include/plant/patch.h index b8c728a19..6396cade7 100644 --- a/inst/include/plant/patch.h +++ b/inst/include/plant/patch.h @@ -244,7 +244,7 @@ size_t Patch::ode_size() const { template size_t Patch::aux_size() const { - // no use for auxillary environment variables (yet) + // no use for auxiliary environment variables (yet) return ode::aux_size(species.begin(), species.end());// + environment.ode_size(); } diff --git a/inst/include/plant/strategy.h b/inst/include/plant/strategy.h index 98e5f45b3..d1f72c912 100644 --- a/inst/include/plant/strategy.h +++ b/inst/include/plant/strategy.h @@ -35,7 +35,7 @@ class Strategy { std::map aux_index; - bool collect_all_auxillary; + bool collect_all_auxiliary; void refresh_indices(); diff --git a/src/ff16_strategy.cpp b/src/ff16_strategy.cpp index f23750312..1409d9a19 100644 --- a/src/ff16_strategy.cpp +++ b/src/ff16_strategy.cpp @@ -7,7 +7,7 @@ namespace plant { // before physiology, before compound things?) // TODO: Consider moving to activating as an initialisation list? FF16_Strategy::FF16_Strategy() { - collect_all_auxillary = false; + collect_all_auxiliary = false; // build the string state/aux name to index map refresh_indices(); name = "FF16"; @@ -86,7 +86,7 @@ double FF16_Strategy::mass_above_ground(double mass_leaf, double mass_bark, return mass_leaf + mass_bark + mass_sapwood + mass_root; } -// for updating auxillary state +// for updating auxiliary state void FF16_Strategy::update_dependent_aux(const int index, Internals& vars) { if (index == HEIGHT_INDEX) { double height = vars.state(HEIGHT_INDEX); @@ -126,7 +126,7 @@ void FF16_Strategy::compute_rates(const FF16_Environment& environment, const double mass_sapwood_ = mass_sapwood(area_sapwood_, height); vars.set_rate(state_index.at("mass_heartwood"), mass_heartwood_dt(mass_sapwood_)); - if (collect_all_auxillary) { + if (collect_all_auxiliary) { vars.set_aux(aux_index.at("area_sapwood"), area_sapwood_); } } else { diff --git a/src/ff16r_strategy.cpp b/src/ff16r_strategy.cpp index 4518c1b9f..a951a2b8e 100644 --- a/src/ff16r_strategy.cpp +++ b/src/ff16r_strategy.cpp @@ -11,7 +11,7 @@ FF16r_Strategy::FF16r_Strategy() { // height above hmat at which allocation to reproduction is half its max value a_f2 = 2; // [dimensionless] - collect_all_auxillary = false; + collect_all_auxiliary = false; // build the string state/aux name to index map refresh_indices(); name = "FF16r"; diff --git a/src/ff16w_strategy.cpp b/src/ff16w_strategy.cpp index d1c3f28ef..c05b8c2ae 100644 --- a/src/ff16w_strategy.cpp +++ b/src/ff16w_strategy.cpp @@ -3,7 +3,7 @@ namespace plant { FF16w_Strategy::FF16w_Strategy() { - collect_all_auxillary = false; + collect_all_auxiliary = false; // build the string state/aux name to index map refresh_indices(); name = "FF16w"; diff --git a/tests/testthat/test-scm-support.R b/tests/testthat/test-scm-support.R index 6611f646e..dbaec21ef 100644 --- a/tests/testthat/test-scm-support.R +++ b/tests/testthat/test-scm-support.R @@ -65,7 +65,7 @@ test_that("expand_parameters", { expect_silent(p2 <- expand_parameters(trait_matrix(0.2, "lma"), p1, mutant=FALSE)) }) -test_that("collect_auxillary_variables", { +test_that("collect_auxiliary_variables", { env <- make_environment("FF16") ctrl <- scm_base_control() p0 <- scm_base_parameters("FF16") @@ -73,7 +73,7 @@ test_that("collect_auxillary_variables", { p0$disturbance_mean_interval <- 30.0 p1 <- expand_parameters(trait_matrix(0.08, "lma"), p0, hyperpar, FALSE) - res <- run_scm_collect(p1, env, ctrl, collect_auxillary_variables=TRUE) + res <- run_scm_collect(p1, env, ctrl, collect_auxiliary_variables=TRUE) state <- res$species[[1]] expect_equal(nrow(state), 9) expect_equal(rownames(state)[8:9], c("competition_effect", "net_mass_production_dt")) diff --git a/tests/testthat/test-strategy-ff16-reference-comparison.R b/tests/testthat/test-strategy-ff16-reference-comparison.R index 030ea05bd..7b22430b4 100644 --- a/tests/testthat/test-strategy-ff16-reference-comparison.R +++ b/tests/testthat/test-strategy-ff16-reference-comparison.R @@ -12,7 +12,7 @@ test_that("FF16_Strategy parameters agree with reference model", { expect_true(all(v %in% names(s))) ## And v.v., except for a few additions: - extra <- c("control", "S_D", "collect_all_auxillary") + extra <- c("control", "S_D", "collect_all_auxiliary") common <- setdiff(names(s), extra) expect_true(all(extra %in% names(s))) expect_true(all(common %in% names(cmp_pars))) @@ -42,7 +42,7 @@ test_that("Reference comparison", { p$set_state("height", h0) expect_identical(p$state("height"), h0) - # testing set auxillary state as well as area_leaf/competition_effect depends on height only + # testing set auxiliary state as well as area_leaf/competition_effect depends on height only expect_equal(p$aux("competition_effect"), cmp$LeafArea(h0)) # expect_equal(p$state("mass_leaf"), cmp$LeafMass(cmp$traits$lma, cmp$LeafArea(h0))) # expect_equal(p$state("mass_sapwood"), cmp$SapwoodMass(cmp$traits$rho, cmp$LeafArea(h0), h0)) diff --git a/tests/testthat/test-strategy-ff16.R b/tests/testthat/test-strategy-ff16.R index aa97bea8b..a63478f16 100644 --- a/tests/testthat/test-strategy-ff16.R +++ b/tests/testthat/test-strategy-ff16.R @@ -34,7 +34,7 @@ test_that("Defaults", { theta = 1.0/4669, k_I = 0.5, control = Control(), - collect_all_auxillary = FALSE) + collect_all_auxiliary = FALSE) keys <- sort(names(expected)) @@ -45,7 +45,7 @@ test_that("Defaults", { expect_identical(unclass(s)[keys], expected[keys]) }) -test_that("FF16 collect_all_auxillary option", { +test_that("FF16 collect_all_auxiliary option", { s <- FF16_Strategy() p <- FF16_Individual(s) @@ -56,8 +56,8 @@ test_that("FF16 collect_all_auxillary option", { "net_mass_production_dt" )) - s <- FF16_Strategy(collect_all_auxillary=TRUE) - expect_true(s$collect_all_auxillary) + s <- FF16_Strategy(collect_all_auxiliary=TRUE) + expect_true(s$collect_all_auxiliary) p <- FF16_Individual(s) expect_equal(p$aux_size, 3) expect_equal(length(p$internals$auxs), 3) diff --git a/tests/testthat/test-strategy-ff16r.R b/tests/testthat/test-strategy-ff16r.R index b60444991..952b20f4b 100644 --- a/tests/testthat/test-strategy-ff16r.R +++ b/tests/testthat/test-strategy-ff16r.R @@ -35,7 +35,7 @@ test_that("Defaults", { theta = 1.0/4669, k_I = 0.5, control = Control(), - collect_all_auxillary = FALSE) + collect_all_auxiliary = FALSE) keys <- sort(names(expected)) @@ -46,7 +46,7 @@ test_that("Defaults", { expect_identical(unclass(s)[keys], expected[keys]) }) -test_that("FF16r collect_all_auxillary option", { +test_that("FF16r collect_all_auxiliary option", { s <- FF16r_Strategy() p <- FF16r_Individual(s) @@ -57,8 +57,8 @@ test_that("FF16r collect_all_auxillary option", { "net_mass_production_dt" )) - s <- FF16r_Strategy(collect_all_auxillary=TRUE) - expect_true(s$collect_all_auxillary) + s <- FF16r_Strategy(collect_all_auxiliary=TRUE) + expect_true(s$collect_all_auxiliary) p <- FF16r_Individual(s) expect_equal(p$aux_size, 3) expect_equal(length(p$internals$auxs), 3) diff --git a/tests/testthat/test-strategy-k93.R b/tests/testthat/test-strategy-k93.R index c4a5c65bc..73a567c62 100644 --- a/tests/testthat/test-strategy-k93.R +++ b/tests/testthat/test-strategy-k93.R @@ -26,7 +26,7 @@ test_that("Defaults", { expect_identical(unclass(s)[keys], expected[keys]) }) -test_that("K93 collect_all_auxillary option", { +test_that("K93 collect_all_auxiliary option", { s <- K93_Strategy() p <- K93_Individual(s) From 1443006fe0794e97884e8d7f12291459c2dca2df Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Mon, 15 Nov 2021 16:37:16 +1100 Subject: [PATCH 05/22] Spline added to FF16 environment with basic API exposed + tests --- R/RcppExports.R | 32 ++++++ R/RcppR6.R | 26 ++++- inst/RcppR6_classes.yml | 24 +++++ inst/include/plant/interpolator.h | 3 + inst/include/plant/models/ff16_environment.h | 52 +++++++++- src/RcppExports.cpp | 100 +++++++++++++++++++ src/RcppR6.cpp | 32 ++++++ src/interpolator.cpp | 24 ++++- tests/testthat/test-environment.R | 32 ++++++ tests/testthat/test-interpolator.R | 1 + 10 files changed, 321 insertions(+), 5 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 93520faaa..479aea22e 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -3293,6 +3293,38 @@ 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__rainfall_init <- function(obj_, x, y) { + invisible(.Call('_plant_FF16_Environment__rainfall_init', PACKAGE = 'plant', obj_, x, y)) +} + +FF16_Environment__rainfall_recompute <- function(obj_) { + invisible(.Call('_plant_FF16_Environment__rainfall_recompute', PACKAGE = 'plant', obj_)) +} + +FF16_Environment__rainfall_eval <- function(obj_, u) { + .Call('_plant_FF16_Environment__rainfall_eval', PACKAGE = 'plant', obj_, u) +} + +FF16_Environment__rainfall_eval_range <- function(obj_, u) { + .Call('_plant_FF16_Environment__rainfall_eval_range', PACKAGE = 'plant', obj_, u) +} + +FF16_Environment__rainfall_add_point <- function(obj_, x, y) { + invisible(.Call('_plant_FF16_Environment__rainfall_add_point', PACKAGE = 'plant', obj_, x, y)) +} + +FF16_Environment__rainfall_add_point_sorted <- function(obj_, x, y) { + invisible(.Call('_plant_FF16_Environment__rainfall_add_point_sorted', PACKAGE = 'plant', obj_, x, y)) +} + +FF16_Environment__rainfall_add_points <- function(obj_, x, y) { + invisible(.Call('_plant_FF16_Environment__rainfall_add_points', PACKAGE = 'plant', obj_, x, y)) +} + +FF16_Environment__rainfall_clear_points <- function(obj_) { + invisible(.Call('_plant_FF16_Environment__rainfall_clear_points', PACKAGE = 'plant', obj_)) +} + FF16_Environment__canopy_openness <- function(obj_, height) { .Call('_plant_FF16_Environment__canopy_openness', PACKAGE = 'plant', obj_, height) } diff --git a/R/RcppR6.R b/R/RcppR6.R index 67d57ccdc..873646c0f 100644 --- 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: 6f2a9dac8f4676a71b9571b0be4de76a ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -4814,6 +4814,30 @@ StochasticPatchRunner <- function(T, E) { initialize = function(ptr) { self$.ptr <- ptr }, + rainfall_init = function(x, y) { + FF16_Environment__rainfall_init(self, x, y) + }, + rainfall_recompute = function() { + FF16_Environment__rainfall_recompute(self) + }, + rainfall_eval = function(u) { + FF16_Environment__rainfall_eval(self, u) + }, + rainfall_eval_range = function(u) { + FF16_Environment__rainfall_eval_range(self, u) + }, + rainfall_add_point = function(x, y) { + FF16_Environment__rainfall_add_point(self, x, y) + }, + rainfall_add_point_sorted = function(x, y) { + FF16_Environment__rainfall_add_point_sorted(self, x, y) + }, + rainfall_add_points = function(x, y) { + FF16_Environment__rainfall_add_points(self, x, y) + }, + rainfall_clear_points = function() { + FF16_Environment__rainfall_clear_points(self) + }, canopy_openness = function(height) { FF16_Environment__canopy_openness(self, height) }, diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index 815517dbe..ac15146b8 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -871,6 +871,30 @@ FF16_Environment: FF16_Environment object @export methods: + rainfall_init: + args: [x: "std::vector", y: "std::vector"] + return_type: void + rainfall_recompute: + args: [] + return_type: void + rainfall_eval: + args: [u: double] + return_type: double + rainfall_eval_range: + args: [u: "std::vector"] + return_type: "std::vector" + rainfall_add_point: + args: [x: double, y: double] + return_type: void + rainfall_add_point_sorted: + args: [x: double, y: double] + return_type: void + rainfall_add_points: + args: [x: "std::vector", y: "std::vector"] + return_type: void + rainfall_clear_points: + args: [] + return_type: void canopy_openness: args: [height: double] return_type: double diff --git a/inst/include/plant/interpolator.h b/inst/include/plant/interpolator.h index 4fca2f521..56d2b4549 100644 --- a/inst/include/plant/interpolator.h +++ b/inst/include/plant/interpolator.h @@ -16,9 +16,11 @@ 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; @@ -26,6 +28,7 @@ class Interpolator { std::vector get_x() const; std::vector get_y() const; + // * R interface SEXP r_get_xy() const; diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h index 0b2b8cd0f..7dfae7f7f 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -4,6 +4,7 @@ #include #include +#include using namespace Rcpp; @@ -22,11 +23,58 @@ class FF16_Environment : public Environment { soil_infiltration_rate = 0.0; vars = Internals(soil_number_of_depths); set_soil_water_state(std::vector(soil_number_of_depths, 0.0)); + + rainfall_spline = interpolator::Interpolator(); }; + // first spline computation + void rainfall_init(std::vector const& x, std::vector const& y) { + rainfall_spline.init(x, y); // set x, y points and compute spline + } + + // if any points are added to the spline, need to recompute spline + void rainfall_recompute() { + rainfall_spline.initialise(); + } + + double rainfall_eval(double u) const { + return rainfall_spline.eval(u); + } + + std::vector rainfall_eval_range(std::vector const& u) const { + return rainfall_spline.r_eval(u); + } + + /* if xi is not larger than every other x in the spline, recomputing the spline + will throw an error and will need to clear_points() and reconstruct spline */ + void rainfall_add_point(double xi, double yi) { + rainfall_spline.add_point(xi, yi); + } + + // safe add point + void rainfall_add_point_sorted(double xi, double yi) { + rainfall_spline.add_point_sorted(xi, yi); + } + + /* not sure if these two are necessary for the rainfall model. could also move these + to be general interpolator.h functions? */ + void rainfall_add_points(std::vector const& x, std::vector const& y) { + util::check_length(x.size(), y.size()); + for (size_t i = 0; i < x.size(); ++i) { + rainfall_add_point(x[i], y[i]); // refactor? + } + } + + // rainfall_add_points_sorted ? + + void rainfall_clear_points() { + rainfall_spline.clear(); + } + // Light interface bool canopy_rescale_usually; + // private? Canopy canopy; // Should this be here or in canopy? @@ -92,7 +140,9 @@ class FF16_Environment : public Environment { void clear_environment() { canopy.clear(); } - + +private: + interpolator::Interpolator rainfall_spline; }; inline Rcpp::NumericMatrix get_state(const FF16_Environment environment) { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 2ce5f7266..602fb875a 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -9235,6 +9235,98 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// FF16_Environment__rainfall_init +void FF16_Environment__rainfall_init(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y); +RcppExport SEXP _plant_FF16_Environment__rainfall_init(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::vector >::type x(xSEXP); + Rcpp::traits::input_parameter< std::vector >::type y(ySEXP); + FF16_Environment__rainfall_init(obj_, x, y); + return R_NilValue; +END_RCPP +} +// FF16_Environment__rainfall_recompute +void FF16_Environment__rainfall_recompute(plant::RcppR6::RcppR6 obj_); +RcppExport SEXP _plant_FF16_Environment__rainfall_recompute(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + FF16_Environment__rainfall_recompute(obj_); + return R_NilValue; +END_RCPP +} +// FF16_Environment__rainfall_eval +double FF16_Environment__rainfall_eval(plant::RcppR6::RcppR6 obj_, double u); +RcppExport SEXP _plant_FF16_Environment__rainfall_eval(SEXP obj_SEXP, SEXP uSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type u(uSEXP); + rcpp_result_gen = Rcpp::wrap(FF16_Environment__rainfall_eval(obj_, u)); + return rcpp_result_gen; +END_RCPP +} +// FF16_Environment__rainfall_eval_range +std::vector FF16_Environment__rainfall_eval_range(plant::RcppR6::RcppR6 obj_, std::vector u); +RcppExport SEXP _plant_FF16_Environment__rainfall_eval_range(SEXP obj_SEXP, SEXP uSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::vector >::type u(uSEXP); + rcpp_result_gen = Rcpp::wrap(FF16_Environment__rainfall_eval_range(obj_, u)); + return rcpp_result_gen; +END_RCPP +} +// FF16_Environment__rainfall_add_point +void FF16_Environment__rainfall_add_point(plant::RcppR6::RcppR6 obj_, double x, double y); +RcppExport SEXP _plant_FF16_Environment__rainfall_add_point(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type x(xSEXP); + Rcpp::traits::input_parameter< double >::type y(ySEXP); + FF16_Environment__rainfall_add_point(obj_, x, y); + return R_NilValue; +END_RCPP +} +// FF16_Environment__rainfall_add_point_sorted +void FF16_Environment__rainfall_add_point_sorted(plant::RcppR6::RcppR6 obj_, double x, double y); +RcppExport SEXP _plant_FF16_Environment__rainfall_add_point_sorted(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type x(xSEXP); + Rcpp::traits::input_parameter< double >::type y(ySEXP); + FF16_Environment__rainfall_add_point_sorted(obj_, x, y); + return R_NilValue; +END_RCPP +} +// FF16_Environment__rainfall_add_points +void FF16_Environment__rainfall_add_points(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y); +RcppExport SEXP _plant_FF16_Environment__rainfall_add_points(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::vector >::type x(xSEXP); + Rcpp::traits::input_parameter< std::vector >::type y(ySEXP); + FF16_Environment__rainfall_add_points(obj_, x, y); + return R_NilValue; +END_RCPP +} +// FF16_Environment__rainfall_clear_points +void FF16_Environment__rainfall_clear_points(plant::RcppR6::RcppR6 obj_); +RcppExport SEXP _plant_FF16_Environment__rainfall_clear_points(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + FF16_Environment__rainfall_clear_points(obj_); + return R_NilValue; +END_RCPP +} // FF16_Environment__canopy_openness double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6 obj_, double height); RcppExport SEXP _plant_FF16_Environment__canopy_openness(SEXP obj_SEXP, SEXP heightSEXP) { @@ -10578,6 +10670,14 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_Weibull_Disturbance_Regime__icdf", (DL_FUNC) &_plant_Weibull_Disturbance_Regime__icdf, 2}, {"_plant_FF16_Strategy__ctor", (DL_FUNC) &_plant_FF16_Strategy__ctor, 0}, {"_plant_FF16_Environment__ctor", (DL_FUNC) &_plant_FF16_Environment__ctor, 2}, + {"_plant_FF16_Environment__rainfall_init", (DL_FUNC) &_plant_FF16_Environment__rainfall_init, 3}, + {"_plant_FF16_Environment__rainfall_recompute", (DL_FUNC) &_plant_FF16_Environment__rainfall_recompute, 1}, + {"_plant_FF16_Environment__rainfall_eval", (DL_FUNC) &_plant_FF16_Environment__rainfall_eval, 2}, + {"_plant_FF16_Environment__rainfall_eval_range", (DL_FUNC) &_plant_FF16_Environment__rainfall_eval_range, 2}, + {"_plant_FF16_Environment__rainfall_add_point", (DL_FUNC) &_plant_FF16_Environment__rainfall_add_point, 3}, + {"_plant_FF16_Environment__rainfall_add_point_sorted", (DL_FUNC) &_plant_FF16_Environment__rainfall_add_point_sorted, 3}, + {"_plant_FF16_Environment__rainfall_add_points", (DL_FUNC) &_plant_FF16_Environment__rainfall_add_points, 3}, + {"_plant_FF16_Environment__rainfall_clear_points", (DL_FUNC) &_plant_FF16_Environment__rainfall_clear_points, 1}, {"_plant_FF16_Environment__canopy_openness", (DL_FUNC) &_plant_FF16_Environment__canopy_openness, 2}, {"_plant_FF16_Environment__clear", (DL_FUNC) &_plant_FF16_Environment__clear, 1}, {"_plant_FF16_Environment__set_fixed_environment", (DL_FUNC) &_plant_FF16_Environment__set_fixed_environment, 3}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index 41566f8ec..a7e1660cd 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -3752,6 +3752,38 @@ plant::FF16_Environment FF16_Environment__ctor(bool canopy_rescale_usually, int return plant::FF16_Environment(canopy_rescale_usually, soil_number_of_depths); } // [[Rcpp::export]] +void FF16_Environment__rainfall_init(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y) { + obj_->rainfall_init(x, y); +} +// [[Rcpp::export]] +void FF16_Environment__rainfall_recompute(plant::RcppR6::RcppR6 obj_) { + obj_->rainfall_recompute(); +} +// [[Rcpp::export]] +double FF16_Environment__rainfall_eval(plant::RcppR6::RcppR6 obj_, double u) { + return obj_->rainfall_eval(u); +} +// [[Rcpp::export]] +std::vector FF16_Environment__rainfall_eval_range(plant::RcppR6::RcppR6 obj_, std::vector u) { + return obj_->rainfall_eval_range(u); +} +// [[Rcpp::export]] +void FF16_Environment__rainfall_add_point(plant::RcppR6::RcppR6 obj_, double x, double y) { + obj_->rainfall_add_point(x, y); +} +// [[Rcpp::export]] +void FF16_Environment__rainfall_add_point_sorted(plant::RcppR6::RcppR6 obj_, double x, double y) { + obj_->rainfall_add_point_sorted(x, y); +} +// [[Rcpp::export]] +void FF16_Environment__rainfall_add_points(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y) { + obj_->rainfall_add_points(x, y); +} +// [[Rcpp::export]] +void FF16_Environment__rainfall_clear_points(plant::RcppR6::RcppR6 obj_) { + obj_->rainfall_clear_points(); +} +// [[Rcpp::export]] double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6 obj_, double height) { return obj_->canopy_openness(height); } diff --git a/src/interpolator.cpp b/src/interpolator.cpp index 6695a27cc..af3e0e568 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -21,6 +21,9 @@ void Interpolator::init(const std::vector& x_, // Compute the interpolated function from the points contained in 'x' // and 'y'. void Interpolator::initialise() { + if (not std::is_sorted(x.begin(), x.end())) { + util::stop("spline control points must be in non-descending order"); + } if (x.size() > 0) { tk_spline.set_points(x, y); active = true; @@ -34,6 +37,14 @@ void Interpolator::add_point(double xi, double yi) { y.push_back(yi); } +// adds point in sorted position (slower than above) +void Interpolator::add_point_sorted(double xi, double yi) { + auto x_upper = std::upper_bound(x.begin(), x.end(), xi); // find smallest number larger than xi + x.insert(x_upper, xi); // add xi below that number + auto y_upper = std::upper_bound(y.begin(), y.end(), yi); + y.insert(y_upper, yi); +} + // Remove all the contents, being ready to be refilled. void Interpolator::clear() { x.clear(); @@ -43,6 +54,12 @@ void Interpolator::clear() { // Compute the value of the interpolated function at point `x=u` double Interpolator::eval(double u) const { + check_active(); + return tk_spline(u); +} + +// faster version +double Interpolator::operator()(double u) const { return tk_spline(u); } @@ -84,9 +101,10 @@ std::vector Interpolator::r_eval(std::vector u) const { check_active(); const size_t n = u.size(); std::vector ret(n); - for (size_t i = 0; i < n; ++i) { - ret[i] = eval(u[i]); - } + // for (size_t i = 0; i < n; ++i) { + // ret[i] = operator()(u[i]); + // } + std::transform(u.begin(), u.end(), ret.begin(), [&](double x){return operator()(x);}); return ret; } diff --git a/tests/testthat/test-environment.R b/tests/testthat/test-environment.R index 4e68dfef4..358e94de4 100644 --- a/tests/testthat/test-environment.R +++ b/tests/testthat/test-environment.R @@ -41,3 +41,35 @@ for (x in names(strategy_types)) { expect_identical(sapply(hmid, env$canopy$canopy_interpolator$eval), sapply(hmid, interplator$eval)) }) } + +test_that("FF16 rainfall spline", { + context(sprintf("Rainfall-FF16",x)) + env <- make_environment("FF16") + + ## simple quadratic + x <- seq(-10, 10, 2) + y <- x^2 + env$rainfall_init(x, y) + + ## discrete points used for spline creation + expect_equal(env$rainfall_eval(2), 4) + expect_equal(env$rainfall_eval(-2), 4) + + ## interpolated points + expect_equal(env$rainfall_eval(3), 9, tolerance=1e-2) + expect_equal(env$rainfall_eval(-3), 9, tolerance=1e-2) + expect_equal(env$rainfall_eval(5.5), 30.25, tolerance=1e-2) + expect_equal(env$rainfall_eval(-5.5), 30.25, tolerance=1e-2) + + ## interpolated range of points + expect_equal(env$rainfall_eval_range(c(-7, 1, 7.8345)), c(49, 1, 61.37939025), tolerance=1e-2) + + ## clear points, set new points and recompute spline + env$rainfall_clear_points() + x <- seq(-10, 10, 1) + y <- 2*x + env$rainfall_add_points(x, y) + env$rainfall_recompute() + expect_equal(env$rainfall_eval(3), 6, tolerance=1e-2) + expect_equal(env$rainfall_eval(3.3), 6.6, tolerance=1e-2) +}) diff --git a/tests/testthat/test-interpolator.R b/tests/testthat/test-interpolator.R index 0a7a587ab..3a789fe33 100644 --- a/tests/testthat/test-interpolator.R +++ b/tests/testthat/test-interpolator.R @@ -26,6 +26,7 @@ test_that("Splines require sensible data", { s <- Interpolator() expect_error(s$init(c(1, 2, 3), 1), "Incorrect length input") expect_error(s$init(numeric(0), numeric(0)), "insufficient number of points") + expect_error(s$init(c(3, 2, 1), c(1, 1, 1)), "spline control points must be in non-descending order") }) test_that("Spline contains correct data", { From 701118933626f8283487a01c0f29329e87a5de57 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Wed, 24 Nov 2021 14:38:55 +1100 Subject: [PATCH 06/22] basic test for rainfall spline over run_scm --- inst/include/plant/models/ff16_environment.h | 3 +- tests/testthat/test-strategy-ff16w.R | 38 ++++++++++++++++++++ 2 files changed, 40 insertions(+), 1 deletion(-) diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h index 7dfae7f7f..fb49fd27b 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -5,6 +5,7 @@ #include #include #include +#include using namespace Rcpp; @@ -104,7 +105,7 @@ class FF16_Environment : public Environment { 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, rainfall_eval(time) / (i+1)); } } diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index 59c428ec1..0f327c96a 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -80,3 +80,41 @@ test_that("Basic run", { expect_equal(out$offspring_production, 16.88946, tolerance=1e-5) expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-5) }) + +test_that("Rainfall spline", { + # one species + p0 <- scm_base_parameters("FF16w") + p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16w_hyperpar,FALSE) + p1$birth_rate <- 20 + + env <- make_environment("FF16w", + soil_number_of_depths = 10, + soil_initial_state = rep(1, 10), + soil_infiltration_rate = 1) + + # init rainfall spline for env + x <- seq(0, 110, 1) + integrand <- function(x) {x^2} + y <- integrand(x) + env$rainfall_init(x, y) + + ctrl <- scm_base_control() + + out <- run_scm(p1, env, ctrl) + + expect_equal(out$patch$environment$ode_size, 10) + + # model is currently 1/z * int(0, t)(f(x)), where z is depth, + # x is time, f is spline, t is number of years + + expect_equal(out$patch$environment$soil$rates, + sapply(seq(1, 10), function(x) { (out$time^2)/x }), + tolerance = .01) + + expect_equal(out$patch$environment$soil$states, + sapply(seq(1, 10), function(i) { (1/i) * integrate(integrand, 0, out$time)$value }), + tolerance = 0.1) + + expect_equal(out$offspring_production, 16.88946, tolerance=1e-5) + expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-5) +}) \ No newline at end of file From da576f1dc11b1d4ced89738cc7529f533b7b82b0 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Wed, 24 Nov 2021 23:35:55 +1100 Subject: [PATCH 07/22] Default rainfall spline of y = 0 --- R/ff16w.R | 7 ++- tests/testthat/test-strategy-ff16w.R | 74 ++++++++++++++++------------ 2 files changed, 47 insertions(+), 34 deletions(-) diff --git a/R/ff16w.R b/R/ff16w.R index fc089ef79..69603281e 100644 --- a/R/ff16w.R +++ b/R/ff16w.R @@ -73,7 +73,9 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4, canopy_rescale_usually = TRUE, soil_number_of_depths = 1, soil_initial_state = 0.0, - soil_infiltration_rate = 0.0) { + soil_infiltration_rate = 0.0, + rainfall_x = seq(0, 9, 1), + rainfall_y = rep(0, 10)) { if(soil_number_of_depths < 1) stop("FF16w Environment must have at least one soil layer") @@ -85,7 +87,7 @@ 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") @@ -94,6 +96,7 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4, } e$set_soil_infiltration_rate(soil_infiltration_rate) + e$rainfall_init(rainfall_x, rainfall_y) return(e) } diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index 0f327c96a..208c4915a 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -17,11 +17,15 @@ test_that("FF16w Environment", { expect_equal(e$soil$rates, 0) # Make it rain - e$set_soil_infiltration_rate(10) + x <- seq(0, 9, 1) + y <- rep(1, 10) + #e$set_soil_infiltration_rate(10) + e$rainfall_init(x, y) expect_equal(e$soil$states, 0) e$compute_rates() - expect_equal(e$soil$rates, 10) + expect_equal(e$soil$rates, 1) + # not needed anymore? should we still be able to set the states manually? # Water logged e$set_soil_water_state(100) expect_equal(e$soil$states, 100) @@ -29,7 +33,9 @@ test_that("FF16w Environment", { # Check construction e <- make_environment("FF16w", soil_number_of_depths = 1, - soil_infiltration_rate = 10) + soil_infiltration_rate = 10, # not needed anymore? + rainfall_x = seq(0, 9, 1), + rainfall_y = rep(10, 10)) expect_equal(e$soil$states, 0) e$compute_rates() @@ -52,36 +58,36 @@ test_that("FF16w Environment", { }) -test_that("Basic run", { - # one species - p0 <- scm_base_parameters("FF16w") - p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16w_hyperpar,FALSE) - p1$birth_rate <- 20 - - env <- make_environment("FF16w", - soil_number_of_depths = 10, - soil_initial_state = rep(1, 10), - soil_infiltration_rate = 1) - - ctrl <- scm_base_control() - - out <- run_scm(p1, env, ctrl) - - expect_equal(out$patch$environment$ode_size, 10) - - expect_equal(out$patch$environment$soil$rates, - c(1.00, 0.50, 0.33, 0.25, 0.20, 0.17, 0.14, 0.13, 0.11, 0.10), - tolerance = .01) - - expect_equal(out$patch$environment$soil$states, - c(105, 52, 35, 26, 21, 17, 15, 13, 11, 10), - tolerance = 0.1) - - expect_equal(out$offspring_production, 16.88946, tolerance=1e-5) - expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-5) -}) +# test_that("Basic run", { +# # one species +# p0 <- scm_base_parameters("FF16w") +# p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16w_hyperpar,FALSE) +# p1$birth_rate <- 20 +# +# env <- make_environment("FF16w", +# soil_number_of_depths = 10, +# soil_initial_state = rep(1, 10), +# soil_infiltration_rate = 1) +# +# ctrl <- scm_base_control() +# +# out <- run_scm(p1, env, ctrl) +# +# expect_equal(out$patch$environment$ode_size, 10) +# +# expect_equal(out$patch$environment$soil$rates, +# c(1.00, 0.50, 0.33, 0.25, 0.20, 0.17, 0.14, 0.13, 0.11, 0.10), +# tolerance = .01) +# +# expect_equal(out$patch$environment$soil$states, +# c(105, 52, 35, 26, 21, 17, 15, 13, 11, 10), +# tolerance = 0.1) +# +# expect_equal(out$offspring_production, 16.88946, tolerance=1e-5) +# expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-5) +# }) -test_that("Rainfall spline", { +test_that("Rainfall spline basic run", { # one species p0 <- scm_base_parameters("FF16w") p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16w_hyperpar,FALSE) @@ -107,10 +113,14 @@ test_that("Rainfall spline", { # model is currently 1/z * int(0, t)(f(x)), where z is depth, # x is time, f is spline, t is number of years + # check the rates are correct, ie 105^2/depth expect_equal(out$patch$environment$soil$rates, sapply(seq(1, 10), function(x) { (out$time^2)/x }), tolerance = .01) + # check the states are correct + # model is currently 1/z * int(0, t)(f(x)), where z is depth, + # x is time, f is spline, t is number of years expect_equal(out$patch$environment$soil$states, sapply(seq(1, 10), function(i) { (1/i) * integrate(integrand, 0, out$time)$value }), tolerance = 0.1) From 96c55e65a9f108022b88c97dcc7ab6c2f6cbd574 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Thu, 25 Nov 2021 15:04:54 +1100 Subject: [PATCH 08/22] Removed depricated infiltration rate settings, updated documentation --- R/RcppExports.R | 4 --- R/RcppR6.R | 5 +-- R/ff16w.R | 8 ++--- inst/RcppR6_classes.yml | 3 -- inst/include/plant/models/ff16_environment.h | 8 ----- src/RcppExports.cpp | 12 ------- src/RcppR6.cpp | 4 --- tests/testthat/test-strategy-ff16w.R | 34 +------------------- 8 files changed, 5 insertions(+), 73 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 479aea22e..b2e42d627 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -3337,10 +3337,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 index 873646c0f..2d4a4e689 100644 --- 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: 6f2a9dac8f4676a71b9571b0be4de76a +## Hash: 2149ebb5397c92102688df8a4c991aa7 ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -4847,9 +4847,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 index 69603281e..2fd266900 100644 --- 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_x list of x (time) positions for rainfall spline +##' @param rainfall_y list of y (rainfall) positions for rainfall spline 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_x = seq(0, 9, 1), rainfall_y = rep(0, 10)) { @@ -95,7 +94,6 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4, e$set_soil_water_state(soil_initial_state) } - e$set_soil_infiltration_rate(soil_infiltration_rate) e$rainfall_init(rainfall_x, rainfall_y) return(e) diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index ac15146b8..44de54a1b 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -903,9 +903,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"] return_type: void diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h index fb49fd27b..b5762cc2d 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -21,10 +21,8 @@ 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(soil_number_of_depths, 0.0)); - rainfall_spline = interpolator::Interpolator(); }; @@ -100,8 +98,6 @@ 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++) { @@ -123,10 +119,6 @@ class FF16_Environment : public Environment { } } - void set_soil_infiltration_rate(double rate) { - soil_infiltration_rate = rate; - } - // Core functions template void compute_environment(Function f_compute_competition, double height_max) { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 602fb875a..ec91688a3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -9361,17 +9361,6 @@ BEGIN_RCPP return R_NilValue; END_RCPP } -// FF16_Environment__set_soil_infiltration_rate -void FF16_Environment__set_soil_infiltration_rate(plant::RcppR6::RcppR6 obj_, double rate); -RcppExport SEXP _plant_FF16_Environment__set_soil_infiltration_rate(SEXP obj_SEXP, SEXP rateSEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); - Rcpp::traits::input_parameter< double >::type rate(rateSEXP); - FF16_Environment__set_soil_infiltration_rate(obj_, rate); - return R_NilValue; -END_RCPP -} // FF16_Environment__set_soil_water_state void FF16_Environment__set_soil_water_state(plant::RcppR6::RcppR6 obj_, std::vector state); RcppExport SEXP _plant_FF16_Environment__set_soil_water_state(SEXP obj_SEXP, SEXP stateSEXP) { @@ -10681,7 +10670,6 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_FF16_Environment__canopy_openness", (DL_FUNC) &_plant_FF16_Environment__canopy_openness, 2}, {"_plant_FF16_Environment__clear", (DL_FUNC) &_plant_FF16_Environment__clear, 1}, {"_plant_FF16_Environment__set_fixed_environment", (DL_FUNC) &_plant_FF16_Environment__set_fixed_environment, 3}, - {"_plant_FF16_Environment__set_soil_infiltration_rate", (DL_FUNC) &_plant_FF16_Environment__set_soil_infiltration_rate, 2}, {"_plant_FF16_Environment__set_soil_water_state", (DL_FUNC) &_plant_FF16_Environment__set_soil_water_state, 2}, {"_plant_FF16_Environment__compute_rates", (DL_FUNC) &_plant_FF16_Environment__compute_rates, 1}, {"_plant_FF16_Environment__time__get", (DL_FUNC) &_plant_FF16_Environment__time__get, 1}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index a7e1660cd..0f90a6557 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -3796,10 +3796,6 @@ void FF16_Environment__set_fixed_environment(plant::RcppR6::RcppR6set_fixed_environment(value, height_max); } // [[Rcpp::export]] -void FF16_Environment__set_soil_infiltration_rate(plant::RcppR6::RcppR6 obj_, double rate) { - obj_->set_soil_infiltration_rate(rate); -} -// [[Rcpp::export]] void FF16_Environment__set_soil_water_state(plant::RcppR6::RcppR6 obj_, std::vector state) { obj_->set_soil_water_state(state); } diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index 208c4915a..87dd302c0 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -19,7 +19,6 @@ test_that("FF16w Environment", { # Make it rain x <- seq(0, 9, 1) y <- rep(1, 10) - #e$set_soil_infiltration_rate(10) e$rainfall_init(x, y) expect_equal(e$soil$states, 0) e$compute_rates() @@ -33,7 +32,6 @@ test_that("FF16w Environment", { # Check construction e <- make_environment("FF16w", soil_number_of_depths = 1, - soil_infiltration_rate = 10, # not needed anymore? rainfall_x = seq(0, 9, 1), rainfall_y = rep(10, 10)) @@ -58,35 +56,6 @@ test_that("FF16w Environment", { }) -# test_that("Basic run", { -# # one species -# p0 <- scm_base_parameters("FF16w") -# p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16w_hyperpar,FALSE) -# p1$birth_rate <- 20 -# -# env <- make_environment("FF16w", -# soil_number_of_depths = 10, -# soil_initial_state = rep(1, 10), -# soil_infiltration_rate = 1) -# -# ctrl <- scm_base_control() -# -# out <- run_scm(p1, env, ctrl) -# -# expect_equal(out$patch$environment$ode_size, 10) -# -# expect_equal(out$patch$environment$soil$rates, -# c(1.00, 0.50, 0.33, 0.25, 0.20, 0.17, 0.14, 0.13, 0.11, 0.10), -# tolerance = .01) -# -# expect_equal(out$patch$environment$soil$states, -# c(105, 52, 35, 26, 21, 17, 15, 13, 11, 10), -# tolerance = 0.1) -# -# expect_equal(out$offspring_production, 16.88946, tolerance=1e-5) -# expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-5) -# }) - test_that("Rainfall spline basic run", { # one species p0 <- scm_base_parameters("FF16w") @@ -95,8 +64,7 @@ test_that("Rainfall spline basic run", { env <- make_environment("FF16w", soil_number_of_depths = 10, - soil_initial_state = rep(1, 10), - soil_infiltration_rate = 1) + soil_initial_state = rep(1, 10)) # init rainfall spline for env x <- seq(0, 110, 1) From 47bac04fd3be763dc905adca097ac0dff99ab07e Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Sun, 28 Nov 2021 16:36:28 +1100 Subject: [PATCH 09/22] some changes to testing with regards to tolerance --- src/interpolator.cpp | 2 +- tests/testthat/test-environment.R | 22 ++++++++++------------ tests/testthat/test-interpolator.R | 2 +- tests/testthat/test-strategy-ff16w.R | 10 +++++----- 4 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/interpolator.cpp b/src/interpolator.cpp index af3e0e568..ce08479c9 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -22,7 +22,7 @@ void Interpolator::init(const std::vector& x_, // and 'y'. void Interpolator::initialise() { if (not std::is_sorted(x.begin(), x.end())) { - util::stop("spline control points must be in non-descending order"); + util::stop("spline control points must be unique and in ascending order"); } if (x.size() > 0) { tk_spline.set_points(x, y); diff --git a/tests/testthat/test-environment.R b/tests/testthat/test-environment.R index 358e94de4..05d4c97db 100644 --- a/tests/testthat/test-environment.R +++ b/tests/testthat/test-environment.R @@ -47,29 +47,27 @@ test_that("FF16 rainfall spline", { env <- make_environment("FF16") ## simple quadratic - x <- seq(-10, 10, 2) + x <- seq(-10, 10, 0.41) y <- x^2 env$rainfall_init(x, y) - ## discrete points used for spline creation + # interpolated points expect_equal(env$rainfall_eval(2), 4) expect_equal(env$rainfall_eval(-2), 4) - - ## interpolated points - expect_equal(env$rainfall_eval(3), 9, tolerance=1e-2) - expect_equal(env$rainfall_eval(-3), 9, tolerance=1e-2) - expect_equal(env$rainfall_eval(5.5), 30.25, tolerance=1e-2) - expect_equal(env$rainfall_eval(-5.5), 30.25, tolerance=1e-2) + expect_equal(env$rainfall_eval(3), 9, tolerance=1e-7) + expect_equal(env$rainfall_eval(-3), 9, tolerance=1e-7) + expect_equal(env$rainfall_eval(5.5), 30.25, tolerance=1e-7) + expect_equal(env$rainfall_eval(-5.5), 30.25, tolerance=1e-7) ## interpolated range of points - expect_equal(env$rainfall_eval_range(c(-7, 1, 7.8345)), c(49, 1, 61.37939025), tolerance=1e-2) + expect_equal(env$rainfall_eval_range(c(-7, 1, 7.8345)), c(49, 1, 61.37939025), tolerance=1e-6) ## clear points, set new points and recompute spline env$rainfall_clear_points() - x <- seq(-10, 10, 1) + x <- seq(-10, 10, 0.23) y <- 2*x env$rainfall_add_points(x, y) env$rainfall_recompute() - expect_equal(env$rainfall_eval(3), 6, tolerance=1e-2) - expect_equal(env$rainfall_eval(3.3), 6.6, tolerance=1e-2) + expect_equal(env$rainfall_eval(3), 6, tolerance=1e-7) + expect_equal(env$rainfall_eval(3.3), 6.6, tolerance=1e-7) }) diff --git a/tests/testthat/test-interpolator.R b/tests/testthat/test-interpolator.R index 3a789fe33..3e154759f 100644 --- a/tests/testthat/test-interpolator.R +++ b/tests/testthat/test-interpolator.R @@ -26,7 +26,7 @@ test_that("Splines require sensible data", { s <- Interpolator() expect_error(s$init(c(1, 2, 3), 1), "Incorrect length input") expect_error(s$init(numeric(0), numeric(0)), "insufficient number of points") - expect_error(s$init(c(3, 2, 1), c(1, 1, 1)), "spline control points must be in non-descending order") + expect_error(s$init(c(3, 2, 1), c(1, 1, 1)), "spline control points must be unique and in ascending order") }) test_that("Spline contains correct data", { diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index 87dd302c0..8fa0bfb88 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -67,7 +67,7 @@ test_that("Rainfall spline basic run", { soil_initial_state = rep(1, 10)) # init rainfall spline for env - x <- seq(0, 110, 1) + x <- seq(0, 110, 0.1) integrand <- function(x) {x^2} y <- integrand(x) env$rainfall_init(x, y) @@ -84,15 +84,15 @@ test_that("Rainfall spline basic run", { # check the rates are correct, ie 105^2/depth expect_equal(out$patch$environment$soil$rates, sapply(seq(1, 10), function(x) { (out$time^2)/x }), - tolerance = .01) + tolerance = 1e-7) # check the states are correct # model is currently 1/z * int(0, t)(f(x)), where z is depth, # x is time, f is spline, t is number of years expect_equal(out$patch$environment$soil$states, sapply(seq(1, 10), function(i) { (1/i) * integrate(integrand, 0, out$time)$value }), - tolerance = 0.1) + tolerance = 1e-6) - expect_equal(out$offspring_production, 16.88946, tolerance=1e-5) - expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-5) + expect_equal(out$offspring_production, 16.88946, tolerance=1e-7) + expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-7) }) \ No newline at end of file From aec5d5b1209b9ae98bf45430eb0ceee83cff10c1 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Sun, 28 Nov 2021 16:39:21 +1100 Subject: [PATCH 10/22] fixed incorrect tolerance --- tests/testthat/test-strategy-ff16w.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index 8fa0bfb88..0ef50c659 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -91,7 +91,7 @@ test_that("Rainfall spline basic run", { # x is time, f is spline, t is number of years expect_equal(out$patch$environment$soil$states, sapply(seq(1, 10), function(i) { (1/i) * integrate(integrand, 0, out$time)$value }), - tolerance = 1e-6) + tolerance = 1e-5) expect_equal(out$offspring_production, 16.88946, tolerance=1e-7) expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-7) From ed4000a1770620f199df7446fb498ca0d7f56986 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Sun, 28 Nov 2021 17:09:42 +1100 Subject: [PATCH 11/22] fixed comp function for std::is_sorted + more specific testsasdasd --- src/interpolator.cpp | 3 ++- tests/testthat/test-interpolator.R | 5 ++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/interpolator.cpp b/src/interpolator.cpp index ce08479c9..aa32fe57f 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -21,7 +21,8 @@ void Interpolator::init(const std::vector& x_, // Compute the interpolated function from the points contained in 'x' // and 'y'. void Interpolator::initialise() { - if (not std::is_sorted(x.begin(), x.end())) { + // https://stackoverflow.com/questions/17769114/stdis-sorted-and-strictly-less-comparison + if (not std::is_sorted(x.begin(), x.end(), std::less_equal())) { util::stop("spline control points must be unique and in ascending order"); } if (x.size() > 0) { diff --git a/tests/testthat/test-interpolator.R b/tests/testthat/test-interpolator.R index 3e154759f..4e621645f 100644 --- a/tests/testthat/test-interpolator.R +++ b/tests/testthat/test-interpolator.R @@ -26,7 +26,10 @@ test_that("Splines require sensible data", { s <- Interpolator() expect_error(s$init(c(1, 2, 3), 1), "Incorrect length input") expect_error(s$init(numeric(0), numeric(0)), "insufficient number of points") - expect_error(s$init(c(3, 2, 1), c(1, 1, 1)), "spline control points must be unique and in ascending order") + # ascending order + expect_error(s$init(c(3, 2, 4), c(1, 1, 1)), "spline control points must be unique and in ascending order") + # unique + expect_error(s$init(c(1, 1, 2), c(1, 1, 1)), "spline control points must be unique and in ascending order") }) test_that("Spline contains correct data", { From e9a421105013efc7255a549dd245ca116d00e294 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Sun, 28 Nov 2021 17:11:49 +1100 Subject: [PATCH 12/22] less lenient tolerance to pass windows tests --- tests/testthat/test-strategy-ff16w.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index 0ef50c659..2a4077615 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -93,6 +93,6 @@ test_that("Rainfall spline basic run", { sapply(seq(1, 10), function(i) { (1/i) * integrate(integrand, 0, out$time)$value }), tolerance = 1e-5) - expect_equal(out$offspring_production, 16.88946, tolerance=1e-7) + expect_equal(out$offspring_production, 16.88946, tolerance=1e-5) expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-7) }) \ No newline at end of file From 23ff3de62d84b948b325d32db1b0bc92b7ca04b1 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Sun, 5 Dec 2021 22:05:29 +1100 Subject: [PATCH 13/22] Refactored FF16w rainfall splines in favour of more generic 'extrinsic drivers' --- R/RcppExports.R | 32 ++----- R/RcppR6.R | 29 ++---- R/ff16w.R | 10 +-- inst/RcppR6_classes.yml | 28 ++---- inst/include/plant/environment.h | 28 ++++++ inst/include/plant/models/ff16_environment.h | 47 +--------- src/RcppExports.cpp | 94 ++++---------------- src/RcppR6.cpp | 32 ++----- src/interpolator.cpp | 16 +--- tests/testthat/test-environment.R | 35 ++++---- tests/testthat/test-strategy-ff16w.R | 13 ++- 11 files changed, 105 insertions(+), 259 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index b2e42d627..01df5ac7c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -3293,36 +3293,16 @@ 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__rainfall_init <- function(obj_, x, y) { - invisible(.Call('_plant_FF16_Environment__rainfall_init', PACKAGE = 'plant', obj_, x, y)) +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__rainfall_recompute <- function(obj_) { - invisible(.Call('_plant_FF16_Environment__rainfall_recompute', PACKAGE = 'plant', obj_)) +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__rainfall_eval <- function(obj_, u) { - .Call('_plant_FF16_Environment__rainfall_eval', PACKAGE = 'plant', obj_, u) -} - -FF16_Environment__rainfall_eval_range <- function(obj_, u) { - .Call('_plant_FF16_Environment__rainfall_eval_range', PACKAGE = 'plant', obj_, u) -} - -FF16_Environment__rainfall_add_point <- function(obj_, x, y) { - invisible(.Call('_plant_FF16_Environment__rainfall_add_point', PACKAGE = 'plant', obj_, x, y)) -} - -FF16_Environment__rainfall_add_point_sorted <- function(obj_, x, y) { - invisible(.Call('_plant_FF16_Environment__rainfall_add_point_sorted', PACKAGE = 'plant', obj_, x, y)) -} - -FF16_Environment__rainfall_add_points <- function(obj_, x, y) { - invisible(.Call('_plant_FF16_Environment__rainfall_add_points', PACKAGE = 'plant', obj_, x, y)) -} - -FF16_Environment__rainfall_clear_points <- function(obj_) { - invisible(.Call('_plant_FF16_Environment__rainfall_clear_points', PACKAGE = 'plant', obj_)) +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__canopy_openness <- function(obj_, height) { diff --git a/R/RcppR6.R b/R/RcppR6.R index 2d4a4e689..9a749375b 100644 --- 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: 2149ebb5397c92102688df8a4c991aa7 +## Hash: d3e496d60c31b1fbfc22f90defa136ea ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -4814,29 +4814,14 @@ StochasticPatchRunner <- function(T, E) { initialize = function(ptr) { self$.ptr <- ptr }, - rainfall_init = function(x, y) { - FF16_Environment__rainfall_init(self, x, y) + set_extrinsic_driver = function(driver_name, x, y) { + FF16_Environment__set_extrinsic_driver(self, driver_name, x, y) }, - rainfall_recompute = function() { - FF16_Environment__rainfall_recompute(self) + extrinsic_driver_evaluate = function(driver_name, u) { + FF16_Environment__extrinsic_driver_evaluate(self, driver_name, u) }, - rainfall_eval = function(u) { - FF16_Environment__rainfall_eval(self, u) - }, - rainfall_eval_range = function(u) { - FF16_Environment__rainfall_eval_range(self, u) - }, - rainfall_add_point = function(x, y) { - FF16_Environment__rainfall_add_point(self, x, y) - }, - rainfall_add_point_sorted = function(x, y) { - FF16_Environment__rainfall_add_point_sorted(self, x, y) - }, - rainfall_add_points = function(x, y) { - FF16_Environment__rainfall_add_points(self, x, y) - }, - rainfall_clear_points = function() { - FF16_Environment__rainfall_clear_points(self) + extrinsic_driver_evaluate_range = function(driver_name, u) { + FF16_Environment__extrinsic_driver_evaluate_range(self, driver_name, u) }, canopy_openness = function(height) { FF16_Environment__canopy_openness(self, height) diff --git a/R/ff16w.R b/R/ff16w.R index 2fd266900..b7b8df402 100644 --- a/R/ff16w.R +++ b/R/ff16w.R @@ -65,16 +65,14 @@ FF16w_StochasticPatchRunner <- function(p) { ##' @export ##' @rdname FF16_Environment ##' @param soil_number_of_depths the number of soil layers -##' @param rainfall_x list of x (time) positions for rainfall spline -##' @param rainfall_y list of y (rainfall) positions for rainfall spline +##' @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, - rainfall_x = seq(0, 9, 1), - rainfall_y = rep(0, 10)) { + rainfall = 1) { if(soil_number_of_depths < 1) stop("FF16w Environment must have at least one soil layer") @@ -94,7 +92,9 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4, e$set_soil_water_state(soil_initial_state) } - e$rainfall_init(rainfall_x, rainfall_y) + x <- seq(0, 10, 1) + y <- rep(rainfall, 11) + e$set_extrinsic_driver("rainfall", x, y) return(e) } diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index 44de54a1b..8e402908a 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -871,30 +871,16 @@ FF16_Environment: FF16_Environment object @export methods: - rainfall_init: - args: [x: "std::vector", y: "std::vector"] + # unforunately need to expose these base Environment functions for every derived class... + set_extrinsic_driver: + args: [driver_name: "std::string", x: "std::vector", y: "std::vector"] return_type: void - rainfall_recompute: - args: [] - return_type: void - rainfall_eval: - args: [u: double] + extrinsic_driver_evaluate: + args: [driver_name: "std::string", u: double] return_type: double - rainfall_eval_range: - args: [u: "std::vector"] + extrinsic_driver_evaluate_range: + args: [driver_name: "std::string", u: "std::vector"] return_type: "std::vector" - rainfall_add_point: - args: [x: double, y: double] - return_type: void - rainfall_add_point_sorted: - args: [x: double, y: double] - return_type: void - rainfall_add_points: - args: [x: "std::vector", y: "std::vector"] - return_type: void - rainfall_clear_points: - args: [] - return_type: void canopy_openness: args: [height: double] return_type: double diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h index 67073917d..bf3880542 100644 --- a/inst/include/plant/environment.h +++ b/inst/include/plant/environment.h @@ -8,9 +8,11 @@ #include #include #include +#include #include using namespace Rcpp; +// interpolator has its own namespace, doing this for cleaner code namespace plant { @@ -69,6 +71,32 @@ 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 const& x, std::vector 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[driver_name].init(x, y); + } + } + + // 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 { + extrinsic_drivers.at(driver_name).eval(u); + } + + // // evaluate/query interpolated spline for driver at vector of points, return vector of values + std::vector extrinsic_driver_evaluate_range(std::string driver_name, std::vector u) const { + return extrinsic_drivers.at(driver_name).r_eval(u); + } + +// private? +protected: + std::unordered_map extrinsic_drivers; }; } diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h index b5762cc2d..6cb312061 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -23,52 +23,9 @@ class FF16_Environment : public Environment { canopy = Canopy(); vars = Internals(soil_number_of_depths); set_soil_water_state(std::vector(soil_number_of_depths, 0.0)); - rainfall_spline = interpolator::Interpolator(); + extrinsic_drivers["rainfall"] = interpolator::Interpolator(); }; - // first spline computation - void rainfall_init(std::vector const& x, std::vector const& y) { - rainfall_spline.init(x, y); // set x, y points and compute spline - } - - // if any points are added to the spline, need to recompute spline - void rainfall_recompute() { - rainfall_spline.initialise(); - } - - double rainfall_eval(double u) const { - return rainfall_spline.eval(u); - } - - std::vector rainfall_eval_range(std::vector const& u) const { - return rainfall_spline.r_eval(u); - } - - /* if xi is not larger than every other x in the spline, recomputing the spline - will throw an error and will need to clear_points() and reconstruct spline */ - void rainfall_add_point(double xi, double yi) { - rainfall_spline.add_point(xi, yi); - } - - // safe add point - void rainfall_add_point_sorted(double xi, double yi) { - rainfall_spline.add_point_sorted(xi, yi); - } - - /* not sure if these two are necessary for the rainfall model. could also move these - to be general interpolator.h functions? */ - void rainfall_add_points(std::vector const& x, std::vector const& y) { - util::check_length(x.size(), y.size()); - for (size_t i = 0; i < x.size(); ++i) { - rainfall_add_point(x[i], y[i]); // refactor? - } - } - - // rainfall_add_points_sorted ? - - void rainfall_clear_points() { - rainfall_spline.clear(); - } // Light interface bool canopy_rescale_usually; @@ -101,7 +58,7 @@ class FF16_Environment : public Environment { virtual void compute_rates() { for (size_t i = 0; i < vars.state_size; i++) { - vars.set_rate(i, rainfall_eval(time) / (i+1)); + vars.set_rate(i, extrinsic_drivers["rainfall"].eval(time) / (i+1)); } } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index ec91688a3..640d830e9 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -9235,98 +9235,45 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// FF16_Environment__rainfall_init -void FF16_Environment__rainfall_init(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y); -RcppExport SEXP _plant_FF16_Environment__rainfall_init(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { +// FF16_Environment__set_extrinsic_driver +void FF16_Environment__set_extrinsic_driver(plant::RcppR6::RcppR6 obj_, std::string driver_name, std::vector x, std::vector y); +RcppExport SEXP _plant_FF16_Environment__set_extrinsic_driver(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP xSEXP, SEXP ySEXP) { BEGIN_RCPP Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP); Rcpp::traits::input_parameter< std::vector >::type x(xSEXP); Rcpp::traits::input_parameter< std::vector >::type y(ySEXP); - FF16_Environment__rainfall_init(obj_, x, y); + FF16_Environment__set_extrinsic_driver(obj_, driver_name, x, y); return R_NilValue; END_RCPP } -// FF16_Environment__rainfall_recompute -void FF16_Environment__rainfall_recompute(plant::RcppR6::RcppR6 obj_); -RcppExport SEXP _plant_FF16_Environment__rainfall_recompute(SEXP obj_SEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); - FF16_Environment__rainfall_recompute(obj_); - return R_NilValue; -END_RCPP -} -// FF16_Environment__rainfall_eval -double FF16_Environment__rainfall_eval(plant::RcppR6::RcppR6 obj_, double u); -RcppExport SEXP _plant_FF16_Environment__rainfall_eval(SEXP obj_SEXP, SEXP uSEXP) { +// FF16_Environment__extrinsic_driver_evaluate +double FF16_Environment__extrinsic_driver_evaluate(plant::RcppR6::RcppR6 obj_, std::string driver_name, double u); +RcppExport SEXP _plant_FF16_Environment__extrinsic_driver_evaluate(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP uSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP); Rcpp::traits::input_parameter< double >::type u(uSEXP); - rcpp_result_gen = Rcpp::wrap(FF16_Environment__rainfall_eval(obj_, u)); + rcpp_result_gen = Rcpp::wrap(FF16_Environment__extrinsic_driver_evaluate(obj_, driver_name, u)); return rcpp_result_gen; END_RCPP } -// FF16_Environment__rainfall_eval_range -std::vector FF16_Environment__rainfall_eval_range(plant::RcppR6::RcppR6 obj_, std::vector u); -RcppExport SEXP _plant_FF16_Environment__rainfall_eval_range(SEXP obj_SEXP, SEXP uSEXP) { +// FF16_Environment__extrinsic_driver_evaluate_range +std::vector FF16_Environment__extrinsic_driver_evaluate_range(plant::RcppR6::RcppR6 obj_, std::string driver_name, std::vector u); +RcppExport SEXP _plant_FF16_Environment__extrinsic_driver_evaluate_range(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP uSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP); Rcpp::traits::input_parameter< std::vector >::type u(uSEXP); - rcpp_result_gen = Rcpp::wrap(FF16_Environment__rainfall_eval_range(obj_, u)); + rcpp_result_gen = Rcpp::wrap(FF16_Environment__extrinsic_driver_evaluate_range(obj_, driver_name, u)); return rcpp_result_gen; END_RCPP } -// FF16_Environment__rainfall_add_point -void FF16_Environment__rainfall_add_point(plant::RcppR6::RcppR6 obj_, double x, double y); -RcppExport SEXP _plant_FF16_Environment__rainfall_add_point(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); - Rcpp::traits::input_parameter< double >::type x(xSEXP); - Rcpp::traits::input_parameter< double >::type y(ySEXP); - FF16_Environment__rainfall_add_point(obj_, x, y); - return R_NilValue; -END_RCPP -} -// FF16_Environment__rainfall_add_point_sorted -void FF16_Environment__rainfall_add_point_sorted(plant::RcppR6::RcppR6 obj_, double x, double y); -RcppExport SEXP _plant_FF16_Environment__rainfall_add_point_sorted(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); - Rcpp::traits::input_parameter< double >::type x(xSEXP); - Rcpp::traits::input_parameter< double >::type y(ySEXP); - FF16_Environment__rainfall_add_point_sorted(obj_, x, y); - return R_NilValue; -END_RCPP -} -// FF16_Environment__rainfall_add_points -void FF16_Environment__rainfall_add_points(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y); -RcppExport SEXP _plant_FF16_Environment__rainfall_add_points(SEXP obj_SEXP, SEXP xSEXP, SEXP ySEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); - Rcpp::traits::input_parameter< std::vector >::type x(xSEXP); - Rcpp::traits::input_parameter< std::vector >::type y(ySEXP); - FF16_Environment__rainfall_add_points(obj_, x, y); - return R_NilValue; -END_RCPP -} -// FF16_Environment__rainfall_clear_points -void FF16_Environment__rainfall_clear_points(plant::RcppR6::RcppR6 obj_); -RcppExport SEXP _plant_FF16_Environment__rainfall_clear_points(SEXP obj_SEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); - FF16_Environment__rainfall_clear_points(obj_); - return R_NilValue; -END_RCPP -} // FF16_Environment__canopy_openness double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6 obj_, double height); RcppExport SEXP _plant_FF16_Environment__canopy_openness(SEXP obj_SEXP, SEXP heightSEXP) { @@ -10659,14 +10606,9 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_Weibull_Disturbance_Regime__icdf", (DL_FUNC) &_plant_Weibull_Disturbance_Regime__icdf, 2}, {"_plant_FF16_Strategy__ctor", (DL_FUNC) &_plant_FF16_Strategy__ctor, 0}, {"_plant_FF16_Environment__ctor", (DL_FUNC) &_plant_FF16_Environment__ctor, 2}, - {"_plant_FF16_Environment__rainfall_init", (DL_FUNC) &_plant_FF16_Environment__rainfall_init, 3}, - {"_plant_FF16_Environment__rainfall_recompute", (DL_FUNC) &_plant_FF16_Environment__rainfall_recompute, 1}, - {"_plant_FF16_Environment__rainfall_eval", (DL_FUNC) &_plant_FF16_Environment__rainfall_eval, 2}, - {"_plant_FF16_Environment__rainfall_eval_range", (DL_FUNC) &_plant_FF16_Environment__rainfall_eval_range, 2}, - {"_plant_FF16_Environment__rainfall_add_point", (DL_FUNC) &_plant_FF16_Environment__rainfall_add_point, 3}, - {"_plant_FF16_Environment__rainfall_add_point_sorted", (DL_FUNC) &_plant_FF16_Environment__rainfall_add_point_sorted, 3}, - {"_plant_FF16_Environment__rainfall_add_points", (DL_FUNC) &_plant_FF16_Environment__rainfall_add_points, 3}, - {"_plant_FF16_Environment__rainfall_clear_points", (DL_FUNC) &_plant_FF16_Environment__rainfall_clear_points, 1}, + {"_plant_FF16_Environment__set_extrinsic_driver", (DL_FUNC) &_plant_FF16_Environment__set_extrinsic_driver, 4}, + {"_plant_FF16_Environment__extrinsic_driver_evaluate", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate, 3}, + {"_plant_FF16_Environment__extrinsic_driver_evaluate_range", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate_range, 3}, {"_plant_FF16_Environment__canopy_openness", (DL_FUNC) &_plant_FF16_Environment__canopy_openness, 2}, {"_plant_FF16_Environment__clear", (DL_FUNC) &_plant_FF16_Environment__clear, 1}, {"_plant_FF16_Environment__set_fixed_environment", (DL_FUNC) &_plant_FF16_Environment__set_fixed_environment, 3}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index 0f90a6557..c8870dad9 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -3752,36 +3752,16 @@ plant::FF16_Environment FF16_Environment__ctor(bool canopy_rescale_usually, int return plant::FF16_Environment(canopy_rescale_usually, soil_number_of_depths); } // [[Rcpp::export]] -void FF16_Environment__rainfall_init(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y) { - obj_->rainfall_init(x, y); +void FF16_Environment__set_extrinsic_driver(plant::RcppR6::RcppR6 obj_, std::string driver_name, std::vector x, std::vector y) { + obj_->set_extrinsic_driver(driver_name, x, y); } // [[Rcpp::export]] -void FF16_Environment__rainfall_recompute(plant::RcppR6::RcppR6 obj_) { - obj_->rainfall_recompute(); +double FF16_Environment__extrinsic_driver_evaluate(plant::RcppR6::RcppR6 obj_, std::string driver_name, double u) { + return obj_->extrinsic_driver_evaluate(driver_name, u); } // [[Rcpp::export]] -double FF16_Environment__rainfall_eval(plant::RcppR6::RcppR6 obj_, double u) { - return obj_->rainfall_eval(u); -} -// [[Rcpp::export]] -std::vector FF16_Environment__rainfall_eval_range(plant::RcppR6::RcppR6 obj_, std::vector u) { - return obj_->rainfall_eval_range(u); -} -// [[Rcpp::export]] -void FF16_Environment__rainfall_add_point(plant::RcppR6::RcppR6 obj_, double x, double y) { - obj_->rainfall_add_point(x, y); -} -// [[Rcpp::export]] -void FF16_Environment__rainfall_add_point_sorted(plant::RcppR6::RcppR6 obj_, double x, double y) { - obj_->rainfall_add_point_sorted(x, y); -} -// [[Rcpp::export]] -void FF16_Environment__rainfall_add_points(plant::RcppR6::RcppR6 obj_, std::vector x, std::vector y) { - obj_->rainfall_add_points(x, y); -} -// [[Rcpp::export]] -void FF16_Environment__rainfall_clear_points(plant::RcppR6::RcppR6 obj_) { - obj_->rainfall_clear_points(); +std::vector FF16_Environment__extrinsic_driver_evaluate_range(plant::RcppR6::RcppR6 obj_, std::string driver_name, std::vector u) { + return obj_->extrinsic_driver_evaluate_range(driver_name, u); } // [[Rcpp::export]] double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6 obj_, double height) { diff --git a/src/interpolator.cpp b/src/interpolator.cpp index aa32fe57f..85e49f7bc 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -55,12 +55,7 @@ void Interpolator::clear() { // Compute the value of the interpolated function at point `x=u` double Interpolator::eval(double u) const { - check_active(); - return tk_spline(u); -} - -// faster version -double Interpolator::operator()(double u) const { + check_active(); // slower to check this every time return tk_spline(u); } @@ -99,13 +94,8 @@ SEXP Interpolator::r_get_xy() const { // Compute the value of the interpolated function at a vector of // points `x=u`, returning a vector of the same length. std::vector Interpolator::r_eval(std::vector u) const { - check_active(); - const size_t n = u.size(); - std::vector ret(n); - // for (size_t i = 0; i < n; ++i) { - // ret[i] = operator()(u[i]); - // } - std::transform(u.begin(), u.end(), ret.begin(), [&](double x){return operator()(x);}); + std::vector ret(u.size()); + std::transform(u.begin(), u.end(), ret.begin(), [&](double x){ return eval(x); }); return ret; } diff --git a/tests/testthat/test-environment.R b/tests/testthat/test-environment.R index 05d4c97db..503fd4d29 100644 --- a/tests/testthat/test-environment.R +++ b/tests/testthat/test-environment.R @@ -44,30 +44,29 @@ for (x in names(strategy_types)) { test_that("FF16 rainfall spline", { context(sprintf("Rainfall-FF16",x)) - env <- make_environment("FF16") + env <- make_environment("FF16w") + # test extrapolation on default spline of y = 1 + expect_equal(env$extrinsic_driver_evaluate("rainfall", 100), 1) + expect_equal(env$extrinsic_driver_evaluate("rainfall", 10000000), 1) + + # test extrapolation on spline of y = 5.613432 + env <- make_environment("FF16w", rainfall=5.613432) + expect_equal(env$extrinsic_driver_evaluate("rainfall", 100), 5.613432) + expect_equal(env$extrinsic_driver_evaluate("rainfall", 10000000), 5.613432) ## simple quadratic x <- seq(-10, 10, 0.41) y <- x^2 - env$rainfall_init(x, y) + env$set_extrinsic_driver("rainfall", x, y) # overwrites previously created spline # interpolated points - expect_equal(env$rainfall_eval(2), 4) - expect_equal(env$rainfall_eval(-2), 4) - expect_equal(env$rainfall_eval(3), 9, tolerance=1e-7) - expect_equal(env$rainfall_eval(-3), 9, tolerance=1e-7) - expect_equal(env$rainfall_eval(5.5), 30.25, tolerance=1e-7) - expect_equal(env$rainfall_eval(-5.5), 30.25, tolerance=1e-7) + expect_equal(env$extrinsic_driver_evaluate("rainfall", 2), 4) + expect_equal(env$extrinsic_driver_evaluate("rainfall", -2), 4) + expect_equal(env$extrinsic_driver_evaluate("rainfall", 3), 9, tolerance=1e-7) + expect_equal(env$extrinsic_driver_evaluate("rainfall", -3), 9, tolerance=1e-7) + expect_equal(env$extrinsic_driver_evaluate("rainfall", 5.5), 30.25, tolerance=1e-7) + expect_equal(env$extrinsic_driver_evaluate("rainfall", -5.5), 30.25, tolerance=1e-7) ## interpolated range of points - expect_equal(env$rainfall_eval_range(c(-7, 1, 7.8345)), c(49, 1, 61.37939025), tolerance=1e-6) - - ## clear points, set new points and recompute spline - env$rainfall_clear_points() - x <- seq(-10, 10, 0.23) - y <- 2*x - env$rainfall_add_points(x, y) - env$rainfall_recompute() - expect_equal(env$rainfall_eval(3), 6, tolerance=1e-7) - expect_equal(env$rainfall_eval(3.3), 6.6, tolerance=1e-7) + expect_equal(env$extrinsic_driver_evaluate_range("rainfall", c(-7, 1, 7.8345)), c(49, 1, 61.37939025), tolerance=1e-6) }) diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index 2a4077615..fef2d8d36 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -14,15 +14,15 @@ test_that("FF16w Environment", { expect_equal(e$soil$states, 0) e$compute_rates() - expect_equal(e$soil$rates, 0) + expect_equal(e$soil$rates, 1) # default rainfall is now y = 1 # Make it rain x <- seq(0, 9, 1) - y <- rep(1, 10) - e$rainfall_init(x, y) + y <- rep(5, 10) + e$set_extrinsic_driver("rainfall", x, y) expect_equal(e$soil$states, 0) e$compute_rates() - expect_equal(e$soil$rates, 1) + expect_equal(e$soil$rates, 5) # not needed anymore? should we still be able to set the states manually? # Water logged @@ -32,8 +32,7 @@ test_that("FF16w Environment", { # Check construction e <- make_environment("FF16w", soil_number_of_depths = 1, - rainfall_x = seq(0, 9, 1), - rainfall_y = rep(10, 10)) + rainfall = 10) expect_equal(e$soil$states, 0) e$compute_rates() @@ -70,7 +69,7 @@ test_that("Rainfall spline basic run", { x <- seq(0, 110, 0.1) integrand <- function(x) {x^2} y <- integrand(x) - env$rainfall_init(x, y) + env$set_extrinsic_driver("rainfall", x, y) ctrl <- scm_base_control() From 1604176a6c5ce027f869bb53032a7bc92f2fcbf5 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Sun, 5 Dec 2021 22:24:16 +1100 Subject: [PATCH 14/22] Ability to get extrinsic driver names in R + small refactoring --- R/RcppExports.R | 4 ++++ R/RcppR6.R | 5 ++++- inst/RcppR6_classes.yml | 3 +++ inst/include/plant/environment.h | 8 ++++++++ src/RcppExports.cpp | 12 ++++++++++++ src/RcppR6.cpp | 4 ++++ src/interpolator.cpp | 7 +++++-- tests/testthat/test-environment.R | 3 +++ 8 files changed, 43 insertions(+), 3 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 01df5ac7c..237fdc21a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -3305,6 +3305,10 @@ FF16_Environment__extrinsic_driver_evaluate_range <- function(obj_, driver_name, .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) } diff --git a/R/RcppR6.R b/R/RcppR6.R index 9a749375b..73a07333c 100644 --- 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: d3e496d60c31b1fbfc22f90defa136ea +## Hash: 845ce882240dceae4ad386bd438b712d ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -4823,6 +4823,9 @@ StochasticPatchRunner <- function(T, E) { 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) }, diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index 8e402908a..c4d28568b 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -881,6 +881,9 @@ FF16_Environment: extrinsic_driver_evaluate_range: args: [driver_name: "std::string", u: "std::vector"] return_type: "std::vector" + get_extrinsic_driver_names: + args: [] + return_type: "std::vector" canopy_openness: args: [height: double] return_type: double diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h index bf3880542..d99aa07d0 100644 --- a/inst/include/plant/environment.h +++ b/inst/include/plant/environment.h @@ -94,6 +94,14 @@ class Environment { return extrinsic_drivers.at(driver_name).r_eval(u); } + std::vector get_extrinsic_driver_names() { + auto ret = std::vector(); + for (auto const& driver : extrinsic_drivers) { + ret.push_back(driver.first); + } + return ret; + } + // private? protected: std::unordered_map extrinsic_drivers; diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 640d830e9..f2919971b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -9274,6 +9274,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// FF16_Environment__get_extrinsic_driver_names +std::vector FF16_Environment__get_extrinsic_driver_names(plant::RcppR6::RcppR6 obj_); +RcppExport SEXP _plant_FF16_Environment__get_extrinsic_driver_names(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(FF16_Environment__get_extrinsic_driver_names(obj_)); + return rcpp_result_gen; +END_RCPP +} // FF16_Environment__canopy_openness double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6 obj_, double height); RcppExport SEXP _plant_FF16_Environment__canopy_openness(SEXP obj_SEXP, SEXP heightSEXP) { @@ -10609,6 +10620,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_FF16_Environment__set_extrinsic_driver", (DL_FUNC) &_plant_FF16_Environment__set_extrinsic_driver, 4}, {"_plant_FF16_Environment__extrinsic_driver_evaluate", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate, 3}, {"_plant_FF16_Environment__extrinsic_driver_evaluate_range", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate_range, 3}, + {"_plant_FF16_Environment__get_extrinsic_driver_names", (DL_FUNC) &_plant_FF16_Environment__get_extrinsic_driver_names, 1}, {"_plant_FF16_Environment__canopy_openness", (DL_FUNC) &_plant_FF16_Environment__canopy_openness, 2}, {"_plant_FF16_Environment__clear", (DL_FUNC) &_plant_FF16_Environment__clear, 1}, {"_plant_FF16_Environment__set_fixed_environment", (DL_FUNC) &_plant_FF16_Environment__set_fixed_environment, 3}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index c8870dad9..66b05f93c 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -3764,6 +3764,10 @@ std::vector FF16_Environment__extrinsic_driver_evaluate_range(plant::Rcp return obj_->extrinsic_driver_evaluate_range(driver_name, u); } // [[Rcpp::export]] +std::vector FF16_Environment__get_extrinsic_driver_names(plant::RcppR6::RcppR6 obj_) { + return obj_->get_extrinsic_driver_names(); +} +// [[Rcpp::export]] double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6 obj_, double height) { return obj_->canopy_openness(height); } diff --git a/src/interpolator.cpp b/src/interpolator.cpp index 85e49f7bc..05e64f78a 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -94,8 +94,11 @@ SEXP Interpolator::r_get_xy() const { // Compute the value of the interpolated function at a vector of // points `x=u`, returning a vector of the same length. std::vector Interpolator::r_eval(std::vector u) const { - std::vector ret(u.size()); - std::transform(u.begin(), u.end(), ret.begin(), [&](double x){ return eval(x); }); + auto ret = std::vector(); + ret.reserve(u.size()); // fast to do this once rather than multiple times with push_back + for (auto const& x : u) { + ret.push_back(eval(x)); + } return ret; } diff --git a/tests/testthat/test-environment.R b/tests/testthat/test-environment.R index 503fd4d29..d1410a716 100644 --- a/tests/testthat/test-environment.R +++ b/tests/testthat/test-environment.R @@ -45,6 +45,9 @@ for (x in names(strategy_types)) { test_that("FF16 rainfall spline", { context(sprintf("Rainfall-FF16",x)) env <- make_environment("FF16w") + # get list of extrinsic drivers for the environment + expect_equal(env$get_extrinsic_driver_names(), c("rainfall")) + # test extrapolation on default spline of y = 1 expect_equal(env$extrinsic_driver_evaluate("rainfall", 100), 1) expect_equal(env$extrinsic_driver_evaluate("rainfall", 10000000), 1) From 518019e91e525e6b9a6ee3c8c95d0b79bbfc4d7f Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Tue, 7 Dec 2021 20:41:56 +1100 Subject: [PATCH 15/22] fixed old typo in patch tests --- tests/testthat/test-patch.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-patch.R b/tests/testthat/test-patch.R index ab3b3379d..e50d7d489 100644 --- a/tests/testthat/test-patch.R +++ b/tests/testthat/test-patch.R @@ -42,8 +42,8 @@ for (x in names(strategy_types)) { 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)) + expect_identical(patch$ode_state, env_state) + expect_identical(patch$ode_rates, env_rates) ## Empty environment: patch$compute_environment() From 333592594bebb8f11367540f79613ce9f6e1563e Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Thu, 9 Dec 2021 11:37:17 +1100 Subject: [PATCH 16/22] extrapolation setting in interpolator --- R/RcppExports.R | 4 ++++ R/RcppR6.R | 5 ++++- R/ff16w.R | 1 + inst/RcppR6_classes.yml | 3 +++ inst/include/plant/environment.h | 5 +++++ inst/include/plant/interpolator.h | 4 +++- inst/include/plant/models/ff16_environment.h | 1 - src/RcppExports.cpp | 13 +++++++++++++ src/RcppR6.cpp | 4 ++++ src/interpolator.cpp | 20 +++++++++++++++++--- 10 files changed, 54 insertions(+), 6 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 237fdc21a..3d6a713c9 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -3297,6 +3297,10 @@ 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) } diff --git a/R/RcppR6.R b/R/RcppR6.R index 73a07333c..7c5a8a115 100644 --- 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: 845ce882240dceae4ad386bd438b712d +## Hash: 2cb4fee1111abff105bb05242a8b85a0 ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -4817,6 +4817,9 @@ StochasticPatchRunner <- function(T, E) { 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) }, diff --git a/R/ff16w.R b/R/ff16w.R index b7b8df402..6e85b1fc8 100644 --- a/R/ff16w.R +++ b/R/ff16w.R @@ -95,6 +95,7 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4, 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 index c4d28568b..2c9a50572 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -875,6 +875,9 @@ FF16_Environment: set_extrinsic_driver: args: [driver_name: "std::string", x: "std::vector", y: "std::vector"] 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 diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h index d99aa07d0..c43689deb 100644 --- a/inst/include/plant/environment.h +++ b/inst/include/plant/environment.h @@ -81,9 +81,14 @@ class Environment { util::stop(driver_name + " doesn't exist in the list of extrinsic_drivers."); } else { extrinsic_drivers[driver_name].init(x, y); + extrinsic_drivers[driver_name].setExtrapolate(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 { extrinsic_drivers.at(driver_name).eval(u); diff --git a/inst/include/plant/interpolator.h b/inst/include/plant/interpolator.h index 56d2b4549..738e2de16 100644 --- a/inst/include/plant/interpolator.h +++ b/inst/include/plant/interpolator.h @@ -25,6 +25,7 @@ class Interpolator { double min() const; double max() const; + void setExtrapolate(bool e); std::vector get_x() const; std::vector get_y() const; @@ -38,7 +39,8 @@ class Interpolator { void check_active() const; std::vector 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 index 6cb312061..8189e8e16 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -55,7 +55,6 @@ class FF16_Environment : public Environment { canopy.r_init_interpolators(state); } - virtual void compute_rates() { for (size_t i = 0; i < vars.state_size; i++) { vars.set_rate(i, extrinsic_drivers["rainfall"].eval(time) / (i+1)); diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f2919971b..ed33a1714 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -9248,6 +9248,18 @@ BEGIN_RCPP return R_NilValue; END_RCPP } +// FF16_Environment__extrinsic_driver_extrapolation +void FF16_Environment__extrinsic_driver_extrapolation(plant::RcppR6::RcppR6 obj_, std::string driver_name, bool extrapolate); +RcppExport SEXP _plant_FF16_Environment__extrinsic_driver_extrapolation(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP extrapolateSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP); + Rcpp::traits::input_parameter< bool >::type extrapolate(extrapolateSEXP); + FF16_Environment__extrinsic_driver_extrapolation(obj_, driver_name, extrapolate); + return R_NilValue; +END_RCPP +} // FF16_Environment__extrinsic_driver_evaluate double FF16_Environment__extrinsic_driver_evaluate(plant::RcppR6::RcppR6 obj_, std::string driver_name, double u); RcppExport SEXP _plant_FF16_Environment__extrinsic_driver_evaluate(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP uSEXP) { @@ -10618,6 +10630,7 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_FF16_Strategy__ctor", (DL_FUNC) &_plant_FF16_Strategy__ctor, 0}, {"_plant_FF16_Environment__ctor", (DL_FUNC) &_plant_FF16_Environment__ctor, 2}, {"_plant_FF16_Environment__set_extrinsic_driver", (DL_FUNC) &_plant_FF16_Environment__set_extrinsic_driver, 4}, + {"_plant_FF16_Environment__extrinsic_driver_extrapolation", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_extrapolation, 3}, {"_plant_FF16_Environment__extrinsic_driver_evaluate", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate, 3}, {"_plant_FF16_Environment__extrinsic_driver_evaluate_range", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate_range, 3}, {"_plant_FF16_Environment__get_extrinsic_driver_names", (DL_FUNC) &_plant_FF16_Environment__get_extrinsic_driver_names, 1}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index 66b05f93c..672bea1c0 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -3756,6 +3756,10 @@ void FF16_Environment__set_extrinsic_driver(plant::RcppR6::RcppR6set_extrinsic_driver(driver_name, x, y); } // [[Rcpp::export]] +void FF16_Environment__extrinsic_driver_extrapolation(plant::RcppR6::RcppR6 obj_, std::string driver_name, bool extrapolate) { + obj_->extrinsic_driver_extrapolation(driver_name, extrapolate); +} +// [[Rcpp::export]] double FF16_Environment__extrinsic_driver_evaluate(plant::RcppR6::RcppR6 obj_, std::string driver_name, double u) { return obj_->extrinsic_driver_evaluate(driver_name, u); } diff --git a/src/interpolator.cpp b/src/interpolator.cpp index 05e64f78a..8864314c8 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -18,8 +18,7 @@ void Interpolator::init(const std::vector& x_, initialise(); } -// Compute the interpolated function from the points contained in 'x' -// and 'y'. +// Compute the interpolated function from the points contained in 'x' and 'y'. void Interpolator::initialise() { // https://stackoverflow.com/questions/17769114/stdis-sorted-and-strictly-less-comparison if (not std::is_sorted(x.begin(), x.end(), std::less_equal())) { @@ -55,7 +54,17 @@ void Interpolator::clear() { // Compute the value of the interpolated function at point `x=u` double Interpolator::eval(double u) const { - check_active(); // slower to check this every time + check_active(); + if (not extrapolate and (u < min() or u > max())) { + auto err = std::ostringstream{}; + err << "Extrapolation disabled and evaluation point of " << u << " outside of interpolated domain."; + util::stop(err.str()); + } + return tk_spline(u); +} + +// faster version of above +double Interpolator::operator()(double u) const { return tk_spline(u); } @@ -75,6 +84,10 @@ double Interpolator::max() const { return size() > 0 ? x.back() : R_NegInf; } +void Interpolator::setExtrapolate(bool e) { + extrapolate = e; +} + std::vector Interpolator::get_x() const { return x; } @@ -94,6 +107,7 @@ SEXP Interpolator::r_get_xy() const { // Compute the value of the interpolated function at a vector of // points `x=u`, returning a vector of the same length. std::vector Interpolator::r_eval(std::vector u) const { + check_active(); auto ret = std::vector(); ret.reserve(u.size()); // fast to do this once rather than multiple times with push_back for (auto const& x : u) { From b1e2a86ca5bdcf26c3a55cb6e9e8bee7ba46e456 Mon Sep 17 00:00:00 2001 From: Mitchell Henry Date: Thu, 9 Dec 2021 12:17:21 +1100 Subject: [PATCH 17/22] trying with .at() rather than operator[] --- inst/include/plant/environment.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h index c43689deb..eb162209c 100644 --- a/inst/include/plant/environment.h +++ b/inst/include/plant/environment.h @@ -80,8 +80,8 @@ class Environment { if (extrinsic_drivers.count(driver_name) == 0) { util::stop(driver_name + " doesn't exist in the list of extrinsic_drivers."); } else { - extrinsic_drivers[driver_name].init(x, y); - extrinsic_drivers[driver_name].setExtrapolate(false); // default no extrapolation + extrinsic_drivers.at(driver_name).init(x, y); + extrinsic_drivers.at(driver_name).setExtrapolate(false); // default no extrapolation } } From 0c3e2661d8a6f746355b82229d0951aef7a311d7 Mon Sep 17 00:00:00 2001 From: aornugent Date: Mon, 13 Dec 2021 09:58:37 +1100 Subject: [PATCH 18/22] Remove rainfall_spline; update Rcpp --- inst/include/plant/models/ff16_environment.h | 3 --- src/RcppExports.cpp | 5 ----- 2 files changed, 8 deletions(-) diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h index 8189e8e16..727d82804 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -89,9 +89,6 @@ class FF16_Environment : public Environment { void clear_environment() { canopy.clear(); } - -private: - interpolator::Interpolator rainfall_spline; }; inline Rcpp::NumericMatrix get_state(const FF16_Environment environment) { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index ed33a1714..5e3ef2683 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -6,11 +6,6 @@ using namespace Rcpp; -#ifdef RCPP_USE_GLOBAL_ROSTREAM -Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); -Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); -#endif - // Lorenz__ctor plant::ode::test::Lorenz Lorenz__ctor(double sigma, double R, double b); RcppExport SEXP _plant_Lorenz__ctor(SEXP sigmaSEXP, SEXP RSEXP, SEXP bSEXP) { From 7facab18f7e5578403176d5c232839ba7f990545 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Mon, 13 Dec 2021 11:38:20 +1100 Subject: [PATCH 19/22] removed stringstream in favour of classic std::string for error handling --- src/interpolator.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/interpolator.cpp b/src/interpolator.cpp index 8864314c8..a0a99ed8d 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -56,9 +56,7 @@ void Interpolator::clear() { double Interpolator::eval(double u) const { check_active(); if (not extrapolate and (u < min() or u > max())) { - auto err = std::ostringstream{}; - err << "Extrapolation disabled and evaluation point of " << u << " outside of interpolated domain."; - util::stop(err.str()); + util::stop("Extrapolation disabled and evaluation point of " + std::to_string(u) + " outside of interpolated domain."); } return tk_spline(u); } From ffda217700084682d73945d67cbe5696f5976d5a Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Mon, 13 Dec 2021 12:44:42 +1100 Subject: [PATCH 20/22] more useful error messages for extrapolation --- src/interpolator.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/interpolator.cpp b/src/interpolator.cpp index a0a99ed8d..19e76931b 100644 --- a/src/interpolator.cpp +++ b/src/interpolator.cpp @@ -56,7 +56,11 @@ void Interpolator::clear() { double Interpolator::eval(double u) const { check_active(); if (not extrapolate and (u < min() or u > max())) { - util::stop("Extrapolation disabled and evaluation point of " + std::to_string(u) + " outside of interpolated domain."); + auto err = std::string{"Extrapolation disabled and evaluation point of "}; + err += std::to_string(u); + err += " outside of interpolated domain: "; + err += "[" + std::to_string(min()) + ", " + std::to_string(min()) + "]"; + util::stop(err); } return tk_spline(u); } From e38a6231c9e3cc9df1bbed10bfcd89a002f8274e Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Mon, 13 Dec 2021 13:05:52 +1100 Subject: [PATCH 21/22] adding return statement in evaluate --- inst/include/plant/environment.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h index eb162209c..73d67995f 100644 --- a/inst/include/plant/environment.h +++ b/inst/include/plant/environment.h @@ -91,7 +91,7 @@ class Environment { // 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 { - extrinsic_drivers.at(driver_name).eval(u); + return extrinsic_drivers.at(driver_name).eval(u); } // // evaluate/query interpolated spline for driver at vector of points, return vector of values From 846c0d8d12b2337c22976a0154371d7eb5dd1af6 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Mon, 13 Dec 2021 13:52:56 +1100 Subject: [PATCH 22/22] Removed comments + minor refactor --- inst/include/plant/environment.h | 4 +--- inst/include/plant/models/ff16_environment.h | 2 +- tests/testthat/test-strategy-ff16w.R | 1 - 3 files changed, 2 insertions(+), 5 deletions(-) diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h index 73d67995f..c5edea29f 100644 --- a/inst/include/plant/environment.h +++ b/inst/include/plant/environment.h @@ -12,7 +12,6 @@ #include using namespace Rcpp; -// interpolator has its own namespace, doing this for cleaner code namespace plant { @@ -81,7 +80,7 @@ class Environment { util::stop(driver_name + " doesn't exist in the list of extrinsic_drivers."); } else { extrinsic_drivers.at(driver_name).init(x, y); - extrinsic_drivers.at(driver_name).setExtrapolate(false); // default no extrapolation + extrinsic_driver_extrapolation(driver_name, false); // default no extrapolation } } @@ -107,7 +106,6 @@ class Environment { return ret; } -// private? protected: std::unordered_map extrinsic_drivers; }; diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h index 727d82804..563e84240 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -57,7 +57,7 @@ class FF16_Environment : public Environment { virtual void compute_rates() { for (size_t i = 0; i < vars.state_size; i++) { - vars.set_rate(i, extrinsic_drivers["rainfall"].eval(time) / (i+1)); + vars.set_rate(i, extrinsic_driver_evaluate("rainfall", time) / (i+1)); } } diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R index fef2d8d36..bd6cd8d90 100644 --- a/tests/testthat/test-strategy-ff16w.R +++ b/tests/testthat/test-strategy-ff16w.R @@ -24,7 +24,6 @@ test_that("FF16w Environment", { e$compute_rates() expect_equal(e$soil$rates, 5) - # not needed anymore? should we still be able to set the states manually? # Water logged e$set_soil_water_state(100) expect_equal(e$soil$states, 100)