Requirements relaxed for simulation of transfers

- Scenarios are not required to start at t=0, negative output times are also possible
- Transfer times no longer need to part of the output times

R/set_transfer.R

@@ -151,3 +151,10 @@ set_transfer2 <- function(x, interval, times, biomass, scaled_comp) {
   }
   x
 }
+
+set_notransfer <- function(x) {
+  if(length(x) > 1)
+    return(sapply(x, set_notransfer))
+
+  set_transfer(x, interval=-1)
+}

R/simulate.R

@@ -180,27 +180,35 @@ simulate_transfer <- function(scenario, times, in_sequence=FALSE, ...) {
 
   times <- scenario@times
   t_min <- min(times)
-  if(t_min < 0)
-    stop("output times must not be negative")
   t_max <- max(times)
 
   # find transfer time points that occur during the simulated period
   tr_points <- vector("numeric")
   ends_on_transfer <- FALSE
   if(has_regular_transfer(scenario)) {
-    # sim has to be long enough for transfers to occur
-    if(t_max >= scenario@transfer.interval) {
-      tr_points <- seq(scenario@transfer.interval, t_max, scenario@transfer.interval)
-    }
+    seq_start <- ceiling(t_min / scenario@transfer.interval) * scenario@transfer.interval
+    seq_end <- t_max
+    if(seq_start <= seq_end) # avoid errors in case no transfer occurs
+      tr_points <- seq(seq_start, seq_end, scenario@transfer.interval)
   } else {
     tr_points <- scenario@transfer.times
   }
+  # does simulation period end on a transfer?
   if(length(tr_points) > 0)
-    ends_on_transfer <- tail(tr_points, 1) == t_max
+    ends_on_transfer <- tr_points[length(tr_points)] == t_max
   # limit vector to transfer time points which occur during the simulated period.
   # but exclude the starting time as a potential transfer
   tr_points <- tr_points[tr_points > t_min & tr_points <= t_max]
 
+  # if no transfers occurs during simulated period -> early exit
+  if(length(tr_points) == 0 & !ends_on_transfer)
+    return(simulate_scenario(set_times(scenario, times), ...))
+
+  # add all transfer time points to the output times vector
+  tr_inserted <- setdiff(tr_points, times)
+  if(length(tr_inserted) > 0)
+    times <- sort(c(times, tr_inserted))
+
   # biomass level after each transfer
   if(length(scenario@transfer.biomass) == 1)
     tr_biomass <- rep(scenario@transfer.biomass, times=length(tr_points))
@@ -209,34 +217,26 @@ simulate_transfer <- function(scenario, times, in_sequence=FALSE, ...) {
   else
     stop("length of biomass and transfer times vectors do not match", call.=FALSE)
 
-  # if no transfers occur -> early exit
-  if(length(tr_points) == 0 & !ends_on_transfer)
-    return(simulate_scenario(scenario, times=times, ...))
-
-  # transfer time points must be contained in output time points
-  if(length(setdiff(tr_points, times)) > 0)
-    stop(paste("transfer time points missing in output times:", paste(setdiff(tr_points, times), sep=",")))
-
   df <- data.frame()
   t_start <- t_min
-  #tr_points <- c(tr_points, t_max)
-  #tr_biomass <- c(tr_biomass, NA_real_) # dummy value
+  # add last time point to also simulate the remainder of the scenario, if it exists
+  if(!ends_on_transfer)
+    tr_points <- c(tr_points, t_max)
   # simulate population until next transfer to new medium
-  for(i in seq_along(tr_points)) {
-    t <- tr_points[i]
-    period <- times[times >= t_start & times <= t]
-    # simulate
+  for(i in seq_along(tr_points))
+  {
+    tr <- tr_points[i]
+    period <- times[times >= t_start & times <= tr]
     out <- solver(set_times(scenario, period), ...)
 
-    # append results to output data.frame
-    if(i==1) {
-      df <- rbind(df, out)
-    } else {
-      # append simulation results, but exclude values for the starting time point
-      # where the transfer occurred. the starting time point `t` could be included more
-      # than once in the `times` or `period` vector
-      df <- rbind(df, out[-seq(sum(period == t)),])
+    # select rows/time points:
+    # 1) that do NOT occur directly after a transfer, i.e. the start of a simulation
+    selector <- i == 1 | period != t_start
+    # 2) that were NOT inserted for technical reasons
+    if(length(tr_inserted) > 0) {
+      selector <- selector & !(period %in% tr_inserted)
     }
+    df <- rbind(df, out[selector, ])
 
     # Transfer a fixed number of biomass to a new medium:
     # use last state as starting point
@@ -248,13 +248,7 @@ simulate_transfer <- function(scenario, times, in_sequence=FALSE, ...) {
     # scale other compartments relative to new biomass
     init[scenario@transfer.comp.scaled] <- init[scenario@transfer.comp.scaled] * BM.fac
     scenario <- set_init(scenario, init)
-    t_start <- t
-  }
-  # simulate the remainder of the scenario if any exists
-  if(t_max > max(tr_points)) {
-    period <- times[times >= t_start]
-    out <- solver(set_times(scenario, period), ...)
-    df <- rbind(df, out[-seq(sum(period == t)),])
+    t_start <- tr
   }
 
   # add metadata at which system state the subsequent simulation should start

tests/testthat/test-simulate_transfer.R

@@ -1,36 +1,159 @@
-test_that("regular intervals", {
-  metsulfuron %>%
+test_that("regular intervals, starting and ending on transfer", {
+  sc <- metsulfuron %>%
     set_init(c(BM=1)) %>%
+    set_times(0:10) %>%
     set_noexposure() %>%
-    set_transfer(interval=2, biomass=1) %>%
-    simulate() -> rs
+    set_transfer(interval=5, biomass=1)
+
+  ref <- sc %>%
+    set_notransfer() %>%
+    set_times(0:5) %>%
+    simulate()
+
+  rs <- simulate(sc)
+  # first batch must equal to ref sim
+  expect_equal(rs[1:6, ], ref, tolerance=1e-6, ignore_attr=TRUE)
+  # second batch must also equal to ref sim, apart from the timestamp and sim start (t==0)
+  expect_equal(rs[7:11, -1], ref[2:6, -1], tolerance=1e-6, ignore_attr=TRUE)
+
+  # negative starting time
+  rs <- sc %>%
+    set_times(-5:5) %>%
+    simulate()
+  expect_equal(rs[1:6, -1], ref[, -1], tolerance=1e-6, ignore_attr=TRUE)
+  expect_equal(rs[7:11, -1], ref[2:6, -1], tolerance=1e-6, ignore_attr=TRUE)
+})
+
+test_that("regular intervals, not starting on transfer", {
+  sc <- metsulfuron %>%
+    set_init(c(BM=1)) %>%
+    set_times(1:10) %>%
+    set_noexposure() %>%
+    set_transfer(interval=5, biomass=1)
+
+  ref <- sc %>%
+    set_notransfer() %>%
+    set_times(0:5) %>%
+    simulate()
 
-  expect_equal(rs[,-1], rs[c(1:3, rep(2:3, 6)),-1], tolerance=1e-6, ignore_attr=T)
+  rs <- simulate(sc)
+  expect_equal(rs[1:5, -1], ref[1:5, -1], tolerance=1e-6, ignore_attr=TRUE)
+  expect_equal(rs[6:10, -1], ref[2:6, -1], tolerance=1e-6, ignore_attr=TRUE)
+
+  # negative starting time
+  rs <- sc %>%
+    set_times(-4:5) %>%
+    simulate()
+  expect_equal(rs[1:5, -1], ref[1:5, -1], tolerance=1e-6, ignore_attr=TRUE)
+  expect_equal(rs[6:10, -1], ref[2:6, -1], tolerance=1e-6, ignore_attr=TRUE)
+})
+
+test_that("regular intervals, not ending on transfer", {
+  sc <- metsulfuron %>%
+    set_init(c(BM=1)) %>%
+    set_times(0:9) %>%
+    set_noexposure() %>%
+    set_transfer(interval=5, biomass=1)
+
+  ref <- sc %>%
+    set_notransfer() %>%
+    set_times(0:5) %>%
+    simulate()
+
+  rs <- simulate(sc)
+  expect_equal(rs[1:6, ], ref[1:6, ], tolerance=1e-6, ignore_attr=TRUE)
+  expect_equal(rs[7:10, -1], ref[2:5, -1], tolerance=1e-6, ignore_attr=TRUE)
+
+  # negative starting time
+  rs <- sc %>%
+    set_times(-5:4) %>%
+    simulate()
+  expect_equal(rs[1:6, -1], ref[1:6, -1], tolerance=1e-6, ignore_attr=TRUE)
+  expect_equal(rs[7:10, ], ref[2:5, ], tolerance=1e-6, ignore_attr=TRUE)
 })
 
+test_that("regular intervals, no transfers relevant", {
+  # starts and ends transfer
+  sc <- metsulfuron %>%
+    set_init(c(BM=1)) %>%
+    set_times(0:5) %>%
+    set_noexposure() %>%
+    set_transfer(interval=5, biomass=1)
+  rs <- simulate(sc)
+  ref <- sc %>%
+    set_notransfer() %>%
+    simulate()
+
+  expect_equal(rs, ref, tolerance=1e-6)
+
+  # negative starting time
+  rs <- sc %>%
+    set_times(-5:0) %>%
+    simulate()
+  expect_equal(rs[, -1], ref[, -1], tolerance=1e-6, ignore_attr=TRUE)
+
+  # starts and ends somewhere in between transfers
+  sc <- sc %>% set_times(1:4)
+  rs <- sc %>% simulate()
+  ref <- sc %>%
+    set_notransfer() %>%
+    simulate()
+
+  expect_equal(rs, ref, tolerance=1e-6)
+})
+
+
 test_that("regular intervals with redundant time-points", {
-  metsulfuron %>%
+  times <- c(0, 0, 1:5, 5, 6, 6, 7:10)
+
+  sc <- metsulfuron %>%
     set_init(c(BM=1)) %>%
     set_noexposure() %>%
-    set_times(rep(0:13, each=2)) %>%
-    set_transfer(interval=2, biomass=1) %>%
-    simulate() -> rs
+    set_times(times) %>%
+    set_transfer(interval=5, biomass=1)
+  rs <- simulate(sc)
 
-  expect_equal(rs$time, rep(0:13, each=2))
-  # starts at 0 with init, then values at t+1, t+2 should repeat until the end
-  expect_equal(rs$BM[1:2], c(1, 1))
-  expect_equal(rs$BM[c(3,4,7,8,27,28)], rep(1.0881825, 6), tolerance=1e-5)
-  expect_equal(rs$BM[c(5,6,9,10)], rep(1.18407, 4), tolerance=1e-5)
+  ref <- sc %>%
+    set_notransfer() %>%
+    set_times(times[times <= 5]) %>%
+    simulate()
+
+  expect_equal(rs$time, times)
+  expect_equal(rs[1:8, ], ref, tolerance=1e-6)
+  expect_equal(rs[9:14, -1], ref[c(3, 3:7), -1], tolerance=1e-6, ignore_attr=TRUE)
+})
+
+test_that("transfer not in output times", {
+  sc <- metsulfuron %>%
+    set_init(c(BM=1)) %>%
+    set_noexposure() %>%
+    set_times(c(1, 10)) %>%
+    set_transfer(interval=5, biomass=1)
+  rs <- simulate(sc)
+
+  ref <- sc %>%
+    set_notransfer() %>%
+    set_times(0:5) %>%
+    simulate()
+
+  expect_equal(rs[, 1], c(1, 10))
+  expect_equal(rs[1, -1], ref[1, -1], tolerance=1e-6, ignore_attr=TRUE)
+  expect_equal(rs[2, -1], ref[6, -1], tolerance=1e-6, ignore_attr=TRUE)
 })
 
 test_that("custom time points", {
-  metsulfuron %>%
+  sc <- metsulfuron %>%
     set_init(c(BM=1)) %>%
     set_noexposure() %>%
-    set_transfer(times=c(2,4,6,8,10,12), biomass=1) %>%
-    simulate() -> rs
+    set_transfer(times=c(3, 8), biomass=1)
+  rs <- simulate(sc)
+  ref <- sc %>%
+    set_notransfer() %>%
+    simulate()
 
-  expect_equal(rs[,-1], rs[c(1:3, rep(2:3, 6)),-1], tolerance=1e-6, ignore_attr=T)
+  expect_equal(rs[1:4, -1], ref[1:4, -1], tolerance=1e-6) # first batch
+  expect_equal(rs[5:9, -1], ref[2:6, -1], tolerance=1e-6, ignore_attr=TRUE) # second batch
+  expect_equal(rs[10:15, -1], ref[2:7, -1], tolerance=1e-6, ignore_attr=TRUE) # third batch
 })
 
 test_that("custom biomass", {