tests/testthat/test_accumulate.R

In hydroloom: Utilities to Weave Hydrologic Fabrics

test_that("accumulate downstream", {
  test_data <- data.frame(id = c(1, 2, 3, 4, 6, 7, 8, 9),
    toid = c(2, 3, 4, 0, 7, 8, 9, 4),
    a = c(1, 2, 3, 4, 1, 2, 3, 4))

  expect_message(
    expect_equal(accumulate_downstream(test_data, "a"),
      c(1, 3, 6, 20, 1, 3, 6, 10)), "Dendritic routing will be applied")

  expect_error(accumulate_downstream(test_data, "b"), "b must be in x")

  expect_error(accumulate_downstream(dplyr::rename(test_data, borked = toid), "a"),
    "requires hy_topo")

  x <- sf::read_sf(system.file("extdata", "walker.gpkg", package = "hydroloom"))

  x <- add_toids(x, return_dendritic = TRUE)

  expect_message(
    y <- accumulate_downstream(x, "AreaSqKM"),
    "Dendritic")

  expect(mean(abs(y - x$TotDASqKM)) < 1e-3, "drainage area not close enough")
  expect(max(abs(y - x$TotDASqKM)) < 1e-2, "drainage area not close enough")

  x$AreaSqKM[1] <- NA

  expect_warning(accumulate_downstream(x, "AreaSqKM", quiet = TRUE), "NA values found")
})

test_that("divergences with total", {
  x <- sf::read_sf(system.file("extdata/new_hope.gpkg", package = "hydroloom"))

  net <- navigate_network_dfs(x, 8893236, "up")

  x <- x[x$COMID %in% unlist(net), ]

  z <- x |>
    dplyr::select(COMID, FromNode, ToNode, Divergence, AreaSqKM, TotDASqKM)

  start <- c(8895516, 8893238, 8893180, 8893238, 8893216)
  remove <- navigate_network_dfs(z, start, "up")

  remove <- unlist(remove)
  remove <- remove[!remove %in% start]

  z <- z[!z$COMID %in% remove, ]
  # nolint start
  # d <- sf::st_transform(sf::st_simplify(sf::st_transform(z, 5070),
  #                                       dTolerance = units::as_units(100, "m")),
  #                       4326)
  # d <- add_toids(d, return_dendritic = FALSE)
  # d <- sort_network(d)
  # d <- select(d, -toid) |> distinct()
  # d$id <- 1:nrow(d)
  #
  # nodes <- 1:length(unique(c(d$FromNode, d$ToNode)))
  # nodes <- data.frame(new_node = nodes,
  #                     old_node = unique(c(d$FromNode, d$ToNode)))
  #
  # d <- left_join(d, rename(nodes, fromnode = "new_node"),
  #                by = c("FromNode" = "old_node"))
  # d <- left_join(d, rename(nodes, tonode = "new_node"),
  #                by = c("ToNode" = "old_node"))
  #
  # d <- select(d, id, fromnode, tonode,
  #             divergence = Divergence,
  #             areasqkm = AreaSqKM)

  # sf::write_sf(d, "tests/testthat/data/simple_diversions.geojson")
  # nolint end

  z$dend_totdasqkm <- accumulate_downstream(z, var = "AreaSqKM", quiet = TRUE)

  # for dendritic, we shoud have the outlet be the sum of everything
  expect_equal(max(z$dend_totdasqkm), sum(z$AreaSqKM), tolerance = 1e-6)
  # a diversion should have its own area as its total
  expect_equal(z$dend_totdasqkm[z$COMID == 8893212], z$AreaSqKM[z$COMID == 8893212], tolerance = 1e-6)

  z$divergence_fraction <- 1
  z$divergence_fraction[z$Divergence == 2] <- 0.25
  z$divergence_fraction[z$Divergence == 1] <- 0.75

  z$div_totdasqkm <- accumulate_downstream(z, var = "AreaSqKM")

  # we shoud have the outlet be the sum of everything
  expect_equal(max(z$div_totdasqkm), sum(z$AreaSqKM), tolerance = 1e-6)
  # a diversion should have its own area plus its diversion fraction from upstream
  expect_equal(z$div_totdasqkm[z$COMID == 8893212],
    z$AreaSqKM[z$COMID == 8893212] + 0.25 * z$div_totdasqkm[z$COMID == 8893184], tolerance = 1e-6)
  # a main should have its own area plus its diversion fraction from upstream
  expect_equal(z$div_totdasqkm[z$COMID == 8893182],
    z$AreaSqKM[z$COMID == 8893182] + 0.75 * z$div_totdasqkm[z$COMID == 8893184], tolerance = 1e-6)

  z$tot_totdasqkm <- accumulate_downstream(z, "AreaSqKM", total = TRUE)

  # a diversion should have its own area plus everything from upstream
  expect_equal(z$tot_totdasqkm[z$COMID == 8893210],
    z$AreaSqKM[z$COMID == 8893210] + z$tot_totdasqkm[z$COMID == 8893190], tolerance = 1e-6)
  # a main should have its own area plus its diversion fraction from upstream
  expect_equal(z$tot_totdasqkm[z$COMID == 8893186],
    z$AreaSqKM[z$COMID == 8893186] + z$tot_totdasqkm[z$COMID == 8893190], tolerance = 1e-6)

  # need all the stuff upstream but only reachable on a diversion in this goofy edge case
  expect_equal(z$tot_totdasqkm[z$COMID == 8893202],
    z$AreaSqKM[z$COMID == 8893202] + sum(z$AreaSqKM[z$COMID == 8893192],
      z$AreaSqKM[z$COMID == 8893172],
      z$AreaSqKM[z$COMID == 8893176],
      z$AreaSqKM[z$COMID == 8893174],
      z$tot_totdasqkm[z$COMID == 8893194]), tolerance = 1e-6)

  expect_equal(z$tot_totdasqkm[z$COMID == 8893218],
    z$AreaSqKM[z$COMID == 8893218] + sum(z$tot_totdasqkm[z$COMID == 8893202],
      z$AreaSqKM[z$COMID == 8893198],
      z$AreaSqKM[z$COMID == 8893208],
      z$AreaSqKM[z$COMID == 8893212]), tolerance = 1e-6)

  # we shoud have the outlet be the sum of everything
  expect_equal(max(z$tot_totdasqkm), sum(z$AreaSqKM), tolerance = 1e-6)

  z <- x |>
    dplyr::select(COMID, LevelPathI, FromNode, ToNode, Divergence, AreaSqKM, TotDASqKM)

  z$tot_totdasqkm <- accumulate_downstream(z, "AreaSqKM", total = TRUE)

  expect_equal(z$tot_totdasqkm, z$TotDASqKM, tolerance = 1e-6)
})

