Basic parsing of dependencies and phases

R/dependencies.R

@@ -23,3 +23,30 @@ find_dependencies <- function(expr) {
        variables = unique(variables$get()))
 
 }
+
+
+topological_order <- function(deps) {
+  if (all(lengths(deps) == 0)) {
+    return(list(success = TRUE, value = seq_along(deps)))
+  }
+
+  m <- matrix(FALSE, length(deps), length(deps))
+  for (i in seq_along(deps)) {
+    m[, i] <- unname(names(deps) %in% deps[[i]])
+  }
+
+  pending <- rep(TRUE, length(deps))
+  ret <- integer(0)
+  while (any(pending)) {
+    i <- which(pending)[colSums(m[, pending, drop = FALSE]) == 0]
+    if (length(i) > 0L) {
+      ret <- c(ret, i)
+      pending[i] <- FALSE
+      m[i, ] <- FALSE
+    } else {
+      return(list(success = FALSE, error = which(pending)))
+    }
+  }
+
+  list(success = TRUE, value = ret)
+}

R/parse.R

@@ -4,10 +4,16 @@ odin_parse <- function(expr, input_type = NULL) {
   exprs <- lapply(dat$exprs, function(x) parse_expr(x$value, x, call = call))
 
   system <- parse_system_overall(exprs, call)
+  equations <- parse_system_depends(system$exprs$equations, system$variables,
+                                    call)
+  phases <- parse_system_phases(system$exprs, equations, system$variables, call)
+
   ret <- list(time = system$time,
               class = "odin",
               variables = system$variables,
               parameters = system$parameters,
+              equations = equations,
+              phases = phases,
               data = system$data)
 
   ret

R/parse_system.R

@@ -7,7 +7,7 @@ parse_system_overall <- function(exprs, call) {
   is_compare <- special == "compare"
   is_data <- special == "data"
   is_parameter <- special == "parameter"
-  is_equation <- special == ""
+  is_equation <- special %in% c("", "parameter", "data")
 
   ## We take initial as the set of variables:
   variables <- vcapply(exprs[is_initial], function(x) x$lhs$name)
@@ -96,3 +96,124 @@ parse_system_overall <- function(exprs, call) {
        data = data,
        exprs = exprs)
 }
+
+
+parse_system_depends <- function(equations, variables, call) {
+  implicit <- c(variables, "time", "dt")
+
+  names(equations) <- vcapply(equations, function(eq) eq$lhs$name)
+  deps <- lapply(equations, function(eq) {
+    ## In an earlier proof-of-concept here we also removed eq$lhs$name
+    ## from the dependencies - we do need to do that for arrays at
+    ## least, so at some point some more effort is required here.
+    setdiff(eq$rhs$depends$variables, implicit)
+  })
+  res <- topological_order(deps)
+  if (!res$success) {
+    nms <- names(deps)[res$error]
+    details <- vcapply(nms, function(x) {
+      sprintf("%s: depends on: %s", x, paste(deps[[x]], collapse = ", "))
+    })
+    src <- unname(lapply(equations[res$error], "[[", "src"))
+    odin_parse_error(
+      c("Cyclic dependency detected within equation{?s} {squote(nms)}",
+        set_names(details, "i")),
+      "E2005", src, call)
+  }
+
+  deps <- deps[res$value]
+  deps_recursive <- list()
+  for (nm in names(deps)) {
+    vars <- deps[[nm]]
+    deps_recursive[[nm]] <- union(
+      vars,
+      unlist(deps_recursive[vars], FALSE, FALSE))
+    equations[[nm]]$rhs$depends$variables_recursive <- deps_recursive[[nm]]
+  }
+
+  equations[names(deps)]
+}
+
+
+parse_system_phases <- function(exprs, equations, variables, call) {
+  ## First compute the 'stage' that things occur in; there are only
+  ## three of these, but "time" covers a multitude of sins and
+  ## includes things like the compare function as well as deriv/update
+  ## (and in the case of mixed models *both* deriv/update are
+  ## considered time).
+  stages <- c(system_create = 1,
+              parameter_update = 2,
+              time = 3)
+  implicit <- c(variables, "time", "dt")
+  stage <- c(
+    set_names(rep(stages[["time"]], length(implicit)), implicit))
+  for (nm in names(equations)) {
+    rhs <- equations[[nm]]$rhs
+    if (identical(rhs$type, "parameter")) {
+      is_constant <- isTRUE(rhs$args$constant)
+      stage[[nm]] <-
+        stages[[if (is_constant) "system_create" else "parameter_update"]]
+    } else {
+      stage[[nm]] <- max(stages[["system_create"]],
+                         stage[rhs$depends$variables])
+    }
+  }
+  stage <- set_names(names(stages)[stage], names(stage))
+
+  ## Now, we try and work out which parts of the graph are needed at
+  ## different "phases".  These roughly correspond to dust functions.
+
+  used <- character()
+
+  deps_recursive <- lapply(equations, function(x) {
+    x$rhs$depends$variables_recursive
+  })
+
+  used <- character()
+  required <- character()
+
+  phases <- set_names(vector("list", 5),
+                      c("update", "deriv", "output", "initial", "compare"))
+
+  for (phase in names(phases)) {
+    e <- exprs[[phase]]
+    if (length(e) > 0) {
+      deps <- unique(unlist(lapply(e, function(x) x$rhs$depends$variables),
+                            FALSE, FALSE))
+      eqs <- intersect(names(equations), deps)
+      eqs <- union(eqs, unlist(deps_recursive[eqs], FALSE, FALSE))
+      used <- union(used, eqs)
+
+      eqs_time <- intersect(names(equations), eqs[stage[eqs] == "time"])
+      unpack <- intersect(deps, variables)
+      required <- union(required, eqs[stage[eqs] != "time"])
+
+      if (phase %in% c("update", "deriv", "output")) {
+        phases[[phase]] <- list(unpack = unpack,
+                                equations = eqs_time,
+                                variables = e)
+      } else if (phase == "initial") {
+        ## I forget what the trick was here, but there's some extra
+        ## effort required.
+        if (length(unpack) > 0) {
+          odin_parse_error(
+            "Dependencies within initial conditions not yet supported",
+            "E0001", NULL, call)
+        }
+        phases[[phase]] <- list(equations = eqs_time,
+                                variables = e)
+      } else if (phase == "compare") {
+        phases[[phase]] <- list(equations = eqs_time,
+                                unpack = unpack,
+                                compare = e)
+      }
+    }
+  }
+
+  eqs_shared <- intersect(names(equations), required)
+  phases$build_shared <- list(equations = eqs_shared)
+  phases$update_shared <- list(
+    equations = eqs_shared[stage[eqs_shared] == "parameter_update"])
+
+  phases
+}

R/util.R

@@ -68,3 +68,15 @@ collector <- function(init = character(0)) {
 data_frame <- function(...) {
   data.frame(..., stringsAsFactors = FALSE, check.names = FALSE)
 }
+
+
+set_names <- function(x, nms) {
+  if (length(nms) == 1 && length(x) != 1) {
+    if (is.null(x)) {
+      return(NULL)
+    }
+    nms <- rep_len(nms, length(x))
+  }
+  names(x) <- nms
+  x
+}

tests/testthat/test-dependencies.R

@@ -0,0 +1,25 @@
+test_that("can solve dependencies for trivial system", {
+  expect_equal(topological_order(list()),
+               list(success = TRUE, value = integer()))
+  expect_equal(topological_order(list(a = character())),
+               list(success = TRUE, value = 1))
+  expect_equal(topological_order(list(a = character(), b = character())),
+               list(success = TRUE, value = 1:2))
+})
+
+
+test_that("can solve dependencies for simple system", {
+  deps <- list(a = character(), b = "a", c = "b")
+  expect_equal(topological_order(deps),
+               list(success = TRUE, value = 1:3))
+  deps <- list(a = character(), b = "a", c = "b")
+  expect_equal(topological_order(rev(deps)),
+               list(success = TRUE, value = 3:1))
+})
+
+
+test_that("can report back on cyclic dependencies", {
+  deps <- list(x = character(), y = "x", a = "c", b = "a", c = "a")
+  expect_equal(topological_order(deps),
+               list(success = FALSE, error = 3:5))
+})

tests/testthat/test-parse-system.R

@@ -120,3 +120,42 @@ test_that("parameters default to constant in face of differentiability", {
                           differentiate = c(FALSE, FALSE, TRUE),
                           constant = c(TRUE, TRUE, FALSE)))
 })
