From acae90c22a6fb5bfda703b6cc9eab5a9f1eb6567 Mon Sep 17 00:00:00 2001 From: Jaime Ashander Date: Mon, 23 Feb 2015 14:58:43 -0800 Subject: [PATCH] Full simulator with Luis Features: simulate autocorrelated environment with trait and population dynamics. Also ceiling density regulation. Expand README and add LICENSE (MIT) Add Luis as package author Add testing: - unit testing at highest (ts simulation) level - env sim - stub for tests - a regression test Ignore build and autosave files - Add rds files to ignore --- .gitignore | 6 + DESCRIPTION | 15 +- LICENSE | 21 ++ NAMESPACE | 20 +- R/RcppExports.R | 102 ++++--- R/envsim.R | 52 ---- R/traitfitpop.R | 194 ------------- README.Rmd | 12 + README.md | 2 - rcppexample.R => inst/rcppexample.R | 8 +- man/Beta.Rd | 26 -- man/Env_shift.Rd | 27 -- man/Env_smooth.Rd | 26 -- man/Log_R_bar.Rd | 21 -- man/Log_W_bar.Rd | 8 +- man/Log_W_bar_lethal.Rd | 19 -- man/Pheno_demo.Rd | 28 -- man/Pheno_demo_econ.Rd | 32 --- man/Pheno_demo_lande.Rd | 28 -- man/Pheno_lande.Rd | 27 -- man/R_bar.Rd | 19 -- man/R_bar_ceiling.Rd | 25 ++ man/{R_bar_dd.Rd => R_bar_thetalog.Rd} | 10 +- man/W_bar.Rd | 8 +- man/make_env.Rd | 19 ++ man/{pdTS.Rd => simulate_pheno_ts.Rd} | 6 +- man/stab_curve.Rd | 19 -- src/RcppExports.cpp | 208 ++++++-------- src/traitfitpop.cpp | 360 ++++++++++++------------- tests/testthat.R | 4 + tests/testthat/test-demography.R | 6 + tests/testthat/test-simulate-env.R | 25 ++ tests/testthat/test-simulate-ts.R | 55 ++++ 33 files changed, 511 insertions(+), 927 deletions(-) create mode 100644 LICENSE delete mode 100644 R/envsim.R create mode 100644 README.Rmd delete mode 100644 README.md rename rcppexample.R => inst/rcppexample.R (93%) delete mode 100644 man/Beta.Rd delete mode 100644 man/Env_shift.Rd delete mode 100644 man/Env_smooth.Rd delete mode 100644 man/Log_R_bar.Rd delete mode 100644 man/Log_W_bar_lethal.Rd delete mode 100644 man/Pheno_demo.Rd delete mode 100644 man/Pheno_demo_econ.Rd delete mode 100644 man/Pheno_demo_lande.Rd delete mode 100644 man/Pheno_lande.Rd delete mode 100644 man/R_bar.Rd create mode 100644 man/R_bar_ceiling.Rd rename man/{R_bar_dd.Rd => R_bar_thetalog.Rd} (71%) create mode 100644 man/make_env.Rd rename man/{pdTS.Rd => simulate_pheno_ts.Rd} (88%) delete mode 100644 man/stab_curve.Rd create mode 100644 tests/testthat.R create mode 100644 tests/testthat/test-demography.R create mode 100644 tests/testthat/test-simulate-env.R create mode 100644 tests/testthat/test-simulate-ts.R diff --git a/.gitignore b/.gitignore index dc528e5..01daad1 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,7 @@ .Rhistory +*.rds +..Rcheck +*.o +*.so +.\#* +\#*\# diff --git a/DESCRIPTION b/DESCRIPTION index 1aef5b2..6db2b3d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,15 +4,16 @@ Description: A set of codes to simulate evolutionary, plastic and ecological constantly moving optima Title: Simulations of eco-evo trait and demographic dynamics under changing environments -Version: 0.0-1 -Date: 2013-10-31 +Version: 0.1.0 +Date: 2015-11-17 Authors@R: c(person("Jaime", "Ashander", role = c("aut", "cre"), email= - 'jashander@ucdavis.edu')) + "jashander@ucdavis.edu"), person("Luis-Miguel", "Chevin", role=c("aut"), + email="luis-miguel.chevin@cefe.cnrs.fr")) Depends: RcppArmadillo (>= 0.4.3) Imports: - Rcpp, - deSolve, - MASS + Rcpp LinkingTo: Rcpp, RcppArmadillo -License: BSD +License: MIT +Suggests: + testthat diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..c47ffbc --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2015 Jaime Ashander and Luis-Miguel Chevin + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/NAMESPACE b/NAMESPACE index a554c3b..ba1586f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,23 +1,11 @@ # Generated by roxygen2 (4.0.2): do not edit by hand -export(Beta) -export(Env_shift) -export(Env_shift_cpp) -export(Env_smooth) export(G) -export(Log_R_bar) export(Log_W_bar) -export(Log_W_bar_lethal) -export(Pheno_demo) -export(Pheno_demo_econ) -export(Pheno_demo_lande) -export(Pheno_lande) -export(R_bar) -export(R_bar_dd) +export(R_bar_ceiling) +export(R_bar_thetalog) export(Va) export(W_bar) -export(pdTS) -export(stab_curve) -import(MASS) -import(deSolve) +export(make_env) +export(simulate_pheno_ts) useDynLib(phenoecosim) diff --git a/R/RcppExports.R b/R/RcppExports.R index cc92efe..ee955b4 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,70 +1,71 @@ # This file was generated by Rcpp::compileAttributes # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +#' Compute a white noise environment with a shift +#' @param t the time point +#' @param env_args other args: +#' t.jump the location for the jump, +#' delta env change that takes place at t.jump +#' sd noise in env +#' sdc noise in cue +#' @details nada +NULL + #' Compute mean fitness under stabilizing selection #' @param zbar mean trait prior to selection #' @param theta enironmental optimum -#' @param Oz2 strength of stabilizing selection -#' @param gamma 1/(Oz2 + Vz2), with Vz2 phenotypic variance -#' @details no details +#' @param oz2 strength of stabilizing selection +#' @param gamma 1/(oz2 + Vz2), with Vz2 phenotypic variance +#' @details IN USE #' @export -Log_W_bar <- function(zbar, theta, Oz2, gamma) { - .Call('phenoecosim_Log_W_bar', PACKAGE = 'phenoecosim', zbar, theta, Oz2, gamma) +Log_W_bar <- function(zbar, theta, oz2, gamma) { + .Call('phenoecosim_Log_W_bar', PACKAGE = 'phenoecosim', zbar, theta, oz2, gamma) } #' Compute mean fitness under stabilizing selection #' @inheritParams Log_W_bar -#' @param LOG (=TRUE) whether to return the log of fitness -#' @details no details +#' @param LOG (=TRUE) whether to return the log of fitness +#' @details IN USE #' @export -W_bar <- function(zbar, theta, Oz2, gamma, LOG) { - .Call('phenoecosim_W_bar', PACKAGE = 'phenoecosim', zbar, theta, Oz2, gamma, LOG) +W_bar <- function(zbar, theta, oz2, gamma, LOG) { + .Call('phenoecosim_W_bar', PACKAGE = 'phenoecosim', zbar, theta, oz2, gamma, LOG) } -#' Compute population growth rate under stabilizing selection +#' Compute population growth rate under stabilizing selection and no regulation #' @param R0 basic reproductive number #' @param Wbar average fitness -#' @details Assumes density-independent -#' @export -R_bar <- function(R0, Wbar) { - .Call('phenoecosim_R_bar', PACKAGE = 'phenoecosim', R0, Wbar) -} - -#' Compute LOG population growth rate under stabilizing selection -#' @inheritParams R_bar -#' @param logWbar log average fitness -#' @details Assumes theta-logistic population regulation +#' @param N number of individuals in this generation +#' @param K carrying capacity +#' @details Assumes ceiling population regulation #' would be good to have separate DD function specified after #' chevin and lande 2010 #' @export -Log_R_bar <- function(R0, logWbar) { - .Call('phenoecosim_Log_R_bar', PACKAGE = 'phenoecosim', R0, logWbar) +R_bar <- function(R0, Wbar, N) { + .Call('phenoecosim_R_bar', PACKAGE = 'phenoecosim', R0, Wbar, N) } -#' Compute population growth rate under stabilizing selection -#' @inheritParams R_bar +#' Compute population growth rate under stabilizing selection and ceiling regulation +#' @param R0 basic reproductive number +#' @param Wbar average fitness #' @param N number of individuals in this generation #' @param K carrying capacity -#' @param thetaL theta-logistic parameter for density dependence -#' @details Assumes theta-logistic population regulation +#' @details Assumes ceiling population regulation #' would be good to have separate DD function specified after #' chevin and lande 2010 #' @export -R_bar_dd <- function(R0, Wbar, N, K, thetaL) { - .Call('phenoecosim_R_bar_dd', PACKAGE = 'phenoecosim', R0, Wbar, N, K, thetaL) +R_bar_ceiling <- function(R0, Wbar, N, K) { + .Call('phenoecosim_R_bar_ceiling', PACKAGE = 'phenoecosim', R0, Wbar, N, K) } -#' Selection on plasticity as as function of environment assuming stabilizing selection -#' @param gamma = 1/(Oz2 + Vz2), with Vz2 phenotypic variance -#' @param A the environmental optimum RN int -#' @param B the environmental optimum RN slope -#' @param a current value of RN int -#' @param b current value of RN slope -#' @param e_t environment now -#' @param e_plast env that cues plasticity +#' Compute population growth rate under stabilizing selection and theta-logistic regulation +#' @inheritParams R_bar_ceiling +#' @param thetaL theta-logistic parameter for density dependence +#' @details Assumes theta-logistic population regulation +#' would be good to have separate DD function specified after +#' chevin and lande 2010 #' @export -Beta <- function(gamma, A, B, a, b, e_t, e_plast) { - .Call('phenoecosim_Beta', PACKAGE = 'phenoecosim', gamma, A, B, a, b, e_t, e_plast) +R_bar_thetalog <- function(R0, Wbar, N, K, thetaL) { + .Call('phenoecosim_R_bar_thetalog', PACKAGE = 'phenoecosim', R0, Wbar, N, K, thetaL) } #' Va additive genetic variance in the phenotype as a function of the environment @@ -75,28 +76,25 @@ Va <- function(env, GG) { .Call('phenoecosim_Va', PACKAGE = 'phenoecosim', env, GG) } -#' Compute a white noise environment with a shift -#' @param t the time point -#' @param env_args other args: -#' t.jump the location for the jump, -#' delta env change that takes place at t.jump -#' sd noise in env -#' sdc noise in cue -#' @details nada +#' Compute correctly-scaled environment of development and selection +#' @param env_args environment args +#' @param T how many gens +#' @details given a pre-generated vector of random environments #' @export -Env_shift_cpp <- function(t, env_args) { - .Call('phenoecosim_Env_shift_cpp', PACKAGE = 'phenoecosim', t, env_args) +make_env <- function(T, env_args) { + .Call('phenoecosim_make_env', PACKAGE = 'phenoecosim', T, env_args) } -#' Compute phenotypic dynamic Time Series of trait + demographic change under stabilizing selection as +#' Compute phenotypic dynamic Time Series of trait + demographic change under stabilizing selection as #' function of environment (after Lande Chevin) #' @param T end time, assuming start time of 1 #' @param X parameters (z, a, b, wbar, logN, theta) -#' @param params a list with (gamma_sh, omegaz, A, B, R0, Va, Vb, delta, sigma_xi, rho_tau, fractgen) +#' @param params a list with (gamma_sh, omegaz, A, B, R0, var_a, Vb, delta, sigma_xi, rho_tau, fractgen) #' @param env_args extra args for env.fn #' @details NB - for now assume Tchange = 0 and demography after CL 2010 +#' @value a long matrix #' @export -pdTS <- function(T, X, params, env_args) { - .Call('phenoecosim_pdTS', PACKAGE = 'phenoecosim', T, X, params, env_args) +simulate_pheno_ts <- function(T, X, params, env_args) { + .Call('phenoecosim_simulate_pheno_ts', PACKAGE = 'phenoecosim', T, X, params, env_args) } diff --git a/R/envsim.R b/R/envsim.R deleted file mode 100644 index 915c58e..0000000 --- a/R/envsim.R +++ /dev/null @@ -1,52 +0,0 @@ -#' Compute a white noise environment with a shift -#' @param env previous year's environment NOT used in current formulation -#' @param t the time point -#' @param rho the correlation between environment of selection and the cue -#' @param env.args other args: -#' t.jump the location for the jump, -#' delta env change that takes place at t.jump -#' sd noise in env -#' sdc noise in cue -#' @details nada -#' @import MASS -#' @export -Env_shift <- function(env, t, rho, env.args){ - ## signature function(environment, time, ...) - t.jump <- env.args[['t.jump']] - k <- env.args[['k']] - sd <- env.args[['sd']] - sdc <- env.args[['sdc']] - ## for red noise do somehting like k * env + sqrt(1 - k^2) * rnorm (n=1, mean=0, sd=sd)) but hard - ## to figure this out with the cue / correlation thing - out <- mvrnorm(1, mu = c(0,0) , Sigma = matrix(c(sd, sd*sdc*rho, sd*sdc*rho, sdc), nrow=2)) - if( t < t.jump) - return(out) - else if(t >= t.jump) - return(env.args[['delta']] + out ) - else - return("You should never see this") -} - -#' Simulate one step of a linearly changing environment, possibly with noise -#' @param env environment in previous time step -#' @param t the time point -#' @param rho the correlation between environment of selection and the cue -#' @param env.args other args: -#' eta rate of environmental change, -#' sd noise in env -#' sdc noise in cue -#' @details other args in env args -#' @import MASS -#' @export -Env_smooth <- function(env, t, rho, env.args){ - ## signature function(environment, time, ...) - eta <- env.args[['eta']] - sdc <- env.args[['sdc']] - sd <- env.args[['sd']] - ## cue <- env.args[['cue']] ## want to use as current cue, but ... argument passing harder, maybe - eta + mvrnorm(1, mu = c(env,env) , Sigma = matrix(c(sd, sd*sdc*rho, sd*sdc*rho, sdc), nrow=2)) -} - - - -#' Env with autocorrelation: need to add this and makesure I can re-obtain results from democonstriants_ms diff --git a/R/traitfitpop.R b/R/traitfitpop.R index c9d6c18..126558d 100644 --- a/R/traitfitpop.R +++ b/R/traitfitpop.R @@ -1,23 +1,3 @@ -#' Mean fitness -- with constraint, need to make as function of environment... -#' @param x the environment -#' @param LL lethal limit (symmetric) -#' @details fitness surf infinitesimal optimum need to revisit -#' @export -stab_curve<- function(x, LL) -(x-LL)*(x+LL) * 1/LL^2 - -#' Mean fitness -- with constraint, need to make as function of environment... -#' @param env the environment -#' @param LL lethal limit -#' @details need to revisit this -#' @export -Log_W_bar_lethal <- function(env, LL=4) { - if(env > -LL & env < LL) - return(stab_curve(env, LL)) - else - return(0) -} - - #' G matrix #' @param Gaa additive genetic variance in RN intercept #' @param Gbb additive genetic variance in RN slope @@ -25,177 +5,3 @@ Log_W_bar_lethal <- function(env, LL=4) { #' @export G <- function(Gaa, Gbb, Gab) matrix( c(Gaa, Gab, Gab, Gbb), nrow=2) #eqn 3a - -#' Compute trait change under stabilizing selection as function of environment -#' @param t timestep -#' @param X parameters (a, b, env) -#' @param p (rho, gamma, A, B) -#' @param G the (constant) G matrix -#' @param env.fn function to compute environment with signature function(environment, time, ...) -#' @param env.args extra args for env.fn -#' @details function signature for use with deSolve and 'iterate' -#' @import deSolve -#' @export -Pheno_lande <- function(t, X, p, G, env.fn, env.args) { - ## state vars - a <- X[1] # elevation - b <- X[2] # slope - env <- X[3] - ## params - rho <- p[1] - gamma <- p[2] - Oz2 <- p[3] - A <- p[4] - B <- p[5] - EE <- env.fn(env, t, rho, env.args) - env <- EE[1] - e.plast <- EE[2] - z.bar <- a + b * e.plast - log.W.bar <- W_bar(z.bar, A + B*env, Oz2, gamma) - beta <- Beta(gamma, A, B, a, b, env, e.plast) - X <- c(a, b) + G %*% beta - va <- Va(e.plast, G) - list(c(X, env, z.bar, va)) -} - -#' Compute trait + demographic change under stabilizing selection as function of environment -#' @param t timestep -#' @param X parameters (a, b, env, N) -#' @param p (rho, gamma, A, B, R0) -#' @param G the (constant) G matrix -#' @param env.fn function to compute environment with signature function(environment, time, ...) -#' @param env.args extra args for env.fn -#' @details function signature for use with deSolve and 'iterate', imposes fitness based on -#' assumed tolerance curve -#' @import deSolve -#' @export -Pheno_demo <- function(t, X, p, G, env.fn, env.args) { - ## state vars - a <- X[1] # elevation - b <- X[2] # slope - env <- X[3] - N = X[4] - - ## params - rho <- p[1] - gamma <- p[2] - Oz2 <- p[3] - - A <- p[4] - B <- p[5] - R0 <- p[6] - - EE <- env.fn(env, t, rho, env.args) - env <- EE[1] - e.plast <- EE[2] - - log.W.bar <- Log_W_bar_lethal(env) - N <- log.W.bar * R0 * N ## LOG demo update based on mean fitness - - beta <- Beta(gamma, A, B, a, b, env, e.plast) - X <- c(a, b) + G %*% beta - va <- Va(e.plast, G) - list(c(X, env, N, va)) -} - - -#' Compute trait + demographic change under stabilizing selection as function of environment via Lande Chevin -#' @param t timestep -#' @param X parameters (a, b, env, N) -#' @param p (rho, gamma, A, B, R0, K, thetaL) -#' @param G the (constant) G matrix -#' @param env.fn function to compute environment with signature function(environment, time, ...) -#' @param env.args extra args for env.fn -#' @details function signature for use with deSolve and 'iterate', imposes fitness based on Lande -#' and demography after CL 2010 -#' @import deSolve -#' @export -Pheno_demo_lande <- function(t, X, p, G, env.fn, env.args) { - ## state vars - a <- X[1] # elevation - b <- X[2] # slope - env <- X[3] - N = X[4] - - ## params - rho <- p[1] - gamma <- p[2] - Oz2 <- p[3] - - A <- p[4] - B <- p[5] - R0 <- p[6] - K <- p[7] - thetaL <- p[8] - - EE <- env.fn(env, t, rho, env.args) - env <- EE[1] - e.plast <- EE[2] - - z.bar <- a + b * e.plast - log.W.bar <- W_bar(z.bar, A + B*env, Oz2, gamma) - N <- R_bar( R0, exp(log.W.bar)) # for _dd N , K, thetaL) * N ## demo update based on mean fitness - - beta <- Beta(gamma, A, B, a, b, env, e.plast) - X <- c(a, b) + G %*% beta - va <- Va(e.plast, G) - list(c(X, env, N, va)) -} - -#' Compute trait + demographic change under stabilizing selection as function of environment via Lande Chevin -#' but with intervention -#' @param t timestep -#' @param X parameters (a, b, env, N) -#' @param p (rho, gamma, A, B, R0, K, thetaL) -#' @param G the (constant) G matrix -#' @param env.fn function to compute environment with signature function(environment, time, ...) -#' @param env.args extra args for env.fn -#' @param interv a function of the time, policy to modify the cue -#' @details function signature for use with deSolve and 'iterate', imposes fitness based on Lande -#' and demography after CL 2010 -#' @import deSolve -#' @export -Pheno_demo_econ <- function(t, X, p, G, env.fn, env.args, interv) { - ## state vars - a <- X[1] # elevation - b <- X[2] # slope - env <- X[3] - N = X[4] - - ## params - rho <- p[1] - gamma <- p[2] - Oz2 <- p[3] - - A <- p[4] - B <- p[5] - R0 <- p[6] - K <- p[7] - thetaL <- p[8] - - EE <- env.fn(env, t, rho, env.args) - env <- EE[1] - e.plast <- EE[2] + interv(t) - - z.bar <- a + b * e.plast - log.W.bar <- W_bar(z.bar, A + B*env, Oz2, gamma) - N <- R_bar(R0, exp(log.W.bar)) # for _dd N, K, thetaL) * N ## demo update based on mean fitness - - beta <- Beta(gamma, A, B, a, b, env, e.plast) - X <- c(a, b) + G %*% beta - va <- Va(e.plast, G) - list(c(X, env, N, va)) -} - - - -## for examples use #' examples \dontrun{} - - -## tolerance-functions -## fitsurf Lande (mean) -- not sure what i'm doing with this -## lande.fit <- function(x) exp( - (x * (E.B-B))^2/(2 * 50)) - - - - diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..80a792f --- /dev/null +++ b/README.Rmd @@ -0,0 +1,12 @@ +phenoecosim +=========== + +## Simulations of eco-evo trait and demographic dynamics under changing environments + +Jaime Ashander and Luis-Miguel Chevin + +Simulate evolutionary, plastic, and ecological responses to +changing environments, including discrete shifts and constantly +moving optima. + +Built with [Armadillo](http://arma.sourceforge.net/) and [Rcpp](http://www.rcpp.org/) (via RcppArmadillo package). diff --git a/README.md b/README.md deleted file mode 100644 index d20242a..0000000 --- a/README.md +++ /dev/null @@ -1,2 +0,0 @@ -phenoecosim -=========== diff --git a/rcppexample.R b/inst/rcppexample.R similarity index 93% rename from rcppexample.R rename to inst/rcppexample.R index 3fbc9ea..c416c50 100644 --- a/rcppexample.R +++ b/inst/rcppexample.R @@ -27,8 +27,10 @@ cppFunction('double testarma(arma::vec x, arma::rowvec y) { ## I get one warning: ## ld: warning: directory '/usr/local/lib/x86_64' following -L not found +## os x 10.10 +## ld: warning: directory not found for option '-L/usr/local/lib/gcc/x86_64-apple-darwin13.0.0/4.8.2' testarma(1:2, 2:3) - +# 8 ## below also works with sessionInfo() as above (this code is example from ?Rcpp::cppFunction) cppFunction(depends = "RcppArmadillo", @@ -112,6 +114,10 @@ Rcpp::NumericVector colNorm(Rcpp::NumericMatrix sM) { return n; // return vector }') + +## fails on 10.10 with +## Error in if (file.exists(pkgHeaderPath)) { : argument is of length zero +## Calls: cppFunction -> .linkingToIncludes code <- ' arma::mat GG = Rcpp::as(G); arma::colvec env = Rcpp::as(e); diff --git a/man/Beta.Rd b/man/Beta.Rd deleted file mode 100644 index 8c3de0e..0000000 --- a/man/Beta.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Beta} -\alias{Beta} -\title{Selection on plasticity as as function of environment assuming stabilizing selection} -\usage{ -Beta(gamma, A, B, a, b, e_t, e_plast) -} -\arguments{ -\item{gamma}{= 1/(Oz2 + Vz2), with Vz2 phenotypic variance} - -\item{A}{the environmental optimum RN int} - -\item{B}{the environmental optimum RN slope} - -\item{a}{current value of RN int} - -\item{b}{current value of RN slope} - -\item{e_t}{environment now} - -\item{e_plast}{env that cues plasticity} -} -\description{ -Selection on plasticity as as function of environment assuming stabilizing selection -} - diff --git a/man/Env_shift.Rd b/man/Env_shift.Rd deleted file mode 100644 index 13f5325..0000000 --- a/man/Env_shift.Rd +++ /dev/null @@ -1,27 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Env_shift} -\alias{Env_shift} -\title{Compute a white noise environment with a shift} -\usage{ -Env_shift(env, t, rho, env.args) -} -\arguments{ -\item{env}{previous year's environment NOT used in current formulation} - -\item{t}{the time point} - -\item{rho}{the correlation between environment of selection and the cue} - -\item{env.args}{other args: -t.jump the location for the jump, -delta env change that takes place at t.jump -sd noise in env -sdc noise in cue} -} -\description{ -Compute a white noise environment with a shift -} -\details{ -nada -} - diff --git a/man/Env_smooth.Rd b/man/Env_smooth.Rd deleted file mode 100644 index 2784c1e..0000000 --- a/man/Env_smooth.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Env_smooth} -\alias{Env_smooth} -\title{Simulate one step of a linearly changing environment, possibly with noise} -\usage{ -Env_smooth(env, t, rho, env.args) -} -\arguments{ -\item{env}{environment in previous time step} - -\item{t}{the time point} - -\item{rho}{the correlation between environment of selection and the cue} - -\item{env.args}{other args: -eta rate of environmental change, -sd noise in env -sdc noise in cue} -} -\description{ -Simulate one step of a linearly changing environment, possibly with noise -} -\details{ -other args in env args -} - diff --git a/man/Log_R_bar.Rd b/man/Log_R_bar.Rd deleted file mode 100644 index a66916b..0000000 --- a/man/Log_R_bar.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Log_R_bar} -\alias{Log_R_bar} -\title{Compute LOG population growth rate under stabilizing selection} -\usage{ -Log_R_bar(R0, logWbar) -} -\arguments{ -\item{R0}{basic reproductive number} - -\item{logWbar}{log average fitness} -} -\description{ -Compute LOG population growth rate under stabilizing selection -} -\details{ -Assumes theta-logistic population regulation -would be good to have separate DD function specified after -chevin and lande 2010 -} - diff --git a/man/Log_W_bar.Rd b/man/Log_W_bar.Rd index 79e8229..6ec8e9e 100644 --- a/man/Log_W_bar.Rd +++ b/man/Log_W_bar.Rd @@ -3,21 +3,21 @@ \alias{Log_W_bar} \title{Compute mean fitness under stabilizing selection} \usage{ -Log_W_bar(zbar, theta, Oz2, gamma) +Log_W_bar(zbar, theta, oz2, gamma) } \arguments{ \item{zbar}{mean trait prior to selection} \item{theta}{enironmental optimum} -\item{Oz2}{strength of stabilizing selection} +\item{oz2}{strength of stabilizing selection} -\item{gamma}{1/(Oz2 + Vz2), with Vz2 phenotypic variance} +\item{gamma}{1/(oz2 + Vz2), with Vz2 phenotypic variance} } \description{ Compute mean fitness under stabilizing selection } \details{ -no details +IN USE } diff --git a/man/Log_W_bar_lethal.Rd b/man/Log_W_bar_lethal.Rd deleted file mode 100644 index b72cb26..0000000 --- a/man/Log_W_bar_lethal.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Log_W_bar_lethal} -\alias{Log_W_bar_lethal} -\title{Mean fitness -- with constraint, need to make as function of environment...} -\usage{ -Log_W_bar_lethal(env, LL = 4) -} -\arguments{ -\item{env}{the environment} - -\item{LL}{lethal limit} -} -\description{ -Mean fitness -- with constraint, need to make as function of environment... -} -\details{ -need to revisit this -} - diff --git a/man/Pheno_demo.Rd b/man/Pheno_demo.Rd deleted file mode 100644 index d053239..0000000 --- a/man/Pheno_demo.Rd +++ /dev/null @@ -1,28 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Pheno_demo} -\alias{Pheno_demo} -\title{Compute trait + demographic change under stabilizing selection as function of environment} -\usage{ -Pheno_demo(t, X, p, G, env.fn, env.args) -} -\arguments{ -\item{t}{timestep} - -\item{X}{parameters (a, b, env, N)} - -\item{p}{(rho, gamma, A, B, R0)} - -\item{G}{the (constant) G matrix} - -\item{env.fn}{function to compute environment with signature function(environment, time, ...)} - -\item{env.args}{extra args for env.fn} -} -\description{ -Compute trait + demographic change under stabilizing selection as function of environment -} -\details{ -function signature for use with deSolve and 'iterate', imposes fitness based on -assumed tolerance curve -} - diff --git a/man/Pheno_demo_econ.Rd b/man/Pheno_demo_econ.Rd deleted file mode 100644 index 8f48585..0000000 --- a/man/Pheno_demo_econ.Rd +++ /dev/null @@ -1,32 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Pheno_demo_econ} -\alias{Pheno_demo_econ} -\title{Compute trait + demographic change under stabilizing selection as function of environment via Lande Chevin -but with intervention} -\usage{ -Pheno_demo_econ(t, X, p, G, env.fn, env.args, interv) -} -\arguments{ -\item{t}{timestep} - -\item{X}{parameters (a, b, env, N)} - -\item{p}{(rho, gamma, A, B, R0, K, thetaL)} - -\item{G}{the (constant) G matrix} - -\item{env.fn}{function to compute environment with signature function(environment, time, ...)} - -\item{env.args}{extra args for env.fn} - -\item{interv}{a function of the time, policy to modify the cue} -} -\description{ -Compute trait + demographic change under stabilizing selection as function of environment via Lande Chevin -but with intervention -} -\details{ -function signature for use with deSolve and 'iterate', imposes fitness based on Lande -and demography after CL 2010 -} - diff --git a/man/Pheno_demo_lande.Rd b/man/Pheno_demo_lande.Rd deleted file mode 100644 index ec5b712..0000000 --- a/man/Pheno_demo_lande.Rd +++ /dev/null @@ -1,28 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Pheno_demo_lande} -\alias{Pheno_demo_lande} -\title{Compute trait + demographic change under stabilizing selection as function of environment via Lande Chevin} -\usage{ -Pheno_demo_lande(t, X, p, G, env.fn, env.args) -} -\arguments{ -\item{t}{timestep} - -\item{X}{parameters (a, b, env, N)} - -\item{p}{(rho, gamma, A, B, R0, K, thetaL)} - -\item{G}{the (constant) G matrix} - -\item{env.fn}{function to compute environment with signature function(environment, time, ...)} - -\item{env.args}{extra args for env.fn} -} -\description{ -Compute trait + demographic change under stabilizing selection as function of environment via Lande Chevin -} -\details{ -function signature for use with deSolve and 'iterate', imposes fitness based on Lande -and demography after CL 2010 -} - diff --git a/man/Pheno_lande.Rd b/man/Pheno_lande.Rd deleted file mode 100644 index 754f752..0000000 --- a/man/Pheno_lande.Rd +++ /dev/null @@ -1,27 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{Pheno_lande} -\alias{Pheno_lande} -\title{Compute trait change under stabilizing selection as function of environment} -\usage{ -Pheno_lande(t, X, p, G, env.fn, env.args) -} -\arguments{ -\item{t}{timestep} - -\item{X}{parameters (a, b, env)} - -\item{p}{(rho, gamma, A, B)} - -\item{G}{the (constant) G matrix} - -\item{env.fn}{function to compute environment with signature function(environment, time, ...)} - -\item{env.args}{extra args for env.fn} -} -\description{ -Compute trait change under stabilizing selection as function of environment -} -\details{ -function signature for use with deSolve and 'iterate' -} - diff --git a/man/R_bar.Rd b/man/R_bar.Rd deleted file mode 100644 index df0f3dc..0000000 --- a/man/R_bar.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{R_bar} -\alias{R_bar} -\title{Compute population growth rate under stabilizing selection} -\usage{ -R_bar(R0, Wbar) -} -\arguments{ -\item{R0}{basic reproductive number} - -\item{Wbar}{average fitness} -} -\description{ -Compute population growth rate under stabilizing selection -} -\details{ -Assumes density-independent -} - diff --git a/man/R_bar_ceiling.Rd b/man/R_bar_ceiling.Rd new file mode 100644 index 0000000..e82f37b --- /dev/null +++ b/man/R_bar_ceiling.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2 (4.0.2): do not edit by hand +\name{R_bar_ceiling} +\alias{R_bar_ceiling} +\title{Compute population growth rate under stabilizing selection and ceiling regulation} +\usage{ +R_bar_ceiling(R0, Wbar, N, K) +} +\arguments{ +\item{R0}{basic reproductive number} + +\item{Wbar}{average fitness} + +\item{N}{number of individuals in this generation} + +\item{K}{carrying capacity} +} +\description{ +Compute population growth rate under stabilizing selection and ceiling regulation +} +\details{ +Assumes ceiling population regulation +would be good to have separate DD function specified after +chevin and lande 2010 +} + diff --git a/man/R_bar_dd.Rd b/man/R_bar_thetalog.Rd similarity index 71% rename from man/R_bar_dd.Rd rename to man/R_bar_thetalog.Rd index 7acec1e..7dfa592 100644 --- a/man/R_bar_dd.Rd +++ b/man/R_bar_thetalog.Rd @@ -1,9 +1,9 @@ % Generated by roxygen2 (4.0.2): do not edit by hand -\name{R_bar_dd} -\alias{R_bar_dd} -\title{Compute population growth rate under stabilizing selection} +\name{R_bar_thetalog} +\alias{R_bar_thetalog} +\title{Compute population growth rate under stabilizing selection and theta-logistic regulation} \usage{ -R_bar_dd(R0, Wbar, N, K, thetaL) +R_bar_thetalog(R0, Wbar, N, K, thetaL) } \arguments{ \item{R0}{basic reproductive number} @@ -17,7 +17,7 @@ R_bar_dd(R0, Wbar, N, K, thetaL) \item{thetaL}{theta-logistic parameter for density dependence} } \description{ -Compute population growth rate under stabilizing selection +Compute population growth rate under stabilizing selection and theta-logistic regulation } \details{ Assumes theta-logistic population regulation diff --git a/man/W_bar.Rd b/man/W_bar.Rd index 9708dd1..f2efce3 100644 --- a/man/W_bar.Rd +++ b/man/W_bar.Rd @@ -3,16 +3,16 @@ \alias{W_bar} \title{Compute mean fitness under stabilizing selection} \usage{ -W_bar(zbar, theta, Oz2, gamma, LOG) +W_bar(zbar, theta, oz2, gamma, LOG) } \arguments{ \item{zbar}{mean trait prior to selection} \item{theta}{enironmental optimum} -\item{Oz2}{strength of stabilizing selection} +\item{oz2}{strength of stabilizing selection} -\item{gamma}{1/(Oz2 + Vz2), with Vz2 phenotypic variance} +\item{gamma}{1/(oz2 + Vz2), with Vz2 phenotypic variance} \item{LOG}{(=TRUE) whether to return the log of fitness} } @@ -20,6 +20,6 @@ W_bar(zbar, theta, Oz2, gamma, LOG) Compute mean fitness under stabilizing selection } \details{ -no details +IN USE } diff --git a/man/make_env.Rd b/man/make_env.Rd new file mode 100644 index 0000000..9be4a09 --- /dev/null +++ b/man/make_env.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2 (4.0.2): do not edit by hand +\name{make_env} +\alias{make_env} +\title{Compute correctly-scaled environment of development and selection} +\usage{ +make_env(T, env_args) +} +\arguments{ +\item{T}{how many gens} + +\item{env_args}{environment args} +} +\description{ +Compute correctly-scaled environment of development and selection +} +\details{ +given a pre-generated vector of random environments +} + diff --git a/man/pdTS.Rd b/man/simulate_pheno_ts.Rd similarity index 88% rename from man/pdTS.Rd rename to man/simulate_pheno_ts.Rd index 1a46c51..cd36eb0 100644 --- a/man/pdTS.Rd +++ b/man/simulate_pheno_ts.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2 (4.0.2): do not edit by hand -\name{pdTS} -\alias{pdTS} +\name{simulate_pheno_ts} +\alias{simulate_pheno_ts} \title{Compute phenotypic dynamic Time Series of trait + demographic change under stabilizing selection as function of environment (after Lande Chevin)} \usage{ -pdTS(T, X, params, env_args) +simulate_pheno_ts(T, X, params, env_args) } \arguments{ \item{T}{end time, assuming start time of 1} diff --git a/man/stab_curve.Rd b/man/stab_curve.Rd deleted file mode 100644 index 373b367..0000000 --- a/man/stab_curve.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2 (4.0.2): do not edit by hand -\name{stab_curve} -\alias{stab_curve} -\title{Mean fitness -- with constraint, need to make as function of environment...} -\usage{ -stab_curve(x, LL) -} -\arguments{ -\item{x}{the environment} - -\item{LL}{lethal limit (symmetric)} -} -\description{ -Mean fitness -- with constraint, need to make as function of environment... -} -\details{ -fitness surf infinitesimal optimum need to revisit -} - diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 558f6c5..fec2a1e 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -7,161 +7,111 @@ using namespace Rcpp; // Log_W_bar -double Log_W_bar(double zbar, double theta, double Oz2, double gamma); -RcppExport SEXP phenoecosim_Log_W_bar(SEXP zbarSEXP, SEXP thetaSEXP, SEXP Oz2SEXP, SEXP gammaSEXP) { +double Log_W_bar(double zbar, double theta, double oz2, double gamma); +RcppExport SEXP phenoecosim_Log_W_bar(SEXP zbarSEXP, SEXP thetaSEXP, SEXP oz2SEXP, SEXP gammaSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< double >::type zbar(zbarSEXP ); - Rcpp::traits::input_parameter< double >::type theta(thetaSEXP ); - Rcpp::traits::input_parameter< double >::type Oz2(Oz2SEXP ); - Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP ); - double __result = Log_W_bar(zbar, theta, Oz2, gamma); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< double >::type zbar(zbarSEXP); + Rcpp::traits::input_parameter< double >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< double >::type oz2(oz2SEXP); + Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP); + __result = Rcpp::wrap(Log_W_bar(zbar, theta, oz2, gamma)); + return __result; END_RCPP } // W_bar -double W_bar(double zbar, double theta, double Oz2, double gamma, bool LOG); -RcppExport SEXP phenoecosim_W_bar(SEXP zbarSEXP, SEXP thetaSEXP, SEXP Oz2SEXP, SEXP gammaSEXP, SEXP LOGSEXP) { +double W_bar(double zbar, double theta, double oz2, double gamma, bool LOG); +RcppExport SEXP phenoecosim_W_bar(SEXP zbarSEXP, SEXP thetaSEXP, SEXP oz2SEXP, SEXP gammaSEXP, SEXP LOGSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< double >::type zbar(zbarSEXP ); - Rcpp::traits::input_parameter< double >::type theta(thetaSEXP ); - Rcpp::traits::input_parameter< double >::type Oz2(Oz2SEXP ); - Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP ); - Rcpp::traits::input_parameter< bool >::type LOG(LOGSEXP ); - double __result = W_bar(zbar, theta, Oz2, gamma, LOG); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< double >::type zbar(zbarSEXP); + Rcpp::traits::input_parameter< double >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< double >::type oz2(oz2SEXP); + Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP); + Rcpp::traits::input_parameter< bool >::type LOG(LOGSEXP); + __result = Rcpp::wrap(W_bar(zbar, theta, oz2, gamma, LOG)); + return __result; END_RCPP } // R_bar -double R_bar(double R0, double Wbar); -RcppExport SEXP phenoecosim_R_bar(SEXP R0SEXP, SEXP WbarSEXP) { +double R_bar(double R0, double Wbar, double N); +RcppExport SEXP phenoecosim_R_bar(SEXP R0SEXP, SEXP WbarSEXP, SEXP NSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< double >::type R0(R0SEXP ); - Rcpp::traits::input_parameter< double >::type Wbar(WbarSEXP ); - double __result = R_bar(R0, Wbar); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< double >::type R0(R0SEXP); + Rcpp::traits::input_parameter< double >::type Wbar(WbarSEXP); + Rcpp::traits::input_parameter< double >::type N(NSEXP); + __result = Rcpp::wrap(R_bar(R0, Wbar, N)); + return __result; END_RCPP } -// Log_R_bar -double Log_R_bar(double R0, double logWbar); -RcppExport SEXP phenoecosim_Log_R_bar(SEXP R0SEXP, SEXP logWbarSEXP) { +// R_bar_ceiling +double R_bar_ceiling(double R0, double Wbar, double N, double K); +RcppExport SEXP phenoecosim_R_bar_ceiling(SEXP R0SEXP, SEXP WbarSEXP, SEXP NSEXP, SEXP KSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< double >::type R0(R0SEXP ); - Rcpp::traits::input_parameter< double >::type logWbar(logWbarSEXP ); - double __result = Log_R_bar(R0, logWbar); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< double >::type R0(R0SEXP); + Rcpp::traits::input_parameter< double >::type Wbar(WbarSEXP); + Rcpp::traits::input_parameter< double >::type N(NSEXP); + Rcpp::traits::input_parameter< double >::type K(KSEXP); + __result = Rcpp::wrap(R_bar_ceiling(R0, Wbar, N, K)); + return __result; END_RCPP } -// R_bar_dd -double R_bar_dd(double R0, double Wbar, double N, double K, double thetaL); -RcppExport SEXP phenoecosim_R_bar_dd(SEXP R0SEXP, SEXP WbarSEXP, SEXP NSEXP, SEXP KSEXP, SEXP thetaLSEXP) { +// R_bar_thetalog +double R_bar_thetalog(double R0, double Wbar, double N, double K, double thetaL); +RcppExport SEXP phenoecosim_R_bar_thetalog(SEXP R0SEXP, SEXP WbarSEXP, SEXP NSEXP, SEXP KSEXP, SEXP thetaLSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< double >::type R0(R0SEXP ); - Rcpp::traits::input_parameter< double >::type Wbar(WbarSEXP ); - Rcpp::traits::input_parameter< double >::type N(NSEXP ); - Rcpp::traits::input_parameter< double >::type K(KSEXP ); - Rcpp::traits::input_parameter< double >::type thetaL(thetaLSEXP ); - double __result = R_bar_dd(R0, Wbar, N, K, thetaL); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; -END_RCPP -} -// Beta -arma::vec Beta(double gamma, double A, double B, double a, double b, double e_t, double e_plast); -RcppExport SEXP phenoecosim_Beta(SEXP gammaSEXP, SEXP ASEXP, SEXP BSEXP, SEXP aSEXP, SEXP bSEXP, SEXP e_tSEXP, SEXP e_plastSEXP) { -BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< double >::type gamma(gammaSEXP ); - Rcpp::traits::input_parameter< double >::type A(ASEXP ); - Rcpp::traits::input_parameter< double >::type B(BSEXP ); - Rcpp::traits::input_parameter< double >::type a(aSEXP ); - Rcpp::traits::input_parameter< double >::type b(bSEXP ); - Rcpp::traits::input_parameter< double >::type e_t(e_tSEXP ); - Rcpp::traits::input_parameter< double >::type e_plast(e_plastSEXP ); - arma::vec __result = Beta(gamma, A, B, a, b, e_t, e_plast); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< double >::type R0(R0SEXP); + Rcpp::traits::input_parameter< double >::type Wbar(WbarSEXP); + Rcpp::traits::input_parameter< double >::type N(NSEXP); + Rcpp::traits::input_parameter< double >::type K(KSEXP); + Rcpp::traits::input_parameter< double >::type thetaL(thetaLSEXP); + __result = Rcpp::wrap(R_bar_thetalog(R0, Wbar, N, K, thetaL)); + return __result; END_RCPP } // Va double Va(arma::vec env, arma::mat GG); RcppExport SEXP phenoecosim_Va(SEXP envSEXP, SEXP GGSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< arma::vec >::type env(envSEXP ); - Rcpp::traits::input_parameter< arma::mat >::type GG(GGSEXP ); - double __result = Va(env, GG); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< arma::vec >::type env(envSEXP); + Rcpp::traits::input_parameter< arma::mat >::type GG(GGSEXP); + __result = Rcpp::wrap(Va(env, GG)); + return __result; END_RCPP } -// Env_shift_cpp -arma::vec Env_shift_cpp(int t, List env_args); -RcppExport SEXP phenoecosim_Env_shift_cpp(SEXP tSEXP, SEXP env_argsSEXP) { +// make_env +arma::mat make_env(int T, List env_args); +RcppExport SEXP phenoecosim_make_env(SEXP TSEXP, SEXP env_argsSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< int >::type t(tSEXP ); - Rcpp::traits::input_parameter< List >::type env_args(env_argsSEXP ); - arma::vec __result = Env_shift_cpp(t, env_args); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< int >::type T(TSEXP); + Rcpp::traits::input_parameter< List >::type env_args(env_argsSEXP); + __result = Rcpp::wrap(make_env(T, env_args)); + return __result; END_RCPP } -// pdTS -arma::mat pdTS(int T, arma::rowvec X, List params, List env_args); -RcppExport SEXP phenoecosim_pdTS(SEXP TSEXP, SEXP XSEXP, SEXP paramsSEXP, SEXP env_argsSEXP) { +// simulate_pheno_ts +arma::mat simulate_pheno_ts(int T, arma::rowvec X, List params, List env_args); +RcppExport SEXP phenoecosim_simulate_pheno_ts(SEXP TSEXP, SEXP XSEXP, SEXP paramsSEXP, SEXP env_argsSEXP) { BEGIN_RCPP - SEXP __sexp_result; - { - Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< int >::type T(TSEXP ); - Rcpp::traits::input_parameter< arma::rowvec >::type X(XSEXP ); - Rcpp::traits::input_parameter< List >::type params(paramsSEXP ); - Rcpp::traits::input_parameter< List >::type env_args(env_argsSEXP ); - arma::mat __result = pdTS(T, X, params, env_args); - PROTECT(__sexp_result = Rcpp::wrap(__result)); - } - UNPROTECT(1); - return __sexp_result; + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< int >::type T(TSEXP); + Rcpp::traits::input_parameter< arma::rowvec >::type X(XSEXP); + Rcpp::traits::input_parameter< List >::type params(paramsSEXP); + Rcpp::traits::input_parameter< List >::type env_args(env_argsSEXP); + __result = Rcpp::wrap(simulate_pheno_ts(T, X, params, env_args)); + return __result; END_RCPP } diff --git a/src/traitfitpop.cpp b/src/traitfitpop.cpp index 564e2dd..1aa0e0c 100644 --- a/src/traitfitpop.cpp +++ b/src/traitfitpop.cpp @@ -1,93 +1,79 @@ -#include "RcppArmadillo.h" // after version 4.3 +#include // after version 4.3 +#include // to suppress linter warnings // [[Rcpp::depends(RcppArmadillo)]] using namespace Rcpp; // done by sourceCpp? - -// ** below at least compile when used in -// cppFunction(code="//code here", depends="RcppArmadillo") - //' Compute mean fitness under stabilizing selection //' @param zbar mean trait prior to selection //' @param theta enironmental optimum -//' @param Oz2 strength of stabilizing selection -//' @param gamma 1/(Oz2 + Vz2), with Vz2 phenotypic variance -//' @details no details +//' @param oz2 strength of stabilizing selection +//' @param gamma 1/(oz2 + Vz2), with Vz2 phenotypic variance +//' @details IN USE //' @export // [[Rcpp::export]] -double Log_W_bar (double zbar, double theta, double Oz2, double gamma){ - return 0.5 * log(gamma * Oz2) - gamma /2 * pow(zbar - theta, 2); //## compare with eqn 2c lande 2009 +double Log_W_bar (double zbar, double theta, double oz2, double gamma){ + return 0.5 * log(gamma * oz2) - gamma / 2 * pow(zbar - theta, 2); //## compare with eqn 2c lande 2009 } //' Compute mean fitness under stabilizing selection //' @inheritParams Log_W_bar -//' @param LOG (=TRUE) whether to return the log of fitness -//' @details no details +//' @param LOG (=TRUE) whether to return the log of fitness +//' @details IN USE //' @export // [[Rcpp::export]] -double W_bar (double zbar, double theta, double Oz2, double gamma, bool LOG){ - if(LOG){ - return Log_W_bar(zbar, theta, Oz2, gamma); - } - else { - return sqrt(gamma * Oz2) * exp(- gamma /2 * pow(zbar - theta,2)); - } +double W_bar (double zbar, double theta, double oz2, double gamma, bool LOG){ + if(LOG){ + return Log_W_bar(zbar, theta, oz2, gamma); + } + else { + return sqrt(gamma * oz2) * exp(- gamma / 2 * pow(zbar - theta, 2)); + } } -//' Compute population growth rate under stabilizing selection +//' Compute population growth rate under stabilizing selection and no regulation //' @param R0 basic reproductive number //' @param Wbar average fitness -//' @details Assumes density-independent +//' @param N number of individuals in this generation +//' @param K carrying capacity +//' @details Assumes ceiling population regulation +//' would be good to have separate DD function specified after +//' chevin and lande 2010 //' @export // [[Rcpp::export]] -double R_bar (double R0, double Wbar){ - return R0 * Wbar; - } +double R_bar(double R0, double Wbar, double N){ + return R0 * Wbar * N; +} -//' Compute LOG population growth rate under stabilizing selection -//' @inheritParams R_bar -//' @param logWbar log average fitness -//' @details Assumes theta-logistic population regulation +//' Compute population growth rate under stabilizing selection and ceiling regulation +//' @param R0 basic reproductive number +//' @param Wbar average fitness +//' @param N number of individuals in this generation +//' @param K carrying capacity +//' @details Assumes ceiling population regulation //' would be good to have separate DD function specified after //' chevin and lande 2010 //' @export // [[Rcpp::export]] -double Log_R_bar (double R0, double logWbar){ - return log(R0) + logWbar; - } +double R_bar_ceiling(double R0, double Wbar, double N, double K){ + N = R0 * Wbar * N; + if(N > K) { + return K; + } else { + return N; + } +} -//' Compute population growth rate under stabilizing selection -//' @inheritParams R_bar -//' @param N number of individuals in this generation -//' @param K carrying capacity +//' Compute population growth rate under stabilizing selection and theta-logistic regulation +//' @inheritParams R_bar_ceiling //' @param thetaL theta-logistic parameter for density dependence //' @details Assumes theta-logistic population regulation //' would be good to have separate DD function specified after //' chevin and lande 2010 //' @export // [[Rcpp::export]] -double R_bar_dd(double R0, double Wbar, double N, double K, double thetaL){ - return pow(R0, (1 -pow(N/K, thetaL))) * Wbar; - } - -//' Selection on plasticity as as function of environment assuming stabilizing selection -//' @param gamma = 1/(Oz2 + Vz2), with Vz2 phenotypic variance -//' @param A the environmental optimum RN int -//' @param B the environmental optimum RN slope -//' @param a current value of RN int -//' @param b current value of RN slope -//' @param e_t environment now -//' @param e_plast env that cues plasticity -//' @export -// [[Rcpp::export]] -arma::vec Beta(double gamma, double A, double B, double a, double b, double e_t, double e_plast){ - //## eqn 3b - double beta1; - beta1 = -gamma * (a - A + b * e_plast - B * e_t); - arma::vec beta(2); - beta(0) = beta1; - beta(1) = beta1 * e_plast; - return beta; +double R_bar_thetalog(double R0, double Wbar, double N, double K, double thetaL){ + return pow(R0, (1 - pow(N / K, thetaL))) * Wbar; } //' Va additive genetic variance in the phenotype as a function of the environment @@ -95,9 +81,9 @@ arma::vec Beta(double gamma, double A, double B, double a, double b, double e_t, //' @param GG additive genetic variance-covariance matrix //' @export // [[Rcpp::export]] -double Va (arma::vec env, arma::mat GG){ - // env should be column vector - return as_scalar(trans(env) * GG * env); +double Va(arma::vec env, arma::mat GG){ + // env should be column vector + return as_scalar(trans(env) * GG * env); } //' Compute a white noise environment with a shift @@ -108,141 +94,143 @@ double Va (arma::vec env, arma::mat GG){ //' sd noise in env //' sdc noise in cue //' @details nada -//' @export -// [[Rcpp::export]] arma::vec Env_shift_cpp(int t, List env_args) { - // returns (selection env, cue env) - double t_jump = as(env_args["t.jump"]); - arma::mat sigma = as(env_args["sigma"]); - //double rho = as(env_args["rho"]); - //double k = as(env_args["k"]); - arma::vec out(2); // maybe make a row vec to avoid trans() in output below? - - if( t < t_jump){ - out.zeros(); - } else { - out.fill(Rcpp::as(env_args["delta"])); - } - // Below equivalent to R code: - // > out <- mvrnorm(1, mu = c(0,0) , Sigma = matrix(c(sd, sd*sdc*rho, sd*sdc*rho, sdc), nrow=2) - // generate mvrandom w/ chol, - // [todo] - break into sep mvrnormarma see http://gallery.rcpp.org/articles/simulate-multivariate-normal - int ncols = sigma.n_cols; - arma::mat Z = arma::randn(1, ncols); - return out + arma::trans(Z * arma::chol(sigma)); // wld use arma::repmat(out, 1, n) for n>1 - // for red noise do somehting like k * env + sqrt(1 - k^2) * rnorm (n=1, mean=0, sd=sd)) but hard - // to figure this out with the cue / correlation thing - + // returns (selection env, cue env) + double t_jump = as(env_args["t.jump"]); + arma::mat sigma = as(env_args["sigma"]); + //double rho = as(env_args["rho"]); + //double k = as(env_args["k"]); + arma::vec out(2); // maybe make a row vec to avoid trans() in output below? + + if( t < t_jump){ + out.zeros(); + } else { + out.fill(Rcpp::as(env_args["delta"])); + } + // Below equivalent to R code: + // > out <- mvrnorm(1, mu = c(0,0) , Sigma = matrix(c(sd, sd*sdc*rho, sd*sdc*rho, sdc), nrow=2) + // generate mvrandom w/ chol, + // [todo] - break into sep mvrnormarma see http://gallery.rcpp.org/articles/simulate-multivariate-normal + int ncols = sigma.n_cols; + arma::mat Z = arma::randn(1, ncols); + return out + arma::trans(Z * arma::chol(sigma)); // wld use arma::repmat(out, 1, n) for n>1 + // for red noise do somehting like k * env + sqrt(1 - k^2) * rnorm (n=1, mean=0, sd=sd)) but hard + // to figure this out with the cue / correlation thing } +//' Compute correctly-scaled environment of development and selection +//' @param env_args environment args +//' @param T how many gens +//' @details given a pre-generated vector of random environments +//' @export +// [[Rcpp::export]] +arma::mat make_env(int T, List env_args) { + double sigma_xi = as(env_args["sigma_xi"]); + double rho_tau = as(env_args["rho_tau"]); + double fractgen = as(env_args["fractgen"]); + double delta = as(env_args["delta"]); + + //temp vars + double rho2 = pow(rho_tau, 2); + double xi_dev; + double xi_sel = 0; + arma::vec all_env = arma::randn((T + 1) * fractgen); // random normal for every timestep + int env_ctr = 0; + all_env = all_env * sigma_xi; // scale noise + + arma::mat xi_all(2, T + 1); + for(int t=0; t <= T; t++) { + xi_dev = xi_sel; + + for(int i=1; i < fractgen; i++) { + xi_dev = rho_tau * xi_dev + sqrt(1 - rho2) * all_env(env_ctr); + env_ctr = env_ctr + 1; + } + xi_sel = rho_tau * xi_dev + sqrt(1 - rho2) * all_env(env_ctr); + + xi_all(0, t) = xi_dev + delta; + xi_all(1, t) = xi_sel + delta; + } + return xi_all; +} -//' Compute phenotypic dynamic Time Series of trait + demographic change under stabilizing selection as +//' Compute phenotypic dynamic Time Series of trait + demographic change under stabilizing selection as //' function of environment (after Lande Chevin) //' @param T end time, assuming start time of 1 //' @param X parameters (z, a, b, wbar, logN, theta) -//' @param params a list with (gamma_sh, omegaz, A, B, R0, Va, Vb, delta, sigma_xi, rho_tau, fractgen) +//' @param params a list with (gamma_sh, omegaz, A, B, R0, var_a, Vb, delta, sigma_xi, rho_tau, fractgen) //' @param env_args extra args for env.fn //' @details NB - for now assume Tchange = 0 and demography after CL 2010 +//' @value a long matrix //' @export // [[Rcpp::export]] - -arma::mat pdTS (int T, arma::rowvec X, List params, - List env_args) { - // TODO - make above functions consistent with this approach - // NB - for now assume Tchange = 0 - // state vars in X (R indexing) - // zbar <- X[1] - // abar <- X[2] // elevation - // bbar <- X[3] // slope - // Wbar <- X[4] - // Npop <- X[5] - // Theta <- X[6] - - //tmp vars for bio - double zbart; - double abart; - double bbart; - double wbart; - double thetat; - double betat; - - - double gamma_sh = as(params["gamma_sh"]); - double omegaz = as(params["omegaz"]); - double A = as(params["A"]); - double B = as(params["B"]); - double R0 = as(params["R0"]); - double Va = as(params["Va"]); - double Vb = as(params["Vb"]); - - double delta = as(env_args["delta"]); - double sigma_xi = as(env_args["sigma_xi"]); - double rho_tau = as(env_args["rho_tau"]); - double fractgen = as(env_args["fractgen"]); - - double gamma = 0; - gamma = gamma_sh; - - int len = X.size(); - arma::mat out(len, T + 1); - out.col(0) = trans(X); - - // tmp vars for env - double xi_dev_temp; - double xi_sel_temp; - double eps_sel; - double eps_dev; - - xi_sel_temp = 0; - - arma::vec all_env = arma::randn((T + 1) * fractgen); - int env_ctr = 0; - - all_env = all_env * sigma_xi; // scale noise - - for(int t=0; t <= T; t++) { - - if (t == 0) - xi_dev_temp = 0; - else - xi_dev_temp = xi_sel_temp; - - for(int i=1; i < fractgen; i++) { - xi_dev_temp = rho_tau * xi_dev_temp + sqrt(1 - pow(rho_tau, 2)) * all_env(env_ctr); // make drawenv - env_ctr = env_ctr + 1; - } - - - xi_sel_temp = rho_tau * xi_dev_temp + sqrt(1 - pow(rho_tau, 2)) * all_env(env_ctr); // make drawenv - env_ctr = env_ctr + 1; - - //## computing the corresponding mean/variance genetic values, and optimal environmen - eps_dev = delta + xi_dev_temp; // ##xi_temp = xi_dev - eps_sel = delta + xi_sel_temp; - - abart = out(1, t); - bbart = out(2, t); - zbart = abart + bbart * eps_dev; // zbar = abar + bbar * eps_dev - out(0, t) = zbart; - // varz <- Va + 2 * Vab * eps_dev + Vb * pow(eps_dev, 2) + Ve; - thetat = A + B * eps_sel; //Theta - out(5, t) = thetat; - // ## evolutionary change with and population growth - wbart = R0 * sqrt(gamma * pow(omegaz, 2)) * exp(-gamma / 2 * pow(zbart - thetat, 2)); - out(3, t) = wbart; - // Wbar= f( zbar - theta) - - //## in R version increment here to update the vars with init conditions already filled in - // t = t + 1; - - if(t < T) { - betat = gamma * (zbart - thetat); - // abar = abar - gamma(zbar - theta) * Va - out(1, t + 1) = abart - betat * Va; - // bbar = bbar - gamma (zbar - theta) eps_dev Vb - out(2, t + 1) = bbart - betat * eps_dev * Vb ; - out(4, t + 1) = out(4, t) * wbart; // N = N * wbar - } - } - return out; +arma::mat simulate_pheno_ts(int T, arma::rowvec X, List params, + List env_args) { + // TODO - make above functions consistent with this approach + // NB - for now assume Tchange = 0 + // state vars in X (R indexing) + // zbar <- X[1] + // abar <- X[2] // elevation + // bbar <- X[3] // slope + // Wbar <- X[4] // Npop <- X[5] // Theta <- X[6] + + //tmp vars for bio + double zbart; + double abart; + double bbart; + double wbart; + double thetat; + double betat; + + double gamma_sh = as(params["gamma_sh"]); + double omegaz = as(params["omegaz"]); + double A = as(params["A"]); + double B = as(params["B"]); + double R0 = as(params["R0"]); + double var_a = as(params["var_a"]); + double Vb = as(params["Vb"]); + + double gamma = gamma_sh; + double oz2 = pow(omegaz, 2); + + int len = X.size(); + arma::mat out(T + 1, len); + out.row(0) = X; + + // tmp vars for env + double eps_sel; + double eps_dev; + //Rfprintf("Making env .......................") + arma::mat eps_all = make_env(T, env_args); + + for(int t=0; t <= T; t++) { + // get environments of dev, sel for this gen + eps_dev = eps_all(0, t); + eps_sel = eps_all(1, t); + + // compute current-generation values for RN params, environmental optimum, and fitness + abart = out( t, 1); + bbart = out( t, 2); + zbart = abart + bbart * eps_dev; + out( t, 0) = zbart; + thetat = A + B * eps_sel; + out( t, 5) = thetat; + + // evolutionary change with and population growth + wbart = W_bar(zbart, thetat, oz2, gamma, FALSE); + out( t, 3) = R0 * wbart; + + if(t < T) { + // compute next gen values for RN params, population size + // \delta (a, b) = \beta (var_a, var_b) = gamma ( z - theta) (var_a, var_b) + // given def betat, bbart line below defines Vb(env) = Vb * env ^2 + // + // N = N * wbar + betat = gamma * (zbart - thetat); + out( t + 1, 1) = abart - betat * var_a; + out( t + 1, 2) = bbart - betat * eps_dev * Vb ; + out( t + 1, 4) = R_bar(R0, wbart, out( t, 4)); + } + } + return out; } diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..410dd43 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(phenoecosim) + +test_check("phenoecosim") diff --git a/tests/testthat/test-demography.R b/tests/testthat/test-demography.R new file mode 100644 index 0000000..16a5257 --- /dev/null +++ b/tests/testthat/test-demography.R @@ -0,0 +1,6 @@ +library(phenoecosim) +context("Demography and related functions testing") + +# function for Rbar_ceiling +# function for rbar thetalog +# function for Wbar, Wbar log diff --git a/tests/testthat/test-simulate-env.R b/tests/testthat/test-simulate-env.R new file mode 100644 index 0000000..a4d3ee2 --- /dev/null +++ b/tests/testthat/test-simulate-env.R @@ -0,0 +1,25 @@ +library(phenoecosim) + +ERROR_TOL <- 0.15 +for(rho in 1:8 * 0.1) { + context(paste("Environmental simulation testing for rho =", rho)) + max_T = 1000 + start_T = 1 + for(delta in 1:8) + for(fractgen in 5:10) + for(sigma_xi in c(0.5, 0.75, 1, 3)) { + env.list <- list(delta=delta, sigma_xi=sigma_xi, rho_tau=rho, fractgen=fractgen) + env_dev_sel = make_env(T=max_T, env_args=env.list) + par_names <- paste0(paste(names(env.list), env.list, sep='='), collapse=',') + test_that(paste("Environmental simulation produces expected correlations, mean, sd for", par_names), { + envcor <- cor(env_dev_sel[1, start_T:max_T], env_dev_sel[2, start_T:max_T]) + envmeans <- rowMeans(env_dev_sel) + envsd <- apply(env_dev_sel, 1, sd) + expect_equal(envcor, rho, tolerance=ERROR_TOL, scale=sigma_xi) + expect_equal(envmeans, rep(delta, 2), tolerance=ERROR_TOL, scale=sigma_xi) + expect_equal(envsd, rep(sigma_xi, 2), tolerance=ERROR_TOL, scale=sigma_xi) + }) + test_that(paste("Environmental simulation produces mean", delta, "for", par_names), { + }) + } +} diff --git a/tests/testthat/test-simulate-ts.R b/tests/testthat/test-simulate-ts.R new file mode 100644 index 0000000..cb24eb6 --- /dev/null +++ b/tests/testthat/test-simulate-ts.R @@ -0,0 +1,55 @@ +library(phenoecosim) +context("Phenotypic timeseries regression test") + +## variance load +Vz_ <- function(siga2, sigb2, delta, sige2) + siga2 + sigb2 * delta ^ 2 + sige2 + +## unvarying parameters +A <- 0 +B <- 2 #environmental cline in the optimum phenotype +fractgen <- 5 +tau <- 1 / fractgen +Ve <- 0.5 +Tchange <- 0 #time for a change in env regime +Npop0 <- 1e4 +rho1 <- 0.5 +sigma_xi <- 1 +omega_z <- sqrt(20) +omega_Wmax <- 1.10 #maximum fitness across environments + +for( delta in c(1.5, 5, 8)) + test_that(paste("Simulate_timeseries produces right object across parameters for delta=",delta), { + ## varying params + for( Tlim in c(250, 500, 1000)) + for( Va in c(0.05, 0.1)) + for( Vb in c(0.025, 0.05)) + for( omegaz in sqrt(c(10, 25))) + for( rho0 in c(0.3, 0.7)) + for( alpha in c(0.3, 0.7)) + { + Vz <- Vz_(Va, Vb, delta, Ve) + gamma_sh <- 1 / (omegaz ^ 2 + Vz) + rmax <- log(omega_Wmax) - 1 / 2 * log(1 + Vz / (omegaz ^ 2)) + abar0 <- A; + bbar0 <- alpha * B + + param.list <- list(Vz=Vz, gamma_sh=gamma_sh, rmax=rmax, omegaz=omegaz, + A=A, B=B, R0=omega_Wmax, var_a=Va, Vb=Vb, Ve=Ve) + + env.list <- list(delta=delta, sigma_xi=sigma_xi, rho_tau=rho1, fractgen=fractgen) + X0 <- c(zbar=NA, abar=abar0, bbar=bbar0, Wbar=NA, Npop=Npop0, theta=NA) + + set.seed(111) + tmp <- simulate_pheno_ts(Tlim, X0, param.list, env.list) + target_parms <- paste(Tlim, Va, Vb, omegaz, omega_Wmax, rho0, rho1, + alpha, delta, sigma_xi, sep='-') + target_filename <- paste0('simulate', target_parms, '.rds') + expect_equal_to_reference(tmp, target_filename, scale=1, tolerance=0.1, + label=paste('Pars:', + paste(paste(names(env.list), env.list, sep='='), + paste(names(param.list), param.list, sep='='), sep=',', + collapse=' ') + )) + } + })