Resolve initial size distribution as well as introduction schedule

R/build_schedule.R

@@ -13,7 +13,7 @@
 ##' \code{birth_rate}.
 ##' @author Rich FitzJohn
 ##' @export
-build_schedule <- function(p) {
+build_schedule <- function(p, state=NULL) {
   p <- validate(p)
 
   n_spp <- length(p$strategies)
@@ -24,7 +24,7 @@ build_schedule <- function(p) {
   eps <- control$schedule_eps
 
   for (i in seq_len(control$schedule_nsteps)) {
-    res <- run_scm_error(p)
+    res <- run_scm_error(p, state)
     net_reproduction_ratios <- res[["net_reproduction_ratios"]]
     split <- lapply(res$err$total, function(x) x > eps)
 
@@ -33,12 +33,24 @@ build_schedule <- function(p) {
     }
 
     ## Prepare for the next iteration:
-    times <- p$cohort_schedule_times
+    times <- res$schedule
+
+    # split cohorts
     for (idx in seq_len(n_spp)) {
+
+      # by introduction time
       times[[idx]] <- split_times(times[[idx]], split[[idx]])
+
+      # or by initial size
+      if(!is.null(state)) {
+        state$species[[idx]] <- split_state(state$species[[idx]], 
+                                            times[[idx]], split[[idx]],
+                                            res$min_heights[[idx]])
+      }
     }
+
+    # set schedule for next patch
     p$cohort_schedule_times <- times
-
     msg <- sprintf("%d: Splitting {%s} times (%s)",
                    i,
                    paste(sapply(split, sum),    collapse=","),
@@ -50,7 +62,7 @@ build_schedule <- function(p) {
   ## Useful to record the last offspring produced:
   attr(p, "net_reproduction_ratios") <- net_reproduction_ratios
 
-  p
+  return(list(parameters = p, state = state))
 }
 
 split_times <- function(times, i) {
@@ -61,5 +73,40 @@ split_times <- function(times, i) {
   ## point twice; one from upstream and one from downstream.
   dt <- diff(times)
   i <- which(i)
+
+  # can't split cohorts introduced in the same time step
+  i <- i[dt[i] > 0]
+
   sort(c(times, times[i] - dt[i-1]/2))
 }
+
+split_state <- function(state, times, i, min_height) {
+  ## oversimplified splitting scheme, introduce a new cohort 
+  # half a step between two existing cohorts, even thought this leaves density
+  # uncorrected and results in a different initial size distribution
+
+  # can only split cohorts introduced at t0
+  i <- which(i)
+
+  # which is tricky when there's only one - we want to bring the 
+  # first introduced (non-initialised) cohort up to meet it
+  if(ncol(state) == 1 & i[1] == 2) 
+    return(cbind(state, state - state / 2))
+
+  # otherwise proceed normally
+  i <- i[times[i] == 0]
+
+  if(length(i) == 0)
+    return(state)
+
+  # append a blank cohort
+  dh <- matrix(0, ncol = ncol(state), nrow = nrow(state), 
+               dimnames = dimnames(state)) 
+
+  dh["height", ] <- diff(c(state["height", ], min_height))
+
+  st <- cbind(state, state[, i-1] + dh[, i-1] / 2)
+
+  return(st[, order(st["height", ], decreasing = T)])
+}
+

R/scm_support.R

@@ -178,24 +178,30 @@ make_scm <- function(p, state=NULL) {
     if(n_spp != n_str)
       stop("State object has more species than strategies defined in Parameters")
 
-    cohorts <- sapply(state$species, ncol)
-
     # need to append cohort times to enable integration of net fecundity
     if(state$time != 0)
       message("Solver must start from 0, resetting initial state time")
-    
+
     times <- scm$cohort_schedule$all_times
-    start_time <- sapply(times, function(t) min(t[-1]))
-    new_times <- mapply(function(i, t) c(rep(0, i), t[-1]), 
-                        cohorts, times, SIMPLIFY = F)
+
+    initial_cohorts <- sapply(times, function(t) sum(t == 0), simplify = F)
+    new_cohorts <- mapply(function(s, i) max(0, ncol(s) - i), 
+                          state$species, initial_cohorts, SIMPLIFY = F)
 
+    new_times <- mapply(function(i, t) c(rep(0, i), t), 
+                        new_cohorts, times, SIMPLIFY = F)
+
     # this introduces one more ind. than necessary, but if we
     # overwrite the oldest cohorts first then we can just start at t1
     scm$set_cohort_schedule_times(new_times)
     scm$run_next()
-
-
-    scm$set_state(min(start_time), unlist(state$species), n = cohorts)
+
+    # next step starts from first non-zero time
+    start_time <- sapply(times, function(t) min(t[t>0]))
+
+    # add as many new cohorts as required to fit `state` object
+    scm$set_state(min(start_time), unlist(state$species), 
+                  n = mapply(`+`, new_cohorts, initial_cohorts))
   }  
 
   return(scm)  
@@ -241,8 +247,12 @@ run_scm_error <- function(p, state=NULL) {
   total <- lapply(seq_len(n_spp), function(idx)
                   f(rbind(lai_error[[idx]], net_reproduction_ratio_errors[[idx]])))
 
+  # schedule is needed to update parameters, not sure why ode_times is carried through
+  # saving min height here is lazy, but I'm not sure where else build_schedule can get it
   list(net_reproduction_ratios=scm$net_reproduction_ratios,
        err=list(lai=lai_error, net_reproduction_ratios=net_reproduction_ratio_errors, total=total),
+       schedule=scm$cohort_schedule$all_times,
+       min_heights = sapply(scm$state$species, function(s) min(s["height", ])),
        ode_times=scm$ode_times)
 }
 

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

@@ -3,18 +3,6 @@ context("Initial size distribution")
 # using SCM collect to generate starting conditions is slow, I should just
 # transcribe a few cohorts with realish initial states.
 
-# p0 <- scm_base_parameters("FF16")
-# p0$patch_type = 'fixed'
-# 
-# p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16_hyperpar,FALSE)
-# p1$birth_rate <- 20
-# 
-# out <- run_scm_collect(p1)
-# 
-# # select a time slice
-# i = 120
-# x <- scm_state(i, out)
-# 
 # test_that("Set patch state", {
 #   # editable patch object
 #   types <- extract_RcppR6_template_types(p1, "Parameters")
@@ -77,15 +65,13 @@ test_that("Multi-species case", {
   # Create some arbitrary state
   vars <- sapply(scm$state$species, rownames, simplify = F)
 
-  # warns that Species #3 has no data
-  expect_warning(
-    state <- mapply(function(names, values, n) 
-      matrix(c(values, rep(0, length(names) - 1)),
-             nrow=length(names), ncol=n, dimnames=list(names)),
-      vars, c(2, 3, 4), c(2, 1, 0), SIMPLIFY=FALSE)
-    )
-
-  expect_equal(sapply(state, dim)[2, ], c(2, 1, 0))
+  # this breaks if not introducing a cohort for one species
+  state <- mapply(function(names, values, n) 
+    matrix(c(values, rep(0, length(names) - 1)),
+           nrow=length(names), ncol=n, dimnames=list(names)),
+    vars, c(2, 3, 4), c(2, 1, 1), SIMPLIFY=FALSE)
+
+  expect_equal(sapply(state, dim)[2, ], c(2, 1, 1))
 
   z <- list(time = 2,
             species = state)
@@ -104,19 +90,19 @@ test_that("Multi-species case", {
   scm$set_state(min(start_time), unlist(z$species), n = cohorts)
 
   # state above + one default cohort (min 2cm height)
-  expect_equal(sapply(scm$state$species, rowSums)["height", ], c(6, 5, 2))
+  expect_equal(sapply(scm$state$species, rowSums)["height", ], c(6, 5, 6))
 
 
   # works through the run_scm api
   x <- run_scm_collect(p2, z)
 
   # check fitness
-  expect_equal(x$net_reproduction_ratios, c(1.176e-05, 1.176e-03, 9.132e-04), tolerance = 0.0001)
+  expect_equal(x$net_reproduction_ratios, c(3.981e-06, 4.585e-04, 3.514e-04), tolerance = 0.0001)
 
   # check # states, time steps, cohorts
-  expect_equal(sapply(x$species, dim)[3, ], c(107, 106, 105))
+  expect_equal(sapply(x$species, dim)[3, ], c(107, 106, 106))
 
-  # This is kinda cool, plot the results
+  # Plot the results
   # t2 <- x$time
   # h1 <- x$species[[1]]["height", , ]
   # h2 <- x$species[[2]]["height", , ]
@@ -131,9 +117,54 @@ test_that("Multi-species case", {
   # # Species 2 - blue
   # matlines(t2, h2, lty=1, col=make_transparent(cols[[2]], .25))
   # 
-  # 
   # # Species 3 - black
   # matlines(t2, h3, lty=1, col=make_transparent(cols[[3]], .25))
 })
 
+# Not actually recommending this as a test but it's a good demo for a hard
+# problem - starting with just a single cohort makes it difficult to resolve
+# whether cohorts should be introduced from the initial size distribution
+# or from the cohort schedule.
+
+# `run_scm_error` is somewhat specialised for the build_schedule function,
+# particularly persisting the schedule with cohorts at t0. 
 
+# Left unattended, cohorts introduced below the minimum viable height
+# of the FF16 strategy trip up the adaptive interpolator. Some ugly patches
+# pass the minimum back to split_densities. split_densities doesn't handle
+# density properly.
+test_that("Build schedule works", {
+  p0 <- scm_base_parameters("FF16")
+  p0$patch_type = 'fixed'
+  p0$control$schedule_nsteps = 5
+
+  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$birth_rate <- 20
+
+  init <- matrix(c(10, 0, 0, 0, 0, 0, -10))
+  rownames(init) <- c("height", "mortality", "fecundity", "area_heartwood",
+                      "mass_heartwood", "offspring_produced_survival_weighted",
+                      "log_density")
+
+  state <- list(time = 10, species = list(init))
+  x <- run_scm_collect(p1, state)
+
+  t2 <- x$time
+  h1 <- x$species[[1]]["height", , ]
+
+  # Species 1 - red
+  matplot(t2, h1, lty=1, type="l", col = "black",
+          las=1, xlab="Time (years)", ylab="Height (m)")
+
+  res <- build_schedule(p1, state)
+
+  x <- run_scm_collect(res$parameters, res$state)
+
+  t2 <- x$time
+  h1 <- x$species[[1]]["height", , ]
+
+  # Species 1 - red
+  matplot(t2, h1, lty=1, type="l", col = "black",
+          las=1, xlab="Time (years)", ylab="Height (m)")
+})

tests/testthat/test-schedule-build.R

@@ -21,7 +21,14 @@ test_that("Schedule building", {
     p$strategies <- list(strategy_types[[x]]())
     p$birth_rate <- 0.1
     p <- build_schedule(p)
-    expect_equal(length(p$cohort_schedule_times_default), 141)
-    expect_equal(length(p$cohort_schedule_times[[1]]), 176)
+
+    # This has changed, perhaps due split_times missing a cohort
+    pars <- p$parameters
+    expect_equal(length(pars$cohort_schedule_times_default), 141)
+    expect_equal(length(pars$cohort_schedule_times[[1]]), 176)
+
+    # not really relevant, just showing it's here
+    state <- p$state
+    expect_null(state)
   }
 })