+
+
+test_that("fail informatively if recursive depenency in equations", {
+  ## This the simplest, it's just impossible:
+  err <- expect_error(
+    odin_parse({
+      initial(x) <- a
+      deriv(x) <- a
+      a <- a + 1
+    }),
+    "Cyclic dependency detected within equation 'a'")
+})
+
+
+test_that("fail informatively if recursive depenency in several equations", {
+  ## This the simplest, it's just impossible:
+  err <- expect_error(
+    odin_parse({
+      initial(x) <- a
+      deriv(x) <- a
+      a <- c + 1
+      b <- a * 2
+      c <- b / 2
+    }),
+    "Cyclic dependency detected within equations 'a', 'b', and 'c'")
+})
+
+
+test_that("prevent dependencies among variables in initial conditions", {
+  ## ...at least for now
+  expect_error(
+    odin_parse({
+      initial(a) <- 1
+      initial(b) <- a + 1
+      deriv(a) <- 0
+      deriv(b) <- 0
+    }),
+    "Dependencies within initial conditions not yet supported")
+})

tests/testthat/test-util.R

@@ -32,3 +32,16 @@ test_that("can match a simple call", {
     suppressWarnings(match_call(quote(f(fo = 2)), fn)),
     list(success = TRUE, value = quote(f(foo = 2, bar = 1))))
 })
