witchtools 0.6.0: fast region-conversion engine + japan/korea/southafrica region sets

.Rbuildignore

@@ -10,3 +10,4 @@
 ^\.github$
 ^CLAUDE\.md$
 ^\.claude$
+^benchmarks$

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: witchtools
 Title: Data Management for Integrated Assessment Models
-Version: 0.5.1.9000
+Version: 0.6.0
 Authors@R: 
     c(person(given = "Laurent",
            family = "Drouet",

NAMESPACE

@@ -19,12 +19,15 @@ export(iamc_native_regions_yml)
 export(iamc_region_mappings_yml)
 export(india_regions)
 export(indonesia_regions)
+export(japan_regions)
+export(korea_regions)
 export(oceania_regions)
 export(oecd_regions)
 export(premise_region_mapping)
 export(region_id)
 export(require_gdxtools)
 export(require_package)
+export(southafrica_regions)
 export(ssa_regions)
 export(usa_regions)
 export(witch_data)

NEWS.md

@@ -1,4 +1,28 @@
-# witchtools (development version)
+# witchtools 0.6.0
+
+## New features
+
+- New region sets `japan_regions()`, `korea_regions()` and
+  `southafrica_regions()`, emitted as `is_japan`, `is_korea` and
+  `is_southafrica` in the generated `regions.inc`. Selection uses the same
+  GDP-majority rule as the other region sets, so under `witch17` the
+  `jpnkor` region belongs to `is_japan` and `is_korea` is empty.
+
+## Performance
+
+- Region-to-region conversion (`convert_region()`, and `convert_table()` /
+  `convert_gdx()` / `convert_duckdb()` / `convert_sqlite()` through it) no
+  longer expands the data to country level. A new engine converts through a
+  precomputed region-pair coefficient table: on 1M-row tables it is ~5x
+  faster with up to ~10x lower peak memory, and it converts the 2.1M-row
+  GLOBIOM reporting tables witch17->witch20 in ~6s within 1 GB where the
+  previous implementation exhausted 6 GB (see `benchmarks/results/`).
+  Results are identical up to floating-point summation order (< 1e-12).
+  The previous engine remains available with
+  `options(witchtools.convert_region_engine = "legacy")`; country-level
+  (iso3) input and the `set1` operator always use it (`set1`'s rounding is
+  discontinuous, so summation reassociation could flip values at .5).
+- `subset()` replaced by direct indexing in the conversion hot paths.
 
 # witchtools 0.5.1
 

R/convert_duckdb.R

@@ -113,7 +113,7 @@ convert_duckdb <- function(duckdb,
 
     # add an additional parameter '_info' when using sumby
     if (!is.null(.info_share)) {
-      indices <- subset(colnames(.data), colnames(.data) != "value")
+      indices <- colnames(.data)[colnames(.data) != "value"]
       data.table::setcolorder(.info_share, data_indices)
       names(.info_share) <- c(indices, "value")
       .i <- list(.info_share)

R/convert_gdx.R

@@ -109,7 +109,7 @@ convert_gdx <- function(gdxfile,
     if (length(colnames(.data)) == 1) {
       names(.data) <- "value"
     } else {
-      indices <- subset(colnames(.data), colnames(.data) != "value")
+      indices <- colnames(.data)[colnames(.data) != "value"]
       indices <- ifelse(indices %in% c(region_name, "t"), indices, "*")
       names(.data) <- c(indices, "value")
     }

R/convert_region.R

@@ -20,6 +20,15 @@
 #' and another one named as the regional mapping (for region name).
 #' The name in the list should also be the regional mapping name.
 #'
+#' Region-to-region conversions use a fast engine that converts through a
+#' precomputed region-pair coefficient table instead of expanding the data to
+#' country level, which drastically reduces memory use and run time on large
+#' tables. Results are identical up to floating-point summation order
+#' (relative differences below 1e-12). The previous implementation remains
+#' available with \code{options(witchtools.convert_region_engine = "legacy")}.
+#' Country-level (iso3) input and the \code{set1} operator always use the
+#' legacy engine.
+#'
 #' @family conversion functions
 #' @seealso \code{\link{convert_table}},
 #' \code{\link{convert_gdx}}.
@@ -129,123 +138,17 @@ convert_region <- function(.x,
     stop(paste0("to_reg == iso3 is not yet implemented."))
   }
 
-  # Add iso3 and data_reg mapping
-  if (rname0 == "iso3") {
-    .x <- merge(.x, rmap1, by = "iso3")
-  } else {
-    .r <- merge(rmap0, rmap1, by = "iso3")
-    .x <- merge(.x, .r, by = rname0, allow.cartesian = TRUE)
-  }
-
-  # Add weight
-  .x <- merge(.x, agg_weight, by = "iso3")
-  .x <- .x[!is.na(get(rname1))]
-
-  dkeys <- function(dd) {
-    return(c(colnames(dd)[!colnames(dd) %in% c(
-      "value",
-      "weight", "sum_weight",
-      "iso3", rname0, rname1
-    )]))
-  }
-
-  # Disaggregation
-  if (rname0 != "iso3") {
-    if (agg_operator %in% c("sum","sumby")) {
-      # total weights are computed because of missing zeros values
-      .w <- merge(rmap0, agg_weight, by = "iso3")
-      .w <- .w[iso3 %in% unique(.x$iso3)]
-      .w <- .w[, list(sum_weight = sum(weight)), by = rname0]
-      .x <- merge(.x, .w, by = rname0)
-      .x <- .x[, list(iso3,
-        rname1 = get(rname1),
-        value = value * weight / sum_weight,
-        weight
-      ),
-      by = c(dkeys(.x), rname0)
-      ]
-      if (agg_operator %in% c("sum")) {
-        .x[, weight := NULL]
-      }
-    } else {
-      if (agg_operator %in% c("mean", "set1", "min", "minw", "max", "maxw")) {
-        .x <- .x[, .(iso3,
-          rname1 = get(rname1),
-          value,
-          weight
-        ),
-        by = c(dkeys(.x), rname0)
-        ]
-      } else {
-        stop(paste("Operator ", agg_operator, "not implemented"))
-      }
-    }
-  } else {
-    data.table::setnames(.x, rname1, "rname1")
-  }
-
-  # informed share
-  .info_share <- NULL
-  if (agg_operator %in% c("sumby")) {
-    .w <- merge(rmap1, agg_weight, by = "iso3")
-    .w <- .w[, .(sum_weight = sum(weight)), by = rname1]
-    data.table::setnames(.w, rname1, "rname1")
-    .x <- merge(.x, .w, by = "rname1")
-    .info_share <- .x[, .(value = sum(weight) / mean(sum_weight)),
-      by = c(dkeys(.x))
-    ]
-    data.table::setnames(.info_share, "rname1", rname1)
-  }
-
-  # Aggregation
-  if (agg_operator %in% c("sum", "sumby")) {
-    .x <- .x[, .(value = sum(value)), by = c(dkeys(.x))]
-  } else {
-    .w <- merge(rmap1, agg_weight, by = "iso3")
-    .w <- .w[, .(sum_weight = sum(weight)), by = rname1]
-    data.table::setnames(.w, rname1, "rname1")
-    .x <- merge(.x, .w, by = "rname1")
-    if (agg_operator == "mean") {
-      if (agg_missing == "zero") {
-        .x <- .x[, .(value = sum(value * weight / sum_weight)),
-          by = c(dkeys(.x))
-        ]
-      }
-      if (agg_missing == "NA") {
-        .x <- .x[!is.na(value), .(value = sum(value * weight / sum(weight))),
-          by = c(dkeys(.x))
-        ]
-      }
-    } else if (agg_operator == "set1") {
-      if (agg_missing == "zero") {
-        .x <- .x[, .(value = round(sum(value * weight / sum_weight))),
-          by = c(dkeys(.x))
-        ]
-      }
-      if (agg_missing == "NA") {
-        .x <- .x[, .(value = round(sum(value * weight / sum(weight)))),
-          by = c(dkeys(.x))
-        ]
-      }
-    } else if (agg_operator %in% c("min", "minw")) {
-      .x <- .x[, .(value = min(value[which(weight == min(weight))])),
-        by = c(dkeys(.x))
-      ]
-    } else if (agg_operator %in% c("max", "maxw")) {
-      .x <- .x[, .(value = max(value[which(weight == max(weight))])),
-        by = c(dkeys(.x))
-      ]
-    } else {
-      stop(paste("Operator", agg_operator, "not implemented"))
-    }
-  }
-
-  # Change the region column name
-  data.table::setnames(.x, "rname1", rname1)
-
-  if (info) {
-    return(list(data = .x, info = .info_share))
-  } else {
-    return(.x)
+  # Engine dispatch. The "fast" engine converts region->region input through
+  # a small pair-coefficient table instead of the iso3 explosion; iso3-level
+  # input is already linear-size and stays on the legacy engine. "set1" also
+  # stays on the legacy engine: its round() is discontinuous, so the summation
+  # reassociation of the fast engine could flip a value sitting exactly on a
+  # .5 boundary, and set1 tables are tiny anyway.
+  engine <- getOption("witchtools.convert_region_engine", "fast")
+  if (identical(engine, "fast") && rname0 != "iso3" && agg_operator != "set1") {
+    return(convert_region_fast(.x, rmap0, rname0, rmap1, rname1,
+                               agg_operator, agg_weight, agg_missing, info))
   }
+  convert_region_via_iso3(.x, rmap0, rname0, rmap1, rname1,
+                          agg_operator, agg_weight, agg_missing, info)
 }

R/convert_region_fast.R

@@ -0,0 +1,162 @@
+# Fast region-conversion engine.
+#
+# The legacy engine downscales the table to iso3 country level (~x15 rows for
+# witch17 input) and re-aggregates. But the downscale->upscale is a linear
+# map: for every operator the result only depends on per-(from, to) pair
+# aggregates of the country weights. This engine precomputes a tiny pair
+# coefficient table (21 rows for witch17->witch20) and converts with a single
+# join (~x1.2 rows) and one grouped aggregation, instead of the iso3 explosion
+# and its four full-size intermediate copies.
+#
+# Numerical results match the legacy engine up to floating-point summation
+# order (tolerance ~1e-12, covered by the equivalence tests in
+# tests/testthat/test-convert-region-equivalence.R).
+
+# Build the per-(from, to) coefficient table.
+#
+# Column semantics (mirroring the legacy engine's intermediate quantities):
+#   w_pair   sum of country weights in (from n to)      [legacy :141 join]
+#   w_min/w_max  min/max country weight in the pair     [for min*/max* ops]
+#   iso3_min smallest iso3 of the pair                  [row-order emulation]
+#   sw_from  sum of w_pair by from-region               [legacy :156-158]
+#   sw_to    sum of weights over the FULL to-mapping    [legacy :190/:204]
+#
+# The legacy sum denominator filters countries to those present in the data
+# (`iso3 %in% unique(.x$iso3)`), but since the explosion expands every present
+# region to all its member countries and denominators are grouped by
+# from-region, the filter is data-independent: sw_from computed here is
+# identical for every from-region that appears in the data at all.
+#' @noRd
+build_region_coeff <- function(rmap0, rname0, rmap1, rname1, agg_weight) {
+  iso3 <- weight <- w_pair <- .from <- sw_from <- NULL # due to NSE notes in R CMD check
+
+  pw <- merge(rmap0, rmap1, by = "iso3")
+  pw <- merge(pw, agg_weight, by = "iso3")
+  pw <- pw[!is.na(get(rname1))]
+  cf <- pw[, list(
+    w_pair = sum(weight),
+    w_min = min(weight),
+    w_max = max(weight),
+    iso3_min = min(iso3)
+  ), by = c(rname0, rname1)]
+  data.table::setnames(cf, c(rname0, rname1), c(".from", ".to"))
+  cf[, sw_from := sum(w_pair), by = .from]
+  sw_to <- merge(rmap1, agg_weight, by = "iso3")
+  sw_to <- sw_to[, list(sw_to = sum(weight)), by = rname1]
+  data.table::setnames(sw_to, rname1, ".to")
+  cf <- merge(cf, sw_to, by = ".to")
+  return(cf)
+}
+
+#' @noRd
+convert_region_fast <- function(.x,
+                                rmap0, rname0,
+                                rmap1, rname1,
+                                agg_operator,
+                                agg_weight,
+                                agg_missing,
+                                info) {
+
+  value <- w_pair <- sw_from <- sw_to <- NULL # due to NSE notes in R CMD check
+  .row0 <- .from <- .to <- iso3_min <- w_min <- w_max <- gw <- NULL
+
+  known_ops <- c("sum", "sumby", "mean", "min", "minw", "max", "maxw")
+  if (!agg_operator %in% known_ops) {
+    if (agg_operator == "set1") {
+      # Dispatched to the legacy engine by convert_region(): round() is
+      # discontinuous, so reassociated summation could flip .5 boundaries.
+      stop("set1 is handled by the legacy engine.")
+    }
+    # Same message and trigger point as the legacy disaggregation branch.
+    stop(paste("Operator ", agg_operator, "not implemented"))
+  }
+
+  cf <- build_region_coeff(rmap0, rname0, rmap1, rname1, agg_weight)
+
+  idcols <- setdiff(names(.x), c(rname0, "value"))
+
+  # Transient input-row tag for row-order emulation; removed on exit so the
+  # caller's table is left untouched even on error.
+  .x[, .row0 := .I]
+  on.exit(
+    if (".row0" %in% names(.x)) .x[, .row0 := NULL],
+    add = TRUE
+  )
+
+  # One small join instead of the iso3 explosion. nomatch=NULL drops input
+  # rows whose region is absent from the crosswalk/weights, as the legacy
+  # inner merges do.
+  xt <- cf[.x, on = c(".from" = rname0), allow.cartesian = TRUE, nomatch = NULL]
+
+  if (agg_operator %in% c("sum", "sumby")) {
+
+    # Legacy row order: groups appear in (from-region, input-row) order for
+    # sum, with an extra leading to-region sort for sumby (its merge at :193).
+    if (agg_operator == "sum") {
+      data.table::setorder(xt, .from, .row0, iso3_min)
+      out <- xt[, list(value = sum(value * w_pair / sw_from)),
+        by = c(idcols, ".to")
+      ]
+      data.table::setnames(out, ".to", rname1)
+      data.table::setcolorder(out, c(idcols, rname1, "value"))
+      if (info) {
+        return(list(data = out, info = NULL))
+      }
+      return(out)
+    }
+
+    data.table::setorder(xt, .to, .from, .row0, iso3_min)
+    out <- xt[, list(value = sum(value * w_pair / sw_from)),
+      by = c(".to", idcols)
+    ]
+    data.table::setnames(out, ".to", rname1)
+    data.table::setcolorder(out, c(rname1, idcols, "value"))
+
+    .info_share <- NULL
+    .info_share <- xt[, list(value = sum(w_pair / sw_to)),
+      by = c(".to", idcols)
+    ]
+    data.table::setnames(.info_share, ".to", rname1)
+    data.table::setcolorder(.info_share, c(rname1, idcols, "value"))
+
+    if (info) {
+      return(list(data = out, info = .info_share))
+    }
+    return(out)
+  }
+
+  # mean / set1 / min* / max* : legacy replicates values to countries and
+  # aggregates with the full to-mapping weight sums (sw_to).
+  data.table::setorder(xt, .to, .from, .row0, iso3_min)
+
+  if (agg_operator == "mean") {
+    if (agg_missing == "zero") {
+      out <- xt[, list(value = sum(value * w_pair) / sw_to[1L]),
+        by = c(".to", idcols)
+      ]
+    } else {
+      # agg_missing == "NA": legacy filters NA rows before aggregating, so
+      # groups whose contributions are all NA disappear from the output.
+      out <- xt[!is.na(value),
+        list(value = sum(value * w_pair) / sum(w_pair)),
+        by = c(".to", idcols)
+      ]
+    }
+  } else if (agg_operator %in% c("min", "minw")) {
+    # Legacy: min(value[weight == min(weight)]) over country rows; replicated
+    # values make this expressible with the per-pair minimum weight.
+    xt[, gw := min(w_min), by = c(".to", idcols)]
+    out <- xt[w_min == gw, list(value = min(value)), by = c(".to", idcols)]
+  } else {
+    xt[, gw := max(w_max), by = c(".to", idcols)]
+    out <- xt[w_max == gw, list(value = max(value)), by = c(".to", idcols)]
+  }
+
+  data.table::setnames(out, ".to", rname1)
+  data.table::setcolorder(out, c(rname1, idcols, "value"))
+
+  if (info) {
+    return(list(data = out, info = NULL))
+  }
+  return(out)
+}