diff --git a/.Rbuildignore b/.Rbuildignore index 7a6244ac..694a0c82 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -11,3 +11,4 @@ ^arf.*\.tgz$ ^doc$ ^Meta$ +^bench$ diff --git a/.gitignore b/.gitignore index 6ad90e01..e470cecf 100644 --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,5 @@ attic/ branch/ /doc/ /Meta/ +bench/logs/ +bench/results/ diff --git a/DESCRIPTION b/DESCRIPTION index ee21c86c..538224cb 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: arf Title: Adversarial Random Forests -Version: 0.2.5 +Version: 0.3.0 Authors@R: c( person("Marvin N.", "Wright", , "cran@wrig.de", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-8542-6291")), @@ -9,7 +9,9 @@ Authors@R: c( person("Kristin", "Blesch", role = "aut", comment = c(ORCID = "0000-0001-6241-3079")), person("Jan", "Kapar", role = "aut", - comment = c(ORCID = "0009-0000-6408-2840")) + comment = c(ORCID = "0009-0000-6408-2840")), + person("Lukas", "Burk", , "cran@lukasburk.de", role = "ctb", + comment = c(ORCID = "0000-0001-7528-3795")) ) Description: Adversarial random forests (ARFs) recursively partition data into fully factorized leaves, where features are jointly independent. @@ -33,19 +35,22 @@ Imports: ranger, stringr, truncnorm -Suggests: +Suggests: doFuture, doParallel, ggplot2, knitr, + mirai, mlbench, + mori, palmerpenguins, + pkgload, rmarkdown, testthat (>= 3.0.0), tibble -VignetteBuilder: +VignetteBuilder: knitr +Config/roxygen2/version: 8.0.0 Config/testthat/edition: 3 Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.2 diff --git a/NEWS.md b/NEWS.md index 3504869b..8e1718a4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,10 @@ +# arf 0.3.0 +* Add mirai/mori parallel backend as an alternative to foreach/doParallel (#62) + * Shares large read-only inputs (training data, forest, learned parameters) across workers via mori, lowering memory use in `adversarial_rf()`, `forde()`, `forge()`, `expct()`, and `lik()` + * Enable with active mirai daemons or `options(arf.backend)`, see `?arf-options` +* Reduce peak memory of `expct()` by processing conditions in bounded blocks (up to ~16x lower in internal benchmarks on large forests with many conditions, tune via `options(arf.block_rows)`, see `?arf-options`) +* Reduce memory and dispatch overhead in `forde()` (fused per-tree parameter pass, per-tree coverage) and `lik()` (per-batch reduction) + # arf 0.2.5 * Export sample_from_leaves() for intra-leaf marginal sampling @@ -27,4 +34,4 @@ * Speed boost for the adversarial resampling step * Early stopping option for adversarial training * alpha parameter for regularizing multinomial distributions in forde -* Unified treatment of colnames with internal semantics (y, obs, tree, leaf) \ No newline at end of file +* Unified treatment of colnames with internal semantics (y, obs, tree, leaf) diff --git a/R/adversarial_rf.R b/R/adversarial_rf.R index ebd88af9..60d28a10 100644 --- a/R/adversarial_rf.R +++ b/R/adversarial_rf.R @@ -16,8 +16,10 @@ #' @param prune Impose \code{min_node_size} by pruning? #' @param verbose Print discriminator accuracy after each round? Will also show #' additional warnings. -#' @param parallel Compute in parallel? Must register backend beforehand, e.g. -#' via \code{doParallel} or \code{doFuture}; see examples. +#' @param parallel Compute in parallel? Enables multithreaded ranger training +#' (no backend needed) and parallelizes the pruning step, which requires a +#' registered \code{foreach} backend (\code{doParallel}, \code{doFuture}) or +#' active \code{mirai} daemons. See \code{\link{arf-options}}. #' @param ... Extra parameters to be passed to \code{ranger}. #' #' @details @@ -92,6 +94,9 @@ #' # ... or with doFuture #' doFuture::registerDoFuture() #' future::plan("multisession", workers = 4) +#' +#' # ... or with mirai (shares large read-only inputs across workers via mori) +#' mirai::daemons(4) #' } #' #' @seealso @@ -183,41 +188,52 @@ adversarial_rf <- function( } } - # Prune leaves to ensure min_node_size w.r.t. real data + # Prune leaves to ensure min_node_size w.r.t. real data. Per-tree work lives in + # arf_prune_tree() (prune_workers.R). This is a meaningful share of runtime once + # ranger training is threaded (the prune loop is pure R, unaffected by + # ranger's num.threads), so it gets the same backend treatment as the rest. if (isTRUE(prune)) { pred <- stats::predict(rf0, x_real, type = 'terminalNodes')$predictions + 1L - prune <- function(tree) { - # Nodes to prune are leaves which contain fewer than min_node_size real samples - out <- rf0$forest$child.nodeIDs[[tree]] - leaves <- which(out[[1]] == 0L) - to_prune <- leaves[!(leaves %in% which(tabulate(pred[, tree]) >= min_node_size))] - while(length(to_prune) > 0) { - if (1 %in% to_prune) { - # Never prune the root - break - } - for (tp in to_prune) { - # Find parent - parent <- which((out[[1]] + 1L) == tp) - if (length(parent) > 0) { - # If node to prune (tp) is the left child of parent, replace left child with right child - out[[1]][parent] <- out[[2]][parent] - } else { - # If node to prune (tp) is the right child of parent, replace right child with left child - parent <- which((out[[2]] + 1L) == tp) - out[[2]][parent] <- out[[1]][parent] - } - } - # If both children of a parent are to be pruned, prune the parent in the next round - # This happens if both children have been pruned - to_prune <- which((out[[1]] + 1L) %in% to_prune) - } - return(out) + prune_one <- function(b) { + arf_prune_tree(b, rf0$forest$child.nodeIDs, pred, min_node_size) } - if (isTRUE(parallel)) { - rf0$forest$child.nodeIDs <- foreach(b = seq_len(num_trees)) %dopar% prune(b) + use_mirai <- FALSE + if (num_trees > 1) { + backend <- arf_select_backend(parallel) + use_mirai <- identical(backend, "mirai") + } + if (use_mirai) { + # arf_prune_tree's body is base-R, so pass it as an object (daemons need no + # arf loaded). Share the two big read-only objects (pred is n x num_trees; + # child.nodeIDs is the forest) once via mori. Chunk contiguously and c() so + # the flat result keeps tree order 1..num_trees (unname: mirai_map names + # chunks, but child.nodeIDs must stay an unnamed list). + pred_shared <- mori::share(pred) + child_shared <- mori::share(rf0$forest$child.nodeIDs) + chunks <- arf_tree_chunks(num_trees, mirai::status()$connections) + chunk_fn <- function(trees, worker, child_nodeIDs, pred, min_node_size) { + lapply(trees, worker, child_nodeIDs = child_nodeIDs, + pred = pred, min_node_size = min_node_size) + } + # Both functions are base-R with explicit args: strip their environments + # before shipping. chunk_fn's would otherwise be THIS frame (rf0, dat, + # x_real: hundreds of MB serialized into every task); arf_prune_tree's + # namespace env would force daemons to load arf. + # See the closure note above arf_mirai_tree_map() in mirai_helpers.R. + environment(chunk_fn) <- globalenv() + prune_worker <- arf_prune_tree + environment(prune_worker) <- globalenv() + res <- mirai::mirai_map(chunks, chunk_fn, + .args = list(worker = prune_worker, + child_nodeIDs = child_shared, + pred = pred_shared, + min_node_size = min_node_size))[] + arf_stop_on_mirai_error(res) + rf0$forest$child.nodeIDs <- unname(do.call(c, res)) + } else if (isTRUE(parallel) && num_trees > 1) { + rf0$forest$child.nodeIDs <- foreach(b = seq_len(num_trees)) %dopar% prune_one(b) } else { - rf0$forest$child.nodeIDs <- foreach(b = seq_len(num_trees)) %do% prune(b) + rf0$forest$child.nodeIDs <- foreach(b = seq_len(num_trees)) %do% prune_one(b) } } diff --git a/R/arf-options.R b/R/arf-options.R new file mode 100644 index 00000000..e3841c7c --- /dev/null +++ b/R/arf-options.R @@ -0,0 +1,103 @@ +#' arf package options +#' +#' Options controlling the parallel backend and its messaging, set via +#' \code{\link{options}}. +#' +#' @details +#' \describe{ +#' \item{\code{arf.backend}}{Parallel backend used when \code{parallel = TRUE}: +#' \code{"foreach"} or \code{"mirai"}. If unset, arf uses \code{"mirai"} when +#' mirai daemons are running and \code{"foreach"} otherwise.} +#' \item{\code{arf.verbose}}{Report the selected backend once per backend +#' configuration per session? Default \code{TRUE}; set \code{FALSE} to +#' silence.} +#' \item{\code{arf.block_rows}}{Cap on rows materialized per block of +#' conditions in \code{\link{expct}}. Default \code{5e6}. Lower it to +#' trade speed for memory on large forests with many conditions.} +#' } +#' +#' \code{arf.block_rows} does not affect results, only peak memory and speed. +#' For memory-constrained hardware, combine a small daemon count with a lower +#' \code{arf.block_rows} and the \code{batch}/\code{stepsize} arguments of +#' \code{\link{lik}}, \code{\link{forge}} and \code{\link{expct}}. +#' +#' The \code{"foreach"} backend uses whatever adapter is registered (e.g. +#' \code{doParallel}, \code{doFuture}). The \code{"mirai"} backend uses +#' \code{mirai} daemons and shares large read-only inputs (training data, +#' forest, learned parameters) across workers via \code{mori}, so workers do +#' not each copy them. Speed is comparable between the backends. The memory +#' benefit is largest for the tree-parallel operations (\code{\link{forde}}, +#' \code{\link{adversarial_rf}}) on large forests with many workers, where +#' \code{foreach} memory grows with the worker count and \code{mirai} stays +#' much flatter (roughly half to a third at 16 workers in internal benchmarks). +#' For small workloads the daemon pool adds a fixed overhead that can outweigh +#' the sharing, so prefer \code{"mirai"} at scale and either backend otherwise. +#' Daemons started via \code{future.mirai} (e.g. +#' \code{plan(future.mirai::mirai_multisession)}) are detected like any other +#' mirai daemons, so futureverse users get the \code{"mirai"} backend +#' automatically. All backend packages are in Suggests; install the ones you +#' use. +#' +#' Reproducibility of stochastic operations (\code{\link{forge}}, categorical +#' \code{\link{expct}}) under parallel execution: \code{\link{set.seed}} only +#' governs the calling process, not the workers. With the \code{"mirai"} +#' backend, seed the daemons instead: \code{mirai::daemons(n, seed = 42)} +#' gives reproducible results provided the daemon count, the seed and the +#' sequence of calls on a fresh daemon pool are kept fixed (changing the +#' daemon count changes how work is chunked and therefore the random stream +#' assignment). This also applies to pools started via \code{future.mirai}: +#' \code{future}'s own seed machinery covers only work dispatched through its +#' API, which arf's backend bypasses, and \code{plan()} does not accept a +#' \code{seed} argument, so seed the pool with \code{mirai::daemons()} +#' directly. For the \code{"foreach"} backend, register \code{doRNG} on top of +#' the adapter (\code{doRNG::registerDoRNG(42)} after +#' \code{registerDoParallel()} or \code{registerDoFuture()}): results are then +#' reproducible, identical across adapters, and independent of the worker +#' count, provided \code{stepsize} is set explicitly (its default depends on +#' the worker count). Without \code{doRNG}, \code{doFuture} flags the +#' stochastic operations with "UNRELIABLE VALUE" warnings. Sequential +#' execution (\code{parallel = FALSE}) with \code{set.seed} is exact as +#' always. +#' +#' @examples +#' \dontrun{ +#' arf <- adversarial_rf(iris) +#' +#' # foreach backend +#' doParallel::registerDoParallel(cores = 4) +#' psi <- forde(arf, iris) +#' +#' # mirai backend: start daemons, then call as usual +#' mirai::daemons(4) +#' psi <- forde(arf, iris) +#' mirai::daemons(0) # shut down when done +#' +#' # futureverse: future.mirai daemons are detected automatically +#' future::plan(future.mirai::mirai_multisession, workers = 4) +#' psi <- forde(arf, iris) +#' +#' # force a backend regardless of what is registered +#' options(arf.backend = "mirai") +#' +#' # silence the backend message +#' options(arf.verbose = FALSE) +#' +#' # reproducible parallel sampling with mirai: seeded daemons +#' # (fixed daemon count, fresh pool) +#' evi <- data.frame(Species = sample(levels(iris$Species), 100, replace = TRUE)) +#' mirai::daemons(4, seed = 42) +#' # stepsize = evidence rows per step; 25 = 100 conditions / 4 workers, +#' # i.e. the default sizing, made explicit +#' x <- forge(psi, n_synth = 1, evidence = evi, stepsize = 25) +#' mirai::daemons(0) +#' +#' # reproducible parallel sampling with foreach: doRNG on top of the +#' # adapter; set stepsize explicitly (its default depends on worker count) +#' doParallel::registerDoParallel(cores = 4) +#' doRNG::registerDoRNG(42) +#' x <- forge(psi, n_synth = 1, evidence = evi, stepsize = 25) +#' } +#' +#' @name arf-options +#' @aliases arf.backend arf.verbose arf.block_rows +NULL diff --git a/R/arf-package.R b/R/arf-package.R index 50933e75..5e76fd76 100644 --- a/R/arf-package.R +++ b/R/arf-package.R @@ -36,6 +36,9 @@ #' # ... or with doFuture #' doFuture::registerDoFuture() #' future::plan("multisession", workers = 4) +#' +#' # ... or with mirai (shares large read-only inputs across workers via mori) +#' mirai::daemons(4) #' } "_PACKAGE" diff --git a/R/cforde_workers.R b/R/cforde_workers.R new file mode 100644 index 00000000..359d3f92 --- /dev/null +++ b/R/cforde_workers.R @@ -0,0 +1,150 @@ +#' @keywords internal +#' @noRd + +# Per-step worker for cforde(), extracted so foreach and mirai backends share one +# definition (params + condition_long mori-shared under mirai). Calls arf/data.table +# internals via bare names, so mirai daemons must have arf loaded. +arf_cforde_step <- function(step_, condition_long, conds_conditioned, + nconds_conditioned, stepsize, cat_cols, cnt_cols, + params, family) { + # To avoid data.table check issues + . <- c_idx <- cvg <- cvg_arf <- cvg_factor <- f_idx <- f_idx_uncond <- i.max <- + i.min <- leaf <- max.x <- max.y <- min.x <- min.y <- mu <- prob <- sigma <- + tree <- V1 <- val <- variable <- min <- max <- NULL + + forest <- params$forest + cat <- params$cat + cnt <- params$cnt + + # Define subset of conditions for step_ + index_start <- conds_conditioned[(step_ - 1) * stepsize + 1] + index_end <- conds_conditioned[min(step_ * stepsize, nconds_conditioned)] + condition_long_step <- condition_long[.(index_start:index_end), nomatch = NULL] + + # Store cat and cnt conditions separately + cat_conds <- condition_long_step[variable %in% cat_cols, c("c_idx", "variable", "val")][, variable := factor(variable)] + cnt_conds <- condition_long_step[variable %in% cnt_cols, c("c_idx", "variable", "min", "max", "val")][, `:=`(variable = factor(variable), + val = as.numeric(val))] + # If cat conditions exist, calculate matching leaves + if (nrow(cat_conds) != 0) { + + # Save leaf indices f_idx cat params in list column grouped by variable and val (value) and merge with cat conditions + cat_relevant <- cat[, .(.(f_idx)), by = .(variable, val)] + setkey(cat_relevant, variable, val) + setkey(cat_conds, variable, val) + cat_relevant <- cat_conds[cat_relevant, on = .(variable, val), nomatch = NULL] + setkey(cat_relevant, c_idx) + + # or-combine different conditions on the same feature + if (uniqueN(cat_relevant, by = c("c_idx", "variable")) != nrow(cat_relevant)) { + cat_relevant <- cat_relevant[, .(.(Reduce(union, V1))), by = .(c_idx, variable)] + } + + # Determine matching leaves for cat conditions + relevant_leaves_changed_cat <- cat_relevant[, Reduce(intersect, V1), by = c_idx][, .(c_idx, f_idx = V1)] + setorder(relevant_leaves_changed_cat) + conditions_unchanged_cat <- setdiff(condition_long_step[, c_idx], cat_conds[, c_idx]) + relevant_leaves_unchanged_cat <- data.table(c_idx = rep(conditions_unchanged_cat, each = nrow(forest)), f_idx = rep(forest[, f_idx], length(conditions_unchanged_cat))) + relevant_leaves_cat <- rbind(relevant_leaves_changed_cat, relevant_leaves_unchanged_cat, fill = TRUE) + relevant_leaves_cat_list <- relevant_leaves_cat[, .(f_idx = .(f_idx)), by = c_idx] + } else { + relevant_leaves_cat <- data.table(c_idx = integer(), f_idx = integer()) + } + + # If cnt conditions exist, calculate matching leaves + if (nrow(cnt_conds) != 0) { + + # Save min, max in cnt params in list columns grouped by variable and merge with cnt conditions + cnt_relevant <- cnt[, .(min = .(min), max = .(max)), by = variable] + cnt_conds_compact <- copy(cnt_conds) + cnt_conds_compact[!is.na(val), `:=`(min = val, max = val)][, val := NULL] + cnt_relevant <- cnt_conds_compact[cnt_relevant, on = .(variable), nomatch = NULL] + setkey(cnt_relevant, c_idx) + + # If cat conds exist, use only matching subset of potentially relevant leaves for cnt conditions + if (nrow(cat_conds) != 0) { + cnt_relevant <- cnt_relevant[relevant_leaves_cat_list, on = .(c_idx)] + } else { + cnt_relevant[, f_idx := NA] + } + + # Determine matching leaves for cnt conditions per row + cnt_relevant <- cnt_relevant[, .( + c_idx, + variable, + f_idx = Map(function(f_idx, min, max, i.min, i.max) { + if (!inherits(f_idx, "logical")) { + rel_cnt_min <- i.min[f_idx] + rel_cnt_max <- i.max[f_idx] + rel_min <- f_idx[which(max > rel_cnt_min)] + rel_max <- f_idx[which(min <= rel_cnt_max)] + } else { + rel_cnt_min <- i.min + rel_cnt_max <- i.max + rel_min <- which(max > rel_cnt_min) + rel_max <- which(min <= rel_cnt_max) + } + intersect(rel_min, rel_max) + }, f_idx = f_idx, min = min, max = max, i.min = i.min, i.max = i.max))] + + # or-combine different conditions on the same feature + or_within_row_cnt <- uniqueN(cnt_relevant, by = c("c_idx", "variable")) != nrow(cnt_relevant) + if (or_within_row_cnt) { + cnt_relevant <- cnt_relevant[, .(f_idx = .(Reduce(union, f_idx))), by = .(c_idx, variable)] + } + + # Determine matching leaves for cnt conditions + relevant_leaves_changed_cnt <- cnt_relevant[, Reduce(intersect, f_idx), by = c_idx][, .(c_idx, f_idx = V1)] + conditions_unchanged_cnt <- setdiff(condition_long_step[, c_idx], cnt_conds[, c_idx]) + relevant_leaves_unchanged_cnt <- data.table(c_idx = rep(conditions_unchanged_cnt, each = nrow(forest)), f_idx = rep(forest[, f_idx], length(conditions_unchanged_cnt))) + relevant_leaves_cnt <- rbind(relevant_leaves_changed_cnt, relevant_leaves_unchanged_cnt) + + # Calculate updates for cnt params matching cnt conditions + cnt_new <- merge(merge(relevant_leaves_cnt, cnt_conds, by = "c_idx", allow.cartesian = TRUE, sort = FALSE), cnt, by = c("f_idx", "variable"), all.x = TRUE, allow.cartesian = TRUE, sort = FALSE) + cnt_new[!is.na(val), `:=`(min = min.y, + max = max.y)] + cnt_new[is.na(val), `:=`(min = pmax(min.x, min.y, na.rm = TRUE), + max = pmin(max.x, max.y, na.rm = TRUE))] + cnt_new <- cnt_new[min <= max & sum(max.x, val, na.rm = TRUE) != min.y, ] + cnt_new[, prob := NA_real_] + if (family == "truncnorm") { + cnt_new[!is.na(val), prob := dtruncnorm(val, a = min.y, b = max.y, mean = mu, sd = sigma) * (val != min.y)] + cnt_new[is.na(val) & (min == min.y) & (max == max.y), prob := 1] + cnt_new[is.na(val) & is.na(prob), prob := ptruncnorm(max, a = min.y, b = max.y, mean = mu, sd = sigma) - ptruncnorm(min, a = min.y, max.y, mean = mu, sd = sigma)] + } else if (family == "unif") { + cnt_new[!is.na(val), prob := dunif(val, min = min.y, max = max.y) * (val != min.y)] + cnt_new[is.na(val) & (min == min.y) & (max == max.y), prob := 1] + cnt_new[is.na(val) & is.na(prob), prob := punif(max, min = min.y, max = max.y) - punif(min, min = min.y, max.y)] + } + cnt_new[, c("min.x", "max.x", "min.y", "max.y") := NULL] + + # If or-combined cnt condition within rows exist, calculate likelihoods for ranges within leaves and norm to 1 + if (or_within_row_cnt) { + cnt_new[, cvg_factor := sum(prob), by = .(f_idx, c_idx, variable)] + cnt_new[, prob := prob / cvg_factor] + } else { + cnt_new[, `:=`(cvg_factor = prob, prob = 1)] + } + + # Calculate final set of matching leaves + if (nrow(cat_conds) > 0) { + relevant_leaves <- merge(relevant_leaves_cnt, relevant_leaves_cat, by = c("c_idx", "f_idx"))[, .(c_idx, f_idx)] + } else { + relevant_leaves <- relevant_leaves_cnt[, .(c_idx, f_idx)] + } + + # If no cnt conditions exist, output empty update table cnt_new for cnt params + } else { + relevant_leaves <- relevant_leaves_cat[, .(c_idx, f_idx)] + cnt_new <- cbind(cnt[FALSE, ], data.table(cvg_factor = numeric(), c_idx = integer(), val = numeric(), prob = numeric())) + } + + # Calculate updates for cat params matching cat conditions + cat_new <- merge(merge(relevant_leaves, cat_conds, by = "c_idx", allow.cartesian = TRUE), cat, by = c("f_idx", "variable", "val")) + + # Ensure probabilities sum to 1 + cat_new[, cvg_factor := sum(prob), by = .(f_idx, c_idx, variable)] + cat_new[, prob := prob / cvg_factor] + + list(cnt_new = cnt_new, cat_new = cat_new, relevant_leaves = relevant_leaves) +} diff --git a/R/expct.R b/R/expct.R index e680531c..970ff9fa 100644 --- a/R/expct.R +++ b/R/expct.R @@ -24,11 +24,17 @@ #' \code{NA} (\code{"na"}). The default is \code{"force"}. #' @param verbose Show warnings, e.g. when no leaf matches a condition? #' @param stepsize How many rows of evidence should be handled at each step? -#' Defaults to \code{nrow(evidence) / num_registered_workers} for +#' Defaults to \code{nrow(evidence)} divided by the number of registered +#' workers or daemons for #' \code{parallel == TRUE}. -#' @param parallel Compute in parallel? Must register backend beforehand, e.g. -#' via \code{doParallel} or \code{doFuture}; see Examples. -#' +#' @param parallel Compute in parallel? Requires a registered \code{foreach} +#' backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +#' daemons. See \code{\link{arf-options}}. With +#' \code{evidence_row_mode = "or"}, parallelization happens inside the +#' conditional circuit computation; in benchmarks this gave little speedup +#' while raising peak memory, so consider \code{parallel = FALSE} for large +#' \code{"or"} queries. +#' #' @details #' This function computes expected values for any subset of features, optionally #' conditioned on some event(s). @@ -90,6 +96,9 @@ #' # ... or with doFuture #' doFuture::registerDoFuture() #' future::plan("multisession", workers = 4) +#' +#' # ... or with mirai (shares large read-only inputs across workers via mori) +#' mirai::daemons(4) #' } #' #' @seealso @@ -120,6 +129,11 @@ expct <- function( variable <- tree <- f_idx <- cvg <- wt <- V1 <- value <- val <- family <- mu <- sigma <- obs <- prob <- f_idx_uncond <- step <- c_idx <- idx <- NA_share <- . <- NULL + + # Defaults so the extracted per-step worker always receives these (set below + # for the conditional cases). + stepsize_cforde <- 0L + parallel_cforde <- FALSE # Prepare evidence and stepsize if (is.null(evidence)) { @@ -128,7 +142,7 @@ expct <- function( evidence <- as.data.table(evidence) if (stepsize == 0) { if (parallel) { - stepsize <- ceiling(nrow(evidence)/foreach::getDoParWorkers()) + stepsize <- ceiling(nrow(evidence)/arf_n_workers()) } else { stepsize <- nrow(evidence) } @@ -164,113 +178,38 @@ expct <- function( } factor_cols <- params$meta[variable %in% query, family == 'multinom'] - # Run in parallel for each step + # Per-step work lives once in arf_expct_step() (expct_workers.R). This closure + # adapts it to foreach's one-argument iteration. par_fun <- function(step_) { - - # Prepare the event space - if (is.null(evidence) || ( ncol(evidence) == 2 && all(colnames(evidence) == c("f_idx", "wt")))) { - cparams <- NULL - } else { - # Call cforde with part of the evidence for this step - index_start <- (step_-1)*stepsize + 1 - index_end <- min(step_*stepsize, nrow(evidence)) - evidence_part <- evidence[index_start:index_end,] - cparams <- cforde(params, evidence_part, evidence_row_mode, nomatch, verbose, - stepsize_cforde, parallel_cforde) - } - - # omega contains the weight (wt) for each leaf (f_idx) for each condition (c_idx) - if (is.null(cparams)) { - if (is.null(evidence)) { - num_trees <- params$forest[, max(tree)] - omega <- params$forest[, .(f_idx, f_idx_uncond = f_idx, cvg)] - omega[, `:=` (c_idx = 1, wt = cvg / num_trees)] - omega[, cvg := NULL] - } else { - omega <- copy(evidence) - omega[, f_idx_uncond := f_idx] - omega[, c_idx := 1] - } - } else { - omega <- cparams$forest[, .(c_idx, f_idx, f_idx_uncond, wt = cvg)] - } - omega <- omega[wt > 0, ] - omega[, idx := .I] - - synth_cnt <- synth_cat <- NULL - # Continuous data - if (any(!factor_cols)) { - if (is.null(cparams) || nrow(cparams$cnt) == 0){ - psi_cond <- data.table() - } else { - psi_cond <- merge(omega, cparams$cnt[variable %in% query, -c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), - sort = FALSE, allow.cartesian = TRUE)[prob > 0,] - # calculate absolute weights for sub-leaf areas (resulting from within-row or-conditions) - if(any(psi_cond[,prob != 1])) { - psi_cond[, wt := wt*prob] - psi_cond[, I := seq_len(.N), by = .(variable, idx)] - } else { - psi_cond[, I := 1] - } - psi_cond[, prob := NULL] - } - psi <- unique(rbind(psi_cond, - merge(omega, params$cnt[variable %in% query, ], by.x = 'f_idx_uncond', by.y = 'f_idx', - sort = FALSE, allow.cartesian = TRUE)[,`:=` (val = NA_real_, I = 1)]), by = c("c_idx", "f_idx", "variable", "I"))[, I := NULL] - psi[NA_share == 1, wt := 0] - cnt <- psi[is.na(val), val := sum(wt * mu)/sum(wt), by = .(c_idx, variable)] - cnt <- unique(cnt[, .(c_idx, variable, val)]) - synth_cnt <- dcast(cnt, c_idx ~ variable, value.var = 'val')[, c_idx := NULL] - } - - - # Categorical data - if (any(factor_cols)) { - if (is.null(cparams) || nrow(cparams$cat) == 0) { - psi <- merge(omega, params$cat[variable %in% query, ], by.x = 'f_idx_uncond', by.y = 'f_idx', sort = FALSE, allow.cartesian = TRUE) - } else { - psi_cond <- merge(omega, cparams$cat[variable %in% query, -c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), - sort = FALSE, allow.cartesian = TRUE) - psi_uncond <- merge(omega, params$cat[variable %in% query, ], by.x = 'f_idx_uncond', by.y = 'f_idx', - sort = FALSE, allow.cartesian = TRUE) - psi_uncond_relevant <- psi_uncond[!psi_cond, on = .(idx, variable)] - psi <- rbind(psi_cond, psi_uncond_relevant) - } - psi[NA_share == 1, wt := 0] - cat <- psi[, sum(wt * prob), by = .(c_idx, variable, val)] - cat <- setDT(cat)[, .SD[which.max.random(V1)], by = .(c_idx, variable)] - synth_cat <- dcast(cat, c_idx ~ variable, value.var = 'val')[, c_idx := NULL] - } - - # Create dataset with expectations - x_synth <- cbind(synth_cnt, synth_cat) - x_synth <- post_x(x_synth, params, round) - - if (evidence_row_mode == "separate" & any(omega[, is.na(f_idx)])) { - setDT(x_synth) - indices_na <- cparams$forest[is.na(f_idx), c_idx] - indices_sampled <- cparams$forest[!is.na(f_idx), unique(c_idx)] - rows_na <- dcast(rbind(data.table(c_idx = 0, variable = params$meta[,variable]), - cparams$evidence_prepped[c_idx %in% indices_na,], - fill = TRUE), - c_idx ~ variable, value.var = "val")[c_idx != 0,] - if (nomatch == "force") { - rows_na_sampled <- expct(params, parallel = parallel, stepsize = stepsize) - rows_na[is.na(rows_na)] <- rows_na_sampled[is.na(rows_na[,-1])] - } - x_synth[, c_idx := indices_sampled] - x_synth <- rbind(x_synth, rows_na, fill = TRUE) - setorder(x_synth, c_idx)[, c_idx := NULL] - x_synth <- post_x(x_synth, params, round) - } - - x_synth - } - if (isTRUE(parallel)) { + arf_expct_step(step_, params, evidence, query, factor_cols, + evidence_row_mode, nomatch, verbose, round, stepsize, + stepsize_cforde, parallel_cforde) + } + # Parallelism is across steps: 1 step is inherently serial for any backend. + # mirai only for step_no > 1; "or" mode already set parallel <- FALSE. + use_mirai <- FALSE + if (step_no > 1) { + backend <- arf_select_backend(parallel) + use_mirai <- identical(backend, "mirai") + } + if (use_mirai) { + arf_load_on_daemons() # daemons need arf: worker calls cforde/post_x/which.max.random + params_shared <- mori::share(params) + evidence_shared <- if (!is.null(evidence)) mori::share(evidence) else NULL + x_synth_ <- arf_mirai_tree_map(step_no, arf_expct_step, list( + params = params_shared, evidence = evidence_shared, query = query, + factor_cols = factor_cols, evidence_row_mode = evidence_row_mode, + nomatch = nomatch, verbose = verbose, round = round, stepsize = stepsize, + stepsize_cforde = stepsize_cforde, parallel_cforde = parallel_cforde), + # package-level combine: an inline closure here would serialize expct's + # whole frame (params included) into every task + # See the closure note in mirai_helpers.R. + combine = arf_rbind_steps) + } else if (isTRUE(parallel) && step_no > 1) { x_synth_ <- foreach(step = 1:step_no, .combine = "rbind") %dopar% par_fun(step) } else { x_synth_ <- foreach(step = 1:step_no, .combine = "rbind") %do% par_fun(step) } return(x_synth_) -} \ No newline at end of file +} diff --git a/R/expct_workers.R b/R/expct_workers.R new file mode 100644 index 00000000..c334a4f9 --- /dev/null +++ b/R/expct_workers.R @@ -0,0 +1,146 @@ +#' @keywords internal +#' @noRd + +# Per-step worker for expct(), extracted so foreach and mirai backends share one +# definition (params/evidence mori-shared under mirai). Calls arf internals +# (cforde, post_x, which.max.random) so mirai daemons must have arf loaded. +arf_expct_step <- function(step_, params, evidence, query, factor_cols, + evidence_row_mode, nomatch, verbose, round, + stepsize, stepsize_cforde, parallel_cforde) { + # To avoid data.table check issues + variable <- tree <- f_idx <- cvg <- wt <- V1 <- value <- val <- family <- + mu <- sigma <- obs <- prob <- f_idx_uncond <- c_idx <- idx <- NA_share <- + . <- I <- NULL + + # Prepare the event space + if (is.null(evidence) || (ncol(evidence) == 2 && all(colnames(evidence) == c("f_idx", "wt")))) { + cparams <- NULL + } else { + # Call cforde with part of the evidence for this step + index_start <- (step_ - 1) * stepsize + 1 + index_end <- min(step_ * stepsize, nrow(evidence)) + evidence_part <- evidence[index_start:index_end, ] + cparams <- cforde(params, evidence_part, evidence_row_mode, nomatch, verbose, + stepsize_cforde, parallel_cforde) + } + + # omega contains the weight (wt) for each leaf (f_idx) for each condition (c_idx) + if (is.null(cparams)) { + if (is.null(evidence)) { + num_trees <- params$forest[, max(tree)] + omega <- params$forest[, .(f_idx, f_idx_uncond = f_idx, cvg)] + omega[, `:=`(c_idx = 1, wt = cvg / num_trees)] + omega[, cvg := NULL] + } else { + omega <- copy(evidence) + omega[, f_idx_uncond := f_idx] + omega[, c_idx := 1] + } + } else { + omega <- cparams$forest[, .(c_idx, f_idx, f_idx_uncond, wt = cvg)] + } + omega <- omega[wt > 0, ] + omega[, idx := .I] + + # Synthesize expectations for one block of conditions (subset of omega + # rows; cparams/params merges below self-restrict via the c_idx join keys). + synth_block <- function(omega_) { + synth_cnt <- synth_cat <- NULL + # Continuous data + if (any(!factor_cols)) { + if (is.null(cparams) || nrow(cparams$cnt) == 0) { + psi_cond <- data.table() + } else { + psi_cond <- merge(omega_, cparams$cnt[variable %in% query, -c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), + sort = FALSE, allow.cartesian = TRUE)[prob > 0, ] + # calculate absolute weights for sub-leaf areas (resulting from within-row or-conditions) + if (any(psi_cond[, prob != 1])) { + psi_cond[, wt := wt * prob] + psi_cond[, I := seq_len(.N), by = .(variable, idx)] + } else { + psi_cond[, I := 1] + } + psi_cond[, prob := NULL] + } + psi <- unique(rbind(psi_cond, + merge(omega_, params$cnt[variable %in% query, ], by.x = 'f_idx_uncond', by.y = 'f_idx', + sort = FALSE, allow.cartesian = TRUE)[, `:=`(val = NA_real_, I = 1)]), by = c("c_idx", "f_idx", "variable", "I"))[, I := NULL] + psi[NA_share == 1, wt := 0] + cnt <- psi[is.na(val), val := sum(wt * mu) / sum(wt), by = .(c_idx, variable)] + cnt <- unique(cnt[, .(c_idx, variable, val)]) + synth_cnt <- dcast(cnt, c_idx ~ variable, value.var = 'val')[, c_idx := NULL] + } + + # Categorical data + if (any(factor_cols)) { + if (is.null(cparams) || nrow(cparams$cat) == 0) { + psi <- merge(omega_, params$cat[variable %in% query, ], by.x = 'f_idx_uncond', by.y = 'f_idx', sort = FALSE, allow.cartesian = TRUE) + } else { + psi_cond <- merge(omega_, cparams$cat[variable %in% query, -c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), + sort = FALSE, allow.cartesian = TRUE) + psi_uncond <- merge(omega_, params$cat[variable %in% query, ], by.x = 'f_idx_uncond', by.y = 'f_idx', + sort = FALSE, allow.cartesian = TRUE) + psi_uncond_relevant <- psi_uncond[!psi_cond, on = .(idx, variable)] + psi <- rbind(psi_cond, psi_uncond_relevant) + } + psi[NA_share == 1, wt := 0] + cat <- psi[, sum(wt * prob), by = .(c_idx, variable, val)] + cat <- setDT(cat)[, .SD[which.max.random(V1)], by = .(c_idx, variable)] + synth_cat <- dcast(cat, c_idx ~ variable, value.var = 'val')[, c_idx := NULL] + } + cbind(synth_cnt, synth_cat) + } + + # The merges in synth_block materialize (#matched leaves x #query variables) + # rows PER CONDITION -- all at once for the step, three times over via + # merge/rbind/unique, on every backend alike; at large forest x condition + # counts this reaches tens of GB. Conditions are independent here + # (every aggregation groups by c_idx, omega arrives sorted by c_idx, and the + # per-group RNG order of which.max.random is ascending c_idx either way), so + # when the estimated join size exceeds block_cap rows, process conditions in + # blocks that keep each materialization bounded. A single condition cannot be + # split; its leaves x variables product is the floor of this algorithm. + # Lower via options(arf.block_rows) for tight-memory runs; see ?arf-options. + block_cap <- max(1, as.numeric(getOption("arf.block_rows", 5e6))) + n_vars <- max(1L, length(query)) + # as.double: nrow * n_vars overflows integer arithmetic exactly in the + # large-forest regime the blocking exists for + if (as.double(nrow(omega)) * n_vars <= block_cap || omega[, uniqueN(c_idx)] == 1L) { + x_synth <- synth_block(omega) + } else { + sizes <- omega[, .N, by = c_idx] # ascending c_idx + g <- integer(nrow(sizes)); gi <- 1L; acc <- 0 + for (i in seq_len(nrow(sizes))) { + r <- sizes$N[i] * n_vars + if (acc > 0 && acc + r > block_cap) { gi <- gi + 1L; acc <- 0 } + g[i] <- gi; acc <- acc + r + } + x_synth <- rbindlist(lapply(split(sizes$c_idx, g), function(cs) { + synth_block(omega[c_idx %in% cs]) + })) + } + + # Create dataset with expectations + x_synth <- post_x(x_synth, params, round) + + if (evidence_row_mode == "separate" & any(omega[, is.na(f_idx)])) { + setDT(x_synth) + indices_na <- cparams$forest[is.na(f_idx), c_idx] + indices_sampled <- cparams$forest[!is.na(f_idx), unique(c_idx)] + rows_na <- dcast(rbind(data.table(c_idx = 0, variable = params$meta[, variable]), + cparams$evidence_prepped[c_idx %in% indices_na, ], + fill = TRUE), + c_idx ~ variable, value.var = "val")[c_idx != 0, ] + if (nomatch == "force") { + # nested recovery runs serial (a worker must not spawn its own backend) + rows_na_sampled <- expct(params, parallel = FALSE) + rows_na[is.na(rows_na)] <- rows_na_sampled[is.na(rows_na[, -1])] + } + x_synth[, c_idx := indices_sampled] + x_synth <- rbind(x_synth, rows_na, fill = TRUE) + setorder(x_synth, c_idx)[, c_idx := NULL] + x_synth <- post_x(x_synth, params, round) + } + + x_synth +} diff --git a/R/forde.R b/R/forde.R index 66ab14b5..16573c67 100644 --- a/R/forde.R +++ b/R/forde.R @@ -25,8 +25,9 @@ #' \code{finite_bounds != "no"}. This avoids zero-density points when test #' data fall outside the support of training data. The gap between lower and #' upper bounds is expanded by a factor of \code{1 + epsilon}. -#' @param parallel Compute in parallel? Must register backend beforehand, e.g. -#' via \code{doParallel} or \code{doFuture}; see examples. +#' @param parallel Compute in parallel? Requires a registered \code{foreach} +#' backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +#' daemons. See \code{\link{arf-options}}. #' #' #' @details @@ -96,6 +97,9 @@ #' # ... or with doFuture #' doFuture::registerDoFuture() #' future::plan("multisession", workers = 4) +#' +#' # ... or with mirai (shares large read-only inputs across workers via mori) +#' mirai::daemons(4) #' } #' #' @@ -110,7 +114,7 @@ #' @importFrom stats predict runif #' @importFrom foreach foreach %do% %dopar% #' - + forde <- function( arf, x, @@ -191,139 +195,80 @@ forde <- function( # Compute leaf bounds and coverage num_trees <- arf$num.trees + # Pick the parallel backend (mirai daemons vs a registered foreach backend) + # and report the choice. Only relevant when parallel = TRUE; see + # arf_select_backend() in mirai_helpers.R. + backend <- arf_select_backend(parallel) + use_mirai <- identical(backend, "mirai") + if (use_mirai) { + # Load arf (and via Imports, data.table with its S3 methods) on the + # daemons; cached once per pool, see arf_load_on_daemons(). + arf_load_on_daemons() + } + # Leaf assignments first: the per-tree bounds worker computes its tree's + # coverage from its pred column. predict() needs prep_x's column names + # (the forest was trained on them). + pred <- stats::predict(arf, x, type = 'terminalNodes')$predictions + 1L + # Restore original column names before x is shared with workers: the psi + # workers match x's melted column names against bnds$variable, which uses + # colnames_x (the bounds worker accesses x by index, never by name). + setnames(x, colnames_x) + # Per-tree workers live once in forde_workers.R, shared by all backends; the + # closures below adapt them to foreach's one-argument iteration. bnd_fn <- function(tree) { - num_nodes <- length(arf$forest$split.varIDs[[tree]]) - lb <- matrix(-Inf, nrow = num_nodes, ncol = d) - ub <- matrix(Inf, nrow = num_nodes, ncol = d) - if (finite_bounds == 'global' & any(!factor_cols)) { - for (j in which(!factor_cols)) { - min_j <- min(x[[j]], na.rm = TRUE) - max_j <- max(x[[j]], na.rm = TRUE) - gap <- max_j - min_j - lb[, j] <- min_j - epsilon / 2 * gap - ub[, j] <- max_j + epsilon / 2 * gap - } - } - for (i in 1:num_nodes) { - left_child <- arf$forest$child.nodeIDs[[tree]][[1]][i] + 1L - right_child <- arf$forest$child.nodeIDs[[tree]][[2]][i] + 1L - splitvarID <- arf$forest$split.varIDs[[tree]][i] + 1L - splitval <- arf$forest$split.values[[tree]][i] - if (left_child > 1) { - ub[left_child, ] <- ub[right_child, ] <- ub[i, ] - lb[left_child, ] <- lb[right_child, ] <- lb[i, ] - if (left_child != right_child) { - # If no pruned node, split changes bounds - ub[left_child, splitvarID] <- lb[right_child, splitvarID] <- splitval - } - } - } - leaves <- which(arf$forest$child.nodeIDs[[tree]][[1]] == 0L) - colnames(lb) <- colnames(ub) <- colnames_x - merge(melt(data.table(tree = tree, leaf = leaves, lb[leaves, , drop = FALSE]), - id.vars = c('tree', 'leaf'), value.name = 'min'), - melt(data.table(tree = tree, leaf = leaves, ub[leaves, , drop = FALSE]), - id.vars = c('tree', 'leaf'), value.name = 'max'), - by = c('tree', 'leaf', 'variable'), sort = FALSE) + arf_bnd_fn(tree, arf$forest, d, finite_bounds, factor_cols, x, epsilon, + colnames_x, pred, arf$inbag.counts, n, oob) } - if (isTRUE(parallel)) { + if (use_mirai) { + forest_slice <- mori::share( + arf$forest[c("split.varIDs", "child.nodeIDs", "split.values")]) + x_shared <- mori::share(x) + pred_shared <- mori::share(pred) + inbag_shared <- if (!is.null(arf$inbag.counts)) { + mori::share(arf$inbag.counts) + } else NULL + bnds <- arf_mirai_tree_map(num_trees, arf_bnd_fn, list( + forest = forest_slice, d = d, finite_bounds = finite_bounds, + factor_cols = factor_cols, x = x_shared, epsilon = epsilon, + colnames_x = colnames_x, pred = pred_shared, + inbag.counts = inbag_shared, n = n, oob = oob)) + } else if (isTRUE(parallel)) { bnds <- foreach(tree = seq_len(num_trees), .combine = rbind) %dopar% bnd_fn(tree) } else { bnds <- foreach(tree = seq_len(num_trees), .combine = rbind) %do% bnd_fn(tree) } - # Compute coverage - pred <- stats::predict(arf, x, type = 'terminalNodes')$predictions + 1L - keep <- data.table('tree' = rep(seq_len(num_trees), each = n), - 'leaf' = as.vector(pred)) - if (isTRUE(oob)) { - keep[, oob := as.vector(sapply(seq_len(num_trees), function(b) { - arf$inbag.counts[[b]][seq_len(n)] == 0L - }))] - keep <- keep[oob == TRUE] - keep <- unique(keep[, cnt := .N, by = .(tree, leaf)]) - keep[, n_oob := sum(oob), by = tree] - keep[, cvg := cnt / n_oob][, c('oob', 'cnt', 'n_oob') := NULL] - keep[, cvg := cvg/sum(cvg), by = tree] - } else if (oob == "inbag") { - keep[, inbag := as.vector(sapply(seq_len(num_trees), function(b) { - arf$inbag.counts[[b]][seq_len(n)] > 0L - }))] - keep <- keep[inbag == TRUE] - keep <- unique(keep[, cnt := .N, by = .(tree, leaf)]) - keep[, n_inbag := sum(inbag), by = tree] - keep[, cvg := cnt / n_inbag][, c('inbag', 'cnt', 'n_inbag') := NULL] - keep[, cvg := cvg/sum(cvg), by = tree] - } else { - keep <- unique(keep[, cnt := .N, by = .(tree, leaf)]) - keep[, cvg := cnt / n][, cnt := NULL] - } - bnds <- merge(bnds, keep, by = c('tree', 'leaf'), sort = FALSE) - rm(keep) # Create forest index setkey(bnds, tree, leaf) bnds[, f_idx := .GRP, by = key(bnds)] - # Calculate distribution parameters for each variable - setnames(x, colnames_x) + # Calculate distribution parameters for each variable: one fused dispatch + # per tree computes continuous and categorical params together (they share + # every input). + psi_fn <- function(tree) { + arf_psi_fn(tree, arf_psi_cnt_fn, arf_psi_cat_fn, x, factor_cols, pred, + arf$inbag.counts, n, oob, bnds, finite_bounds, epsilon, family, + lvl_df_rf, alpha) + } + if (use_mirai) { + bnds_shared <- mori::share(bnds) + psi_pair <- arf_mirai_tree_map(num_trees, arf_psi_fn, list( + cnt_fn = arf_psi_cnt_fn, cat_fn = arf_psi_cat_fn, + x = x_shared, factor_cols = factor_cols, pred = pred_shared, + inbag.counts = inbag_shared, n = n, oob = oob, + bnds = bnds_shared, finite_bounds = finite_bounds, + epsilon = epsilon, family = family, + lvl_df_rf = mori::share(lvl_df_rf), alpha = alpha), + combine = arf_combine_psi) + } else if (isTRUE(parallel)) { + psi_pair <- arf_combine_psi( + foreach(tree = seq_len(num_trees)) %dopar% psi_fn(tree)) + } else { + psi_pair <- arf_combine_psi( + foreach(tree = seq_len(num_trees)) %do% psi_fn(tree)) + } # Continuous case if (any(!factor_cols)) { - psi_cnt_fn <- function(tree) { - dt <- data.table(x[, !factor_cols, drop = FALSE], leaf = pred[, tree]) - if (isTRUE(oob)) { - dt <- dt[arf$inbag.counts[[tree]][1:n] == 0L, ] - dt <- dt[!is.na(leaf)] - } else if (oob == "inbag") { - dt <- dt[arf$inbag.counts[[tree]][1:n] > 0L, ] - dt <- dt[!is.na(leaf)] - } - dt <- melt(dt, id.vars = 'leaf', variable.factor = FALSE)[, tree := tree] - dt <- merge(dt, bnds[, .(tree, leaf, variable, min, max, f_idx)], - by = c('tree', 'leaf', 'variable'), sort = FALSE) - # Caculate bounds for finite_bounds == 'local' - if (finite_bounds == 'local') { - dt[, c('min_emp', 'max_emp') := .(min(value, na.rm = TRUE), max(value, na.rm = TRUE)), by = .(leaf, variable)] - dt[, length_emp := max_emp - min_emp] - # Calculate bounds if min_emp == max_emp in order to be able to sample from cont. distribution - length_emp_0_replace <- min(dt[length_emp > 0, min(length_emp, na.rm = TRUE)], max(epsilon, 1e-12)) - dt[length_emp == 0, c('min_emp', 'max_emp', 'length_emp') := .(min_emp - length_emp_0_replace/2, max_emp + length_emp_0_replace/2, length_emp_0_replace)] - dt[, c('min', 'max', 'min_emp', 'max_emp', 'length_emp') := .(fifelse(!is.finite(min) & !is.na(min_emp), min_emp - length_emp*(epsilon/2), min), - fifelse(!is.finite(max) & !is.na(max_emp), max_emp + length_emp*(epsilon/2), max), - NULL, NULL, NULL)] - } - if (family == 'truncnorm') { - dt[, c('mu', 'sigma', 'NA_share') := .(mean(value, na.rm = TRUE), sd(value, na.rm = TRUE), sum(is.na(value))/.N), - by = .(leaf, variable)] - dt[, c('min_emp', 'max_emp') := .(min(value, na.rm = TRUE), max(value, na.rm = TRUE)), by = variable] - dt[NA_share == 1, c('min', 'max') := .(fifelse(is.infinite(min), min_emp, min), - fifelse(is.infinite(max), max_emp, max))] - dt[, c("min_emp", "max_emp") := NULL] - dt[NA_share == 1, mu := (max + min) / 2] - dt[is.na(sigma), sigma := 0] - if (any(dt[, sigma == 0])) { - dt[, new_min := fifelse(!is.finite(min), min(value, na.rm = TRUE), min), by = variable] - dt[, new_max := fifelse(!is.finite(max), max(value, na.rm = TRUE), max), by = variable] - dt[, mid := (new_min + new_max) / 2] - dt[, sigma0 := (new_max - mid) / stats::qnorm(0.975)] - # This prior places 95% of the density within the bounding box. - # In addition, we set the prior degrees of freedom at nu0 = 2. - # Since the mode of a chisq is max(df-2, 0), this means that - # (1) with a single observation, the posterior reduces to the prior; and - # (2) with more invariant observations, the posterior tends toward zero. - dt[sigma == 0, sigma := sqrt(2 / .N * sigma0^2), by = .(variable, leaf)] - dt[, c('new_min', 'new_max', 'mid', 'sigma0') := NULL] - } - } else if (family == 'unif') { - dt[, NA_share := sum(is.na(value))/.N, by = .(leaf, variable)] - } - return(unique(dt[, c('tree', 'leaf', 'value') := NULL])) - } - if (isTRUE(parallel)) { - psi_cnt <- foreach(tree = seq_len(num_trees), .combine = rbind) %dopar% - psi_cnt_fn(tree) - } else { - psi_cnt <- foreach(tree = seq_len(num_trees), .combine = rbind) %do% - psi_cnt_fn(tree) - } + psi_cnt <- psi_pair$cnt setkey(psi_cnt, f_idx, variable) setcolorder(psi_cnt, c('f_idx', 'variable')) } else { @@ -333,75 +278,7 @@ forde <- function( # Categorical case if (any(factor_cols)) { - psi_cat_fn <- function(tree) { - dt <- data.table(x[, factor_cols, drop = FALSE], leaf = pred[, tree]) - if (isTRUE(oob)) { - dt <- dt[!is.na(leaf)] - } - dt <- melt(dt, id.vars = 'leaf', variable.factor = FALSE, - value.factor = FALSE, value.name = 'val')[, tree := tree] - dt[, NA_share := sum(is.na(val))/.N, by = .(leaf, variable)] - dt <- dt[!(is.na(val) & NA_share != 1)] - if (dt[, any(NA_share == 1)]) { - # Handle leaves where all values for a categorical variable are NA - all_na <- unique(dt[NA_share == 1, ]) - dt <- dt[NA_share != 1, ] - all_na <- merge(all_na, bnds[, .(tree, leaf, variable, min, max, f_idx)], - by = c('tree', 'leaf', 'variable'), sort = FALSE) - all_na[!is.finite(min), min := 0.5] - for (j in names(which(factor_cols))) { - all_na[!is.finite(max) & variable == j, max := lvl_df_rf[variable == j, max(level)]] - } - all_na[!grepl('\\.5', min), min := min + 0.5] - all_na[!grepl('\\.5', max), max := max + 0.5] - all_na[, min := min + 0.5][, max := max - 0.5] - all_na <- all_na[, .(level = seq(min, max), NA_share), by = .(leaf, variable)] - all_na <- merge(all_na, lvl_df_rf, by = c('variable', 'level')) - all_na[, level := NULL][, tree := tree] - setcolorder(all_na, colnames(dt)) - dt <- rbind(dt, all_na) - } - dt[, count := .N, by = .(leaf, variable)] - dt <- merge(dt, bnds[, .(tree, leaf, variable, min, max, f_idx)], - by = c('tree', 'leaf', 'variable'), sort = FALSE) - dt[, c('tree', 'leaf') := NULL] - if (alpha == 0) { - dt <- unique(dt[, prob := .N / count, by = .(f_idx, variable, val)]) - } else { - # Define the range of each variable in each leaf - dt <- unique(dt[, val_count := .N, by = .(f_idx, variable, val)]) - dt <- merge(dt, lvl_df_rf[, .(k = .N), by = variable], by = "variable") - dt[!is.finite(min), min := 0.5][!is.finite(max), max := k + 0.5] - dt[!grepl('\\.5', min), min := min + 0.5][!grepl('\\.5', max), max := max + 0.5] - dt[, k := max - min] - # Enumerate each possible leaf-variable-value combo - tmp <- dt[, seq(min[1] + 0.5, max[1] - 0.5), by = .(f_idx, variable)] - setnames(tmp, 'V1', 'level') - tmp <- merge(tmp, lvl_df_rf, by = c('variable', 'level'), - sort = FALSE)[, level := NULL] - # Populate count, k - tmp <- merge(tmp, unique(dt[, .(f_idx, variable, count, k)]), - by = c('f_idx', 'variable'), sort = FALSE) - # Merge with dt, set val_count = 0 for possible but unobserved levels - dt <- merge(tmp, dt, by = c('f_idx', 'variable', 'val', 'count', 'k'), - all.x = TRUE, sort = FALSE) - dt[is.na(val_count), val_count := 0] - dt[, NA_share := mean(NA_share, na.rm = TRUE), by = .(f_idx, variable)] - # Compute posterior probabilities - dt[, prob := (val_count + alpha) / (count + alpha * k), by = .(f_idx, variable, val)] - dt[, c('val_count', 'k') := NULL] - } - dt[, c('count', 'min', 'max') := NULL] - setcolorder(dt, c("f_idx", "variable", "val", "prob", "NA_share")) - dt - } - if (isTRUE(parallel)) { - psi_cat <- foreach(tree = seq_len(num_trees), .combine = rbind) %dopar% - psi_cat_fn(tree) - } else { - psi_cat <- foreach(tree = seq_len(num_trees), .combine = rbind) %do% - psi_cat_fn(tree) - } + psi_cat <- psi_pair$cat lvl_df_rf[, level := NULL] setkey(psi_cat, f_idx, variable) setcolorder(psi_cat, c('f_idx', 'variable')) diff --git a/R/forde_workers.R b/R/forde_workers.R new file mode 100644 index 00000000..7023789d --- /dev/null +++ b/R/forde_workers.R @@ -0,0 +1,248 @@ +#' @keywords internal +#' @noRd + +# Per-tree worker bodies extracted from forde() so foreach and mirai +# backends share one definition. All dependencies are explicit +# arguments; nothing captured from enclosing scope. This is what makes +# the mirai backend possible — daemons run in clean environments. + +# Worker 1: compute leaf bounds AND coverage for one tree. Coverage lives here +# (not on the calling process) because it is per-tree: tabulating this tree's +# pred column and joining locally keeps the n x num_trees leaf-assignment +# table off the caller, which dominated its memory at large n * trees. +arf_bnd_fn <- function(tree, forest, d, finite_bounds, factor_cols, + x, epsilon, colnames_x, pred, inbag.counts, n, oob) { + # To avoid data.table check issues + variable <- leaf <- cvg <- cnt <- NULL + + num_nodes <- length(forest$split.varIDs[[tree]]) + lb <- matrix(-Inf, nrow = num_nodes, ncol = d) + ub <- matrix(Inf, nrow = num_nodes, ncol = d) + if (finite_bounds == "global" && any(!factor_cols)) { + for (j in which(!factor_cols)) { + min_j <- min(x[[j]], na.rm = TRUE) + max_j <- max(x[[j]], na.rm = TRUE) + gap <- max_j - min_j + lb[, j] <- min_j - epsilon / 2 * gap + ub[, j] <- max_j + epsilon / 2 * gap + } + } + for (i in seq_len(num_nodes)) { + left_child <- forest$child.nodeIDs[[tree]][[1]][i] + 1L + right_child <- forest$child.nodeIDs[[tree]][[2]][i] + 1L + splitvarID <- forest$split.varIDs[[tree]][i] + 1L + splitval <- forest$split.values[[tree]][i] + if (left_child > 1) { + ub[left_child, ] <- ub[right_child, ] <- ub[i, ] + lb[left_child, ] <- lb[right_child, ] <- lb[i, ] + if (left_child != right_child) { + # If no pruned node, split changes bounds + ub[left_child, splitvarID] <- lb[right_child, splitvarID] <- splitval + } + } + } + leaves <- which(forest$child.nodeIDs[[tree]][[1]] == 0L) + colnames(lb) <- colnames(ub) <- colnames_x + bnd <- merge(melt(data.table(tree = tree, leaf = leaves, + lb[leaves, , drop = FALSE]), + id.vars = c("tree", "leaf"), value.name = "min"), + melt(data.table(tree = tree, leaf = leaves, + ub[leaves, , drop = FALSE]), + id.vars = c("tree", "leaf"), value.name = "max"), + by = c("tree", "leaf", "variable"), sort = FALSE) + + # Coverage of this tree's observed leaves; unobserved leaves drop out via + # the inner join. Replicates the original op sequence exactly (including + # the cvg/sum(cvg) renormalization) so results stay bit-identical. + leaf_pred <- pred[, tree] + if (isTRUE(oob)) { + leaf_pred <- leaf_pred[inbag.counts[[tree]][seq_len(n)] == 0L] + } else if (identical(oob, "inbag")) { + leaf_pred <- leaf_pred[inbag.counts[[tree]][seq_len(n)] > 0L] + } + keep <- data.table(tree = tree, leaf = leaf_pred) + keep <- unique(keep[, cnt := .N, by = leaf]) + if (isTRUE(oob) || identical(oob, "inbag")) { + keep[, cvg := cnt / .N][, cnt := NULL] + keep[, cvg := cvg / sum(cvg)] + } else { + keep[, cvg := cnt / n][, cnt := NULL] + } + merge(bnd, keep, by = c("tree", "leaf"), sort = FALSE) +} + +# Workers 2+3 fused: one dispatch per tree computes both continuous and +# categorical params (they share every input: x, pred, bnds). +# Combined across trees by arf_combine_psi() (mirai_helpers.R). +# The sub-workers arrive as ARGUMENTS (cnt_fn/cat_fn), not by name: a bare +# `arf_psi_cnt_fn(...)` call would need the arf namespace resolved on the +# daemon, and deserialization loads the INSTALLED arf there -- stale relative +# to a dev tree. Passing the functions keeps forde's workers fully +# self-contained. +arf_psi_fn <- function(tree, cnt_fn, cat_fn, x, factor_cols, pred, + inbag.counts, n, oob, bnds, finite_bounds, epsilon, + family, lvl_df_rf, alpha) { + list( + cnt = if (any(!factor_cols)) { + cnt_fn(tree, x, factor_cols, pred, inbag.counts, n, oob, bnds, + finite_bounds, epsilon, family) + }, + cat = if (any(factor_cols)) { + cat_fn(tree, x, factor_cols, pred, bnds, oob, lvl_df_rf, alpha) + } + ) +} + +# Worker 2: compute continuous-variable distribution params for one tree. +arf_psi_cnt_fn <- function(tree, x, factor_cols, pred, inbag.counts, n, + oob, bnds, finite_bounds, epsilon, family) { + # To avoid data.table check issues + leaf <- variable <- value <- min_emp <- max_emp <- length_emp <- mu <- + sigma <- NA_share <- new_min <- new_max <- mid <- sigma0 <- f_idx <- . <- + NULL + + dt <- data.table(x[, !factor_cols, drop = FALSE], leaf = pred[, tree]) + if (isTRUE(oob)) { + dt <- dt[inbag.counts[[tree]][1:n] == 0L, ] + dt <- dt[!is.na(leaf)] + } else if (identical(oob, "inbag")) { + dt <- dt[inbag.counts[[tree]][1:n] > 0L, ] + dt <- dt[!is.na(leaf)] + } + dt <- melt(dt, id.vars = "leaf", variable.factor = FALSE)[, tree := tree] + dt <- merge(dt, bnds[, .(tree, leaf, variable, min, max, f_idx)], + by = c("tree", "leaf", "variable"), sort = FALSE) + # Caculate bounds for finite_bounds == 'local' + if (finite_bounds == "local") { + dt[, c("min_emp", "max_emp") := .(min(value, na.rm = TRUE), + max(value, na.rm = TRUE)), + by = .(leaf, variable)] + dt[, length_emp := max_emp - min_emp] + # Calculate bounds if min_emp == max_emp in order to be able to sample from cont. distribution + length_emp_0_replace <- min( + dt[length_emp > 0, min(length_emp, na.rm = TRUE)], + max(epsilon, 1e-12)) + dt[length_emp == 0, + c("min_emp", "max_emp", "length_emp") := + .(min_emp - length_emp_0_replace / 2, + max_emp + length_emp_0_replace / 2, + length_emp_0_replace)] + dt[, c("min", "max", "min_emp", "max_emp", "length_emp") := + .(fifelse(!is.finite(min) & !is.na(min_emp), + min_emp - length_emp * (epsilon / 2), min), + fifelse(!is.finite(max) & !is.na(max_emp), + max_emp + length_emp * (epsilon / 2), max), + NULL, NULL, NULL)] + } + if (family == "truncnorm") { + dt[, c("mu", "sigma", "NA_share") := + .(mean(value, na.rm = TRUE), + stats::sd(value, na.rm = TRUE), + sum(is.na(value)) / .N), + by = .(leaf, variable)] + dt[, c("min_emp", "max_emp") := + .(min(value, na.rm = TRUE), max(value, na.rm = TRUE)), + by = variable] + dt[NA_share == 1, + c("min", "max") := + .(fifelse(is.infinite(min), min_emp, min), + fifelse(is.infinite(max), max_emp, max))] + dt[, c("min_emp", "max_emp") := NULL] + dt[NA_share == 1, mu := (max + min) / 2] + dt[is.na(sigma), sigma := 0] + if (any(dt[, sigma == 0])) { + dt[, new_min := fifelse(!is.finite(min), min(value, na.rm = TRUE), min), + by = variable] + dt[, new_max := fifelse(!is.finite(max), max(value, na.rm = TRUE), max), + by = variable] + dt[, mid := (new_min + new_max) / 2] + dt[, sigma0 := (new_max - mid) / stats::qnorm(0.975)] + # This prior places 95% of the density within the bounding box. + # In addition, we set the prior degrees of freedom at nu0 = 2. + # Since the mode of a chisq is max(df-2, 0), this means that + # (1) with a single observation, the posterior reduces to the prior; and + # (2) with more invariant observations, the posterior tends toward zero. + dt[sigma == 0, sigma := sqrt(2 / .N * sigma0^2), + by = .(variable, leaf)] + dt[, c("new_min", "new_max", "mid", "sigma0") := NULL] + } + } else if (family == "unif") { + dt[, NA_share := sum(is.na(value)) / .N, by = .(leaf, variable)] + } + unique(dt[, c("tree", "leaf", "value") := NULL]) +} + +# Worker 3: compute categorical-variable distribution params for one tree. +arf_psi_cat_fn <- function(tree, x, factor_cols, pred, bnds, oob, + lvl_df_rf, alpha) { + # To avoid data.table check issues + leaf <- variable <- val <- NA_share <- count <- val_count <- k <- + level <- prob <- f_idx <- . <- NULL + + dt <- data.table(x[, factor_cols, drop = FALSE], leaf = pred[, tree]) + if (isTRUE(oob)) { + dt <- dt[!is.na(leaf)] + } + dt <- melt(dt, id.vars = "leaf", variable.factor = FALSE, + value.factor = FALSE, value.name = "val")[, tree := tree] + dt[, NA_share := sum(is.na(val)) / .N, by = .(leaf, variable)] + dt <- dt[!(is.na(val) & NA_share != 1)] + if (dt[, any(NA_share == 1)]) { + # Handle leaves where all values for a categorical variable are NA + all_na <- unique(dt[NA_share == 1, ]) + dt <- dt[NA_share != 1, ] + all_na <- merge(all_na, bnds[, .(tree, leaf, variable, min, max, f_idx)], + by = c("tree", "leaf", "variable"), sort = FALSE) + all_na[!is.finite(min), min := 0.5] + for (j in names(which(factor_cols))) { + all_na[!is.finite(max) & variable == j, + max := lvl_df_rf[variable == j, max(level)]] + } + all_na[!grepl("\\.5", min), min := min + 0.5] + all_na[!grepl("\\.5", max), max := max + 0.5] + all_na[, min := min + 0.5][, max := max - 0.5] + all_na <- all_na[, .(level = seq(min, max), NA_share), + by = .(leaf, variable)] + all_na <- merge(all_na, lvl_df_rf, by = c("variable", "level")) + all_na[, level := NULL][, tree := tree] + setcolorder(all_na, colnames(dt)) + dt <- rbind(dt, all_na) + } + dt[, count := .N, by = .(leaf, variable)] + dt <- merge(dt, bnds[, .(tree, leaf, variable, min, max, f_idx)], + by = c("tree", "leaf", "variable"), sort = FALSE) + dt[, c("tree", "leaf") := NULL] + if (alpha == 0) { + dt <- unique(dt[, prob := .N / count, by = .(f_idx, variable, val)]) + } else { + # Define the range of each variable in each leaf + dt <- unique(dt[, val_count := .N, by = .(f_idx, variable, val)]) + dt <- merge(dt, lvl_df_rf[, .(k = .N), by = variable], by = "variable") + dt[!is.finite(min), min := 0.5][!is.finite(max), max := k + 0.5] + dt[!grepl("\\.5", min), min := min + 0.5][ + !grepl("\\.5", max), max := max + 0.5] + dt[, k := max - min] + # Enumerate each possible leaf-variable-value combo + tmp <- dt[, seq(min[1] + 0.5, max[1] - 0.5), + by = .(f_idx, variable)] + setnames(tmp, "V1", "level") + tmp <- merge(tmp, lvl_df_rf, by = c("variable", "level"), + sort = FALSE)[, level := NULL] + # Populate count, k + tmp <- merge(tmp, unique(dt[, .(f_idx, variable, count, k)]), + by = c("f_idx", "variable"), sort = FALSE) + # Merge with dt, set val_count = 0 for possible but unobserved levels + dt <- merge(tmp, dt, by = c("f_idx", "variable", "val", "count", "k"), + all.x = TRUE, sort = FALSE) + dt[is.na(val_count), val_count := 0] + dt[, NA_share := mean(NA_share, na.rm = TRUE), + by = .(f_idx, variable)] + # Compute posterior probabilities + dt[, prob := (val_count + alpha) / (count + alpha * k), + by = .(f_idx, variable, val)] + dt[, c("val_count", "k") := NULL] + } + dt[, c("count", "min", "max") := NULL] + setcolorder(dt, c("f_idx", "variable", "val", "prob", "NA_share")) + dt +} diff --git a/R/forge.R b/R/forge.R index d42c6043..e6cd6cb8 100644 --- a/R/forge.R +++ b/R/forge.R @@ -20,10 +20,16 @@ #' \code{NA} (\code{"na"}). The default is \code{"force"}. #' @param verbose Show warnings, e.g. when no leaf matches a condition? #' @param stepsize How many rows of evidence should be handled at each step? -#' Defaults to \code{nrow(evidence) / num_registered_workers} for +#' Defaults to \code{nrow(evidence)} divided by the number of registered +#' workers or daemons for #' \code{parallel == TRUE}. -#' @param parallel Compute in parallel? Must register backend beforehand, e.g. -#' via \code{doParallel} or \code{doFuture}; see examples. +#' @param parallel Compute in parallel? Requires a registered \code{foreach} +#' backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +#' daemons. See \code{\link{arf-options}}. With +#' \code{evidence_row_mode = "or"}, parallelization happens inside the +#' conditional circuit computation; in benchmarks this gave little speedup +#' while raising peak memory, so consider \code{parallel = FALSE} for large +#' \code{"or"} queries. #' @param n_synth Number of synthetic samples to generate. #' #' @details @@ -99,6 +105,9 @@ #' # ... or with doFuture #' doFuture::registerDoFuture() #' future::plan("multisession", workers = 4) +#' +#' # ... or with mirai (shares large read-only inputs across workers via mori) +#' mirai::daemons(4) #' } #' #' @seealso @@ -131,7 +140,12 @@ forge <- function( tree <- cvg <- leaf <- idx <- family <- mu <- sigma <- prob <- dat <- variable <- relation <- wt <- j <- f_idx <- val <- . <- step_ <- c_idx <- f_idx_uncond <- N <- step <- V1 <- NULL - + + # Defaults so the extracted per-step worker always receives these (the + # evidence branches below set them for the conditional cases). + stepsize_cforde <- 0L + parallel_cforde <- FALSE + factor_cols <- params$meta[, family == 'multinom'] # Prepare evidence and stepsize @@ -141,7 +155,7 @@ forge <- function( evidence <- as.data.table(evidence) if (stepsize == 0) { if (parallel) { - stepsize <- ceiling(nrow(evidence)/foreach::getDoParWorkers()) + stepsize <- ceiling(nrow(evidence)/arf_n_workers()) } else { stepsize <- nrow(evidence) } @@ -162,143 +176,37 @@ forge <- function( step_no <- ceiling(nrow(evidence)/stepsize) } - # Run in parallel for each step + # Per-step work lives once in arf_forge_step() (forge_workers.R), shared by all + # backends; this closure adapts it to foreach's one-argument iteration. par_fun <- function(step_) { - - # Prepare the event space - if (is.null(evidence) || ( ncol(evidence) == 2 && all(colnames(evidence) == c("f_idx", "wt")))) { - cparams <- NULL - } else { - # Call cforde with part of the evidence for this step - index_start <- (step_-1)*stepsize + 1 - index_end <- min(step_*stepsize, nrow(evidence)) - evidence_part <- evidence[index_start:index_end,] - cparams <- cforde(params, evidence_part, evidence_row_mode, nomatch, verbose, - stepsize_cforde, parallel_cforde) - if (is.null(cparams)) { - n_synth <- n_synth * nrow(evidence_part) - } - } - - # omega contains the weight (wt) for each leaf (f_idx) for each condition (c_idx) - if (is.null(cparams)) { - if (is.null(evidence)) { - num_trees <- params$forest[, max(tree)] - omega <- params$forest[, .(f_idx, f_idx_uncond = f_idx, cvg)] - omega[, `:=` (c_idx = 1, wt = cvg / num_trees)] - omega[, cvg := NULL] - } else { - omega <- copy(evidence) - omega[, f_idx_uncond := f_idx] - omega[, c_idx := 1] - } - } else { - omega <- cparams$forest[, .(c_idx, f_idx, f_idx_uncond, wt = cvg)] - } - omega <- omega[wt > 0, ] - - # For each synthetic sample and condition, draw a leaf according to the leaf weights - if (nrow(omega) == 1) { - omega <- omega[rep(1, n_synth),][, idx := .I] - } else { - if (evidence_row_mode == "or") { - draws <- omega[, .(f_idx = resample(f_idx, size = n_synth, replace = TRUE, prob = wt))] - omega <- merge(draws, omega, by = "f_idx", sort = FALSE)[, idx := .I] - } else { - draws <- omega[, .(f_idx = resample(f_idx, size = n_synth, replace = TRUE, prob = wt)), by = c_idx] - omega <- merge(draws, omega, by = c("c_idx", "f_idx"), sort = FALSE)[, idx := .I] - } - setcolorder(omega, "idx") - } - - # Simulate continuous data - synth_cnt <- synth_cat <- NULL - if (any(!factor_cols)) { - fam <- params$meta[family != 'multinom', unique(family)] - if (is.null(cparams)) { - psi_cond <- data.table() - } else { - psi_cond <- merge(omega, cparams$cnt[,-c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), - sort = FALSE, allow.cartesian = TRUE)[prob > 0,] - # draw sub-leaf areas (resulting from within-row or-conditions) - if(any(psi_cond[,prob != 1])) { - psi_cond[, I := .I] - psi_cond <- psi_cond[sort(c(psi_cond[prob == 1, I], - psi_cond[prob > 0 & prob < 1, fifelse(.N > 1, resample(I, 1, prob = prob), 0), by = .(variable, idx)][,V1])), -"I"] - } - psi_cond[, prob := NULL] - } - psi <- unique(rbind(psi_cond, - merge(omega, params$cnt, by.x = 'f_idx_uncond', by.y = 'f_idx', - sort = FALSE, allow.cartesian = TRUE)[,val := NA_real_]), - by = c("idx", "variable")) - if (fam == 'truncnorm') { - psi[is.na(val), val := truncnorm::rtruncnorm(.N, a = min, b = max, mean = mu, sd = sigma)] - psi[is.na(val), val := mu] - } else if (fam == 'unif') { - psi[is.na(val), val := stats::runif(.N, min = min, max = max)] - } - NA_share_cnt <- psi[,.(idx, variable, NA_share)] - synth_cnt <- dcast(psi, idx ~ variable, value.var = 'val')[, idx := NULL] - } - - # Simulate categorical data - if (any(factor_cols)) { - if (is.null(cparams)) { - psi <- merge(omega, params$cat, by.x = 'f_idx_uncond', by.y = 'f_idx', sort = FALSE, allow.cartesian = TRUE) - } else { - psi_cond <- merge(omega, cparams$cat[,-c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), - sort = FALSE, allow.cartesian = TRUE) - psi_uncond <- merge(omega, params$cat, by.x = 'f_idx_uncond', by.y = 'f_idx', - sort = FALSE, allow.cartesian = TRUE) - psi_uncond_relevant <- psi_uncond[!psi_cond, on = .(idx, variable)] - psi <- rbind(psi_cond, psi_uncond_relevant) - } - psi[prob < 1, val := sample(val, 1, prob = prob), by = .(variable, idx)] - - psi <- unique(psi[, .(idx, variable, val, NA_share)]) - NA_share_cat <- psi[,.(idx, variable, NA_share)] - synth_cat <- dcast(psi, idx ~ variable, value.var = 'val')[, idx := NULL] - } - - # Combine, optionally impose constraint(s) - x_synth <- cbind(synth_cnt, synth_cat) - if (length(x_synth) == 0) { - x_synth <- evidence_part[FALSE,] - } - - # Clean up, export - x_synth <- post_x(x_synth, params, round) - - if (sample_NAs) { - setDT(x_synth) - NA_share <- rbind(NA_share_cnt, NA_share_cat) - setorder(NA_share[,variable := factor(variable, levels = params$meta[,variable])], variable, idx) - NA_share[,dat := rbinom(.N, 1, prob = NA_share)] - x_synth[dcast(NA_share, formula = idx ~ variable, value.var = "dat")[,-"idx"] == 1] <- NA - x_synth <- post_x(x_synth, params, round) - } - if (evidence_row_mode == "separate" & any(omega[, is.na(f_idx)])) { - setDT(x_synth) - indices_na <- cparams$forest[is.na(f_idx), c_idx] - indices_sampled <- cparams$forest[!is.na(f_idx), unique(c_idx)] - rows_na <- dcast(rbind(data.table(c_idx = 0, variable = params$meta[,variable]), - cparams$evidence_prepped[c_idx %in% indices_na,], - fill = TRUE), - c_idx ~ variable, value.var = "val")[c_idx != 0,] - rows_na <- rbindlist(replicate(n_synth, rows_na, simplify = FALSE)) - if (nomatch == "force") { - rows_na_sampled <- forge(params, n_synth = nrow(rows_na), sample_NAs = sample_NAs, parallel = parallel, stepsize = stepsize) - rows_na[is.na(rows_na)] <- rows_na_sampled[is.na(rows_na[,-1])] - } - x_synth[, c_idx := rep(indices_sampled, each = n_synth)] - x_synth <- rbind(x_synth, rows_na, fill = TRUE) - setorder(x_synth, c_idx)[, c_idx := NULL] - x_synth <- post_x(x_synth, params, round) - } - x_synth - } - if (isTRUE(parallel)) { + arf_forge_step(step_, params, evidence, n_synth, factor_cols, + evidence_row_mode, nomatch, verbose, round, sample_NAs, + stepsize, stepsize_cforde, parallel_cforde) + } + # Parallelism is across steps, so a single step is inherently serial regardless + # of backend (1-task %dopar% is pure overhead). Only pick a backend when + # step_no > 1; the "or" branch already set parallel <- FALSE (cforde + # parallelizes there instead). + use_mirai <- FALSE + if (step_no > 1) { + backend <- arf_select_backend(parallel) + use_mirai <- identical(backend, "mirai") + } + if (use_mirai) { + arf_load_on_daemons() # daemons need arf: worker calls cforde/post_x/resample + params_shared <- mori::share(params) + evidence_shared <- if (!is.null(evidence)) mori::share(evidence) else NULL + x_synth_ <- arf_mirai_tree_map(step_no, arf_forge_step, list( + params = params_shared, evidence = evidence_shared, n_synth = n_synth, + factor_cols = factor_cols, evidence_row_mode = evidence_row_mode, + nomatch = nomatch, verbose = verbose, round = round, + sample_NAs = sample_NAs, stepsize = stepsize, + stepsize_cforde = stepsize_cforde, parallel_cforde = parallel_cforde), + # package-level combine: an inline closure here would serialize forge's + # whole frame (params included) into every task + # See the closure note in mirai_helpers.R. + combine = arf_rbind_steps) + } else if (isTRUE(parallel) && step_no > 1) { x_synth_ <- foreach(step = 1:step_no, .combine = "rbind") %dopar% par_fun(step) } else { x_synth_ <- foreach(step = 1:step_no, .combine = "rbind") %do% par_fun(step) diff --git a/R/forge_workers.R b/R/forge_workers.R new file mode 100644 index 00000000..bf21a06a --- /dev/null +++ b/R/forge_workers.R @@ -0,0 +1,150 @@ +#' @keywords internal +#' @noRd + +# Per-step worker for forge(), extracted so the foreach and mirai backends share +# one definition. All dependencies are explicit arguments (params/evidence are +# mori-shared read-only under the mirai backend). Unlike the forde workers this +# calls arf internals (cforde, resample, post_x), so mirai daemons must have arf +# loaded (see arf_load_on_daemons() in forge()). +arf_forge_step <- function(step_, params, evidence, n_synth, factor_cols, + evidence_row_mode, nomatch, verbose, round, + sample_NAs, stepsize, stepsize_cforde, + parallel_cforde) { + # To avoid data.table check issues + tree <- cvg <- leaf <- idx <- family <- mu <- sigma <- prob <- dat <- + variable <- relation <- wt <- j <- f_idx <- val <- . <- step_x <- c_idx <- + f_idx_uncond <- N <- I <- V1 <- min <- max <- NA_share <- NULL + + # Prepare the event space + if (is.null(evidence) || (ncol(evidence) == 2 && all(colnames(evidence) == c("f_idx", "wt")))) { + cparams <- NULL + } else { + # Call cforde with part of the evidence for this step + index_start <- (step_ - 1) * stepsize + 1 + index_end <- min(step_ * stepsize, nrow(evidence)) + evidence_part <- evidence[index_start:index_end, ] + cparams <- cforde(params, evidence_part, evidence_row_mode, nomatch, verbose, + stepsize_cforde, parallel_cforde) + if (is.null(cparams)) { + n_synth <- n_synth * nrow(evidence_part) + } + } + + # omega contains the weight (wt) for each leaf (f_idx) for each condition (c_idx) + if (is.null(cparams)) { + if (is.null(evidence)) { + num_trees <- params$forest[, max(tree)] + omega <- params$forest[, .(f_idx, f_idx_uncond = f_idx, cvg)] + omega[, `:=`(c_idx = 1, wt = cvg / num_trees)] + omega[, cvg := NULL] + } else { + omega <- copy(evidence) + omega[, f_idx_uncond := f_idx] + omega[, c_idx := 1] + } + } else { + omega <- cparams$forest[, .(c_idx, f_idx, f_idx_uncond, wt = cvg)] + } + omega <- omega[wt > 0, ] + + # For each synthetic sample and condition, draw a leaf according to the leaf weights + if (nrow(omega) == 1) { + omega <- omega[rep(1, n_synth), ][, idx := .I] + } else { + if (evidence_row_mode == "or") { + draws <- omega[, .(f_idx = resample(f_idx, size = n_synth, replace = TRUE, prob = wt))] + omega <- merge(draws, omega, by = "f_idx", sort = FALSE)[, idx := .I] + } else { + draws <- omega[, .(f_idx = resample(f_idx, size = n_synth, replace = TRUE, prob = wt)), by = c_idx] + omega <- merge(draws, omega, by = c("c_idx", "f_idx"), sort = FALSE)[, idx := .I] + } + setcolorder(omega, "idx") + } + + # Simulate continuous data + synth_cnt <- synth_cat <- NULL + if (any(!factor_cols)) { + fam <- params$meta[family != 'multinom', unique(family)] + if (is.null(cparams)) { + psi_cond <- data.table() + } else { + psi_cond <- merge(omega, cparams$cnt[, -c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), + sort = FALSE, allow.cartesian = TRUE)[prob > 0, ] + # draw sub-leaf areas (resulting from within-row or-conditions) + if (any(psi_cond[, prob != 1])) { + psi_cond[, I := .I] + psi_cond <- psi_cond[sort(c(psi_cond[prob == 1, I], + psi_cond[prob > 0 & prob < 1, fifelse(.N > 1, resample(I, 1, prob = prob), 0), by = .(variable, idx)][, V1])), -"I"] + } + psi_cond[, prob := NULL] + } + psi <- unique(rbind(psi_cond, + merge(omega, params$cnt, by.x = 'f_idx_uncond', by.y = 'f_idx', + sort = FALSE, allow.cartesian = TRUE)[, val := NA_real_]), + by = c("idx", "variable")) + if (fam == 'truncnorm') { + psi[is.na(val), val := truncnorm::rtruncnorm(.N, a = min, b = max, mean = mu, sd = sigma)] + psi[is.na(val), val := mu] + } else if (fam == 'unif') { + psi[is.na(val), val := stats::runif(.N, min = min, max = max)] + } + NA_share_cnt <- psi[, .(idx, variable, NA_share)] + synth_cnt <- dcast(psi, idx ~ variable, value.var = 'val')[, idx := NULL] + } + + # Simulate categorical data + if (any(factor_cols)) { + if (is.null(cparams)) { + psi <- merge(omega, params$cat, by.x = 'f_idx_uncond', by.y = 'f_idx', sort = FALSE, allow.cartesian = TRUE) + } else { + psi_cond <- merge(omega, cparams$cat[, -c("cvg_factor", "f_idx_uncond")], by = c('c_idx', 'f_idx'), + sort = FALSE, allow.cartesian = TRUE) + psi_uncond <- merge(omega, params$cat, by.x = 'f_idx_uncond', by.y = 'f_idx', + sort = FALSE, allow.cartesian = TRUE) + psi_uncond_relevant <- psi_uncond[!psi_cond, on = .(idx, variable)] + psi <- rbind(psi_cond, psi_uncond_relevant) + } + psi[prob < 1, val := sample(val, 1, prob = prob), by = .(variable, idx)] + psi <- unique(psi[, .(idx, variable, val, NA_share)]) + NA_share_cat <- psi[, .(idx, variable, NA_share)] + synth_cat <- dcast(psi, idx ~ variable, value.var = 'val')[, idx := NULL] + } + + # Combine, optionally impose constraint(s) + x_synth <- cbind(synth_cnt, synth_cat) + if (length(x_synth) == 0) { + x_synth <- evidence_part[FALSE, ] + } + + # Clean up, export + x_synth <- post_x(x_synth, params, round) + + if (sample_NAs) { + setDT(x_synth) + NA_share <- rbind(NA_share_cnt, NA_share_cat) + setorder(NA_share[, variable := factor(variable, levels = params$meta[, variable])], variable, idx) + NA_share[, dat := stats::rbinom(.N, 1, prob = NA_share)] + x_synth[dcast(NA_share, formula = idx ~ variable, value.var = "dat")[, -"idx"] == 1] <- NA + x_synth <- post_x(x_synth, params, round) + } + if (evidence_row_mode == "separate" & any(omega[, is.na(f_idx)])) { + setDT(x_synth) + indices_na <- cparams$forest[is.na(f_idx), c_idx] + indices_sampled <- cparams$forest[!is.na(f_idx), unique(c_idx)] + rows_na <- dcast(rbind(data.table(c_idx = 0, variable = params$meta[, variable]), + cparams$evidence_prepped[c_idx %in% indices_na, ], + fill = TRUE), + c_idx ~ variable, value.var = "val")[c_idx != 0, ] + rows_na <- rbindlist(replicate(n_synth, rows_na, simplify = FALSE)) + if (nomatch == "force") { + # nested recovery draw runs serial (a worker must not spawn its own backend) + rows_na_sampled <- forge(params, n_synth = nrow(rows_na), sample_NAs = sample_NAs, parallel = FALSE) + rows_na[is.na(rows_na)] <- rows_na_sampled[is.na(rows_na[, -1])] + } + x_synth[, c_idx := rep(indices_sampled, each = n_synth)] + x_synth <- rbind(x_synth, rows_na, fill = TRUE) + setorder(x_synth, c_idx)[, c_idx := NULL] + x_synth <- post_x(x_synth, params, round) + } + x_synth +} diff --git a/R/impute.R b/R/impute.R index bf2f15f6..683a6bb5 100644 --- a/R/impute.R +++ b/R/impute.R @@ -43,6 +43,9 @@ #' # ... or with doFuture #' doFuture::registerDoFuture() #' future::plan("multisession", workers = 4) +#' +#' # ... or with mirai (shares large read-only inputs across workers via mori) +#' mirai::daemons(4) #' } #' #' @seealso diff --git a/R/lik.R b/R/lik.R index af2ecce6..547aee01 100644 --- a/R/lik.R +++ b/R/lik.R @@ -21,8 +21,9 @@ #' queries in one round, which is always the fastest option if memory allows. #' However, with large samples or many trees, it can be more memory efficient #' to split the data into batches. This has no impact on results. -#' @param parallel Compute in parallel? Must register backend beforehand, e.g. -#' via \code{doParallel} or \code{doFuture}; see examples. +#' @param parallel Compute in parallel? Requires a registered \code{foreach} +#' backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +#' daemons. See \code{\link{arf-options}}. #' #' #' @details @@ -83,6 +84,9 @@ #' # ... or with doFuture #' doFuture::registerDoFuture() #' future::plan("multisession", workers = 4) +#' +#' # ... or with mirai (shares large read-only inputs across workers via mori) +#' mirai::daemons(4) #' } #' #' @seealso @@ -127,11 +131,13 @@ lik <- function( } x <- suppressWarnings(prep_x(x)) factor_cols <- sapply(x, is.factor) + pure <- all(factor_cols) | all(!factor_cols) # worker arg (was computed inside lik_fn) # Prep evidence conj <- !is.null(evidence) && !(ncol(evidence) == 2 && all(c("f_idx", "wt") %in% colnames(evidence))) # Check ARF + preds <- NULL # set below iff arf-based leaf assignment is used if (d == params$meta[, .N] & !is.null(arf)) { num_trees <- arf$num.trees preds <- stats::predict(arf, x, type = 'terminalNodes')$predictions + 1L @@ -173,145 +179,38 @@ lik <- function( } batch_idx <- suppressWarnings(split(seq_len(n), seq_len(k))) - # Likelihood function + # Per-fold work lives once in arf_lik_fold() (lik_workers.R). This closure + # adapts it to foreach; `arf` is used only via is.null() (preds carries the + # arf-derived leaf assignments), passed to the worker as a boolean. lik_fn <- function(fold, arf) { - - # Prep work - psi_cnt <- psi_cat <- NULL - pure <- all(factor_cols) | all(!factor_cols) - if (is.null(arf) & !isTRUE(pure)) { - omega_tmp <- rbindlist(lapply(batch_idx[[fold]], function(i) { - omega$obs <- i - omega$wt <- NULL - return(omega) - })) - } - - # Continuous data - if (any(!factor_cols)) { - fam <- params$meta[class == 'numeric', unique(family)] - x_long <- melt( - data.table(obs = batch_idx[[fold]], - x[batch_idx[[fold]], !factor_cols, drop = FALSE]), - id.vars = 'obs', variable.factor = FALSE - ) - if (is.null(arf)) { - psi_cnt <- merge(params$cnt[f_idx %in% leaves], x_long, by = 'variable', - sort = FALSE, allow.cartesian = TRUE) - rm(x_long) - } else { - preds_cnt <- merge(preds[f_idx %in% leaves], x_long, by = 'obs', - sort = FALSE, allow.cartesian = TRUE) - rm(x_long) - psi_cnt <- merge(params$cnt[f_idx %in% leaves], preds_cnt, - by = c('f_idx', 'variable'), sort = FALSE) - rm(preds_cnt) - } - if (fam == 'truncnorm') { - psi_cnt[, lik := truncnorm::dtruncnorm(value, a = min, b = max, - mean = mu, sd = sigma)] - } else if (fam == 'unif') { - psi_cnt[, lik := stats::dunif(value, min = min, max = max)] - } - psi_cnt[value == min, lik := 0] - psi_cnt[, lik := prod(lik), by = .(f_idx, obs)] - psi_cnt <- unique(psi_cnt[lik > 0, .(f_idx, obs, lik)]) - if (is.null(arf) & !isTRUE(pure)) { - omega_tmp <- merge(omega_tmp, psi_cnt[, .(f_idx, obs)], - by = c('f_idx', 'obs'), sort = FALSE) - leaves <- omega_tmp[, unique(f_idx)] - } - } - - # Categorical data - if (any(factor_cols)) { - x_tmp <- x[batch_idx[[fold]], factor_cols, drop = FALSE] - n_tmp <- nrow(x_tmp) - x_long <- melt( - data.table(obs = batch_idx[[fold]], x_tmp), - id.vars = 'obs', value.name = 'val', variable.factor = FALSE - ) - # Speedups are possible if there are many duplicates - is_unique <- !duplicated(x_tmp) - if (all(is_unique)) { - x_unique <- x_long - colnames(x_unique)[1] <- 's_idx' - } else { - x_unique <- unique(x_tmp) - x_unique <- melt( - data.table(s_idx = seq_len(nrow(x_unique)), x_unique), - id.vars = 's_idx', value.name = 'val', variable.factor = FALSE - ) - s_idx <- integer(length = n_tmp) - s_idx[is_unique] <- seq_len(sum(is_unique)) - for (i in 2:n_tmp) { - if (s_idx[i] == 0L) { - s_idx[i] <- s_idx[i - 1L] - } - } - idx_dt <- data.table(obs = batch_idx[[fold]], s_idx = s_idx) - } - if (is.null(arf)) { - grd <- rbindlist(lapply(which(factor_cols), function(j) { - expand.grid('f_idx' = leaves, 'variable' = colnames(x)[j], - 'val' = x_long[variable == colnames(x)[j], unique(val)], - stringsAsFactors = FALSE) - })) - rm(x_long) - psi_cat <- merge(params$cat[f_idx %in% leaves], grd, - by = c('f_idx', 'variable', 'val'), - sort = FALSE, all.y = TRUE) - rm(grd) - psi_cat[is.na(prob), prob := 0] - psi_cat <- merge(psi_cat, x_unique, by = c('variable', 'val'), - sort = FALSE, allow.cartesian = TRUE) - psi_cat[, lik := prod(prob), by = .(f_idx, s_idx)] - psi_cat <- unique(psi_cat[lik > 0, .(f_idx, s_idx, lik)]) - if (all(is_unique)) { - setnames(psi_cat, 's_idx', 'obs') - } else { - if (!isTRUE(pure)) { - omega_tmp <- merge(idx_dt, omega_tmp, by = 'obs', sort = FALSE) - psi_cat <- merge(psi_cat, omega_tmp, by = c('f_idx', 's_idx'), - sort = FALSE)[, s_idx := NULL] - rm(omega_tmp) - setcolorder(psi_cat, c('f_idx', 'obs', 'lik')) - psi_cnt <- merge(psi_cnt, psi_cat[, .(f_idx, obs)], - by = c('f_idx', 'obs'), sort = FALSE) - } - } - } else { - preds_cat <- merge(preds[f_idx %in% leaves], x_long, by = 'obs', - sort = FALSE, allow.cartesian = TRUE) - rm(x_long) - psi_cat <- merge(params$cat, preds_cat, by = c('f_idx', 'variable', 'val'), - sort = FALSE, allow.cartesian = TRUE, all.y = TRUE) - rm(preds_cat) - psi_cat[is.na(prob), prob := 0] - psi_cat[, lik := prod(prob), by = .(f_idx, obs)] - psi_cat <- unique(psi_cat[lik > 0, .(f_idx, obs, lik)]) - } - } - - # Put it together - psi_x <- rbind(psi_cnt, psi_cat) - if (!isTRUE(pure)) { - psi_x <- psi_x[, prod(lik), by = .(f_idx, obs)] - setnames(psi_x, 'V1', 'lik') - } - return(psi_x) + arf_lik_fold(fold, params, x, factor_cols, leaves, omega, preds, + batch_idx, pure, !is.null(arf)) + } + # Parallelism is across folds: a single fold is inherently serial. mirai only + # for k > 1. + use_mirai <- FALSE + if (k > 1) { + backend <- arf_select_backend(parallel) + use_mirai <- identical(backend, "mirai") } - if (isTRUE(parallel)) { + if (use_mirai) { + arf_load_on_daemons() # daemons need arf (worker uses bare data.table verbs) + out <- arf_mirai_tree_map(k, arf_lik_fold, list( + params = mori::share(params), x = mori::share(x), + factor_cols = factor_cols, leaves = mori::share(leaves), + omega = mori::share(omega), + preds = if (!is.null(preds)) mori::share(preds) else NULL, + batch_idx = mori::share(batch_idx), pure = pure, + has_arf = !is.null(arf))) + } else if (isTRUE(parallel) && k > 1) { out <- foreach(fold = seq_len(k), .combine = rbind) %dopar% lik_fn(fold, arf) } else { out <- foreach(fold = seq_len(k), .combine = rbind) %do% lik_fn(fold, arf) } - # Compute per-sample likelihoods - out <- merge(out, omega, by = 'f_idx', sort = FALSE) - out <- out[, log(crossprod(wt, lik)), by = obs] - setnames(out, 'V1', 'lik') - + # Folds return per-obs log-likelihoods, already reduced against omega inside + # the worker (folds cover disjoint obs); nothing left to aggregate here. + # Anybody missing? zeros <- setdiff(seq_len(n), out[, obs]) if (length(zeros) > 0L) { diff --git a/R/lik_workers.R b/R/lik_workers.R new file mode 100644 index 00000000..d352493d --- /dev/null +++ b/R/lik_workers.R @@ -0,0 +1,148 @@ +#' @keywords internal +#' @noRd + +# Per-fold worker for lik(), extracted so foreach and mirai backends share one +# definition (params/x/preds/omega mori-shared under mirai). `has_arf` replaces +# the closed-over arf object (only used for is.null() branching; preds carries +# the arf-derived leaf assignments). Uses bare data.table verbs, so mirai +# daemons must have arf loaded (arf imports data.table). +arf_lik_fold <- function(fold, params, x, factor_cols, leaves, omega, preds, + batch_idx, pure, has_arf) { + # To avoid data.table check issues + tree <- cvg <- leaf <- variable <- mu <- sigma <- value <- obs <- prob <- + V1 <- relation <- f_idx <- wt <- val <- family <- f_idx_uncond <- . <- + lik <- s_idx <- min <- max <- NULL + + # Prep work + psi_cnt <- psi_cat <- NULL + if (!has_arf & !isTRUE(pure)) { + omega_tmp <- rbindlist(lapply(batch_idx[[fold]], function(i) { + omega$obs <- i + omega$wt <- NULL + return(omega) + })) + } + + # Continuous data + if (any(!factor_cols)) { + fam <- params$meta[class == 'numeric', unique(family)] + x_long <- melt( + data.table(obs = batch_idx[[fold]], + x[batch_idx[[fold]], !factor_cols, drop = FALSE]), + id.vars = 'obs', variable.factor = FALSE + ) + if (!has_arf) { + psi_cnt <- merge(params$cnt[f_idx %in% leaves], x_long, by = 'variable', + sort = FALSE, allow.cartesian = TRUE) + rm(x_long) + } else { + preds_cnt <- merge(preds[f_idx %in% leaves], x_long, by = 'obs', + sort = FALSE, allow.cartesian = TRUE) + rm(x_long) + psi_cnt <- merge(params$cnt[f_idx %in% leaves], preds_cnt, + by = c('f_idx', 'variable'), sort = FALSE) + rm(preds_cnt) + } + if (fam == 'truncnorm') { + psi_cnt[, lik := truncnorm::dtruncnorm(value, a = min, b = max, + mean = mu, sd = sigma)] + } else if (fam == 'unif') { + psi_cnt[, lik := stats::dunif(value, min = min, max = max)] + } + psi_cnt[value == min, lik := 0] + psi_cnt[, lik := prod(lik), by = .(f_idx, obs)] + psi_cnt <- unique(psi_cnt[lik > 0, .(f_idx, obs, lik)]) + if (!has_arf & !isTRUE(pure)) { + omega_tmp <- merge(omega_tmp, psi_cnt[, .(f_idx, obs)], + by = c('f_idx', 'obs'), sort = FALSE) + leaves <- omega_tmp[, unique(f_idx)] + } + } + + # Categorical data + if (any(factor_cols)) { + x_tmp <- x[batch_idx[[fold]], factor_cols, drop = FALSE] + n_tmp <- nrow(x_tmp) + x_long <- melt( + data.table(obs = batch_idx[[fold]], x_tmp), + id.vars = 'obs', value.name = 'val', variable.factor = FALSE + ) + # Speedups are possible if there are many duplicates + is_unique <- !duplicated(x_tmp) + if (all(is_unique)) { + x_unique <- x_long + colnames(x_unique)[1] <- 's_idx' + } else { + x_unique <- unique(x_tmp) + x_unique <- melt( + data.table(s_idx = seq_len(nrow(x_unique)), x_unique), + id.vars = 's_idx', value.name = 'val', variable.factor = FALSE + ) + s_idx <- integer(length = n_tmp) + s_idx[is_unique] <- seq_len(sum(is_unique)) + for (i in 2:n_tmp) { + if (s_idx[i] == 0L) { + s_idx[i] <- s_idx[i - 1L] + } + } + idx_dt <- data.table(obs = batch_idx[[fold]], s_idx = s_idx) + } + if (!has_arf) { + grd <- rbindlist(lapply(which(factor_cols), function(j) { + expand.grid('f_idx' = leaves, 'variable' = colnames(x)[j], + 'val' = x_long[variable == colnames(x)[j], unique(val)], + stringsAsFactors = FALSE) + })) + rm(x_long) + psi_cat <- merge(params$cat[f_idx %in% leaves], grd, + by = c('f_idx', 'variable', 'val'), + sort = FALSE, all.y = TRUE) + rm(grd) + psi_cat[is.na(prob), prob := 0] + psi_cat <- merge(psi_cat, x_unique, by = c('variable', 'val'), + sort = FALSE, allow.cartesian = TRUE) + psi_cat[, lik := prod(prob), by = .(f_idx, s_idx)] + psi_cat <- unique(psi_cat[lik > 0, .(f_idx, s_idx, lik)]) + if (all(is_unique)) { + setnames(psi_cat, 's_idx', 'obs') + } else { + if (!isTRUE(pure)) { + omega_tmp <- merge(idx_dt, omega_tmp, by = 'obs', sort = FALSE) + psi_cat <- merge(psi_cat, omega_tmp, by = c('f_idx', 's_idx'), + sort = FALSE)[, s_idx := NULL] + rm(omega_tmp) + setcolorder(psi_cat, c('f_idx', 'obs', 'lik')) + psi_cnt <- merge(psi_cnt, psi_cat[, .(f_idx, obs)], + by = c('f_idx', 'obs'), sort = FALSE) + } + } + } else { + preds_cat <- merge(preds[f_idx %in% leaves], x_long, by = 'obs', + sort = FALSE, allow.cartesian = TRUE) + rm(x_long) + psi_cat <- merge(params$cat, preds_cat, by = c('f_idx', 'variable', 'val'), + sort = FALSE, allow.cartesian = TRUE, all.y = TRUE) + rm(preds_cat) + psi_cat[is.na(prob), prob := 0] + psi_cat[, lik := prod(prob), by = .(f_idx, obs)] + psi_cat <- unique(psi_cat[lik > 0, .(f_idx, obs, lik)]) + } + } + + # Put it together + psi_x <- rbind(psi_cnt, psi_cat) + if (!isTRUE(pure)) { + psi_x <- psi_x[, prod(lik), by = .(f_idx, obs)] + setnames(psi_x, 'V1', 'lik') + } + + # Reduce to per-observation log-likelihoods here rather than on the calling + # process: folds cover disjoint obs and omega is a worker argument, so the + # reduction is fold-local. This shrinks the returned object from one row per + # (obs, leaf) to one per obs (less to serialize back under mirai) and + # parallelizes what used to be a serial post-pass over every obs x leaf pair. + psi_x <- merge(psi_x, omega, by = 'f_idx', sort = FALSE) + psi_x <- psi_x[, log(crossprod(wt, lik)), by = obs] + setnames(psi_x, 'V1', 'lik') + psi_x +} diff --git a/R/mirai_helpers.R b/R/mirai_helpers.R new file mode 100644 index 00000000..4b9e6011 --- /dev/null +++ b/R/mirai_helpers.R @@ -0,0 +1,240 @@ +#' @keywords internal +#' @noRd + +# Internal helpers for the mirai+mori backend. +# Activated by options(arf.backend = "mirai"), or automatically when mirai +# daemons are running and the option is unset; see arf_select_backend(). +# Hidden from public API; mirai and mori in Suggests only. + +# Session-scoped state for once-per-run notifications. +.arf_env <- new.env(parent = emptyenv()) + +# Emit a backend notification at most once per distinct state per session, +# and only when getOption("arf.verbose", TRUE) is TRUE (opt-out switch). +arf_backend_inform <- function(msg, key) { + if (!isTRUE(getOption("arf.verbose", TRUE))) { + return(invisible(FALSE)) + } + shown <- get0("backend_shown", envir = .arf_env, ifnotfound = character(0)) + if (key %in% shown) { + return(invisible(FALSE)) + } + assign("backend_shown", c(shown, key), envir = .arf_env) + message(msg) + invisible(TRUE) +} + +# Decide which parallel backend to use (called from all parallelized +# operations: forde, forge, expct, lik, cforde, adversarial_rf), and tell +# the user. +# Only meaningful when parallel = TRUE. Precedence: +# 1. explicit options(arf.backend = "foreach" | "mirai") (validated) +# 2. active mirai daemons -> "mirai" +# 3. otherwise -> "foreach" +# For the foreach path we additionally report whether a real parallel backend +# is registered (>1 worker) or whether it will fall back to sequential. +# Returns one of "sequential", "foreach", "mirai". +arf_select_backend <- function(parallel) { + if (!isTRUE(parallel)) { + return("sequential") + } + mirai_ready <- requireNamespace("mirai", quietly = TRUE) && + requireNamespace("mori", quietly = TRUE) && + { + st <- mirai::status() + !is.null(st$connections) && st$connections >= 1L + } + dopar_workers <- if (requireNamespace("foreach", quietly = TRUE)) { + foreach::getDoParWorkers() + } else { + 1L + } + + opt <- getOption("arf.backend", NULL) + if (!is.null(opt)) { + # Explicit user choice wins; validate and hard-check mirai readiness. + backend <- match.arg(opt, c("foreach", "mirai")) + if (backend == "mirai") { + arf_check_mirai_ready() + } + } else if (mirai_ready) { + backend <- "mirai" + } else { + backend <- "foreach" + } + + # Report the effective backend: parallel = TRUE only, at most once per + # distinct state per session, suppressible via options(arf.verbose = FALSE). + if (backend == "mirai") { + n <- mirai::status()$connections + arf_backend_inform( + paste0("arf: using 'mirai' backend (", n, " daemon", + if (n != 1L) "s" else "", ")."), + key = paste0("mirai:", n)) + } else if (dopar_workers > 1L) { + arf_backend_inform( + paste0("arf: using 'foreach' backend (", foreach::getDoParName(), + ", ", dopar_workers, " workers)."), + key = paste0("foreach:", dopar_workers)) + } else { + arf_backend_inform( + paste0("arf: parallel = TRUE but no parallel backend is registered; ", + "computing sequentially. Register a foreach backend (e.g. ", + "doParallel) or start mirai daemons via mirai::daemons()."), + key = "sequential-fallback") + } + backend +} + +# Load arf on all mirai daemons so workers can call its internals (e.g. +# cforde/resample/post_x in the forge/expct/lik workers) and data.table's S3 +# methods are registered. No-op on daemons where arf is already loaded (e.g. +# dev-loaded in tests; see tests/testthat/helper-mirai.R). +# +# everywhere() round-trips all daemons on every call, which is dead cost for a +# workflow doing many small forge/expct/lik calls against one pool. Register the +# load once per pool and skip on repeat. The pool key is the dispatcher URL +# (status()$daemons), which is minted fresh by every daemons() call: a +# teardown+rebuild yields a new key so the cache self-invalidates, and daemons +# joining an existing pool auto-run the registered everywhere() expression, so +# one call per pool suffices. If status() is unavailable (key NULL) we fall back +# to the old always-load behavior. +arf_load_on_daemons <- function() { + key <- tryCatch(mirai::status()$daemons, error = function(e) NULL) + cached <- get0("arf_loaded_key", envir = .arf_env, ifnotfound = NULL) + if (!is.null(key) && identical(key, cached)) { + return(invisible(FALSE)) + } + mirai::everywhere(suppressMessages(loadNamespace("arf"))) + assign("arf_loaded_key", key, envir = .arf_env) + invisible(TRUE) +} + +# Worker count for sizing step/fold chunks. Whichever parallel backend is active +# reports >1; the inactive one reports 1, so max() picks the right pool without +# re-deriving backend-selection precedence. Needed because stepsize was sized via +# foreach::getDoParWorkers() alone, which is 1 under a pure mirai backend (daemons +# set, no foreach registered) -> step_no 1 -> mirai never engages. +arf_n_workers <- function() { + mirai_conns <- if (requireNamespace("mirai", quietly = TRUE)) { + st <- tryCatch(mirai::status(), error = function(e) NULL) + if (!is.null(st$connections)) as.integer(st$connections) else 0L + } else { + 0L + } + dopar <- if (requireNamespace("foreach", quietly = TRUE)) { + foreach::getDoParWorkers() + } else { + 1L + } + max(1L, mirai_conns, dopar) +} + +arf_check_mirai_ready <- function() { + if (!requireNamespace("mirai", quietly = TRUE)) { + stop("arf.backend = 'mirai' requires the 'mirai' package. ", + "Install it or set options(arf.backend = 'foreach').", + call. = FALSE) + } + if (!requireNamespace("mori", quietly = TRUE)) { + stop("arf.backend = 'mirai' requires the 'mori' package. ", + "Install it or set options(arf.backend = 'foreach').", + call. = FALSE) + } + status <- mirai::status() + if (is.null(status$connections) || status$connections < 1L) { + stop("arf.backend = 'mirai' requires daemons() to be set. ", + "Call mirai::daemons(n) first.", + call. = FALSE) + } + invisible(TRUE) +} + +# Note: data.table objects come back from mori with truelength == 0 +# (selfref dropped by the ALTREP round-trip), but read-only operations +# in the worker bodies tolerate this. No setalloccol() repair is needed +# because workers only read the shared tables and build fresh data.tables +# from the results. + +# Split trees into one contiguous block per worker. Contiguous (sort()) is +# load-bearing: rbindlist/c concatenate blocks in chunk order, so interleaved +# chunks would scramble positional output (forge/expct rows are +# per-evidence-row). Shared by arf_mirai_tree_map() and the prune dispatch in +# adversarial_rf(); the invariant lives here once. A finer-chunks knob (chunks +# per worker > 1) was benchmarked and removed: it increased memory (more +# result objects in flight) without speed gains; step/fold-parallel ops +# control granularity via their stepsize/batch arguments instead. +arf_tree_chunks <- function(num_trees, n_workers) { + n_chunks <- max(1L, min(as.integer(n_workers), num_trees)) + split(seq_len(num_trees), + sort(rep(seq_len(n_chunks), length.out = num_trees))) +} + +# Dispatch a per-tree worker over mirai daemons: one contiguous chunk per +# daemon (arf_tree_chunks); each task loops the worker over its block and +# combines locally, so at most n_workers (not num_trees) result objects are +# serialized back. `combine` stacks per-tree results within a chunk and again +# across chunks; the default rbindlist returns a data.table (forde re-keys), +# forge/expct pass an rbind-based combine to match the class of their serial +# foreach .combine = "rbind" output. +# CAUTION on closures: serializing a function serializes its enclosing +# environment, so a closure defined in a caller's frame drags that whole frame +# into every task, defeating mori sharing. Ship package-level functions, or +# strip base-R-only closures to globalenv() first. See the mirai FAQ: +# https://mirai.r-lib.org/articles/v07-questions.html +arf_mirai_tree_map <- function(num_trees, worker_fn, shared_args, + combine = data.table::rbindlist) { + st <- mirai::status() + n_workers <- max(1L, as.integer(st$connections)) + chunks <- arf_tree_chunks(num_trees, n_workers) + chunk_runner <- function(trees, worker_fn, shared_args, combine) { + parts <- lapply(trees, function(tr) { + do.call(worker_fn, c(list(tr), shared_args)) + }) + combine(parts) + } + # base-R body, all inputs explicit args: strip so this frame (chunks, + # shared_args, combine, ...) is not serialized into every task + environment(chunk_runner) <- globalenv() + res <- mirai::mirai_map( + chunks, + chunk_runner, + .args = list(worker_fn = worker_fn, shared_args = shared_args, + combine = combine) + )[] + arf_stop_on_mirai_error(res) + combine(res) +} + +# A failed task comes back from mirai_map as an error object inside the +# results list; without this check it would be rbind/c-ed into the output and +# corrupt it silently. Covers errors, interrupts and timeouts. +arf_stop_on_mirai_error <- function(res) { + failed <- which(vapply(res, mirai::is_error_value, logical(1))) + if (length(failed)) { + err <- res[[failed[1]]] + msg <- if (mirai::is_mirai_error(err)) conditionMessage(err) else as.character(err) + stop("arf: mirai worker error in chunk ", failed[1], ": ", msg, call. = FALSE) + } + invisible(TRUE) +} + +# Combine for forge/expct step results: matches serial foreach .combine="rbind" +# (same class, clean 1..n row.names). Package-level on purpose: an inline +# closure in forge()/expct() would serialize their whole frame, params +# included, into every task (see note above). +arf_rbind_steps <- function(parts) { + r <- do.call(rbind, parts) + rownames(r) <- NULL + r +} + +# Combine for forde's fused psi worker (arf_psi_fn returns list(cnt, cat) per +# tree): stack each component across parts. Safe under arf_mirai_tree_map's +# two-level application because the output has the same shape as each input. +# rbindlist drops the NULL component of the branch not taken. Package-level +# for the same closure-hygiene reason as arf_rbind_steps. +arf_combine_psi <- function(parts) { + list(cnt = data.table::rbindlist(lapply(parts, `[[`, "cnt")), + cat = data.table::rbindlist(lapply(parts, `[[`, "cat"))) +} diff --git a/R/prune_workers.R b/R/prune_workers.R new file mode 100644 index 00000000..d2470ef3 --- /dev/null +++ b/R/prune_workers.R @@ -0,0 +1,35 @@ +#' @keywords internal +#' @noRd + +# Per-tree worker for adversarial_rf()'s post-hoc leaf pruning (enforces +# min_node_size w.r.t. real data). Base-R only (no arf/data.table internals), so +# mirai daemons need no package loaded. Returns the pruned child.nodeIDs for one +# tree: a list of two integer vectors (left, right children). +arf_prune_tree <- function(tree, child_nodeIDs, pred, min_node_size) { + # Nodes to prune are leaves which contain fewer than min_node_size real samples + out <- child_nodeIDs[[tree]] + leaves <- which(out[[1]] == 0L) + to_prune <- leaves[!(leaves %in% which(tabulate(pred[, tree]) >= min_node_size))] + while (length(to_prune) > 0) { + if (1 %in% to_prune) { + # Never prune the root + break + } + for (tp in to_prune) { + # Find parent + parent <- which((out[[1]] + 1L) == tp) + if (length(parent) > 0) { + # If node to prune (tp) is the left child of parent, replace left child with right child + out[[1]][parent] <- out[[2]][parent] + } else { + # If node to prune (tp) is the right child of parent, replace right child with left child + parent <- which((out[[2]] + 1L) == tp) + out[[2]][parent] <- out[[1]][parent] + } + } + # If both children of a parent are to be pruned, prune the parent in the next round + # This happens if both children have been pruned + to_prune <- which((out[[1]] + 1L) %in% to_prune) + } + out +} diff --git a/R/utils.R b/R/utils.R index 4e0cfd89..990dce4d 100644 --- a/R/utils.R +++ b/R/utils.R @@ -210,8 +210,9 @@ post_x <- function(x, params, round = TRUE) { #' \code{NA} (\code{"na"}). The default is \code{"force"}. #' @param verbose Show warnings, e.g. when no leaf matches a condition? #' @param stepsize Stepsize defining number of condition rows handled in one for each step. -#' @param parallel Compute in parallel? Must register backend beforehand, e.g. -#' via \code{doParallel} or \code{doFuture}; see examples. +#' @param parallel Compute in parallel? Requires a registered \code{foreach} +#' backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +#' daemons. See \code{\link{arf-options}}. #' #' @return List with conditions (\code{evidence_input}), prepared conditions (\code{evidence_prepped}) #' and leaves that match the conditions in evidence with continuous data (\code{cnt}) @@ -271,152 +272,50 @@ cforde <- function(params, # Calculate stepsize for parallelization depending on number of conditions and registered workers if (stepsize == 0) { if (parallel) { - stepsize <- ceiling(nconds_conditioned/getDoParWorkers()) + stepsize <- ceiling(nconds_conditioned/arf_n_workers()) } else { stepsize <- nconds_conditioned } } step_no <- ceiling(nconds_conditioned/stepsize) - # Loop through conditions with defined stepsize to determined matching leaves and updates for cat and cnt params - update_fun <-function(step_) { - - # Define subset of conditions for step_ - index_start <- conds_conditioned[(step_ - 1)*stepsize + 1] - index_end <- conds_conditioned[min(step_ * stepsize, nconds_conditioned)] - condition_long_step <- condition_long[.(index_start:index_end), nomatch = NULL] - - # Store cat and cnt conditions separately - cat_conds <- condition_long_step[variable %in% cat_cols,c("c_idx","variable","val")][, variable := factor(variable)] - cnt_conds <- condition_long_step[variable %in% cnt_cols,c("c_idx","variable","min", "max","val")][,`:=` (variable = factor(variable), - val = as.numeric(val))] - # If cat conditions exist, calculate matching leaves - if (nrow(cat_conds) != 0) { - - # Save leaf indices f_idx cat params in list column grouped by variable and val (value) and merge with cat conditions - cat_relevant <- cat[, .(.(f_idx)), by=.(variable,val)] - setkey(cat_relevant, variable, val) - setkey(cat_conds, variable, val) - cat_relevant <- cat_conds[cat_relevant, on = .(variable, val), nomatch = NULL] - setkey(cat_relevant, c_idx) - - # or-combine different conditions on the same feature - if (uniqueN(cat_relevant, by = c("c_idx", "variable")) != nrow(cat_relevant)) { - cat_relevant <- cat_relevant[, .(.(Reduce(union, V1))), by = .(c_idx, variable)] - } - - # Determine matching leaves for cat conditions - relevant_leaves_changed_cat <- cat_relevant[, Reduce(intersect, V1), by = c_idx][, .(c_idx, f_idx = V1)] - setorder(relevant_leaves_changed_cat) - conditions_unchanged_cat <- setdiff(condition_long_step[, c_idx], cat_conds[, c_idx]) - relevant_leaves_unchanged_cat <- data.table(c_idx = rep(conditions_unchanged_cat, each = nrow(forest) ), f_idx = rep(forest[,f_idx],length(conditions_unchanged_cat))) - relevant_leaves_cat <- rbind(relevant_leaves_changed_cat, relevant_leaves_unchanged_cat, fill = TRUE) - relevant_leaves_cat_list <- relevant_leaves_cat[,.(f_idx = .(f_idx)), by=c_idx] - } else { - relevant_leaves_cat <- data.table(c_idx = integer(), f_idx = integer()) - } - - # If cnt conditions exist, calculate matching leaves - if (nrow(cnt_conds) != 0) { - - # Save min, max in cnt params in list columns grouped by variable and merge with cnt conditions - cnt_relevant <- cnt[, .(min = .(min), max = .(max)), by = variable] - cnt_conds_compact <- copy(cnt_conds) - cnt_conds_compact[!is.na(val), `:=`(min = val, max = val)][, val := NULL] - cnt_relevant <- cnt_conds_compact[cnt_relevant, on = .(variable), nomatch = NULL] - setkey(cnt_relevant, c_idx) - - # If cat conds exist, use only matching subset of potentially relevant leaves for cnt conditions - if (nrow(cat_conds) != 0) { - cnt_relevant <- cnt_relevant[relevant_leaves_cat_list, on = .(c_idx)] - } else { - cnt_relevant[, f_idx := NA] - } - - # Determine matching leaves for cnt conditions per row - cnt_relevant <- cnt_relevant[, .( - c_idx, - variable, - f_idx = Map(function(f_idx, min, max, i.min, i.max) { - if (!inherits(f_idx, "logical")) { - rel_cnt_min <- i.min[f_idx] - rel_cnt_max <- i.max[f_idx] - rel_min <- f_idx[which(max > rel_cnt_min)] - rel_max <- f_idx[which(min <= rel_cnt_max)] - } else { - rel_cnt_min <- i.min - rel_cnt_max <- i.max - rel_min <- which(max > rel_cnt_min) - rel_max <- which(min <= rel_cnt_max) - } - intersect(rel_min,rel_max) - }, f_idx = f_idx, min = min, max = max, i.min = i.min, i.max = i.max))] - - # or-combine different conditions on the same feature - or_within_row_cnt <- uniqueN(cnt_relevant, by = c("c_idx", "variable")) != nrow(cnt_relevant) - if (or_within_row_cnt) { - cnt_relevant <- cnt_relevant[, .(f_idx = .(Reduce(union, f_idx))), by = .(c_idx, variable)] - } - - # Determine matching leaves for cnt conditions - relevant_leaves_changed_cnt <- cnt_relevant[, Reduce(intersect, f_idx),by = c_idx][,.(c_idx, f_idx = V1)] - conditions_unchanged_cnt <- setdiff(condition_long_step[, c_idx], cnt_conds[, c_idx]) - relevant_leaves_unchanged_cnt <- data.table(c_idx = rep(conditions_unchanged_cnt, each = nrow(forest)), f_idx = rep(forest[,f_idx],length(conditions_unchanged_cnt))) - relevant_leaves_cnt <- rbind(relevant_leaves_changed_cnt, relevant_leaves_unchanged_cnt) - - # Calculate updates for cnt params matching cnt conditions - cnt_new <- merge(merge(relevant_leaves_cnt, cnt_conds, by = "c_idx", allow.cartesian = TRUE, sort = FALSE), cnt, by = c("f_idx", "variable"), all.x = TRUE, allow.cartesian = TRUE, sort = FALSE) - cnt_new[!is.na(val),`:=` (min = min.y, - max = max.y)] - cnt_new[is.na(val),`:=` (min = pmax(min.x, min.y, na.rm = TRUE), - max = pmin(max.x, max.y, na.rm = TRUE))] - cnt_new <- cnt_new[min <= max & sum(max.x, val, na.rm = TRUE) != min.y, ] - cnt_new[, prob := NA_real_] - if (family == "truncnorm") { - cnt_new[!is.na(val), prob := dtruncnorm(val, a=min.y, b=max.y, mean=mu, sd=sigma)*(val != min.y)] - cnt_new[is.na(val) & (min == min.y) & (max == max.y), prob := 1] - cnt_new[is.na(val) & is.na(prob), prob := ptruncnorm(max, a=min.y, b=max.y, mean=mu, sd=sigma) - ptruncnorm(min, a=min.y, max.y, mean=mu,sd=sigma)] - } else if (family == "unif") { - cnt_new[!is.na(val), prob := dunif(val, min=min.y, max=max.y)*(val != min.y)] - cnt_new[is.na(val) & (min == min.y) & (max == max.y), prob := 1] - cnt_new[is.na(val) & is.na(prob), prob := punif(max, min=min.y, max=max.y) - punif(min, min=min.y, max.y)] - } - cnt_new[, c("min.x","max.x","min.y","max.y") := NULL] - - # If or-combined cnt condition within rows exist, calculate likelihoods for ranges within leaves and norm to 1 - if (or_within_row_cnt) { - cnt_new[, cvg_factor := sum(prob), by = .(f_idx, c_idx, variable)] - cnt_new[, prob := prob/cvg_factor] - } else { - cnt_new[, `:=` (cvg_factor = prob, prob = 1)] - } - - # Calculate final set of matching leaves - if (nrow(cat_conds) > 0) { - relevant_leaves <- merge(relevant_leaves_cnt, relevant_leaves_cat, by = c("c_idx", "f_idx"))[,.(c_idx, f_idx)] - } else { - relevant_leaves <- relevant_leaves_cnt[,.(c_idx, f_idx)] - } - - # If no cnt conditions exist, output empty update table cnt_new for cnt params - } else { - relevant_leaves <- relevant_leaves_cat[,.(c_idx, f_idx)] - cnt_new <- cbind(cnt[FALSE,], data.table(cvg_factor = numeric(), c_idx = integer(), val = numeric(), prob = numeric())) - } - - # Calculate updates for cat params matching cat conditions - cat_new <- merge(merge(relevant_leaves, cat_conds, by = "c_idx", allow.cartesian = TRUE), cat, by = c("f_idx","variable", "val")) - - # Ensure probabilities sum to 1 - cat_new[, cvg_factor := sum(prob), by = .(f_idx, c_idx, variable)] - cat_new[, prob := prob/cvg_factor] - - list(cnt_new = cnt_new, cat_new = cat_new, relevant_leaves = relevant_leaves) - } - if (isTRUE(parallel)) { - updates_relevant_leaves <- foreach(step = 1:step_no, .combine = "rbind") %dopar% update_fun(step) + # Per-step work lives once in arf_cforde_step() (cforde_workers.R); this closure + # adapts it to foreach's one-argument iteration. + par_fun <- function(step_) { + arf_cforde_step(step_, condition_long, conds_conditioned, nconds_conditioned, + stepsize, cat_cols, cnt_cols, params, family) + } + # Parallelism is across condition steps: 1 step is inherently serial. mirai only + # for step_no > 1. Under "separate" evidence_row_mode the caller (forge/expct) + # already parallelizes and passes parallel = FALSE, so no daemon-in-daemon. + use_mirai <- FALSE + if (step_no > 1) { + backend <- arf_select_backend(parallel) + use_mirai <- identical(backend, "mirai") + } + if (use_mirai) { + arf_load_on_daemons() # daemons call arf/data.table internals via bare names + params_shared <- mori::share(params) + condition_long_shared <- mori::share(condition_long) + # arf_cforde_step returns list(cnt_new, cat_new, relevant_leaves), so combine + # each component across steps (arf_mirai_tree_map assumes data.table results). + res <- mirai::mirai_map( + seq_len(step_no), + arf_cforde_step, + .args = list(condition_long = condition_long_shared, + conds_conditioned = mori::share(conds_conditioned), + nconds_conditioned = nconds_conditioned, stepsize = stepsize, + cat_cols = cat_cols, cnt_cols = cnt_cols, + params = params_shared, family = family))[] + arf_stop_on_mirai_error(res) + updates_relevant_leaves <- list( + cnt_new = rbindlist(lapply(res, `[[`, "cnt_new")), + cat_new = rbindlist(lapply(res, `[[`, "cat_new")), + relevant_leaves = rbindlist(lapply(res, `[[`, "relevant_leaves"))) + } else if (isTRUE(parallel) && step_no > 1) { + updates_relevant_leaves <- foreach(step = 1:step_no, .combine = "rbind") %dopar% par_fun(step) } else { - updates_relevant_leaves <- foreach(step = 1:step_no, .combine = "rbind") %do% update_fun(step) + updates_relevant_leaves <- foreach(step = 1:step_no, .combine = "rbind") %do% par_fun(step) } # Combine results diff --git a/man/adversarial_rf.Rd b/man/adversarial_rf.Rd index c29943c7..710edc9e 100644 --- a/man/adversarial_rf.Rd +++ b/man/adversarial_rf.Rd @@ -40,8 +40,10 @@ round to the next?} \item{verbose}{Print discriminator accuracy after each round? Will also show additional warnings.} -\item{parallel}{Compute in parallel? Must register backend beforehand, e.g. -via \code{doParallel} or \code{doFuture}; see examples.} +\item{parallel}{Compute in parallel? Enables multithreaded ranger training +(no backend needed) and parallelizes the pruning step, which requires a +registered \code{foreach} backend (\code{doParallel}, \code{doFuture}) or +active \code{mirai} daemons. See \code{\link{arf-options}}.} \item{...}{Extra parameters to be passed to \code{ranger}.} } @@ -111,6 +113,9 @@ doParallel::registerDoParallel(cores = 4) # ... or with doFuture doFuture::registerDoFuture() future::plan("multisession", workers = 4) + +# ... or with mirai (shares large read-only inputs across workers via mori) +mirai::daemons(4) } } diff --git a/man/arf-options.Rd b/man/arf-options.Rd new file mode 100644 index 00000000..b8775085 --- /dev/null +++ b/man/arf-options.Rd @@ -0,0 +1,108 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/arf-options.R +\name{arf-options} +\alias{arf-options} +\alias{arf.backend} +\alias{arf.verbose} +\alias{arf.block_rows} +\title{arf package options} +\description{ +Options controlling the parallel backend and its messaging, set via +\code{\link{options}}. +} +\details{ +\describe{ +\item{\code{arf.backend}}{Parallel backend used when \code{parallel = TRUE}: +\code{"foreach"} or \code{"mirai"}. If unset, arf uses \code{"mirai"} when +mirai daemons are running and \code{"foreach"} otherwise.} +\item{\code{arf.verbose}}{Report the selected backend once per backend +configuration per session? Default \code{TRUE}; set \code{FALSE} to +silence.} +\item{\code{arf.block_rows}}{Cap on rows materialized per block of +conditions in \code{\link{expct}}. Default \code{5e6}. Lower it to +trade speed for memory on large forests with many conditions.} +} + +\code{arf.block_rows} does not affect results, only peak memory and speed. +For memory-constrained hardware, combine a small daemon count with a lower +\code{arf.block_rows} and the \code{batch}/\code{stepsize} arguments of +\code{\link{lik}}, \code{\link{forge}} and \code{\link{expct}}. + +The \code{"foreach"} backend uses whatever adapter is registered (e.g. +\code{doParallel}, \code{doFuture}). The \code{"mirai"} backend uses +\code{mirai} daemons and shares large read-only inputs (training data, +forest, learned parameters) across workers via \code{mori}, so workers do +not each copy them. Speed is comparable between the backends. The memory +benefit is largest for the tree-parallel operations (\code{\link{forde}}, +\code{\link{adversarial_rf}}) on large forests with many workers, where +\code{foreach} memory grows with the worker count and \code{mirai} stays +much flatter (roughly half to a third at 16 workers in internal benchmarks). +For small workloads the daemon pool adds a fixed overhead that can outweigh +the sharing, so prefer \code{"mirai"} at scale and either backend otherwise. +Daemons started via \code{future.mirai} (e.g. +\code{plan(future.mirai::mirai_multisession)}) are detected like any other +mirai daemons, so futureverse users get the \code{"mirai"} backend +automatically. All backend packages are in Suggests; install the ones you +use. + +Reproducibility of stochastic operations (\code{\link{forge}}, categorical +\code{\link{expct}}) under parallel execution: \code{\link{set.seed}} only +governs the calling process, not the workers. With the \code{"mirai"} +backend, seed the daemons instead: \code{mirai::daemons(n, seed = 42)} +gives reproducible results provided the daemon count, the seed and the +sequence of calls on a fresh daemon pool are kept fixed (changing the +daemon count changes how work is chunked and therefore the random stream +assignment). This also applies to pools started via \code{future.mirai}: +\code{future}'s own seed machinery covers only work dispatched through its +API, which arf's backend bypasses, and \code{plan()} does not accept a +\code{seed} argument, so seed the pool with \code{mirai::daemons()} +directly. For the \code{"foreach"} backend, register \code{doRNG} on top of +the adapter (\code{doRNG::registerDoRNG(42)} after +\code{registerDoParallel()} or \code{registerDoFuture()}): results are then +reproducible, identical across adapters, and independent of the worker +count, provided \code{stepsize} is set explicitly (its default depends on +the worker count). Without \code{doRNG}, \code{doFuture} flags the +stochastic operations with "UNRELIABLE VALUE" warnings. Sequential +execution (\code{parallel = FALSE}) with \code{set.seed} is exact as +always. +} +\examples{ +\dontrun{ +arf <- adversarial_rf(iris) + +# foreach backend +doParallel::registerDoParallel(cores = 4) +psi <- forde(arf, iris) + +# mirai backend: start daemons, then call as usual +mirai::daemons(4) +psi <- forde(arf, iris) +mirai::daemons(0) # shut down when done + +# futureverse: future.mirai daemons are detected automatically +future::plan(future.mirai::mirai_multisession, workers = 4) +psi <- forde(arf, iris) + +# force a backend regardless of what is registered +options(arf.backend = "mirai") + +# silence the backend message +options(arf.verbose = FALSE) + +# reproducible parallel sampling with mirai: seeded daemons +# (fixed daemon count, fresh pool) +evi <- data.frame(Species = sample(levels(iris$Species), 100, replace = TRUE)) +mirai::daemons(4, seed = 42) +# stepsize = evidence rows per step; 25 = 100 conditions / 4 workers, +# i.e. the default sizing, made explicit +x <- forge(psi, n_synth = 1, evidence = evi, stepsize = 25) +mirai::daemons(0) + +# reproducible parallel sampling with foreach: doRNG on top of the +# adapter; set stepsize explicitly (its default depends on worker count) +doParallel::registerDoParallel(cores = 4) +doRNG::registerDoRNG(42) +x <- forge(psi, n_synth = 1, evidence = evi, stepsize = 25) +} + +} diff --git a/man/arf-package.Rd b/man/arf-package.Rd index be56c58b..37f8487e 100644 --- a/man/arf-package.Rd +++ b/man/arf-package.Rd @@ -38,6 +38,9 @@ doParallel::registerDoParallel(cores = 4) # ... or with doFuture doFuture::registerDoFuture() future::plan("multisession", workers = 4) + +# ... or with mirai (shares large read-only inputs across workers via mori) +mirai::daemons(4) } } \seealso{ @@ -56,9 +59,15 @@ Useful links: Authors: \itemize{ + \item Marvin N. Wright \email{cran@wrig.de} (\href{https://orcid.org/0000-0002-8542-6291}{ORCID}) \item David S. Watson \email{david.s.watson11@gmail.com} (\href{https://orcid.org/0000-0001-9632-2159}{ORCID}) \item Kristin Blesch (\href{https://orcid.org/0000-0001-6241-3079}{ORCID}) \item Jan Kapar (\href{https://orcid.org/0009-0000-6408-2840}{ORCID}) } +Other contributors: +\itemize{ + \item Lukas Burk \email{cran@lukasburk.de} (\href{https://orcid.org/0000-0001-7528-3795}{ORCID}) [contributor] +} + } diff --git a/man/cforde.Rd b/man/cforde.Rd index ea856e7e..970b2b1b 100644 --- a/man/cforde.Rd +++ b/man/cforde.Rd @@ -29,8 +29,9 @@ Options are to force sampling from a random leaf (\code{"force"}) or return \item{stepsize}{Stepsize defining number of condition rows handled in one for each step.} -\item{parallel}{Compute in parallel? Must register backend beforehand, e.g. -via \code{doParallel} or \code{doFuture}; see examples.} +\item{parallel}{Compute in parallel? Requires a registered \code{foreach} +backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +daemons. See \code{\link{arf-options}}.} } \value{ List with conditions (\code{evidence_input}), prepared conditions (\code{evidence_prepped}) diff --git a/man/expct.Rd b/man/expct.Rd index 40759a54..393788d3 100644 --- a/man/expct.Rd +++ b/man/expct.Rd @@ -45,11 +45,17 @@ Options are to force sampling from a random leaf (\code{"force"}) or return \item{verbose}{Show warnings, e.g. when no leaf matches a condition?} \item{stepsize}{How many rows of evidence should be handled at each step? -Defaults to \code{nrow(evidence) / num_registered_workers} for +Defaults to \code{nrow(evidence)} divided by the number of registered +workers or daemons for \code{parallel == TRUE}.} -\item{parallel}{Compute in parallel? Must register backend beforehand, e.g. -via \code{doParallel} or \code{doFuture}; see Examples.} +\item{parallel}{Compute in parallel? Requires a registered \code{foreach} +backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +daemons. See \code{\link{arf-options}}. With +\code{evidence_row_mode = "or"}, parallelization happens inside the +conditional circuit computation; in benchmarks this gave little speedup +while raising peak memory, so consider \code{parallel = FALSE} for large +\code{"or"} queries.} } \value{ A one row data frame with values for all query variables. @@ -108,6 +114,9 @@ doParallel::registerDoParallel(cores = 4) # ... or with doFuture doFuture::registerDoFuture() future::plan("multisession", workers = 4) + +# ... or with mirai (shares large read-only inputs across workers via mori) +mirai::daemons(4) } } diff --git a/man/forde.Rd b/man/forde.Rd index ef58eaf8..621c1820 100644 --- a/man/forde.Rd +++ b/man/forde.Rd @@ -46,8 +46,9 @@ on multinomial likelihoods.} data fall outside the support of training data. The gap between lower and upper bounds is expanded by a factor of \code{1 + epsilon}.} -\item{parallel}{Compute in parallel? Must register backend beforehand, e.g. -via \code{doParallel} or \code{doFuture}; see examples.} +\item{parallel}{Compute in parallel? Requires a registered \code{foreach} +backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +daemons. See \code{\link{arf-options}}.} } \value{ A \code{list} with 5 elements: (1) parameters for continuous data; (2) @@ -106,6 +107,9 @@ doParallel::registerDoParallel(cores = 4) # ... or with doFuture doFuture::registerDoFuture() future::plan("multisession", workers = 4) + +# ... or with mirai (shares large read-only inputs across workers via mori) +mirai::daemons(4) } diff --git a/man/forge.Rd b/man/forge.Rd index d9c3a479..907d20a5 100644 --- a/man/forge.Rd +++ b/man/forge.Rd @@ -45,11 +45,17 @@ Options are to force sampling from a random leaf (\code{"force"}) or return \item{verbose}{Show warnings, e.g. when no leaf matches a condition?} \item{stepsize}{How many rows of evidence should be handled at each step? -Defaults to \code{nrow(evidence) / num_registered_workers} for +Defaults to \code{nrow(evidence)} divided by the number of registered +workers or daemons for \code{parallel == TRUE}.} -\item{parallel}{Compute in parallel? Must register backend beforehand, e.g. -via \code{doParallel} or \code{doFuture}; see examples.} +\item{parallel}{Compute in parallel? Requires a registered \code{foreach} +backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +daemons. See \code{\link{arf-options}}. With +\code{evidence_row_mode = "or"}, parallelization happens inside the +conditional circuit computation; in benchmarks this gave little speedup +while raising peak memory, so consider \code{parallel = FALSE} for large +\code{"or"} queries.} } \value{ A dataset of \code{n_synth} synthetic samples. @@ -119,6 +125,9 @@ doParallel::registerDoParallel(cores = 4) # ... or with doFuture doFuture::registerDoFuture() future::plan("multisession", workers = 4) + +# ... or with mirai (shares large read-only inputs across workers via mori) +mirai::daemons(4) } } diff --git a/man/impute.Rd b/man/impute.Rd index 87897085..9883cc4a 100644 --- a/man/impute.Rd +++ b/man/impute.Rd @@ -72,6 +72,9 @@ doParallel::registerDoParallel(cores = 4) # ... or with doFuture doFuture::registerDoFuture() future::plan("multisession", workers = 4) + +# ... or with mirai (shares large read-only inputs across workers via mori) +mirai::daemons(4) } } diff --git a/man/lik.Rd b/man/lik.Rd index 70233f93..b2122270 100644 --- a/man/lik.Rd +++ b/man/lik.Rd @@ -42,8 +42,9 @@ queries in one round, which is always the fastest option if memory allows. However, with large samples or many trees, it can be more memory efficient to split the data into batches. This has no impact on results.} -\item{parallel}{Compute in parallel? Must register backend beforehand, e.g. -via \code{doParallel} or \code{doFuture}; see examples.} +\item{parallel}{Compute in parallel? Requires a registered \code{foreach} +backend (\code{doParallel}, \code{doFuture}) or active \code{mirai} +daemons. See \code{\link{arf-options}}.} } \value{ A vector of likelihoods, optionally on the log scale. @@ -97,6 +98,9 @@ doParallel::registerDoParallel(cores = 4) # ... or with doFuture doFuture::registerDoFuture() future::plan("multisession", workers = 4) + +# ... or with mirai (shares large read-only inputs across workers via mori) +mirai::daemons(4) } } diff --git a/tests/testthat/helper-mirai.R b/tests/testthat/helper-mirai.R new file mode 100644 index 00000000..562b34ce --- /dev/null +++ b/tests/testthat/helper-mirai.R @@ -0,0 +1,20 @@ +# Set up mirai daemons for backend tests. In CI the package is installed, so the +# package's own arf_load_on_daemons() (loadNamespace) suffices. Under load_all +# (devtools::test()) daemons can't see the source build, so dev-load arf on them +# here; the package's loadNamespace call is then a no-op. +setup_mirai_daemons <- function(n = 2, seed = NULL) { + # mirai/mori are in Suggests, so skip_if_not_installed() does NOT skip on + # CRAN; skip here instead -- daemon processes would brush against CRAN's + # 2-core policy and are a known source of flakiness on its builders. + testthat::skip_on_cran() + if (is.null(seed)) { + mirai::daemons(n) + } else { + mirai::daemons(n, seed = seed) + } + if (requireNamespace("pkgload", quietly = TRUE) && + isTRUE(tryCatch(pkgload::is_dev_package("arf"), error = function(e) FALSE))) { + pdir <- pkgload::pkg_path() + mirai::everywhere(suppressMessages(pkgload::load_all(pdir, quiet = TRUE)), pdir = pdir) + } +} diff --git a/tests/testthat/test-mirai-backend.R b/tests/testthat/test-mirai-backend.R new file mode 100644 index 00000000..5123970b --- /dev/null +++ b/tests/testthat/test-mirai-backend.R @@ -0,0 +1,341 @@ +# Correctness gate for the mirai+mori backend: +# equality vs the sequential/foreach paths on one dataset. + +test_that("mirai backend produces equal forde output on iris", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + + old <- options(arf.backend = "foreach") + on.exit(options(old), add = TRUE) + psi_foreach <- forde(arf, iris, parallel = FALSE) + + options(arf.backend = "mirai") + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + psi_mirai <- forde(arf, iris, parallel = TRUE) + + expect_equal(psi_mirai$cnt, psi_foreach$cnt, ignore_attr = TRUE) + expect_equal(psi_mirai$cat, psi_foreach$cat, ignore_attr = TRUE) + expect_equal(psi_mirai$forest, psi_foreach$forest, ignore_attr = TRUE) +}) + +test_that("mirai backend errors clearly when daemons are not set", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + mirai::daemons(0) + old <- options(arf.backend = "mirai") + on.exit(options(old), add = TRUE) + + expect_error(forde(arf, iris, parallel = TRUE), "daemons") +}) + +test_that("mirai backend produces structurally consistent forge() output", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + # forge() is stochastic, so compare structure (not values) across backends. + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + psi <- forde(arf, iris, parallel = FALSE) + # 21 separate exact conditions (mixed species) -> multi-step + evi <- iris[c(1:7, 51:57, 101:107), "Species", drop = FALSE] + + old <- options(arf.backend = "foreach") + on.exit(options(old), add = TRUE) + x_foreach <- forge(psi, n_synth = 3, evidence = evi, parallel = FALSE, + stepsize = 5, verbose = FALSE) + + options(arf.backend = "mirai") + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + x_mirai <- forge(psi, n_synth = 3, evidence = evi, parallel = TRUE, + stepsize = 5, verbose = FALSE) + + expect_equal(nrow(x_mirai), nrow(x_foreach)) + expect_equal(colnames(x_mirai), colnames(x_foreach)) + expect_equal(sapply(x_mirai, class), sapply(x_foreach, class)) + # exact conditions must map onto their output rows in evidence order: + # catches scrambled step-to-row assembly that structural checks miss + expect_equal(as.character(x_mirai$Species), + as.character(rep(evi$Species, each = 3))) +}) + +test_that("mirai backend gives identical lik() (deterministic)", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + psi <- forde(arf, iris, parallel = FALSE) + + old <- options(arf.backend = "foreach") + on.exit(options(old), add = TRUE) + ll_foreach <- lik(psi, iris, arf = arf, batch = 30, parallel = FALSE) + + options(arf.backend = "mirai") + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + ll_mirai <- lik(psi, iris, arf = arf, batch = 30, parallel = TRUE) + + expect_equal(ll_mirai, ll_foreach) # lik is deterministic +}) + +test_that("mirai backend gives structurally consistent expct()", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + psi <- forde(arf, iris, parallel = FALSE) + evi <- iris[1:20, "Species", drop = FALSE] # 20 separate conditions -> multi-step + + old <- options(arf.backend = "foreach") + on.exit(options(old), add = TRUE) + x_foreach <- expct(psi, evidence = evi, parallel = FALSE, stepsize = 5, verbose = FALSE) + + options(arf.backend = "mirai") + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + x_mirai <- expct(psi, evidence = evi, parallel = TRUE, stepsize = 5, verbose = FALSE) + + expect_equal(dim(x_mirai), dim(x_foreach)) + expect_equal(colnames(x_mirai), colnames(x_foreach)) +}) + +test_that("arf_n_workers reflects active mirai daemons (stepsize sizing)", { + skip_if_not_installed("mirai") + skip_on_cran() + + mirai::daemons(0) + n_idle <- arf_n_workers() # no mirai, no foreach -> 1 + expect_equal(n_idle, 1L) + + mirai::daemons(3) + on.exit(mirai::daemons(0), add = TRUE) + expect_gte(arf_n_workers(), 3L) # must see daemons, else step_no==1 kills mirai +}) + +test_that("mirai backend gives identical cforde() (deterministic)", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + psi <- forde(arf, iris, parallel = FALSE) + set.seed(1) + evi <- data.frame(Sepal.Length = runif(20, 4.5, 7)) # 20 conditions -> multi-step + + old <- options(arf.backend = "foreach") + on.exit(options(old), add = TRUE) + cf_foreach <- arf:::cforde(psi, evi, stepsize = 5, parallel = FALSE, verbose = FALSE) + + options(arf.backend = "mirai") + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + cf_mirai <- arf:::cforde(psi, evi, stepsize = 5, parallel = TRUE, verbose = FALSE) + + expect_equal(cf_mirai$forest, cf_foreach$forest) # cforde is deterministic + expect_equal(cf_mirai$cnt, cf_foreach$cnt) + expect_equal(cf_mirai$cat, cf_foreach$cat) +}) + +test_that("mirai backend preserves row order and class in expct() (regression)", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + # step_no (8) > n_workers (2): guards against interleaved-chunk row scrambling + # and data.table-vs-data.frame class drift in arf_mirai_tree_map. + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + psi <- forde(arf, iris, parallel = FALSE) + evi <- data.frame(Sepal.Length = seq(4.5, 7, length.out = 8)) + + old <- options(arf.backend = "foreach") + on.exit(options(old), add = TRUE) + x_foreach <- expct(psi, query = "Petal.Length", evidence = evi, + parallel = FALSE, stepsize = 1, verbose = FALSE) + + options(arf.backend = "mirai") + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + x_mirai <- expct(psi, query = "Petal.Length", evidence = evi, + parallel = TRUE, stepsize = 1, verbose = FALSE) + + expect_equal(x_mirai, x_foreach) # exact: order + values + class + row.names +}) + +test_that("mirai backend gives identical adversarial_rf() pruning (deterministic)", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + # ranger draws per-tree seeds from the R RNG, so training is reproducible + # under set.seed() regardless of threading; prune is deterministic given a + # forest. That makes the full run comparable end to end, exercising the + # actual mirai prune dispatch in adversarial_rf(). + set.seed(7) + a_serial <- adversarial_rf(iris, num_trees = 50, parallel = FALSE, + verbose = FALSE) + + old <- options(arf.backend = "mirai") + on.exit(options(old), add = TRUE) + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + set.seed(7) + a_mirai <- adversarial_rf(iris, num_trees = 50, parallel = TRUE, + verbose = FALSE) + + expect_identical(a_mirai$forest$child.nodeIDs, a_serial$forest$child.nodeIDs) +}) + +test_that("mirai backend runs adversarial_rf() end to end", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + old <- options(arf.backend = "mirai") + on.exit(options(old), add = TRUE) + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + + a <- adversarial_rf(iris, num_trees = 30, parallel = TRUE, verbose = FALSE) + expect_s3_class(a, "ranger") + expect_length(a$forest$child.nodeIDs, a$num.trees) +}) + +test_that("arf.block_rows caps expct blocks without changing results", { + a <- adversarial_rf(iris, num_trees = 10, parallel = FALSE, verbose = FALSE) + psi <- forde(a, iris, parallel = FALSE) + evi <- data.frame(Species = sample(levels(iris$Species), 6, replace = TRUE)) + evi1 <- evi[1, , drop = FALSE] + + set.seed(1) + ref <- expct(psi, evidence = evi, parallel = FALSE) + set.seed(2) + ref1 <- expct(psi, evidence = evi1, parallel = FALSE) + + old <- options(arf.block_rows = 1) # force one condition per block + on.exit(options(old), add = TRUE) + set.seed(1) + blocked <- expct(psi, evidence = evi, parallel = FALSE) + expect_identical(blocked, ref) + + # single condition: cannot be split, must take the unblocked path + set.seed(2) + blocked1 <- expct(psi, evidence = evi1, parallel = FALSE) + expect_identical(blocked1, ref1) +}) + +test_that("arf_tree_chunks yields contiguous in-order blocks covering all trees", { + for (nt in c(1L, 5L, 7L)) { + for (nw in c(1L, 3L, 10L)) { + ch <- arf_tree_chunks(nt, nw) + # concatenating chunks in chunk order must reproduce 1..nt exactly: + # this is the contiguity invariant positional ops rely on + expect_identical(unlist(ch, use.names = FALSE), seq_len(nt)) + expect_lte(length(ch), max(1L, min(nw, nt))) + } + } +}) + +test_that("arf_select_backend applies the documented precedence", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + skip_on_cran() + + old <- options(arf.backend = NULL, arf.verbose = FALSE) + on.exit(options(old), add = TRUE) + + expect_identical(arf_select_backend(FALSE), "sequential") + + mirai::daemons(0) + expect_identical(arf_select_backend(TRUE), "foreach") # no daemons, option unset + + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + expect_identical(arf_select_backend(TRUE), "mirai") # daemons auto-detected + + options(arf.backend = "foreach") + expect_identical(arf_select_backend(TRUE), "foreach") # explicit option wins + + options(arf.backend = "bogus") + expect_error(arf_select_backend(TRUE)) +}) + +test_that("arf_load_on_daemons caches per daemon pool and self-invalidates", { + skip_if_not_installed("mirai") + + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + assign("arf_loaded_key", NULL, envir = arf:::.arf_env) + + expect_true(arf_load_on_daemons()) # first call loads + expect_false(arf_load_on_daemons()) # same pool -> cached + + mirai::daemons(0) + setup_mirai_daemons(2) # rebuilt pool mints a new key + expect_true(arf_load_on_daemons()) +}) + +test_that("mirai worker errors propagate instead of corrupting results", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + setup_mirai_daemons(2) + on.exit(mirai::daemons(0), add = TRUE) + + boom <- function(tree) stop("boom") + environment(boom) <- globalenv() + expect_error(arf_mirai_tree_map(4, boom, list()), "boom") +}) + +test_that("seeded daemons make stochastic mirai results reproducible", { + skip_if_not_installed("mirai") + skip_if_not_installed("mori") + + # set.seed() only governs the calling process; the documented recipe for + # reproducible forge()/expct() under mirai is daemons(n, seed = ) with a + # fixed daemon count and call sequence on a fresh pool (see ?arf-options). + arf <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + psi <- forde(arf, iris, parallel = FALSE) + evi <- iris[c(1:7, 51:57, 101:107), "Species", drop = FALSE] + + old <- options(arf.backend = "mirai") + on.exit(options(old), add = TRUE) + on.exit(mirai::daemons(0), add = TRUE) + + run_once <- function() { + setup_mirai_daemons(2, seed = 42) + x <- forge(psi, n_synth = 3, evidence = evi, parallel = TRUE, + stepsize = 5, verbose = FALSE) + mirai::daemons(0) + x + } + expect_identical(run_once(), run_once()) +}) + +test_that("foreach doParallel backend gives equal forde output", { + skip_if_not_installed("doParallel") + skip_on_cran() + + a <- adversarial_rf(iris, verbose = FALSE, parallel = FALSE) + psi_seq <- forde(a, iris, parallel = FALSE) + + # PSOCK, not fork: forking after mirai/nanonext threads exist is unsafe + cl <- parallel::makeCluster(2) + on.exit(parallel::stopCluster(cl), add = TRUE) + if (requireNamespace("pkgload", quietly = TRUE) && + isTRUE(tryCatch(pkgload::is_dev_package("arf"), error = function(e) FALSE))) { + pdir <- pkgload::pkg_path() + parallel::clusterCall(cl, function(p) { + suppressMessages(pkgload::load_all(p, quiet = TRUE)) + }, pdir) + } + doParallel::registerDoParallel(cl) + on.exit(foreach::registerDoSEQ(), add = TRUE) + old <- options(arf.backend = "foreach") + on.exit(options(old), add = TRUE) + + psi_par <- forde(a, iris, parallel = TRUE) + + expect_equal(psi_par$cnt, psi_seq$cnt, ignore_attr = TRUE) + expect_equal(psi_par$cat, psi_seq$cat, ignore_attr = TRUE) + expect_equal(psi_par$forest, psi_seq$forest, ignore_attr = TRUE) +}) diff --git a/vignettes/arf.Rmd b/vignettes/arf.Rmd index 738ac8b9..f38210f6 100644 --- a/vignettes/arf.Rmd +++ b/vignettes/arf.Rmd @@ -76,7 +76,41 @@ In either case, we can now execute in parallel. arf_iris <- adversarial_rf(iris, num_trees = 100) ``` -The result is an object of class `ranger`, which we can input to downstream functions. +As an alternative to `foreach`, arf can use the `mirai` backend. Start daemons and arf uses them automatically, with no `foreach` registration needed. Large read-only inputs (training data, forest, learned parameters) are shared across workers via `mori`, so workers do not each copy them. Speed is comparable between the backends. The memory benefit is largest for the tree-parallel operations (`forde`, `adversarial_rf`) on large forests with many workers, where `foreach` memory grows with the worker count and `mirai` stays much flatter. For small workloads the daemon pool adds a fixed overhead that can outweigh the sharing, so prefer `mirai` at scale and either backend otherwise. + +```{r mirai, eval=FALSE} +# mirai backend: start daemons, then call arf as usual +library(mirai) +daemons(2) +arf_iris <- adversarial_rf(iris, num_trees = 100) +daemons(0) # shut down when done +``` + +If both are available, set `options(arf.backend = "foreach")` or `"mirai"` to choose explicitly. See `?"arf-options"` for details. The same backend is used by `forde`, `forge`, `expct`, and `lik`. Daemons started via `future.mirai` (e.g. `plan(future.mirai::mirai_multisession)`) are detected automatically, so futureverse users get the `mirai` backend without extra setup. + +Stochastic operations (`forge`, categorical `expct`) are not governed by `set.seed()` under parallel execution. Both backends offer a reproducible setup: + +```{r repro, eval=FALSE} +# params_iris is the forde() output introduced in the next section; +# stepsize is the number of evidence rows per parallel step (default: +# evidence rows / workers), made explicit here for reproducibility +evi <- data.frame(Species = sample(levels(iris$Species), 100, replace = TRUE)) + +# mirai: seeded daemons; keep the daemon count fixed +daemons(2, seed = 42) +x_synth <- forge(params_iris, n_synth = 1, evidence = evi, stepsize = 25) +daemons(0) + +# foreach: doRNG on top of the adapter; set stepsize explicitly +# (its default depends on the worker count) +doParallel::registerDoParallel(2) +doRNG::registerDoRNG(42) +x_synth <- forge(params_iris, n_synth = 1, evidence = evi, stepsize = 25) +``` + +**Nested parallelism**: ARF splits work across workers, and `data.table` also threads within each worker. These multiply (`workers × threads`), so a product above your core count oversubscribes, which slows things down and inflates load average, since `data.table`'s idle OpenMP threads busy-wait by default. If parallel performance disappoints, cap per-worker threads with `data.table::setDTthreads()` and set `OMP_WAIT_POLICY=passive` before starting R so idle threads sleep. The best split of workers and threads depends on the data and hardware. + +The result is an object of class `ranger`, which we can input to downstream functions. # Parameter Learning