test_that("simple diversions total", {
  x <- sf::read_sf(list.files(pattern = "simple_diversions.geojson", full.names = TRUE, recursive = TRUE))

  x$tot_totareasqkm <- accumulate_downstream(x, "areasqkm", total = TRUE)

  check_fun <- function(check_id) {
    up_net <- navigate_network_dfs(x, check_id, "up")

    expect_equal(x$tot_totareasqkm[x$id == check_id],
      sum(x$areasqkm[x$id %in% unique(unlist(up_net))]), tolerance = 1e-6)
  }

  # need all the stuff upstream but only reachable on a diversion in this goofy edge case
  check_fun(25)
  check_fun(27)
  check_fun(21)
  check_fun(29)

  # we shoud have the outlet be the sum of everything
  expect_equal(max(x$tot_totareasqkm), sum(x$areasqkm), tolerance = 1e-6)

})

test_that("complex diversions", {

  net <- read.csv(list.files(pattern = "diversions.csv", full.names = TRUE, recursive = TRUE))

  net$tot_totareasqkm <- accumulate_downstream(net, "areasqkm", TRUE)

  # these compare to NHDPlus with the addition of 184.3632 that was removed for testing
  # this only holds on the path downstream of the main set of diversions in this example
  net$diff <- net$tot_totareasqkm - (net$totdasqkm - 184.3632)

  expect_equal(net$tot_totareasqkm[net$comid == 14702428], sum(
    net$areasqkm[net$comid == 14702428],
    net$areasqkm[net$comid == 14702926],
    net$areasqkm[net$comid == 14702436],
    net$areasqkm[net$comid == 14703168],
    net$areasqkm[net$comid == 14702438]
  ), tolerance = 1e-6)

  check_fun <- function(check_comid) {
    up_net <- navigate_network_dfs(net, check_comid, "up")

    expect_equal(net$tot_totareasqkm[net$comid == check_comid],
      sum(net$areasqkm[net$comid %in% unique(unlist(up_net))]), tolerance = 1e-6)
  }

  check_fun(14702384)
  check_fun(14702360)
  check_fun(14702406)
  check_fun(14702402)
  check_fun(14702400)
  check_fun(14702376)
  check_fun(14702374)
  check_fun(14702378)
  check_fun(14702336)
  check_fun(14702328)
  check_fun(14702352)

})

