Refactor initialisation spline

R/build_schedule.R

@@ -84,6 +84,72 @@ build_schedule <- function(p, env = make_environment(parameters = p),
   return(list(parameters = p, state = state, n_steps = i, complete = complete))
 }
 
+##' @rdname initialise_scm
+##' @param i Index to extract from \code{x}
+##' @param state List object with initial cohorts for each species
+##' @param size_idx Index of the strategy size characteristic
+##' @export
+init_spline <- function(state, size_idx = 1) {
+  state_to_spline <- function(x) {
+    vars = rownames(x)
+
+    s <- x[vars[size_idx], ]
+
+    if(length(s) < 2)
+      stop("At least two cohorts needed to define a spline")
+
+    m <- c(min(s), max(s))
+
+    # identity for size, zero for zeros, log for everything else
+    f_ <- lapply(vars,
+                 function(v) {
+                   y = x[v, ]
+                   if(v == vars[size_idx])
+                     splinefun(s, y)
+                   else if(v == "log_density")
+                     splinefun_log(s, y)
+                   else if(length(y[y>0]) == 0)
+                     splinefun(s, rep(0, length(s)))
+                   else 
+                     splinefun_loglog(s[y>0], y[y>0])
+                 })
+
+    # set bounds
+    f_ <- lapply(f_, clamp_domain, m)
+
+    # rename to: size_variable
+    spline_names <- paste(vars[size_idx], vars, sep = "_")
+    names(f_) <- spline_names
+
+    # useful attributes
+    attr(f_, 'size_idx') <- size_idx
+    attr(f_, 'domain') <- m
+
+    return(f_)
+  }
+
+  splines <- lapply(state, state_to_spline)
+}
+
+##' @rdname initialise_scm
+##' @param splines Output of \code{\link{init_spline}}
+##' @param sizes Size of initial cohorts
+##' @param n Number of initial cohorts, if sizes is missing (default = 10)
+##' @export
+partition_spline <- function(splines, sizes=NULL, n=10) {
+  if(is.null(sizes)) {
+    m <- attr(splines, 'domain')
+    sizes = seq(m[1], m[2], len = n)
+  }
+
+  state <- t(sapply(splines, function(f) f(rev(sizes))))
+
+  # regex magic to remove first word and underscore
+  rownames(state) <- gsub("^[^_]+(?=_)_", "", names(splines), perl=T)
+
+  return(state)
+}
+
 split_times <- function(times, i) {
   ## Upwind splitting scheme only, which means that we will never
   ## split the last interval [assuming OK for now].  Inefficiently

R/scm_support.R

@@ -158,6 +158,7 @@ run_scm_collect <- function(p, env = make_environment(parameters = p),
 }
 
 ##' Functions for reconstructing a Patch from an SCM
+##' @rdname make_patch
 ##' @title Reconstruct a patch
 ##' @param state State object created by \code{scm_state}
 ##' @param p Parameters object
@@ -174,9 +175,11 @@ make_patch <- function(state, p, env = make_environment(parameters = p),
 }
 
 ##' Functions for reconstructing a Patch from an SCM
-##' @title Reconstruct a patch
+##' @title Construct an SCM
 ##' @param p Parameters object
-##' @param state An optional State object matching the strategies in \code{p}
+##' @param env Environment object (defaults to FF16_Environment)
+##' @param ctrl Control object
+##' @param state List object with initial cohorts for each species
 ##' @export
 make_scm <- function(p, env, ctrl, state=NULL) {
   types <- extract_RcppR6_template_types(p, "Parameters")
@@ -209,74 +212,13 @@ make_scm <- function(p, env, ctrl, state=NULL) {
     # add as many new cohorts as required to fit `state` object
     scm$set_state(min(start_time), unlist(state),
                   n = mapply(`+`, new_cohorts, initial_cohorts))
-  }
-
-  return(scm)
-}
-
-##' @rdname make_patch
-##' @param i Index to extract from \code{x}
-##' @param state List object with initial cohorts for each species
-##' @param size_idx Index of the strategy size characteristic
-##' @export
-init_spline <- function(i, state, size_idx = 1) {
-  state_to_spline <- function(x) {
-    vars = rownames(x)
 
-    s <- x[vars[size_idx], ]
-    m <- c(min(s), max(s))
-
-    # identity for size, zero for zeros, log for everything else
-    f_ <- lapply(vars,
-                 function(v) {
-                   y = x[v, ]
-                   if(v == vars[size_idx])
-                     splinefun(s, y)
-                   else if(v == "log_density")
-                     splinefun_log(s, y)
-                   else if(length(y[y>0]) == 0)
-                     splinefun(s, 0)
-                   else 
-                     splinefun_loglog(s[y>0], y[y>0])
-                   })
-
-    # set bounds
-    f_ <- lapply(f_, clamp_domain, m)
-
-    # rename to: size_variable
-    spline_names <- paste(vars[size_idx], vars, sep = "_")
-    names(f_) <- spline_names
-
-    # useful attributes
-    attr(f_, 'size_idx') <- size_idx
-    attr(f_, 'domain') <- m
-
-    return(f_)
-  }
-
-  splines <- lapply(state, state_to_spline)
-}
-
-##' @rdname make_patch
-##' @param splines Output of \code{\link{scm_spline}}
-##' @param sizes Size of initial cohorts
-##' @param n Number of initial cohorts, if sizes is missing (default = 10)
-##' @export
-partition_spline <- function(splines, sizes=NULL, n=10) {
-  if(is.null(sizes)) {
-    m <- attr(splines, 'domain')
-    sizes = seq(m[1], m[2], len = n)
   }
 
-  state <- t(sapply(splines, function(f) f(rev(sizes))))
-
-  # regex magic to remove first word and underscore
-  rownames(state) <- gsub("^[^_]+(?=_)_", "", names(splines), perl=T)
-
-  return(state)
+  return(scm)
 }
 
-##' @rdname make_patch
+##' @rdname initialise_scm
 ##' @param i Index to extract from \code{x}
 ##' @param x Result of running \code{\link{run_scm_collect}}
 ##' @export

tests/testthat/test-initial-size-distribution.R

@@ -114,44 +114,46 @@ test_that("Build schedule works", {
   ctrl <- scm_base_control()
   p0 <- scm_base_parameters("FF16")
   p0$patch_type = 'fixed'
-  ctrl$schedule_nsteps = 5
+  ctrl$schedule_nsteps = 20
 
   p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16_hyperpar,FALSE)
-  p1$cohort_schedule_times[[1]] <- seq(0.01, p1$max_patch_lifetime, length = 10)
+  p1$cohort_schedule_times[[1]] <- seq(0, p1$max_patch_lifetime, length = 10)
   p1$birth_rate <- 20
 
-  init <- matrix(c(10, 0, 0, 0, 0, 0, -2))
+  # manual construction (see also: scm_state)
+  init <- matrix(c(10, 0, 0, 0, 0, 0, 1,
+                   0.3920458, 0, 0, 0, 0, 0, -4), ncol = 2)
   rownames(init) <- c("height", "mortality", "fecundity", "area_heartwood",
                       "mass_heartwood", "offspring_produced_survival_weighted",
                       "log_density")
 
   # interpolate initial size distribution
-  splines <- init_spline(1, list(init), size_idx = 1)
+  splines <- init_spline(list(init), size_idx = 1)
 
-  res <- build_schedule(p1, env, ctrl, splines, n_init = 2)
+  res <- build_schedule(p1, env, ctrl, splines, n_init = 20)
   expect_false(res$complete)
   expect_equal(res$n_steps, 5)
-  expect_equal(attr(res$parameters, "net_reproduction_ratios"), 4997.393, tolerance = 1e-6)
-
+  expect_equal(attr(res$parameters, "net_reproduction_ratios"), 
+               4997.393, tolerance = 1e-6)
 
   # Plot for fun
-  # par(mfrow=c(2, 1))
-  # x <- run_scm_collect(p1, env, ctrl, state = list(init))
-  # 
-  # t2 <- x$time
-  # h1 <- x$species[[1]]["height", , ]
-  # 
-  # matplot(t2, h1, lty=1, type="l", col = "black",
-  #         las=1, xlab="Time (years)", ylab="Height (m)")
-  # title("Unresolved schedule")
-  # 
-  # x <- run_scm_collect(res$parameters, env, ctrl, res$state)
-  # 
-  # t2 <- x$time
-  # h1 <- x$species[[1]]["height", , ]
-  # 
-  # # Species 1
-  # matplot(t2, h1, lty=1, type="l", col = "black",
-  #         las=1, xlab="Time (years)", ylab="Height (m)")
-  # title("Resolved schedule")
-})
+  par(mfrow=c(2, 1))
+  x <- run_scm_collect(p1, env, ctrl, list(init))
+
+  t2 <- x$time
+  h1 <- x$species[[1]]["height", , ]
+
+  matplot(t2, h1, lty=1, type="l", col = "black",
+          las=1, xlab="Time (years)", ylab="Height (m)")
+  title("Unresolved schedule")
+
+  x <- run_scm_collect(res$parameters, env, ctrl, res$state)
+
+  t2 <- x$time
+  h1 <- x$species[[1]]["height", , ]
+
+  # Species 1
+  matplot(t2, h1, lty=1, type="l", col = "black",
+          las=1, xlab="Time (years)", ylab="Height (m)")
+  title("Resolved schedule")
+})