From 9365267f001273a8c93923ecb39a21e812cdddfe Mon Sep 17 00:00:00 2001 From: Theresa O'Brien Date: Thu, 4 Apr 2024 15:25:21 +1100 Subject: [PATCH] Constant (#10) Finalised the testing structure for models with single and multiple individuals, using the constant growth model as the test case. - Generating a dataset to be loaded with the model in the fixtures. - Running the chosen model on the first individual in the simulated data and test that the parameter estimate converges to a reasonable value (with 0.1 of the true value) - Running the chosen model on multiple individuals and tests whether the output samples are of the right size for the number of chains, iterations, and parameters. * Notation (#5) * Updated notation and function structure in the multi individual constant model file for consistency. Updated rmot_run, rmot_models to reflect the changes. * Removed old constant stan file and cleared out references to it in stanmodels.R, rmot.rmd. * Updated ignore to exclude compiled binaries. * Cleared out .o and .so files. * Added develop branch as a trigger for PR --------- Co-authored-by: Fonti Kar * Added single ind const model, updated rmot_models, rmot_run, rmot vignette. Need to do test files yet. * Updated the single individual constant model stan file and the testing code. * Added some comments to the model testing script. * Expending tests: testing run * Reorganise if statements * Rename file to linear * Added unit testing of model output for linear data. * Re-added testing structure. * Found problem with unit testing for const model outputs: object y_single/y_multi in global environment not being seen within rmot_assign_data function. * Added internal testing data so no data is created within test_that() functions, y_single and y_multi are not evaluated globally #1 * Added internal testing data so no data is created within test_that() functions, y_single and y_multi are not evaluated globally #1 * Skip snapshots on CI * Updated snaps * Removed snaps as a testing framework and using expect_equal #1 * Removed skip_on_ci #1 * Updating roxygen version and package doc * Updated code to generate testing data * Moved helper code in generating fake data, removed internal testing data * Updated rmot_assign and tests passing I think * Completed testing workflow * Added filepath to saving rds and updated TO DO in assign_data * Set tolerance * Updated version for checkout * Added constant data generation and tests based on model output summary for single and multiple individuals. Tests all complete dynamically and when run in the testing console. * Init comit for constant branch * Tess helper funcs * Regenerate test data using stan seed method * using rstan::extract instead of summary * Added constant as trigger * Moved unique code back to make_constant.R * Changed test data generation and single- and multi-individual testing structures for constant model. * Fixed linear regression model testing. * Resolving issues raised in pull request review: changed the set-up of the constant model and removed reference to the constant branch from R-CMD-checn.yaml. --------- Co-authored-by: Fonti Kar Co-authored-by: Daniel Falster --- .Rbuildignore | 3 + .github/workflows/R-CMD-check.yaml | 4 +- .github/workflows/test-coverage.yaml | 2 +- .gitignore | 7 + DESCRIPTION | 5 +- R/rmot-package.R | 1 + R/rmot_assign_data.R | 7 +- R/rmot_models.R | 34 +- R/rmot_run.R | 3 +- R/stanmodels.R | 5 +- inst/stan/constant_multi_ind.stan | 96 +++ inst/stan/constant_single.stan | 83 --- inst/stan/constant_single_ind.stan | 79 +++ man/rmot-package.Rd | 15 + rmot.Rproj | 1 - src/RcppExports.cpp | 6 +- ...e.cc => stanExports_constant_multi_ind.cc} | 6 +- ...gle.h => stanExports_constant_multi_ind.h} | 568 +++++++-------- src/stanExports_constant_single_ind.cc | 32 + src/stanExports_constant_single_ind.h | 657 ++++++++++++++++++ .../constant/constant_data_multi_ind.rds | Bin 0 -> 814 bytes .../constant/constant_data_single_ind.rds | Bin 0 -> 471 bytes .../fixtures/constant/make_constant.R | 39 ++ .../fixtures/linear/lm_baseline_output.rds | Bin 0 -> 3187 bytes tests/testthat/fixtures/linear/lm_data.rds | Bin 0 -> 547 bytes tests/testthat/fixtures/linear/make_lm.R | 22 + tests/testthat/helper-data_generation.R | 99 +++ tests/testthat/helper-testing.R | 72 ++ tests/testthat/test-rmot_models.R | 19 - tests/testthat/test-rmot_models_constant.R | 36 + tests/testthat/test-rmot_models_linearreg.R | 23 + vignettes/rmot.Rmd | 49 +- 32 files changed, 1555 insertions(+), 418 deletions(-) create mode 100644 inst/stan/constant_multi_ind.stan delete mode 100644 inst/stan/constant_single.stan create mode 100644 inst/stan/constant_single_ind.stan rename src/{stanExports_constant_single.cc => stanExports_constant_multi_ind.cc} (93%) rename src/{stanExports_constant_single.h => stanExports_constant_multi_ind.h} (63%) create mode 100644 src/stanExports_constant_single_ind.cc create mode 100644 src/stanExports_constant_single_ind.h create mode 100644 tests/testthat/fixtures/constant/constant_data_multi_ind.rds create mode 100644 tests/testthat/fixtures/constant/constant_data_single_ind.rds create mode 100644 tests/testthat/fixtures/constant/make_constant.R create mode 100644 tests/testthat/fixtures/linear/lm_baseline_output.rds create mode 100644 tests/testthat/fixtures/linear/lm_data.rds create mode 100644 tests/testthat/fixtures/linear/make_lm.R create mode 100644 tests/testthat/helper-data_generation.R create mode 100644 tests/testthat/helper-testing.R delete mode 100644 tests/testthat/test-rmot_models.R create mode 100644 tests/testthat/test-rmot_models_constant.R create mode 100644 tests/testthat/test-rmot_models_linearreg.R diff --git a/.Rbuildignore b/.Rbuildignore index 74b0c8d..3c8ccfe 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -6,3 +6,6 @@ ^LICENSE\.md$ ^\.github$ ^codecov\.yml$ +^data-raw$ +^tests/testthat/_snaps/$ + diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index a3ac618..c60649c 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -4,7 +4,7 @@ on: push: branches: [main, master] pull_request: - branches: [main, master] + branches: [main, master, develop] name: R-CMD-check @@ -29,7 +29,7 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 13955f1..707f4ad 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -4,7 +4,7 @@ on: push: branches: [main, master] pull_request: - branches: [main, master] + branches: [main, master, develop] name: test-coverage diff --git a/.gitignore b/.gitignore index 0042aee..92d552e 100644 --- a/.gitignore +++ b/.gitignore @@ -51,3 +51,10 @@ rsconnect/ inst/doc /doc/ /Meta/ + +# Pre-compiled stan files +*.o +*.so + +# Bugged directory from snapshot package +tests/testthat/_snaps diff --git a/DESCRIPTION b/DESCRIPTION index f238155..8d3253a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -12,7 +12,7 @@ Description: What the package does (one paragraph). License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 Biarch: true Depends: R (>= 3.4.0) @@ -36,6 +36,7 @@ SystemRequirements: GNU make Suggests: knitr, rmarkdown, - testthat (>= 3.0.0) + testthat (>= 3.0.0), + withr VignetteBuilder: knitr Config/testthat/edition: 3 diff --git a/R/rmot-package.R b/R/rmot-package.R index 24722e2..e0619c3 100644 --- a/R/rmot-package.R +++ b/R/rmot-package.R @@ -15,4 +15,5 @@ #' @references #' Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.26.23. https://mc-stan.org #' +"_PACKAGE" NULL diff --git a/R/rmot_assign_data.R b/R/rmot_assign_data.R index 9e62408..9a0ae81 100644 --- a/R/rmot_assign_data.R +++ b/R/rmot_assign_data.R @@ -10,7 +10,7 @@ #' rmot_model("linear") |> rmot_assign_data(X = Loblolly$age, Y = Loblolly$height) rmot_assign_data <- function(model_template, ...){ # Grab user expressions - user_code <- rlang::enexprs(..., .check_assign = TRUE) + user_code <- rlang::enquos(..., .check_assign = TRUE) # Grab the names fields <- names(user_code) @@ -19,13 +19,14 @@ rmot_assign_data <- function(model_template, ...){ # Evaluate the RHS of expressions (the values) data <- purrr::map(user_code, - eval) + ~rlang::eval_tidy(.x, env = rlang::caller_env()) + ) for(i in fields){ model_template <- purrr::list_modify(model_template, !!!data[i]) } - #TODO: Check if N is supplied, if not, assign by default to length(X) and give warning + #TODO: Check if N is supplied, if not, give error return(model_template) } diff --git a/R/rmot_models.R b/R/rmot_models.R index 512fe5b..13e0758 100644 --- a/R/rmot_models.R +++ b/R/rmot_models.R @@ -14,7 +14,8 @@ rmot_model <- function(model=NULL){ output <- switch(model, linear = rmot_lm(), - constant_single = rmot_cgs()) + constant_single_ind = rmot_const_single_ind(), + constant_multi_ind = rmot_const_multi_ind()) class(output) <- "rmot_object" @@ -32,19 +33,32 @@ rmot_lm <- function(){ model = "linear") } +#' Data configuration template for constant growth single individual model +#' @keywords internal +#' @noRd + +rmot_const_single_ind <- function(){ + list(n_obs = NULL, + y_obs = NULL, + obs_index = NULL, + time = NULL, + y_0_obs = NULL, + model = "constant_single_ind") +} + #' Data configuration template for constant growth single species model #' @keywords internal #' @noRd -rmot_cgs <- function(){ - list(N_obs = NULL, - N_ind = NULL, - S_obs = NULL, - census = NULL, - census_interval = NULL, - id_factor = NULL, - S_0_obs = NULL, - model = "constant_single") +rmot_const_multi_ind <- function(){ + list(n_obs = NULL, + n_ind = NULL, + y_obs = NULL, + obs_index = NULL, + time = NULL, + ind_id = NULL, + y_0_obs = NULL, + model = "constant_multi_ind") } diff --git a/R/rmot_run.R b/R/rmot_run.R index 2e0cfa1..193f379 100644 --- a/R/rmot_run.R +++ b/R/rmot_run.R @@ -17,7 +17,8 @@ rmot_run <- function(model_template, ...) { # Detect model out <- switch(model_template$model, linear = rstan::sampling(stanmodels$linear, data = model_template, ...), - constant_single = rstan::sampling(stanmodels$constant_single, data = model_template, ...)) + constant_single_ind = rstan::sampling(stanmodels$constant_single_ind, data = model_template, ...), + constant_multi_ind = rstan::sampling(stanmodels$constant_multi_ind, data = model_template, ...)) return(out) } diff --git a/R/stanmodels.R b/R/stanmodels.R index 5a58cd8..e5e4f0b 100644 --- a/R/stanmodels.R +++ b/R/stanmodels.R @@ -1,10 +1,11 @@ # Generated by rstantools. Do not edit by hand. # names of stan models -stanmodels <- c("constant_single", "linear") +stanmodels <- c("constant_multi_ind", "constant_single_ind", "linear") # load each stan module -Rcpp::loadModule("stan_fit4constant_single_mod", what = TRUE) +Rcpp::loadModule("stan_fit4constant_multi_ind_mod", what = TRUE) +Rcpp::loadModule("stan_fit4constant_single_ind_mod", what = TRUE) Rcpp::loadModule("stan_fit4linear_mod", what = TRUE) # instantiate each stanmodel object diff --git a/inst/stan/constant_multi_ind.stan b/inst/stan/constant_multi_ind.stan new file mode 100644 index 0000000..b84903a --- /dev/null +++ b/inst/stan/constant_multi_ind.stan @@ -0,0 +1,96 @@ +//Constant DE - Single species +functions{ + //DE function + real DE(real beta){ + return beta; + } + + real size_step(real y, real beta, real time){ + return y + DE(beta) * time; + } +} + +// Data structure +data { + int n_obs; + int n_ind; + real y_obs[n_obs]; + int obs_index[n_obs]; + real time[n_obs]; + int ind_id[n_obs]; + real y_0_obs[n_ind]; +} + +// The parameters accepted by the model. +parameters { + //Individual level + real ind_y_0[n_ind]; + real ind_beta[n_ind]; + + real species_beta_mu; + real species_beta_sigma; + + //Global level + real global_error_sigma; +} + +// The model to be estimated. +model { + real y_hat[n_obs]; + + for(i in 1:n_obs){ + //Fits the first size + if(obs_index[i]==1){ + y_hat[i] = ind_y_0[ind_id[i]]; + } + + // Estimate next size + if(i < n_obs){ + if(ind_id[i+1]==ind_id[i]){ + y_hat[i+1] = size_step(y_hat[i], ind_beta[ind_id[i]], (time[i+1]-time[i])); + } + } + } + + //Likelihood + y_obs ~ normal(y_hat, global_error_sigma); + + //Priors + //Individual level + ind_y_0 ~ normal(y_0_obs, global_error_sigma); + ind_beta ~ lognormal(species_beta_mu, + species_beta_sigma); + + //Species level + species_beta_mu ~ normal(0.1, 1); + species_beta_sigma ~cauchy(0.1, 1); + + //Global level + global_error_sigma ~cauchy(0.1, 1); +} + +// The output +generated quantities { + real y_hat[n_obs]; + real Delta_hat[n_obs]; + + for(i in 1:n_obs){ + + //Fits the first size + if(obs_index[i]==1){ + y_hat[i] = ind_y_0[ind_id[i]]; + } + + // Estimate next size + if(i < n_obs){ + if(ind_id[i+1]==ind_id[i]){ + y_hat[i+1] = size_step(y_hat[i], ind_beta[ind_id[i]], (time[i+1]-time[i])); + Delta_hat[i] = y_hat[i+1] - y_hat[i]; + } else { + Delta_hat[i] = DE(ind_beta[ind_id[i]]) * (time[i]-time[i-1]); + } + } else { + Delta_hat[i] = DE(ind_beta[ind_id[i]]) * (time[i]-time[i-1]); + } + } +} diff --git a/inst/stan/constant_single.stan b/inst/stan/constant_single.stan deleted file mode 100644 index 24a984c..0000000 --- a/inst/stan/constant_single.stan +++ /dev/null @@ -1,83 +0,0 @@ -//Constant Growth - Single species - -// Data structure -data { - int N_obs; - int N_ind; - real S_obs[N_obs]; - int census[N_obs]; - real census_interval[N_obs]; - int id_factor[N_obs]; - real S_0_obs[N_ind]; -} - -// The parameters accepted by the model. -parameters { - //Individual level - real ind_S_0[N_ind]; - real ind_beta[N_ind]; - - real species_beta_mu; - real species_beta_sigma; - - //Global level - real global_error_sigma; -} - -// The model to be estimated. -model { - real S_hat[N_obs]; - real G_hat[N_obs]; - - for(i in 1:N_obs){ - if(id_factor[i+1]==id_factor[i]){ - if(census[i]==1){//Fits the first size - S_hat[i] = ind_S_0[id_factor[i]]; - } - - if(i < N_obs){ //Analytic solution - G_hat[i] = ind_beta[id_factor[i]]; - S_hat[i+1] = S_hat[i] + G_hat[i]*census_interval[i+1]; - } - } else { - G_hat[i] = 0; //Gives 0 as the growth estimate for the last data point. - } - } - - //Likelihood - S_obs ~ normal(S_hat, global_error_sigma); - - //Priors - //Individual level - ind_S_0 ~ normal(S_0_obs, global_error_sigma); - ind_beta ~ lognormal(species_beta_mu, - species_beta_sigma); - - //Species level - species_beta_mu ~ normal(0.1, 1); - species_beta_sigma ~cauchy(0.1, 1); - - //Global level - global_error_sigma ~cauchy(0.1, 1); -} - -// The output -generated quantities { - real S_hat[N_obs]; - real G_hat[N_obs]; - - for(i in 1:N_obs){ - if(id_factor[i+1]==id_factor[i]){ - if(census[i]==1){//Fits the first size - S_hat[i] = ind_S_0[id_factor[i]]; - } - - if(i < N_obs){ //Analytic solution - G_hat[i] = ind_beta[id_factor[i]]; - S_hat[i+1] = S_hat[i] + G_hat[i]*census_interval[i+1]; - } - } else { - G_hat[i] = 0; //Gives 0 as the growth estimate for the last data point. - } - } -} diff --git a/inst/stan/constant_single_ind.stan b/inst/stan/constant_single_ind.stan new file mode 100644 index 0000000..cd7570a --- /dev/null +++ b/inst/stan/constant_single_ind.stan @@ -0,0 +1,79 @@ +//Constant DE - Single species +functions{ + //DE function + real DE(real beta){ + return beta; + } + + real size_step(real y, real beta, real time){ + return y + DE(beta) * time; + } +} + +// Data structure +data { + int n_obs; + real y_obs[n_obs]; + int obs_index[n_obs]; + real time[n_obs]; + real y_0_obs; +} + +// The parameters accepted by the model. +parameters { + //Individual level + real ind_y_0; + real ind_beta; + + //Global level + real global_error_sigma; +} + +// The model to be estimated. +model { + real y_hat[n_obs]; + + for(i in 1:n_obs){ + //Fits the first size + if(obs_index[i]==1){ + y_hat[i] = ind_y_0; + } + + // Estimate next size + if(i < n_obs){ + y_hat[i+1] = size_step(y_hat[i], ind_beta, (time[i+1]-time[i])); + } + } + + //Likelihood + y_obs ~ normal(y_hat, global_error_sigma); + + //Priors + //Individual level + ind_y_0 ~ normal(y_0_obs, global_error_sigma); + ind_beta ~ lognormal(0.1, 1); + + //Global level + global_error_sigma ~cauchy(0.1, 1); +} + +// The output +generated quantities { + real y_hat[n_obs]; + real Delta_hat[n_obs]; + + for(i in 1:n_obs){ + //Fits the first size + if(obs_index[i]==1){ + y_hat[i] = ind_y_0; + } + + // Estimate next size + if(i < n_obs){ + y_hat[i+1] = size_step(y_hat[i], ind_beta, (time[i+1]-time[i])); + Delta_hat[i] = y_hat[i+1] - y_hat[i]; + } else { + Delta_hat[i] = DE(ind_beta) * (time[i]-time[i-1]); + } + } +} diff --git a/man/rmot-package.Rd b/man/rmot-package.Rd index 724b6b7..fcd3549 100644 --- a/man/rmot-package.Rd +++ b/man/rmot-package.Rd @@ -11,3 +11,18 @@ A DESCRIPTION OF THE PACKAGE \references{ Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.26.23. https://mc-stan.org } +\author{ +\strong{Maintainer}: Tess O'Brien \email{theresa.obrien@unsw.edu.au} (\href{https://orcid.org/XXXX-XXXX-XXXX-XXXX}{ORCID}) [copyright holder] + +Authors: +\itemize{ + \item Daniel Falster \email{daniel.falster@unsw.edu.au} (\href{https://orcid.org/0000-0002-9814-092X}{ORCID}) [contributor] + \item David Warton \email{david.warton@unsw.edu.au} (\href{https://orcid.org/0000-0002-1642-628X}{ORCID}) [contributor] +} + +Other contributors: +\itemize{ + \item Fonti Kar \email{f.kar@unsw.edu.au} (\href{https://orcid.org/0000-0002-2760-3974}{ORCID}) [contributor] +} + +} diff --git a/rmot.Rproj b/rmot.Rproj index 69fafd4..6a3ede2 100644 --- a/rmot.Rproj +++ b/rmot.Rproj @@ -19,4 +19,3 @@ LineEndingConversion: Posix BuildType: Package PackageUseDevtools: Yes PackageInstallArgs: --no-multiarch --with-keep.source -PackageRoxygenize: rd,collate,namespace diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index c924a64..6a109cd 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,11 +12,13 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif -RcppExport SEXP _rcpp_module_boot_stan_fit4constant_single_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4constant_multi_ind_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4constant_single_ind_mod(); RcppExport SEXP _rcpp_module_boot_stan_fit4linear_mod(); static const R_CallMethodDef CallEntries[] = { - {"_rcpp_module_boot_stan_fit4constant_single_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4constant_single_mod, 0}, + {"_rcpp_module_boot_stan_fit4constant_multi_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4constant_multi_ind_mod, 0}, + {"_rcpp_module_boot_stan_fit4constant_single_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4constant_single_ind_mod, 0}, {"_rcpp_module_boot_stan_fit4linear_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4linear_mod, 0}, {NULL, NULL, 0} }; diff --git a/src/stanExports_constant_single.cc b/src/stanExports_constant_multi_ind.cc similarity index 93% rename from src/stanExports_constant_single.cc rename to src/stanExports_constant_multi_ind.cc index 961bd38..a3ad589 100644 --- a/src/stanExports_constant_single.cc +++ b/src/stanExports_constant_multi_ind.cc @@ -2,12 +2,12 @@ #include using namespace Rcpp ; -#include "stanExports_constant_single.h" +#include "stanExports_constant_multi_ind.h" -RCPP_MODULE(stan_fit4constant_single_mod) { +RCPP_MODULE(stan_fit4constant_multi_ind_mod) { - class_ >("rstantools_model_constant_single") + class_ >("rstantools_model_constant_multi_ind") .constructor() diff --git a/src/stanExports_constant_single.h b/src/stanExports_constant_multi_ind.h similarity index 63% rename from src/stanExports_constant_single.h rename to src/stanExports_constant_multi_ind.h index a617591..d5dede6 100644 --- a/src/stanExports_constant_single.h +++ b/src/stanExports_constant_multi_ind.h @@ -23,7 +23,7 @@ #include // Code generated by stanc v2.26.1-4-gd72b68b7-dirty #include -namespace model_constant_single_namespace { +namespace model_constant_multi_ind_namespace { inline void validate_positive_index(const char* var_name, const char* expr, int val) { if (val < 1) { @@ -75,93 +75,118 @@ using stan::math::pow; stan::math::profile_map profiles__; static int current_statement__= 0; static const std::vector locations_array__ = {" (found before start of program)", - " (in 'constant_single', line 15, column 2 to column 31)", - " (in 'constant_single', line 16, column 2 to column 32)", - " (in 'constant_single', line 17, column 2 to column 23)", - " (in 'constant_single', line 18, column 2 to column 35)", - " (in 'constant_single', line 20, column 2 to column 35)", - " (in 'constant_single', line 54, column 2 to column 20)", - " (in 'constant_single', line 55, column 2 to column 20)", - " (in 'constant_single', line 66, column 6 to column 19)", - " (in 'constant_single', line 65, column 11 to line 67, column 5)", - " (in 'constant_single', line 59, column 8 to column 41)", - " (in 'constant_single', line 58, column 22 to line 60, column 7)", - " (in 'constant_single', line 58, column 6 to line 60, column 7)", - " (in 'constant_single', line 62, column 8 to column 42)", - " (in 'constant_single', line 63, column 8 to column 62)", - " (in 'constant_single', line 61, column 19 to line 64, column 7)", - " (in 'constant_single', line 61, column 6 to line 64, column 7)", - " (in 'constant_single', line 57, column 36 to line 65, column 5)", - " (in 'constant_single', line 57, column 4 to line 67, column 5)", - " (in 'constant_single', line 56, column 19 to line 68, column 3)", - " (in 'constant_single', line 56, column 2 to line 68, column 3)", - " (in 'constant_single', line 24, column 13 to column 18)", - " (in 'constant_single', line 24, column 2 to column 20)", - " (in 'constant_single', line 25, column 13 to column 18)", - " (in 'constant_single', line 25, column 2 to column 20)", - " (in 'constant_single', line 36, column 6 to column 19)", - " (in 'constant_single', line 35, column 11 to line 37, column 5)", - " (in 'constant_single', line 29, column 8 to column 41)", - " (in 'constant_single', line 28, column 22 to line 30, column 7)", - " (in 'constant_single', line 28, column 6 to line 30, column 7)", - " (in 'constant_single', line 32, column 8 to column 42)", - " (in 'constant_single', line 33, column 8 to column 62)", - " (in 'constant_single', line 31, column 19 to line 34, column 7)", - " (in 'constant_single', line 31, column 6 to line 34, column 7)", - " (in 'constant_single', line 27, column 36 to line 35, column 5)", - " (in 'constant_single', line 27, column 4 to line 37, column 5)", - " (in 'constant_single', line 26, column 19 to line 38, column 3)", - " (in 'constant_single', line 26, column 2 to line 38, column 3)", - " (in 'constant_single', line 40, column 2 to column 44)", - " (in 'constant_single', line 43, column 2 to column 48)", - " (in 'constant_single', line 44, column 2 to line 45, column 40)", - " (in 'constant_single', line 47, column 2 to column 35)", - " (in 'constant_single', line 48, column 2 to column 37)", - " (in 'constant_single', line 50, column 2 to column 37)", - " (in 'constant_single', line 4, column 2 to column 12)", - " (in 'constant_single', line 5, column 2 to column 12)", - " (in 'constant_single', line 6, column 13 to column 18)", - " (in 'constant_single', line 6, column 2 to column 20)", - " (in 'constant_single', line 7, column 13 to column 18)", - " (in 'constant_single', line 7, column 2 to column 20)", - " (in 'constant_single', line 8, column 23 to column 28)", - " (in 'constant_single', line 8, column 2 to column 30)", - " (in 'constant_single', line 9, column 16 to column 21)", - " (in 'constant_single', line 9, column 2 to column 23)", - " (in 'constant_single', line 10, column 15 to column 20)", - " (in 'constant_single', line 10, column 2 to column 22)", - " (in 'constant_single', line 15, column 24 to column 29)", - " (in 'constant_single', line 16, column 25 to column 30)", - " (in 'constant_single', line 54, column 13 to column 18)", - " (in 'constant_single', line 55, column 13 to column 18)"}; + " (in 'constant_multi_ind', line 21, column 2 to column 31)", + " (in 'constant_multi_ind', line 22, column 2 to column 32)", + " (in 'constant_multi_ind', line 23, column 2 to column 23)", + " (in 'constant_multi_ind', line 24, column 2 to column 35)", + " (in 'constant_multi_ind', line 26, column 2 to column 35)", + " (in 'constant_multi_ind', line 62, column 2 to column 20)", + " (in 'constant_multi_ind', line 63, column 2 to column 24)", + " (in 'constant_multi_ind', line 67, column 6 to column 36)", + " (in 'constant_multi_ind', line 66, column 23 to line 68, column 5)", + " (in 'constant_multi_ind', line 66, column 4 to line 68, column 5)", + " (in 'constant_multi_ind', line 70, column 4 to column 47)", + " (in 'constant_multi_ind', line 75, column 8 to column 65)", + " (in 'constant_multi_ind', line 73, column 32 to line 76, column 7)", + " (in 'constant_multi_ind', line 73, column 6 to line 76, column 7)", + " (in 'constant_multi_ind', line 72, column 17 to line 77, column 5)", + " (in 'constant_multi_ind', line 72, column 4 to line 77, column 5)", + " (in 'constant_multi_ind', line 64, column 19 to line 78, column 3)", + " (in 'constant_multi_ind', line 64, column 2 to line 78, column 3)", + " (in 'constant_multi_ind', line 30, column 13 to column 18)", + " (in 'constant_multi_ind', line 30, column 2 to column 20)", + " (in 'constant_multi_ind', line 31, column 17 to column 22)", + " (in 'constant_multi_ind', line 31, column 2 to column 24)", + " (in 'constant_multi_ind', line 35, column 6 to column 36)", + " (in 'constant_multi_ind', line 34, column 23 to line 36, column 5)", + " (in 'constant_multi_ind', line 34, column 4 to line 36, column 5)", + " (in 'constant_multi_ind', line 38, column 4 to column 47)", + " (in 'constant_multi_ind', line 43, column 8 to column 65)", + " (in 'constant_multi_ind', line 41, column 32 to line 44, column 7)", + " (in 'constant_multi_ind', line 41, column 6 to line 44, column 7)", + " (in 'constant_multi_ind', line 40, column 17 to line 45, column 5)", + " (in 'constant_multi_ind', line 40, column 4 to line 45, column 5)", + " (in 'constant_multi_ind', line 32, column 19 to line 46, column 3)", + " (in 'constant_multi_ind', line 32, column 2 to line 46, column 3)", + " (in 'constant_multi_ind', line 48, column 2 to column 44)", + " (in 'constant_multi_ind', line 51, column 2 to column 48)", + " (in 'constant_multi_ind', line 52, column 2 to line 53, column 40)", + " (in 'constant_multi_ind', line 55, column 2 to column 35)", + " (in 'constant_multi_ind', line 56, column 2 to column 37)", + " (in 'constant_multi_ind', line 58, column 2 to column 37)", + " (in 'constant_multi_ind', line 10, column 2 to column 12)", + " (in 'constant_multi_ind', line 11, column 2 to column 12)", + " (in 'constant_multi_ind', line 12, column 13 to column 18)", + " (in 'constant_multi_ind', line 12, column 2 to column 20)", + " (in 'constant_multi_ind', line 13, column 16 to column 21)", + " (in 'constant_multi_ind', line 13, column 2 to column 23)", + " (in 'constant_multi_ind', line 14, column 12 to column 17)", + " (in 'constant_multi_ind', line 14, column 2 to column 19)", + " (in 'constant_multi_ind', line 15, column 13 to column 18)", + " (in 'constant_multi_ind', line 15, column 2 to column 20)", + " (in 'constant_multi_ind', line 16, column 15 to column 20)", + " (in 'constant_multi_ind', line 16, column 2 to column 22)", + " (in 'constant_multi_ind', line 21, column 24 to column 29)", + " (in 'constant_multi_ind', line 22, column 25 to column 30)", + " (in 'constant_multi_ind', line 62, column 13 to column 18)", + " (in 'constant_multi_ind', line 63, column 17 to column 22)", + " (in 'constant_multi_ind', line 5, column 4 to column 16)", + " (in 'constant_multi_ind', line 4, column 24 to line 6, column 3)"}; +template +stan::promote_args_t +growth(const T0__& beta, std::ostream* pstream__) { + using local_scalar_t__ = stan::promote_args_t; + const static bool propto__ = true; + (void) propto__; + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + current_statement__ = 56; + return beta; + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); + } + +} +struct growth_functor__ { +template +stan::promote_args_t +operator()(const T0__& beta, std::ostream* pstream__) const +{ +return growth(beta, pstream__); +} +}; #include -class model_constant_single final : public model_base_crtp { +class model_constant_multi_ind final : public model_base_crtp { private: - int N_obs; - int N_ind; - std::vector S_obs; - std::vector census; - std::vector census_interval; - std::vector id_factor; - std::vector S_0_obs; + int n_obs; + int n_ind; + std::vector y_obs; + std::vector obs_index; + std::vector time; + std::vector ind_id; + std::vector y_0_obs; public: - ~model_constant_single() { } + ~model_constant_multi_ind() { } - inline std::string model_name() const final { return "model_constant_single"; } + inline std::string model_name() const final { return "model_constant_multi_ind"; } inline std::vector model_compile_info() const noexcept { return std::vector{"stanc_version = stanc3 v2.26.1-4-gd72b68b7-dirty", "stancflags = "}; } - model_constant_single(stan::io::var_context& context__, - unsigned int random_seed__ = 0, - std::ostream* pstream__ = nullptr) : model_base_crtp(0) { + model_constant_multi_ind(stan::io::var_context& context__, + unsigned int random_seed__ = 0, + std::ostream* pstream__ = nullptr) : model_base_crtp(0) { using local_scalar_t__ = double ; boost::ecuyer1988 base_rng__ = stan::services::util::create_rng(random_seed__, 0); (void) base_rng__; // suppress unused var warning - static const char* function__ = "model_constant_single_namespace::model_constant_single"; + static const char* function__ = "model_constant_multi_ind_namespace::model_constant_multi_ind"; (void) function__; // suppress unused var warning local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); (void) DUMMY_VAR__; // suppress unused var warning @@ -171,79 +196,78 @@ class model_constant_single final : public model_base_crtp::min(); pos__ = 1; - current_statement__ = 44; - context__.validate_dims("data initialization","N_obs","int", + current_statement__ = 40; + context__.validate_dims("data initialization","n_obs","int", context__.to_vec()); - N_obs = std::numeric_limits::min(); + n_obs = std::numeric_limits::min(); + current_statement__ = 40; + n_obs = context__.vals_i("n_obs")[(1 - 1)]; + current_statement__ = 41; + context__.validate_dims("data initialization","n_ind","int", + context__.to_vec()); + n_ind = std::numeric_limits::min(); + + current_statement__ = 41; + n_ind = context__.vals_i("n_ind")[(1 - 1)]; + current_statement__ = 42; + validate_non_negative_index("y_obs", "n_obs", n_obs); + current_statement__ = 43; + context__.validate_dims("data initialization","y_obs","double", + context__.to_vec(n_obs)); + y_obs = std::vector(n_obs, std::numeric_limits::quiet_NaN()); + + current_statement__ = 43; + assign(y_obs, nil_index_list(), context__.vals_r("y_obs"), + "assigning variable y_obs"); current_statement__ = 44; - N_obs = context__.vals_i("N_obs")[(1 - 1)]; + validate_non_negative_index("obs_index", "n_obs", n_obs); current_statement__ = 45; - context__.validate_dims("data initialization","N_ind","int", - context__.to_vec()); - N_ind = std::numeric_limits::min(); + context__.validate_dims("data initialization","obs_index","int", + context__.to_vec(n_obs)); + obs_index = std::vector(n_obs, std::numeric_limits::min()); current_statement__ = 45; - N_ind = context__.vals_i("N_ind")[(1 - 1)]; + assign(obs_index, nil_index_list(), context__.vals_i("obs_index"), + "assigning variable obs_index"); current_statement__ = 46; - validate_non_negative_index("S_obs", "N_obs", N_obs); + validate_non_negative_index("time", "n_obs", n_obs); current_statement__ = 47; - context__.validate_dims("data initialization","S_obs","double", - context__.to_vec(N_obs)); - S_obs = std::vector(N_obs, std::numeric_limits::quiet_NaN()); + context__.validate_dims("data initialization","time","double", + context__.to_vec(n_obs)); + time = std::vector(n_obs, std::numeric_limits::quiet_NaN()); current_statement__ = 47; - assign(S_obs, nil_index_list(), context__.vals_r("S_obs"), - "assigning variable S_obs"); + assign(time, nil_index_list(), context__.vals_r("time"), + "assigning variable time"); current_statement__ = 48; - validate_non_negative_index("census", "N_obs", N_obs); + validate_non_negative_index("ind_id", "n_obs", n_obs); current_statement__ = 49; - context__.validate_dims("data initialization","census","int", - context__.to_vec(N_obs)); - census = std::vector(N_obs, std::numeric_limits::min()); + context__.validate_dims("data initialization","ind_id","int", + context__.to_vec(n_obs)); + ind_id = std::vector(n_obs, std::numeric_limits::min()); current_statement__ = 49; - assign(census, nil_index_list(), context__.vals_i("census"), - "assigning variable census"); + assign(ind_id, nil_index_list(), context__.vals_i("ind_id"), + "assigning variable ind_id"); current_statement__ = 50; - validate_non_negative_index("census_interval", "N_obs", N_obs); + validate_non_negative_index("y_0_obs", "n_ind", n_ind); current_statement__ = 51; - context__.validate_dims("data initialization","census_interval", - "double",context__.to_vec(N_obs)); - census_interval = std::vector(N_obs, std::numeric_limits::quiet_NaN()); + context__.validate_dims("data initialization","y_0_obs","double", + context__.to_vec(n_ind)); + y_0_obs = std::vector(n_ind, std::numeric_limits::quiet_NaN()); current_statement__ = 51; - assign(census_interval, nil_index_list(), - context__.vals_r("census_interval"), - "assigning variable census_interval"); + assign(y_0_obs, nil_index_list(), context__.vals_r("y_0_obs"), + "assigning variable y_0_obs"); current_statement__ = 52; - validate_non_negative_index("id_factor", "N_obs", N_obs); - current_statement__ = 53; - context__.validate_dims("data initialization","id_factor","int", - context__.to_vec(N_obs)); - id_factor = std::vector(N_obs, std::numeric_limits::min()); - + validate_non_negative_index("ind_y_0", "n_ind", n_ind); current_statement__ = 53; - assign(id_factor, nil_index_list(), context__.vals_i("id_factor"), - "assigning variable id_factor"); + validate_non_negative_index("ind_beta", "n_ind", n_ind); current_statement__ = 54; - validate_non_negative_index("S_0_obs", "N_ind", N_ind); - current_statement__ = 55; - context__.validate_dims("data initialization","S_0_obs","double", - context__.to_vec(N_ind)); - S_0_obs = std::vector(N_ind, std::numeric_limits::quiet_NaN()); - + validate_non_negative_index("y_hat", "n_obs", n_obs); current_statement__ = 55; - assign(S_0_obs, nil_index_list(), context__.vals_r("S_0_obs"), - "assigning variable S_0_obs"); - current_statement__ = 56; - validate_non_negative_index("ind_S_0", "N_ind", N_ind); - current_statement__ = 57; - validate_non_negative_index("ind_beta", "N_ind", N_ind); - current_statement__ = 58; - validate_non_negative_index("S_hat", "N_obs", N_obs); - current_statement__ = 59; - validate_non_negative_index("G_hat", "N_obs", N_obs); + validate_non_negative_index("Delta_hat", "n_obs", n_obs); } catch (const std::exception& e) { stan::lang::rethrow_located(e, locations_array__[current_statement__]); // Next line prevents compiler griping about no return @@ -252,8 +276,8 @@ class model_constant_single final : public model_base_crtp lp_accum__; - static const char* function__ = "model_constant_single_namespace::log_prob"; + static const char* function__ = "model_constant_multi_ind_namespace::log_prob"; (void) function__; // suppress unused var warning stan::io::reader in__(params_r__, params_i__); local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); (void) DUMMY_VAR__; // suppress unused var warning try { - std::vector ind_S_0; - ind_S_0 = std::vector(N_ind, DUMMY_VAR__); + std::vector ind_y_0; + ind_y_0 = std::vector(n_ind, DUMMY_VAR__); current_statement__ = 1; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 1; - assign(ind_S_0, cons_list(index_uni(sym1__), nil_index_list()), - in__.scalar(), "assigning variable ind_S_0");} + assign(ind_y_0, cons_list(index_uni(sym1__), nil_index_list()), + in__.scalar(), "assigning variable ind_y_0");} current_statement__ = 1; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 1; if (jacobian__) { current_statement__ = 1; - assign(ind_S_0, cons_list(index_uni(sym1__), nil_index_list()), - stan::math::lb_constrain(ind_S_0[(sym1__ - 1)], 0, lp__), - "assigning variable ind_S_0"); + assign(ind_y_0, cons_list(index_uni(sym1__), nil_index_list()), + stan::math::lb_constrain(ind_y_0[(sym1__ - 1)], 0, lp__), + "assigning variable ind_y_0"); } else { current_statement__ = 1; - assign(ind_S_0, cons_list(index_uni(sym1__), nil_index_list()), - stan::math::lb_constrain(ind_S_0[(sym1__ - 1)], 0), - "assigning variable ind_S_0"); + assign(ind_y_0, cons_list(index_uni(sym1__), nil_index_list()), + stan::math::lb_constrain(ind_y_0[(sym1__ - 1)], 0), + "assigning variable ind_y_0"); }} std::vector ind_beta; - ind_beta = std::vector(N_ind, DUMMY_VAR__); + ind_beta = std::vector(n_ind, DUMMY_VAR__); current_statement__ = 2; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 2; assign(ind_beta, cons_list(index_uni(sym1__), nil_index_list()), in__.scalar(), "assigning variable ind_beta");} current_statement__ = 2; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 2; if (jacobian__) { current_statement__ = 2; @@ -356,58 +380,55 @@ class model_constant_single final : public model_base_crtp S_hat; - S_hat = std::vector(N_obs, DUMMY_VAR__); + current_statement__ = 19; + validate_non_negative_index("y_hat", "n_obs", n_obs); + std::vector y_hat; + y_hat = std::vector(n_obs, DUMMY_VAR__); - current_statement__ = 23; - validate_non_negative_index("G_hat", "N_obs", N_obs); - std::vector G_hat; - G_hat = std::vector(N_obs, DUMMY_VAR__); + current_statement__ = 21; + validate_non_negative_index("Delta_hat", "n_obs", n_obs); + std::vector Delta_hat; + Delta_hat = std::vector(n_obs, DUMMY_VAR__); - current_statement__ = 37; - for (int i = 1; i <= N_obs; ++i) { - current_statement__ = 35; - if (logical_eq(id_factor[((i + 1) - 1)], id_factor[(i - 1)])) { + current_statement__ = 33; + for (int i = 1; i <= n_obs; ++i) { + current_statement__ = 25; + if (logical_eq(obs_index[(i - 1)], 1)) { + current_statement__ = 23; + assign(y_hat, cons_list(index_uni(i), nil_index_list()), + ind_y_0[(ind_id[(i - 1)] - 1)], "assigning variable y_hat"); + } + current_statement__ = 26; + assign(Delta_hat, cons_list(index_uni(i), nil_index_list()), + growth(ind_beta[(ind_id[(i - 1)] - 1)], pstream__), + "assigning variable Delta_hat"); + current_statement__ = 31; + if (logical_lt(i, n_obs)) { current_statement__ = 29; - if (logical_eq(census[(i - 1)], 1)) { + if (logical_eq(ind_id[((i + 1) - 1)], ind_id[(i - 1)])) { current_statement__ = 27; - assign(S_hat, cons_list(index_uni(i), nil_index_list()), - ind_S_0[(id_factor[(i - 1)] - 1)], "assigning variable S_hat"); - } - current_statement__ = 33; - if (logical_lt(i, N_obs)) { - current_statement__ = 30; - assign(G_hat, cons_list(index_uni(i), nil_index_list()), - ind_beta[(id_factor[(i - 1)] - 1)], - "assigning variable G_hat"); - current_statement__ = 31; - assign(S_hat, cons_list(index_uni((i + 1)), nil_index_list()), - (S_hat[(i - 1)] + - (G_hat[(i - 1)] * census_interval[((i + 1) - 1)])), - "assigning variable S_hat"); + assign(y_hat, cons_list(index_uni((i + 1)), nil_index_list()), + (y_hat[(i - 1)] + + (Delta_hat[(i - 1)] * + (time[((i + 1) - 1)] - time[(i - 1)]))), + "assigning variable y_hat"); } - } else { - current_statement__ = 25; - assign(G_hat, cons_list(index_uni(i), nil_index_list()), 0, - "assigning variable G_hat"); - }} - current_statement__ = 38; + } } + current_statement__ = 34; lp_accum__.add( - normal_lpdf(S_obs, S_hat, global_error_sigma)); - current_statement__ = 39; + normal_lpdf(y_obs, y_hat, global_error_sigma)); + current_statement__ = 35; lp_accum__.add( - normal_lpdf(ind_S_0, S_0_obs, global_error_sigma)); - current_statement__ = 40; + normal_lpdf(ind_y_0, y_0_obs, global_error_sigma)); + current_statement__ = 36; lp_accum__.add( lognormal_lpdf(ind_beta, species_beta_mu, species_beta_sigma)); - current_statement__ = 41; + current_statement__ = 37; lp_accum__.add(normal_lpdf(species_beta_mu, 0.1, 1)); - current_statement__ = 42; + current_statement__ = 38; lp_accum__.add(cauchy_lpdf(species_beta_sigma, 0.1, 1)); - current_statement__ = 43; + current_statement__ = 39; lp_accum__.add(cauchy_lpdf(global_error_sigma, 0.1, 1)); } } catch (const std::exception& e) { @@ -428,7 +449,7 @@ class model_constant_single final : public model_base_crtp in__(params_r__, params_i__); - static const char* function__ = "model_constant_single_namespace::write_array"; + static const char* function__ = "model_constant_multi_ind_namespace::write_array"; (void) function__; // suppress unused var warning (void) function__; // suppress unused var warning double lp__ = 0.0; @@ -438,30 +459,30 @@ class model_constant_single final : public model_base_crtp ind_S_0; - ind_S_0 = std::vector(N_ind, std::numeric_limits::quiet_NaN()); + std::vector ind_y_0; + ind_y_0 = std::vector(n_ind, std::numeric_limits::quiet_NaN()); current_statement__ = 1; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 1; - assign(ind_S_0, cons_list(index_uni(sym1__), nil_index_list()), - in__.scalar(), "assigning variable ind_S_0");} + assign(ind_y_0, cons_list(index_uni(sym1__), nil_index_list()), + in__.scalar(), "assigning variable ind_y_0");} current_statement__ = 1; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 1; - assign(ind_S_0, cons_list(index_uni(sym1__), nil_index_list()), - stan::math::lb_constrain(ind_S_0[(sym1__ - 1)], 0), - "assigning variable ind_S_0");} + assign(ind_y_0, cons_list(index_uni(sym1__), nil_index_list()), + stan::math::lb_constrain(ind_y_0[(sym1__ - 1)], 0), + "assigning variable ind_y_0");} std::vector ind_beta; - ind_beta = std::vector(N_ind, std::numeric_limits::quiet_NaN()); + ind_beta = std::vector(n_ind, std::numeric_limits::quiet_NaN()); current_statement__ = 2; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 2; assign(ind_beta, cons_list(index_uni(sym1__), nil_index_list()), in__.scalar(), "assigning variable ind_beta");} current_statement__ = 2; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 2; assign(ind_beta, cons_list(index_uni(sym1__), nil_index_list()), stan::math::lb_constrain(ind_beta[(sym1__ - 1)], 0), @@ -485,9 +506,9 @@ class model_constant_single final : public model_base_crtp S_hat; - S_hat = std::vector(N_obs, std::numeric_limits::quiet_NaN()); - - std::vector G_hat; - G_hat = std::vector(N_obs, std::numeric_limits::quiet_NaN()); - - current_statement__ = 20; - for (int i = 1; i <= N_obs; ++i) { - current_statement__ = 18; - if (logical_eq(id_factor[((i + 1) - 1)], id_factor[(i - 1)])) { - current_statement__ = 12; - if (logical_eq(census[(i - 1)], 1)) { - current_statement__ = 10; - assign(S_hat, cons_list(index_uni(i), nil_index_list()), - ind_S_0[(id_factor[(i - 1)] - 1)], "assigning variable S_hat"); - } - current_statement__ = 16; - if (logical_lt(i, N_obs)) { - current_statement__ = 13; - assign(G_hat, cons_list(index_uni(i), nil_index_list()), - ind_beta[(id_factor[(i - 1)] - 1)], "assigning variable G_hat"); - current_statement__ = 14; - assign(S_hat, cons_list(index_uni((i + 1)), nil_index_list()), - (S_hat[(i - 1)] + - (G_hat[(i - 1)] * census_interval[((i + 1) - 1)])), - "assigning variable S_hat"); - } - } else { + std::vector y_hat; + y_hat = std::vector(n_obs, std::numeric_limits::quiet_NaN()); + + std::vector Delta_hat; + Delta_hat = std::vector(n_obs, std::numeric_limits::quiet_NaN()); + + current_statement__ = 18; + for (int i = 1; i <= n_obs; ++i) { + current_statement__ = 10; + if (logical_eq(obs_index[(i - 1)], 1)) { current_statement__ = 8; - assign(G_hat, cons_list(index_uni(i), nil_index_list()), 0, - "assigning variable G_hat"); - }} - for (int sym1__ = 1; sym1__ <= N_obs; ++sym1__) { - vars__.emplace_back(S_hat[(sym1__ - 1)]);} - for (int sym1__ = 1; sym1__ <= N_obs; ++sym1__) { - vars__.emplace_back(G_hat[(sym1__ - 1)]);} + assign(y_hat, cons_list(index_uni(i), nil_index_list()), + ind_y_0[(ind_id[(i - 1)] - 1)], "assigning variable y_hat"); + } + current_statement__ = 11; + assign(Delta_hat, cons_list(index_uni(i), nil_index_list()), + growth(ind_beta[(ind_id[(i - 1)] - 1)], pstream__), + "assigning variable Delta_hat"); + current_statement__ = 16; + if (logical_lt(i, n_obs)) { + current_statement__ = 14; + if (logical_eq(ind_id[((i + 1) - 1)], ind_id[(i - 1)])) { + current_statement__ = 12; + assign(y_hat, cons_list(index_uni((i + 1)), nil_index_list()), + (y_hat[(i - 1)] + + (Delta_hat[(i - 1)] * (time[((i + 1) - 1)] - time[(i - 1)]))), + "assigning variable y_hat"); + } + } } + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + vars__.emplace_back(y_hat[(sym1__ - 1)]);} + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + vars__.emplace_back(Delta_hat[(sym1__ - 1)]);} } catch (const std::exception& e) { stan::lang::rethrow_located(e, locations_array__[current_statement__]); // Next line prevents compiler griping about no return @@ -555,33 +573,33 @@ class model_constant_single final : public model_base_crtp::min(); pos__ = 1; - std::vector ind_S_0; - ind_S_0 = std::vector(N_ind, std::numeric_limits::quiet_NaN()); + std::vector ind_y_0; + ind_y_0 = std::vector(n_ind, std::numeric_limits::quiet_NaN()); current_statement__ = 1; - assign(ind_S_0, nil_index_list(), context__.vals_r("ind_S_0"), - "assigning variable ind_S_0"); - std::vector ind_S_0_free__; - ind_S_0_free__ = std::vector(N_ind, std::numeric_limits::quiet_NaN()); + assign(ind_y_0, nil_index_list(), context__.vals_r("ind_y_0"), + "assigning variable ind_y_0"); + std::vector ind_y_0_free__; + ind_y_0_free__ = std::vector(n_ind, std::numeric_limits::quiet_NaN()); current_statement__ = 1; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 1; - assign(ind_S_0_free__, + assign(ind_y_0_free__, cons_list(index_uni(sym1__), nil_index_list()), - stan::math::lb_free(ind_S_0[(sym1__ - 1)], 0), - "assigning variable ind_S_0_free__");} + stan::math::lb_free(ind_y_0[(sym1__ - 1)], 0), + "assigning variable ind_y_0_free__");} std::vector ind_beta; - ind_beta = std::vector(N_ind, std::numeric_limits::quiet_NaN()); + ind_beta = std::vector(n_ind, std::numeric_limits::quiet_NaN()); current_statement__ = 2; assign(ind_beta, nil_index_list(), context__.vals_r("ind_beta"), "assigning variable ind_beta"); std::vector ind_beta_free__; - ind_beta_free__ = std::vector(N_ind, std::numeric_limits::quiet_NaN()); + ind_beta_free__ = std::vector(n_ind, std::numeric_limits::quiet_NaN()); current_statement__ = 2; - for (int sym1__ = 1; sym1__ <= N_ind; ++sym1__) { + for (int sym1__ = 1; sym1__ <= n_ind; ++sym1__) { current_statement__ = 2; assign(ind_beta_free__, cons_list(index_uni(sym1__), nil_index_list()), @@ -612,9 +630,9 @@ class model_constant_single final : public model_base_crtp& names__) const { names__.clear(); - names__.emplace_back("ind_S_0"); + names__.emplace_back("ind_y_0"); names__.emplace_back("ind_beta"); names__.emplace_back("species_beta_mu"); names__.emplace_back("species_beta_sigma"); names__.emplace_back("global_error_sigma"); - names__.emplace_back("S_hat"); - names__.emplace_back("G_hat"); + names__.emplace_back("y_hat"); + names__.emplace_back("Delta_hat"); } // get_param_names() inline void get_dims(std::vector>& dimss__) const { dimss__.clear(); - dimss__.emplace_back(std::vector{static_cast(N_ind)}); + dimss__.emplace_back(std::vector{static_cast(n_ind)}); - dimss__.emplace_back(std::vector{static_cast(N_ind)}); + dimss__.emplace_back(std::vector{static_cast(n_ind)}); dimss__.emplace_back(std::vector{}); @@ -650,9 +668,9 @@ class model_constant_single final : public model_base_crtp{}); - dimss__.emplace_back(std::vector{static_cast(N_obs)}); + dimss__.emplace_back(std::vector{static_cast(n_obs)}); - dimss__.emplace_back(std::vector{static_cast(N_obs)}); + dimss__.emplace_back(std::vector{static_cast(n_obs)}); } // get_dims() @@ -662,11 +680,11 @@ class model_constant_single final : public model_base_crtp +using namespace Rcpp ; +#include "stanExports_constant_single_ind.h" + +RCPP_MODULE(stan_fit4constant_single_ind_mod) { + + + class_ >("rstantools_model_constant_single_ind") + + .constructor() + + + .method("call_sampler", &rstan::stan_fit ::call_sampler) + .method("param_names", &rstan::stan_fit ::param_names) + .method("param_names_oi", &rstan::stan_fit ::param_names_oi) + .method("param_fnames_oi", &rstan::stan_fit ::param_fnames_oi) + .method("param_dims", &rstan::stan_fit ::param_dims) + .method("param_dims_oi", &rstan::stan_fit ::param_dims_oi) + .method("update_param_oi", &rstan::stan_fit ::update_param_oi) + .method("param_oi_tidx", &rstan::stan_fit ::param_oi_tidx) + .method("grad_log_prob", &rstan::stan_fit ::grad_log_prob) + .method("log_prob", &rstan::stan_fit ::log_prob) + .method("unconstrain_pars", &rstan::stan_fit ::unconstrain_pars) + .method("constrain_pars", &rstan::stan_fit ::constrain_pars) + .method("num_pars_unconstrained", &rstan::stan_fit ::num_pars_unconstrained) + .method("unconstrained_param_names", &rstan::stan_fit ::unconstrained_param_names) + .method("constrained_param_names", &rstan::stan_fit ::constrained_param_names) + .method("standalone_gqs", &rstan::stan_fit ::standalone_gqs) + ; +} diff --git a/src/stanExports_constant_single_ind.h b/src/stanExports_constant_single_ind.h new file mode 100644 index 0000000..a88757d --- /dev/null +++ b/src/stanExports_constant_single_ind.h @@ -0,0 +1,657 @@ +// Generated by rstantools. Do not edit by hand. + +/* + rmot is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + rmot is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with rmot. If not, see . +*/ +#ifndef MODELS_HPP +#define MODELS_HPP +#define STAN__SERVICES__COMMAND_HPP +#ifndef USE_STANC3 +#define USE_STANC3 +#endif +#include +// Code generated by stanc v2.26.1-4-gd72b68b7-dirty +#include +namespace model_constant_single_ind_namespace { +inline void validate_positive_index(const char* var_name, const char* expr, + int val) { + if (val < 1) { + std::stringstream msg; + msg << "Found dimension size less than one in simplex declaration" + << "; variable=" << var_name << "; dimension size expression=" << expr + << "; expression value=" << val; + std::string msg_str(msg.str()); + throw std::invalid_argument(msg_str.c_str()); + } +} +inline void validate_unit_vector_index(const char* var_name, const char* expr, + int val) { + if (val <= 1) { + std::stringstream msg; + if (val == 1) { + msg << "Found dimension size one in unit vector declaration." + << " One-dimensional unit vector is discrete" + << " but the target distribution must be continuous." + << " variable=" << var_name << "; dimension size expression=" << expr; + } else { + msg << "Found dimension size less than one in unit vector declaration" + << "; variable=" << var_name << "; dimension size expression=" << expr + << "; expression value=" << val; + } + std::string msg_str(msg.str()); + throw std::invalid_argument(msg_str.c_str()); + } +} +using std::istream; +using std::string; +using std::stringstream; +using std::vector; +using std::pow; +using stan::io::dump; +using stan::math::lgamma; +using stan::model::model_base_crtp; +using stan::model::rvalue; +using stan::model::cons_list; +using stan::model::index_uni; +using stan::model::index_max; +using stan::model::index_min; +using stan::model::index_min_max; +using stan::model::index_multi; +using stan::model::index_omni; +using stan::model::nil_index_list; +using namespace stan::math; +using stan::math::pow; +stan::math::profile_map profiles__; +static int current_statement__= 0; +static const std::vector locations_array__ = {" (found before start of program)", + " (in 'constant_single_ind', line 19, column 2 to column 24)", + " (in 'constant_single_ind', line 20, column 2 to column 25)", + " (in 'constant_single_ind', line 22, column 2 to column 35)", + " (in 'constant_single_ind', line 48, column 2 to column 20)", + " (in 'constant_single_ind', line 49, column 2 to column 24)", + " (in 'constant_single_ind', line 52, column 6 to column 25)", + " (in 'constant_single_ind', line 51, column 23 to line 53, column 5)", + " (in 'constant_single_ind', line 51, column 4 to line 53, column 5)", + " (in 'constant_single_ind', line 54, column 4 to column 46)", + " (in 'constant_single_ind', line 56, column 6 to column 63)", + " (in 'constant_single_ind', line 55, column 17 to line 57, column 5)", + " (in 'constant_single_ind', line 55, column 4 to line 57, column 5)", + " (in 'constant_single_ind', line 50, column 19 to line 58, column 3)", + " (in 'constant_single_ind', line 50, column 2 to line 58, column 3)", + " (in 'constant_single_ind', line 26, column 13 to column 18)", + " (in 'constant_single_ind', line 26, column 2 to column 20)", + " (in 'constant_single_ind', line 27, column 17 to column 22)", + " (in 'constant_single_ind', line 27, column 2 to column 24)", + " (in 'constant_single_ind', line 30, column 6 to column 25)", + " (in 'constant_single_ind', line 29, column 23 to line 31, column 5)", + " (in 'constant_single_ind', line 29, column 4 to line 31, column 5)", + " (in 'constant_single_ind', line 32, column 4 to column 46)", + " (in 'constant_single_ind', line 34, column 6 to column 63)", + " (in 'constant_single_ind', line 33, column 17 to line 35, column 5)", + " (in 'constant_single_ind', line 33, column 4 to line 35, column 5)", + " (in 'constant_single_ind', line 28, column 19 to line 36, column 3)", + " (in 'constant_single_ind', line 28, column 2 to line 36, column 3)", + " (in 'constant_single_ind', line 38, column 2 to column 44)", + " (in 'constant_single_ind', line 41, column 2 to column 48)", + " (in 'constant_single_ind', line 42, column 2 to column 31)", + " (in 'constant_single_ind', line 44, column 2 to column 37)", + " (in 'constant_single_ind', line 10, column 2 to column 12)", + " (in 'constant_single_ind', line 11, column 13 to column 18)", + " (in 'constant_single_ind', line 11, column 2 to column 20)", + " (in 'constant_single_ind', line 12, column 16 to column 21)", + " (in 'constant_single_ind', line 12, column 2 to column 23)", + " (in 'constant_single_ind', line 13, column 12 to column 17)", + " (in 'constant_single_ind', line 13, column 2 to column 19)", + " (in 'constant_single_ind', line 14, column 2 to column 15)", + " (in 'constant_single_ind', line 48, column 13 to column 18)", + " (in 'constant_single_ind', line 49, column 17 to column 22)", + " (in 'constant_single_ind', line 5, column 4 to column 16)", + " (in 'constant_single_ind', line 4, column 32 to line 6, column 3)"}; +template +stan::promote_args_t +growth(const T0__& y, const T1__& beta, std::ostream* pstream__) { + using local_scalar_t__ = stan::promote_args_t; + const static bool propto__ = true; + (void) propto__; + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + current_statement__ = 42; + return beta; + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); + } + +} +struct growth_functor__ { +template +stan::promote_args_t +operator()(const T0__& y, const T1__& beta, std::ostream* pstream__) const +{ +return growth(y, beta, pstream__); +} +}; +#include +class model_constant_single_ind final : public model_base_crtp { +private: + int n_obs; + std::vector y_obs; + std::vector obs_index; + std::vector time; + double y_0_obs; + +public: + ~model_constant_single_ind() { } + + inline std::string model_name() const final { return "model_constant_single_ind"; } + inline std::vector model_compile_info() const noexcept { + return std::vector{"stanc_version = stanc3 v2.26.1-4-gd72b68b7-dirty", "stancflags = "}; + } + + + model_constant_single_ind(stan::io::var_context& context__, + unsigned int random_seed__ = 0, + std::ostream* pstream__ = nullptr) : model_base_crtp(0) { + using local_scalar_t__ = double ; + boost::ecuyer1988 base_rng__ = + stan::services::util::create_rng(random_seed__, 0); + (void) base_rng__; // suppress unused var warning + static const char* function__ = "model_constant_single_ind_namespace::model_constant_single_ind"; + (void) function__; // suppress unused var warning + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + int pos__; + pos__ = std::numeric_limits::min(); + + pos__ = 1; + current_statement__ = 32; + context__.validate_dims("data initialization","n_obs","int", + context__.to_vec()); + n_obs = std::numeric_limits::min(); + + current_statement__ = 32; + n_obs = context__.vals_i("n_obs")[(1 - 1)]; + current_statement__ = 33; + validate_non_negative_index("y_obs", "n_obs", n_obs); + current_statement__ = 34; + context__.validate_dims("data initialization","y_obs","double", + context__.to_vec(n_obs)); + y_obs = std::vector(n_obs, std::numeric_limits::quiet_NaN()); + + current_statement__ = 34; + assign(y_obs, nil_index_list(), context__.vals_r("y_obs"), + "assigning variable y_obs"); + current_statement__ = 35; + validate_non_negative_index("obs_index", "n_obs", n_obs); + current_statement__ = 36; + context__.validate_dims("data initialization","obs_index","int", + context__.to_vec(n_obs)); + obs_index = std::vector(n_obs, std::numeric_limits::min()); + + current_statement__ = 36; + assign(obs_index, nil_index_list(), context__.vals_i("obs_index"), + "assigning variable obs_index"); + current_statement__ = 37; + validate_non_negative_index("time", "n_obs", n_obs); + current_statement__ = 38; + context__.validate_dims("data initialization","time","double", + context__.to_vec(n_obs)); + time = std::vector(n_obs, std::numeric_limits::quiet_NaN()); + + current_statement__ = 38; + assign(time, nil_index_list(), context__.vals_r("time"), + "assigning variable time"); + current_statement__ = 39; + context__.validate_dims("data initialization","y_0_obs","double", + context__.to_vec()); + y_0_obs = std::numeric_limits::quiet_NaN(); + + current_statement__ = 39; + y_0_obs = context__.vals_r("y_0_obs")[(1 - 1)]; + current_statement__ = 40; + validate_non_negative_index("y_hat", "n_obs", n_obs); + current_statement__ = 41; + validate_non_negative_index("Delta_hat", "n_obs", n_obs); + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); + } + num_params_r__ = 0U; + + try { + num_params_r__ += 1; + num_params_r__ += 1; + num_params_r__ += 1; + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); + } + } + template * = nullptr, stan::require_vector_like_vt* = nullptr> + inline stan::scalar_type_t log_prob_impl(VecR& params_r__, + VecI& params_i__, + std::ostream* pstream__ = nullptr) const { + using T__ = stan::scalar_type_t; + using local_scalar_t__ = T__; + T__ lp__(0.0); + stan::math::accumulator lp_accum__; + static const char* function__ = "model_constant_single_ind_namespace::log_prob"; +(void) function__; // suppress unused var warning + stan::io::reader in__(params_r__, params_i__); + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + local_scalar_t__ ind_y_0; + ind_y_0 = DUMMY_VAR__; + + current_statement__ = 1; + ind_y_0 = in__.scalar(); + current_statement__ = 1; + if (jacobian__) { + current_statement__ = 1; + ind_y_0 = stan::math::lb_constrain(ind_y_0, 0, lp__); + } else { + current_statement__ = 1; + ind_y_0 = stan::math::lb_constrain(ind_y_0, 0); + } + local_scalar_t__ ind_beta; + ind_beta = DUMMY_VAR__; + + current_statement__ = 2; + ind_beta = in__.scalar(); + current_statement__ = 2; + if (jacobian__) { + current_statement__ = 2; + ind_beta = stan::math::lb_constrain(ind_beta, 0, lp__); + } else { + current_statement__ = 2; + ind_beta = stan::math::lb_constrain(ind_beta, 0); + } + local_scalar_t__ global_error_sigma; + global_error_sigma = DUMMY_VAR__; + + current_statement__ = 3; + global_error_sigma = in__.scalar(); + current_statement__ = 3; + if (jacobian__) { + current_statement__ = 3; + global_error_sigma = stan::math::lb_constrain(global_error_sigma, 0, + lp__); + } else { + current_statement__ = 3; + global_error_sigma = stan::math::lb_constrain(global_error_sigma, 0); + } + { + current_statement__ = 15; + validate_non_negative_index("y_hat", "n_obs", n_obs); + std::vector y_hat; + y_hat = std::vector(n_obs, DUMMY_VAR__); + + current_statement__ = 17; + validate_non_negative_index("Delta_hat", "n_obs", n_obs); + std::vector Delta_hat; + Delta_hat = std::vector(n_obs, DUMMY_VAR__); + + current_statement__ = 27; + for (int i = 1; i <= n_obs; ++i) { + current_statement__ = 21; + if (logical_eq(obs_index[(i - 1)], 1)) { + current_statement__ = 19; + assign(y_hat, cons_list(index_uni(i), nil_index_list()), ind_y_0, + "assigning variable y_hat"); + } + current_statement__ = 22; + assign(Delta_hat, cons_list(index_uni(i), nil_index_list()), + growth(y_hat[(i - 1)], ind_beta, pstream__), + "assigning variable Delta_hat"); + current_statement__ = 25; + if (logical_lt(i, n_obs)) { + current_statement__ = 23; + assign(y_hat, cons_list(index_uni((i + 1)), nil_index_list()), + (y_hat[(i - 1)] + + (Delta_hat[(i - 1)] * (time[((i + 1) - 1)] - time[(i - 1)]))), + "assigning variable y_hat"); + } } + current_statement__ = 28; + lp_accum__.add( + normal_lpdf(y_obs, y_hat, global_error_sigma)); + current_statement__ = 29; + lp_accum__.add( + normal_lpdf(ind_y_0, y_0_obs, global_error_sigma)); + current_statement__ = 30; + lp_accum__.add(lognormal_lpdf(ind_beta, 0.1, 1)); + current_statement__ = 31; + lp_accum__.add(cauchy_lpdf(global_error_sigma, 0.1, 1)); + } + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); + } + lp_accum__.add(lp__); + return lp_accum__.sum(); + } // log_prob_impl() + + template * = nullptr, stan::require_vector_like_vt* = nullptr, stan::require_std_vector_vt* = nullptr> + inline void write_array_impl(RNG& base_rng__, VecR& params_r__, + VecI& params_i__, VecVar& vars__, + const bool emit_transformed_parameters__ = true, + const bool emit_generated_quantities__ = true, + std::ostream* pstream__ = nullptr) const { + using local_scalar_t__ = double; + vars__.resize(0); + stan::io::reader in__(params_r__, params_i__); + static const char* function__ = "model_constant_single_ind_namespace::write_array"; +(void) function__; // suppress unused var warning + (void) function__; // suppress unused var warning + double lp__ = 0.0; + (void) lp__; // dummy to suppress unused var warning + stan::math::accumulator lp_accum__; + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + double ind_y_0; + ind_y_0 = std::numeric_limits::quiet_NaN(); + + current_statement__ = 1; + ind_y_0 = in__.scalar(); + current_statement__ = 1; + ind_y_0 = stan::math::lb_constrain(ind_y_0, 0); + double ind_beta; + ind_beta = std::numeric_limits::quiet_NaN(); + + current_statement__ = 2; + ind_beta = in__.scalar(); + current_statement__ = 2; + ind_beta = stan::math::lb_constrain(ind_beta, 0); + double global_error_sigma; + global_error_sigma = std::numeric_limits::quiet_NaN(); + + current_statement__ = 3; + global_error_sigma = in__.scalar(); + current_statement__ = 3; + global_error_sigma = stan::math::lb_constrain(global_error_sigma, 0); + vars__.emplace_back(ind_y_0); + vars__.emplace_back(ind_beta); + vars__.emplace_back(global_error_sigma); + if (logical_negation((primitive_value(emit_transformed_parameters__) || + primitive_value(emit_generated_quantities__)))) { + return ; + } + if (logical_negation(emit_generated_quantities__)) { + return ; + } + std::vector y_hat; + y_hat = std::vector(n_obs, std::numeric_limits::quiet_NaN()); + + std::vector Delta_hat; + Delta_hat = std::vector(n_obs, std::numeric_limits::quiet_NaN()); + + current_statement__ = 14; + for (int i = 1; i <= n_obs; ++i) { + current_statement__ = 8; + if (logical_eq(obs_index[(i - 1)], 1)) { + current_statement__ = 6; + assign(y_hat, cons_list(index_uni(i), nil_index_list()), ind_y_0, + "assigning variable y_hat"); + } + current_statement__ = 9; + assign(Delta_hat, cons_list(index_uni(i), nil_index_list()), + growth(y_hat[(i - 1)], ind_beta, pstream__), + "assigning variable Delta_hat"); + current_statement__ = 12; + if (logical_lt(i, n_obs)) { + current_statement__ = 10; + assign(y_hat, cons_list(index_uni((i + 1)), nil_index_list()), + (y_hat[(i - 1)] + + (Delta_hat[(i - 1)] * (time[((i + 1) - 1)] - time[(i - 1)]))), + "assigning variable y_hat"); + } } + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + vars__.emplace_back(y_hat[(sym1__ - 1)]);} + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + vars__.emplace_back(Delta_hat[(sym1__ - 1)]);} + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); + } + } // write_array_impl() + + template * = nullptr, stan::require_vector_like_vt* = nullptr> + inline void transform_inits_impl(const stan::io::var_context& context__, + VecI& params_i__, VecVar& vars__, + std::ostream* pstream__ = nullptr) const { + using local_scalar_t__ = double; + vars__.clear(); + vars__.reserve(num_params_r__); + + try { + int pos__; + pos__ = std::numeric_limits::min(); + + pos__ = 1; + double ind_y_0; + ind_y_0 = std::numeric_limits::quiet_NaN(); + + current_statement__ = 1; + ind_y_0 = context__.vals_r("ind_y_0")[(1 - 1)]; + double ind_y_0_free__; + ind_y_0_free__ = std::numeric_limits::quiet_NaN(); + + current_statement__ = 1; + ind_y_0_free__ = stan::math::lb_free(ind_y_0, 0); + double ind_beta; + ind_beta = std::numeric_limits::quiet_NaN(); + + current_statement__ = 2; + ind_beta = context__.vals_r("ind_beta")[(1 - 1)]; + double ind_beta_free__; + ind_beta_free__ = std::numeric_limits::quiet_NaN(); + + current_statement__ = 2; + ind_beta_free__ = stan::math::lb_free(ind_beta, 0); + double global_error_sigma; + global_error_sigma = std::numeric_limits::quiet_NaN(); + + current_statement__ = 3; + global_error_sigma = context__.vals_r("global_error_sigma")[(1 - 1)]; + double global_error_sigma_free__; + global_error_sigma_free__ = std::numeric_limits::quiet_NaN(); + + current_statement__ = 3; + global_error_sigma_free__ = stan::math::lb_free(global_error_sigma, 0); + vars__.emplace_back(ind_y_0_free__); + vars__.emplace_back(ind_beta_free__); + vars__.emplace_back(global_error_sigma_free__); + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); + } + } // transform_inits_impl() + + inline void get_param_names(std::vector& names__) const { + + names__.clear(); + names__.emplace_back("ind_y_0"); + names__.emplace_back("ind_beta"); + names__.emplace_back("global_error_sigma"); + names__.emplace_back("y_hat"); + names__.emplace_back("Delta_hat"); + } // get_param_names() + + inline void get_dims(std::vector>& dimss__) const { + dimss__.clear(); + dimss__.emplace_back(std::vector{}); + + dimss__.emplace_back(std::vector{}); + + dimss__.emplace_back(std::vector{}); + + dimss__.emplace_back(std::vector{static_cast(n_obs)}); + + dimss__.emplace_back(std::vector{static_cast(n_obs)}); + + } // get_dims() + + inline void constrained_param_names( + std::vector& param_names__, + bool emit_transformed_parameters__ = true, + bool emit_generated_quantities__ = true) const + final { + + param_names__.emplace_back(std::string() + "ind_y_0"); + param_names__.emplace_back(std::string() + "ind_beta"); + param_names__.emplace_back(std::string() + "global_error_sigma"); + if (emit_transformed_parameters__) { + + } + + if (emit_generated_quantities__) { + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + { + param_names__.emplace_back(std::string() + "y_hat" + '.' + std::to_string(sym1__)); + }} + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + { + param_names__.emplace_back(std::string() + "Delta_hat" + '.' + std::to_string(sym1__)); + }} + } + + } // constrained_param_names() + + inline void unconstrained_param_names( + std::vector& param_names__, + bool emit_transformed_parameters__ = true, + bool emit_generated_quantities__ = true) const + final { + + param_names__.emplace_back(std::string() + "ind_y_0"); + param_names__.emplace_back(std::string() + "ind_beta"); + param_names__.emplace_back(std::string() + "global_error_sigma"); + if (emit_transformed_parameters__) { + + } + + if (emit_generated_quantities__) { + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + { + param_names__.emplace_back(std::string() + "y_hat" + '.' + std::to_string(sym1__)); + }} + for (int sym1__ = 1; sym1__ <= n_obs; ++sym1__) { + { + param_names__.emplace_back(std::string() + "Delta_hat" + '.' + std::to_string(sym1__)); + }} + } + + } // unconstrained_param_names() + + inline std::string get_constrained_sizedtypes() const { + stringstream s__; + s__ << "[{\"name\":\"ind_y_0\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"ind_beta\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"global_error_sigma\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"y_hat\",\"type\":{\"name\":\"array\",\"length\":" << n_obs << ",\"element_type\":{\"name\":\"real\"}},\"block\":\"generated_quantities\"},{\"name\":\"Delta_hat\",\"type\":{\"name\":\"array\",\"length\":" << n_obs << ",\"element_type\":{\"name\":\"real\"}},\"block\":\"generated_quantities\"}]"; + return s__.str(); + } // get_constrained_sizedtypes() + + inline std::string get_unconstrained_sizedtypes() const { + stringstream s__; + s__ << "[{\"name\":\"ind_y_0\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"ind_beta\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"global_error_sigma\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"y_hat\",\"type\":{\"name\":\"array\",\"length\":" << n_obs << ",\"element_type\":{\"name\":\"real\"}},\"block\":\"generated_quantities\"},{\"name\":\"Delta_hat\",\"type\":{\"name\":\"array\",\"length\":" << n_obs << ",\"element_type\":{\"name\":\"real\"}},\"block\":\"generated_quantities\"}]"; + return s__.str(); + } // get_unconstrained_sizedtypes() + + + // Begin method overload boilerplate + template + inline void write_array(RNG& base_rng, + Eigen::Matrix& params_r, + Eigen::Matrix& vars, + const bool emit_transformed_parameters = true, + const bool emit_generated_quantities = true, + std::ostream* pstream = nullptr) const { + std::vector vars_vec(vars.size()); + std::vector params_i; + write_array_impl(base_rng, params_r, params_i, vars_vec, + emit_transformed_parameters, emit_generated_quantities, pstream); + vars.resize(vars_vec.size()); + for (int i = 0; i < vars.size(); ++i) { + vars.coeffRef(i) = vars_vec[i]; + } + } + template + inline void write_array(RNG& base_rng, std::vector& params_r, + std::vector& params_i, + std::vector& vars, + bool emit_transformed_parameters = true, + bool emit_generated_quantities = true, + std::ostream* pstream = nullptr) const { + write_array_impl(base_rng, params_r, params_i, vars, emit_transformed_parameters, emit_generated_quantities, pstream); + } + template + inline T_ log_prob(Eigen::Matrix& params_r, + std::ostream* pstream = nullptr) const { + Eigen::Matrix params_i; + return log_prob_impl(params_r, params_i, pstream); + } + template + inline T__ log_prob(std::vector& params_r, + std::vector& params_i, + std::ostream* pstream = nullptr) const { + return log_prob_impl(params_r, params_i, pstream); + } + + inline void transform_inits(const stan::io::var_context& context, + Eigen::Matrix& params_r, + std::ostream* pstream = nullptr) const final { + std::vector params_r_vec(params_r.size()); + std::vector params_i; + transform_inits_impl(context, params_i, params_r_vec, pstream); + params_r.resize(params_r_vec.size()); + for (int i = 0; i < params_r.size(); ++i) { + params_r.coeffRef(i) = params_r_vec[i]; + } + } + inline void transform_inits(const stan::io::var_context& context, + std::vector& params_i, + std::vector& vars, + std::ostream* pstream = nullptr) const final { + transform_inits_impl(context, params_i, vars, pstream); + } +}; +} +using stan_model = model_constant_single_ind_namespace::model_constant_single_ind; +#ifndef USING_R +// Boilerplate +stan::model::model_base& new_model( + stan::io::var_context& data_context, + unsigned int seed, + std::ostream* msg_stream) { + stan_model* m = new stan_model(data_context, seed, msg_stream); + return *m; +} +stan::math::profile_map& get_stan_profile_data() { + return model_constant_single_ind_namespace::profiles__; +} +#endif +#endif diff --git a/tests/testthat/fixtures/constant/constant_data_multi_ind.rds b/tests/testthat/fixtures/constant/constant_data_multi_ind.rds new file mode 100644 index 0000000000000000000000000000000000000000..082c76072c2e3878aea31a11a49c1fca0ef656e2 GIT binary patch literal 814 zcmV+}1JV2+iwFP!000001B>8dU|?WoU}0uvU}gm}8CXL@+;lA%7?^~C9FQa*5Hs3; zV1NN$5Fdy|VR8-}a5)fP)WJFFoLOFvmV+1T`l3I6-y8x$>{vu}mpBA3-C_B5d96d( zx0jc$RC_r@d8GGwR{n5^ow;l28IK?bbD@>D9K9R{f0)1l!N9-*3Irf# z17dax2%q}~1-^V#3lDTjedk^Rk@mzbMc3}gtP1V(B}BD%ul zoWx>ugW zb>LD9FxMN;-(CM0nCqua5}TF_%=JH83fFuG=K7ejGlv_2xqg49`Ymp_CoEq@nBToT z5t!?ri&kIW4b1hWpCvCp1?KvujObUOT<>{Z>ybsTjcl8VTc6 sBFS(9nIO5;3Os=Ubv`hcA?akxL-q`RZfQ8dU|?WoU}0uvU}gm}8CXL@+;lA%7?^~C91bAn17b$| z4-7B><~zus@Yx-llg^pt^=LVGv92%rSr{Azh>!%U1BT?nr~{W;fFWr-e|P<3U`S4zBsMJ<7?M9* z3fFuGhGb0HnZu31kldfCev6w5e&dDeXT*(JkmG_iFEKZ@7|0L+djKNHRvBMXRGJEl zu9D1L7@I9KFC{)R1|%l9%)IoR)OZvJ N0RZ60Tu!?K006}a%VYom literal 0 HcmV?d00001 diff --git a/tests/testthat/fixtures/constant/make_constant.R b/tests/testthat/fixtures/constant/make_constant.R new file mode 100644 index 0000000..aab703d --- /dev/null +++ b/tests/testthat/fixtures/constant/make_constant.R @@ -0,0 +1,39 @@ +#Generates data for constant growth model testing +#Differential equation governing dynamics. +DE <- function(y, pars){ + return(pars[1]) +} + +#Function to generate distribution of DE parameters +DE_par_generation <- function(n_ind, pars=list(mean=0, sd=1)){ + par_sample <- data.frame(beta=exp(rnorm(n_ind, mean = pars[[1]], sd=pars[[2]]))) + return(par_sample) +} + +#Set required values +model_name <- "constant" +n_ind <- 3 #Number of individuals for multi-individual data. Single individual takes the first. +n_obs_per_ind <- 7 #How many observations per individual. +interval <- 5 #Time interval between observations +time = seq(from = 0, by = interval, length.out = n_obs_per_ind) + +#Produce parameters +DE_pars <- DE_par_generation(n_ind) +initial_conditions <- data.frame(y_0=exp(rnorm(n_ind, mean = 2, sd=1))) + +#Generate true values +true_data <- rmot_build_true_test_data(n_ind, n_obs_per_ind, interval, + DE_pars, initial_conditions, DE) + +#Generate observations +y_obs <- rmot_add_normal_error(true_data$y_true) + +#Export datasets +rmot_export_test_data(n_obs_per_ind, + n_ind, + y_obs, + time, + DE_pars, + initial_conditions, + true_data, + model_name) diff --git a/tests/testthat/fixtures/linear/lm_baseline_output.rds b/tests/testthat/fixtures/linear/lm_baseline_output.rds new file mode 100644 index 0000000000000000000000000000000000000000..77cfbc4eeee835f4f7ee0ef62b63a9828ba66191 GIT binary patch literal 3187 zcmV-(42<(1iwFP!0000016@}MG?v@f7w35gG*Hs0Nx3D8BKaD&l+@KgNQ$1P& zA}LdmAtgg9v(O+)DU|A*8%>H!0~Mi{L`9}d8NT!Vzkl7e?z7fe`#JkOXa9b`z4v+6 zIu_Cr5)x7pNOF{f6hmJn7&B7?H67|3PbY|uF?87BI(x8jrR*pyJd{U2Oke~F)c(QJVKXSUW-TeG{mJI+J)cfdbj4)U&noO zKMeovXpZ~RG<~ERcHyq`H*U@m{tw~w>1*0Z=Mr9hOz@|OTZGByj62{|hI@)0Mo7;X z#4lgDr!QLOjfXrHFKbEI;9>qr@4B}oxc^n*(`JGV8Fh=-D4yIz7@uwIr1R|r zHGc2oqAc8*)DvtPF^In`8uZne=Ye}U=O^XgNcn9WE0v0r&7$yEoqSfap%8Z!UkkFC z+lo7`o`~35&d0Cuwwl{tFXAuV7Q)lXvvB{W@K3YdIJjTo&_1KM96VAKJo}tJ zhl~n4l`%nS78zv&7N(^txOeaSsrJvF5sA&*uiu8Hepx(llUe0+w+0WmXZ)3AwF3`K zO%A{5IDq@iBWvHh5#q*>paTwh6L6F2^g7wVBHaIUuf~K|2}I)cyQbvwCwREqBBjr2 z74G`()M!ZB@eip*`lE-==bVc&2>AoKLt*!%aaL?udW&_c zPu?1yX^mS(?uhyV-s4u&&vg}bD!6@xe}1vC0q&@Jy+6hI9`0=HYa<2Qac|c`wNqnW z;KWZUwtgf8f4CUFyjEa{D;m}xem3wKzsOLIdgN1&zYS!06kRF91D~fhx?XvShgL_= zDHsOa`JiTOxw`eEPI!+?yhlCknWAocEuNn4dFNaYvE?`n+( zBql~rEXM*;l}TFHJqu91$bZT?HDKhFrS+=P0qHKdmE>v+2>-1xCc_udEf4gxSOlot zIC^-!4lst-YGY=s1ZLT#uMt;;z{I|M_+rySKy`!W<0BIR?e*2rnEnT#$NN9tKM@NI zFOL6yJ%{4#+pCkN52#Y_Zbq3qAV0a$dX>unb=}Wi;_(sCGnb^x`ksJFb!)hyHb9;c z*HvmLe$~n+#$^G(jANdvdz%1bPMcy?0WglTiBf8jfG#N5ciVeWT;67}8O13He-(ez z?H?;tvC?>M4an$WcvhY+)xW_!St%5dLsT*eJ_8I#9MG4M4`^tzv#fJHpn|eB7njQd zIxlKmapW;DCeL+zc!j|1*{m4nU=7Sp=d86SIlveT+l}+e0Ua20))P4cicpbJ%Tb`Z zI(IG;O_!eOwb+AzM$Ci@{{k8(FG%r?8lb~+r!`Yg1G-YF9L<|X{hQCXpVFY; z1={vk&^qfywq`!F0kkJr8m*-|)|Zya-bezJoz1A1G|_%Yn3u7S=5^nth`_iGKnDc| z{O7%Ze96r*x&lCn!TL!a#ei;9=M}$_26X&;hA5Wyb4TplY{fx9%c^VU^Hl&%v(60` zlmXJ3|Iy1>3(y}=1$)XVFGr<``hmlMUhC(d-%sO&<-z}t41*CHtg($@}{o`_TfX-(2tUlj+8%};A@ zHx=NN?)rMlBMW$k%4=RangYjG+qJoLB5>7yDwb4j1=*K&*T?c0;6m?^yp;oRmmc{0 zdf8pzp@;M3zNZ1FCyQTHu^6~7w^&Zns|U{ao=ktKSU@>>9SvzB;5=*5cF(u~+=r*F zZc9G`Zt#TtDRUBmqZ41FZ+8cnthZlga4rFF)s|Y^kpsNgV^_vsR|WnC-i44#SCC$D zyk^|U1K=EOzx~5-A;=ufUi!p$BcL*^?&Ae*z*oNC-YuI_D_|X7J&z@^@ zwR8dSxz$l`)~N$`_ehSDwJtFCdUVY1!~$>V@Y;#tF2K2Lacs8NsNa?!WqO&tqXzsl zwuve`?*e*X>mMWV2c%Ybwm?Y#XXn8aU-kzL^X~)4F}0E$x7STjQ9tr_KPl71g%yUfTgT zU#F)ay9IbEDFJtRuYnij>!Z|;f&1Vd<9A~Pa5pXg6q~yhxXlNH39JInYO_Ay68gWp zmfpW303 zk4-Nz@D&2jD1GZEG#_Mq>T;`J%YbxxXPHlA9L;l({F35i;9K2Ez3N5l|D?$)&B73P zNhiZRGb@3ou8-U=jr=+%g!S9HtFZALVMQCynw$|5wp^Au;!{He^FFUAR*WL-ZFiZ) zxR0=j;waA3GQwU+@O-@VJz-s)lbpP}2pf>xSh;ZxVObqN#Rv()3d~+io9j$ivc)4= zTu<0%hvL@+$q>O5^BDW|$%HL33w)%sjIdUvVI2pi$L*k_{=VZEFcH$2@$*w{^P zv^VDyc2~%K%YX%hwOpd1q?JTi>oq$X$5VU<&djPSUkK}+`C^M*5n)}Lym5wm*ff}@0uY7O0(`h>7K>ml?k-Cr^^WYpP2SpUGl6`o3j6$T%YcA)#rf1aGE z5dV*LwiGL9Q`|$VpF}I)qVa;ptEO!rtj|{|V{HSFTOJhxHX zqDNRyZ_lY0W)SvZ!Q$pvcN(|arm~pkvrtR>rezF`=M?o0a|pZF^5nGH{e-o3>5)9q zNZ2D@iLo8UgiToBoqK|=hc%_j6e&55xl(fQrNE%R3r)?Iw~=ASffYR28s>Xo!FRA*LCilBObSX|-UB&?nB8IeBCm%uSt z%ZB!m$F`i_bjmNJKBsp>FXb^ag`G=vADk?&^7AZVx4#MeJcaVw&0k-8R*vTVSZjsJ zYFd{I6Oa5^O8IZD@bGY?dYtpb2}@|)1nqcjinHfyo|U~8VQpNJkDWS3_1sZxYt|%e z?WDR~u4}S<%J;cnQ@ShVdE$x06{^!}?}_RF9>tG1 z@@@KkdPR|>=VTPUsqp{r@*))MwHMKsasRoz{J&Gq)=}Uj)YVXT61v;F+1l@S|L=nA zuLT=>cfo&sqks8CP7Zqo|NhW*|L)!Yg5R<-#qF0< Z+D_muP~YQ5_tV+%{{qcOU;&N|007YkSE2v_ literal 0 HcmV?d00001 diff --git a/tests/testthat/fixtures/linear/lm_data.rds b/tests/testthat/fixtures/linear/lm_data.rds new file mode 100644 index 0000000000000000000000000000000000000000..086c352d42789d6239719f22a3ded3b716ca8af2 GIT binary patch literal 547 zcmV+;0^I!{iwFP!000001B>8dU|?WoU}0ipU}gm}8CXL@+;lB~V!~hv1_nML4sqaM z00RdRD6Ime^`Nu~l(vM@qv}W1(@{MGES(B**%mCm&gGzbe4R+uLk|a=9f7Y~<#sr@ z+ygP>9lRoE0>LZ?ztThfF&f_-gp9$!LG?OVql2v=NO_QhYZk;DuMa@dw$Qf$@FCg*(Y-u3i!0!%mhtgMwI?ELhdFKuY?U@X+c!`5A zlrL-xF;^`OqTdFp-sMWwL!eO(Udtio`L{vMhx%9P3<#Wcu>1?r?|cj*@3{-4t=hr& zBEPn8|dH%HBT@Q;x5%G5dYcyfSBWI0inI7Lc-1OJj6Uf zpde7jK?NF~*6&0>8Xa6DA@S+4hG5pZ=CFdj-!}P%fxl$5K67|xG lfJUMz=Pb%E*MnIFidUuvAo2hI|207B9{}w{FELvO008l}_ap!S literal 0 HcmV?d00001 diff --git a/tests/testthat/fixtures/linear/make_lm.R b/tests/testthat/fixtures/linear/make_lm.R new file mode 100644 index 0000000..f57a657 --- /dev/null +++ b/tests/testthat/fixtures/linear/make_lm.R @@ -0,0 +1,22 @@ +# Make some data +lm_data <- data.frame(X = Loblolly$age, + Y = Loblolly$height, + N = nrow(Loblolly)) + +saveRDS(lm_data, file = "tests/testthat/fixtures/linear/lm_data.rds") + +# Run a baseline model for the data +suppressWarnings( #Suppresses stan warnings + lm_baseline <- rmot_model("linear") |> + rmot_assign_data(X = lm_data$X, + Y = lm_data$Y, + N = nrow(lm_data)) |> + rmot_run(chains = 1, iter = 300, verbose = FALSE, show_messages = FALSE, seed=1) +) + +# Extract 100 rows of the samples +lm_baseline_output <- rstan::extract(lm_baseline, permuted = FALSE, inc_warmup = FALSE) |> + as.data.frame() |> head(n=100) + +# Save the output to compare againest +saveRDS(lm_baseline_output, "tests/testthat/fixtures/linear/lm_baseline_output.rds") diff --git a/tests/testthat/helper-data_generation.R b/tests/testthat/helper-data_generation.R new file mode 100644 index 0000000..23f51e9 --- /dev/null +++ b/tests/testthat/helper-data_generation.R @@ -0,0 +1,99 @@ +#General functions for data generation +#Runge-Kutta 4th order +rmot_rk4_est <- function(y_0, DE, pars, step_size, n_step){ + runge_kutta_int <- c(y_0) + for(i in 2:n_step){ + k1 <- DE(runge_kutta_int[i-1], pars) + k2 <- DE((runge_kutta_int[i-1] + step_size*k1/2), pars) + k3 <- DE((runge_kutta_int[i-1] + step_size*k2/2), pars) + k4 <- DE((runge_kutta_int[i-1] + step_size*k3), pars) + + runge_kutta_int[i] <- runge_kutta_int[i-1] + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*step_size + } + return(runge_kutta_int) +} + +rmot_build_true_test_data <- function(n_ind, n_obs, interval, + DE_pars, initial_conditions, DE){ + time <- seq(from = 0, by = interval, length.out = n_obs) + + true_data <- data.frame() + for(i in 1:n_ind){ + #Use RK4 with small step size to project forward from initial size + runge_kutta_int <- rmot_rk4_est(y_0 = initial_conditions[i,], + DE = DE, + pars = DE_pars[i,], + step_size = 0.1, + n_step = (1 + n_obs*interval/0.1)) + + #Take a subset of the estimates which are in line with survey structure + runge_kutta_survey <- runge_kutta_int[seq(from=1, + to = length(runge_kutta_int), + by = (length(runge_kutta_int)/n_obs))] + + data_temp <- data.frame( #Build data frame + y_true = runge_kutta_survey, + time = time, + ind_id = rep(i, times=n_obs) + ) + + #Concatenate data + true_data <- rbind(true_data, data_temp) + } + + return(true_data) +} + +rmot_add_normal_error <- function(y, sigma_e=0.001){ + return(y + rnorm(length(y), mean=0, sd=sigma_e)) +} + +#Save data to files +rmot_export_test_data <- function(n_obs_per_ind, + n_ind, + y_obs, + time, + DE_pars, + initial_conditions, + true_data, + model_name){ + + single_ind_data <- list( + step_size = 1, #Number for model RK4 alg + n_obs = n_obs_per_ind, #Number + y_obs = y_obs[1:n_obs_per_ind], #Vector indexed by n_obs + obs_index = 1:n_obs_per_ind, #Vector indexed by n_obs + time = time, #Vector indexed by n_obs + y_0_obs = y_obs[1], #Number + n_pars = ncol(DE_pars), #Number + single_true_data = list( + DE_pars = DE_pars[1,], + initial_conditions = initial_conditions[1,], + ind_id = 1, + true_data = true_data[1:n_obs_per_ind,] + ) + ) + + multi_ind_data <- list( + step_size = 1, #Number + n_obs = length(y_obs), #Number + n_ind = n_ind, #Number + y_obs = y_obs, #Vector indexed by n_obs + obs_index = rep(1:n_obs_per_ind, times = n_ind), #Vector indexed by n_obs + time = rep(time, times=n_ind), #Vector indexed by n_obs + ind_id = sort(rep(1:n_ind, times = n_obs_per_ind)), #Vector indexed by n_obs + y_0_obs = y_obs[seq(from = 1, to=n_ind*n_obs_per_ind, by=n_obs_per_ind)], #Vector indexed by n_ind + n_pars = ncol(DE_pars), + multi_true_data = list( + DE_pars = DE_pars, + initial_conditions = initial_conditions, + ind_id = c(1:n_ind), + true_data = true_data + ) + ) + + if(! dir.exists(test_path("fixtures", model_name))) dir.create(test_path("fixtures", model_name)) + filename <- paste("tests/testthat/fixtures", "/", model_name, "/", model_name, "_data", sep="") + saveRDS(single_ind_data, file=paste(filename, "single_ind.rds", sep="_")) + saveRDS(multi_ind_data, file=paste(filename, "multi_ind.rds", sep="_")) +} diff --git a/tests/testthat/helper-testing.R b/tests/testthat/helper-testing.R new file mode 100644 index 0000000..1e105f6 --- /dev/null +++ b/tests/testthat/helper-testing.R @@ -0,0 +1,72 @@ +rmot_test_model_functions <- function(model_name){ + single_ind_name <- paste0(model_name, "_single_ind") + multi_ind_name <- paste0(model_name, "_multi_ind") + + # Single individual + expect_named(rmot_model(single_ind_name)) + expect_type(rmot_model(single_ind_name), "list") + expect_visible(rmot_model(single_ind_name)) + #Multiple individuals + expect_named(rmot_model(multi_ind_name)) + expect_type(rmot_model(multi_ind_name), "list") + expect_visible(rmot_model(multi_ind_name)) +} + +rmot_test_single_individual <- function(model_name, + par_names){ + data <- readRDS(test_path("fixtures", model_name, + paste0(model_name, "_data_single_ind.rds"))) + + # Test constant single individual + suppressWarnings( #Suppresses stan warnings + single_ind_test <- rmot_model(paste0(model_name, "_single_ind")) |> + rmot_assign_data(n_obs = data$n_obs, #integer + y_obs = data$y_obs, + obs_index = data$obs_index, #vector length N_obs + time = data$time, #Vector length N_obs + y_0_obs = data$y_0_obs #vector length N_ind + ) |> + rmot_run(chains = 1, iter = 1000, verbose = FALSE, show_messages = FALSE) + ) + + # Extract samples and check if parameter estimates are reasonable. + ind_samples <- rstan::extract(single_ind_test, permuted = TRUE, + inc_warmup = FALSE) + par_ests <- c() + for(i in length(par_names)){ + par_ests[i] <- mean(ind_samples[[par_names[i]]]) + } + + initial_condition <- mean(ind_samples$ind_y_0) + expect_equal(par_ests, data$single_true_data$DE_pars, + tolerance = 1e-1) + expect_equal(initial_condition, data$single_true_data$initial_conditions, + tolerance = 1e-1) + + # hecks for output existence and type + expect_visible(single_ind_test) + expect_s4_class(single_ind_test, "stanfit") +} + +rmot_test_multi_individual <- function(model_name, data, est_dim){ + # Test constant multi-individual + suppressWarnings( #Suppresses stan warnings + multi_ind_test <- rmot_model(paste0(model_name, "_multi_ind")) |> + rmot_assign_data(n_obs = data$n_obs, #integer + n_ind = data$n_ind, #integer + y_obs = data$y_obs, #vector length N_obs + obs_index = data$obs_index, #vector length N_obs + time = data$time, #Vector length N_obs + ind_id = data$ind_id, #Vector length N_obs + y_0_obs = data$y_0_obs #vector length N_ind + ) |> + rmot_run(chains = 2, iter = 100, verbose = FALSE, show_messages = FALSE) + ) + + # Extract samples + multi_samples <- rstan::extract(multi_ind_test, permuted = FALSE, inc_warmup = TRUE) + expect_equal(dim(multi_samples), c(100, 2, est_dim)) + + expect_visible(multi_ind_test) + expect_s4_class(multi_ind_test, "stanfit") +} diff --git a/tests/testthat/test-rmot_models.R b/tests/testthat/test-rmot_models.R deleted file mode 100644 index 182fc5a..0000000 --- a/tests/testthat/test-rmot_models.R +++ /dev/null @@ -1,19 +0,0 @@ -test_that("Execution and output", { - expect_named(rmot_model("linear")) - expect_type(rmot_model("linear"), "list") - expect_visible(rmot_model("linear")) - - expect_named(rmot_model("constant_single")) - expect_type(rmot_model("constant_single"), "list") - expect_visible(rmot_model("constant_single")) - - lm_test <- rmot_model("linear") |> - rmot_assign_data(X = Loblolly$age, - Y = Loblolly$height, - N = nrow(Loblolly)) |> - rmot_run(chains = 2, iter = 300, verbose = FALSE, show_messages = FALSE) - - expect_visible(lm_test) - - expect_s4_class(lm_test, "stanfit") -}) diff --git a/tests/testthat/test-rmot_models_constant.R b/tests/testthat/test-rmot_models_constant.R new file mode 100644 index 0000000..56ba31f --- /dev/null +++ b/tests/testthat/test-rmot_models_constant.R @@ -0,0 +1,36 @@ +#Testing for constant model +test_that("Model structures: Constant", { + # Single individual + expect_named(rmot_model("constant_single_ind")) + expect_type(rmot_model("constant_single_ind"), "list") + expect_visible(rmot_model("constant_single_ind")) + #Multiple individuals + expect_named(rmot_model("constant_multi_ind")) + expect_type(rmot_model("constant_multi_ind"), "list") + expect_visible(rmot_model("constant_multi_ind")) +}) + +test_that("Execution: Constant single individual", { + model_name <- "constant" + par_names <- "ind_beta" + + rmot_test_single_individual(model_name, par_names) +}) + +test_that("Execution: Constant multiple individuals", { + model_name <- "constant" + + data <- readRDS(test_path("fixtures", "constant", + "constant_data_multi_ind.rds")) + + #Dimension is: + est_dim <- data$n_ind + #Initial condition + data$n_pars * data$n_ind + #Individual parameters + data$n_pars * 2 + #Population parameters + 1 + #Global error + data$n_obs + #y_ij + data$n_obs + #Delta y_ij + 1 #lp__ + + rmot_test_multi_individual(model_name, data, est_dim) +}) diff --git a/tests/testthat/test-rmot_models_linearreg.R b/tests/testthat/test-rmot_models_linearreg.R new file mode 100644 index 0000000..8abad90 --- /dev/null +++ b/tests/testthat/test-rmot_models_linearreg.R @@ -0,0 +1,23 @@ +test_that("Model structures: Linear", { + expect_named(rmot_model("linear")) + expect_type(rmot_model("linear"), "list") + expect_visible(rmot_model("linear")) +}) + +test_that("Execution and output: Linear", { + lm_data <- readRDS(test_path("fixtures", "linear", "lm_data.rds")) + lm_baseline_output <- readRDS(test_path("fixtures", "linear", "lm_baseline_output.rds")) + + # Test linear model + suppressWarnings( #Suppresses stan warnings + lm_test <- rmot_model("linear") |> + rmot_assign_data(X = lm_data$X, + Y = lm_data$Y, + N = nrow(lm_data)) |> + rmot_run(chains = 1, iter = 300, verbose = FALSE, show_messages = FALSE, + seed = 1) + ) + + expect_visible(lm_test) + expect_s4_class(lm_test, "stanfit") +}) diff --git a/vignettes/rmot.Rmd b/vignettes/rmot.Rmd index 6b708ff..378eee2 100644 --- a/vignettes/rmot.Rmd +++ b/vignettes/rmot.Rmd @@ -3,8 +3,10 @@ title: "rmot" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{rmot} - %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +editor_options: + chunk_output_type: console --- ```{r, include = FALSE} @@ -41,28 +43,47 @@ rmot_model("linear") |> rmot_run() ``` -### Constant Growth - Single Species +### Constant Growth - Single Individual + +```{r, eval=FALSE} +# Create data +true <- c(seq(from=1, by=2.5, length.out=7)) +y <- true + true*rnorm(n=7, 0, 0.02) + rnorm(n=7, 0, 0.2) + +# Configure model template +const_single_ind_output <- rmot_model("constant_single_ind") |> + rmot_assign_data(n_obs = length(y), #integer + y_obs = y, #vector length N_obs + obs_index = 1:7, #vector length N_obs + time = rep(5, times=length(y)), #Vector length N_obs + y_0_obs = y[1] #vector length N_ind + ) |> + rmot_run() + +# Look at output +const_single_ind_output +``` + +### Constant Growth - Multiple Individuals -```{r} +```{r, eval=FALSE} # Create data true <- c(seq(from=1, by=2.5, length.out=7), seq(from=2, by=2, length.out=7)) y <- true + true*rnorm(n=14, 0, 0.02) + rnorm(n=14, 0, 0.2) # Configure model template -cgs_output <- rmot_model("constant_single") |> - rmot_assign_data(N_obs = length(y), - N_ind = 2, - S_obs = y, - census = rep(seq(from=1, to=7, by=1), times=2), #vector length N_obs - census_interval = rep(5, times=length(y)), #Vector length N_obs - id_factor = c(rep(1, times=7), rep(2, times=7)), #Vector length N_obs - S_0_obs = y[c(1, 8)] #vector length N_ind +const_multi_ind_output <- rmot_model("constant_multi_ind") |> + rmot_assign_data(n_obs = length(y), #integer + n_ind = 2, #integer + y_obs = y, #vector length N_obs + obs_index = rep(seq(from=1, to=7, by=1), times=2), #vector length N_obs + time = rep(5, times=length(y)), #Vector length N_obs + ind_id = c(rep(1, times=7), rep(2, times=7)), #Vector length N_obs + y_0_obs = y[c(1, 8)] #vector length N_ind ) |> rmot_run() # Look at output -cgs_output +const_multi_ind_output ``` - -