+
+
+test_that("set_names copes with common pathologies", {
+  expect_equal(set_names(character(), "x"),
+               structure(character(), names = character()))
+  expect_equal(set_names("a", "x"),
+               c("x" = "a"))
+  expect_equal(set_names(c("a", "b"), "x"),
+               c("x" = "a", x = "b"))
+  expect_equal(set_names(c("a", "b"), c("x", "y")),
+               c("x" = "a", y = "b"))
+  expect_null(set_names(NULL, "x"))
+})

vignettes/errors.Rmd

@@ -256,3 +256,25 @@ All variables used in `deriv()` or `update()` require a corresponding entry in `
 Variables are missing calls to `deriv()` or `update()`
 
 You have a system where you use different equations for `deriv()`/`update()` to the variables defined in `initial()`.  This is an error if there are equations in `deriv()`/`update()` that don't have a corresponding equation using `initial()`, or if you have equations in `initial()` that don't have a corresponding `deriv()` or `update()` equation.  The error will highlight all lines that might be involved in the error.
+
+# `E2005`
+
+Cyclic dependency detected within equations.  There are a few ways this can happen.  The simplest is that your equation references itself, for example:
+
+```r
+a <- a + 1
+```
+
+Unlike in R, this is disallowed, as each variable may only be assigned to once within the target function of your system.  Each assignment is much more like mathematical equation than usual programming statements.
+
+You can get more complicated cycles, for example:
+
+```r
+a <- c / 2
+b <- sqrt(a)
+c <- a + 1
+```
+
+Here `a` depends on `c`, `c` depends on `b` and `b` depends on `a`.  The error will reference all the variables involved in cycle.
+
+It is possible that there is more than one cycle within the reported expressions.