test_that("part closed test", {

  net <- read.csv(text =
    "id,toid,divergence
  1,2,0
  1,3,0
  3,4,2
  3,5,2
  4,6,2
  5,6,1
  2,7,1
  6,7,0
  7,0,0")

  net$val <- 1

  expect_equal(accumulate_downstream(net, "val", TRUE), c(1, 1, 2, 2, 3, 3, 2, 5, 7))

  net <- read.csv(text =
    "id,toid,divergence
  1,2,0
  1,3,0
  3,4,1
  3,5,1
  4,6,2
  5,6,1
  2,7,2
  6,7,0
  7,0,0")

  net$val <- 1

  expect_equal(accumulate_downstream(net, "val", TRUE), c(1, 1, 2, 2, 3, 3, 2, 5, 7))

  net <- read.csv(text =
    "id, toid, divergence
                  1,2,0
                  1,3,0
                  1,4,0
                  5,6,0
                  2,6,1
                  3,7,2
                  4,7,2
                  6,8,0
                  7,8,0
                  8,0,0")

  net$val <- 1

  expect_equal(accumulate_downstream(net, "val", total = TRUE), c(1, 1, 1, 1, 2, 2, 2, 4, 4, 8))

})

test_that("accumulation scales linearly in the number of rows", {

  skip_on_cran()

  # The accumulation loop used to write `out[[var]][i] <- ...` into a
  # tibble column on every iteration, which re-assigns the whole column
  # each time and makes the function quadratic in nrow(x). Measured before
  # the fix: 1.6 s at 25k rows, 4.3 s at 50k, 11.6 s at 100k, 31.8 s at
  # 200k. That is invisible on the small fixtures every other test in this
  # file uses, and it is the difference between seconds and hours on a
  # continental network.
  #
  # This guards the shape of the curve, not an absolute time, so it does
  # not fail on a slow machine. A quadratic implementation quadruples when
  # n doubles; a linear one roughly doubles. The threshold is loose enough
  # to absorb ordinary noise and still catch a regression.

  make_net <- function(n) {
    hy(data.frame(id = as.character(seq_len(n)),
                  toid = as.character(c("", floor(seq_len(n)[-1] / 2))),
                  v = 1,
                  stringsAsFactors = FALSE))
  }

  time_it <- function(n) {
    net <- make_net(n)
    t0 <- Sys.time()
    invisible(suppressMessages(accumulate_downstream(net, "v", quiet = TRUE)))
    as.numeric(difftime(Sys.time(), t0, units = "secs"))
  }

  t1 <- time_it(40000L)
  t2 <- time_it(80000L)

  # Guard against a near-zero baseline making the ratio meaningless.
  skip_if(t1 < 0.02, "baseline too fast to measure a ratio reliably")

  expect_lt(t2 / t1, 3)
})

test_that("accumulation is unchanged by hoisting the accumulator", {

  # Pins the values the hoisted loop produces, for both branches, against
  # results computed by hand from the fixture topology.

  # Dendritic: 1 -> 2 -> 3 -> 4, and 6 -> 7 -> 8 -> 9 -> 4.
  test_data <- data.frame(id = c(1, 2, 3, 4, 6, 7, 8, 9),
                          toid = c(2, 3, 4, 0, 7, 8, 9, 4),
                          a = c(1, 2, 3, 4, 1, 2, 3, 4))

  expect_equal(suppressMessages(accumulate_downstream(test_data, "a")),
               c(1, 3, 6, 20, 1, 3, 6, 10))

  # The `total = TRUE` branch shares the same hoisting change.
  net <- hy(data.frame(id = as.character(1:7),
                       toid = as.character(c(3, 3, 5, 5, 7, 7, "")),
                       val = 1, stringsAsFactors = FALSE))
  expect_equal(accumulate_downstream(net, "val", total = TRUE),
               c(1, 1, 3, 1, 5, 1, 7))

  # And a real network still reproduces its own published totals.
  x <- add_toids(
    sf::read_sf(system.file("extdata", "walker.gpkg", package = "hydroloom")),
    return_dendritic = TRUE)
  y <- suppressMessages(accumulate_downstream(x, "AreaSqKM"))
  expect_lt(max(abs(y - x$TotDASqKM)), 1e-2)
})