Each prediction replicate uses one whole parameter set (new parameter .rep)

NEWS.md

@@ -1,5 +1,7 @@
 # fireSense_spreadPredict (development version)
 
+- Predictions no longer average parameter sets (`spreadProbOneELF()` took `rowMeans()` of one map per row, and `yearSpreadSDOf()` the mean `yearSpreadSD`, which broke the pairing of each set's `yearSpreadSD` with its own coefficients). New parameter `.rep` (integer, default 1): replicate `.rep` uses the one whole parameter set `((rep - 1) %% number of sets) + 1` of each ELF, and `fireSense_SpreadSD` is that set's `yearSpreadSD`. `.rep` is a SpaDES-aware parameter: `SpaDES.project::setupProject()` sets `.globals$.rep` from the experiment's `.rep` (PredictiveEcology/SpaDES.project#190), so nothing needs wiring by hand.
+
 - reqdPkgs now lists `data.table`, which the module calls (`copy`, `data.table`, `setDT`) but did not list. Version 1.1.2.
 
 - `spreadProbOneELF()` (fireSense_spreadPredict.R:266-323 pre-fix) went on when a fitted coefficient had no covariate column, and died in `rowMeans(spreadProbMat)` with "'x' must be an array of at least two dimensions" (a predict-only run whose fit had `nfLCC_100_60` and `nfLCC_40_50_80` but whose covariates had a single `nf`). It now stops at once, naming the coefficients without a covariate and the covariates available, and says the non-forest groups / fuel classes differ from the fit's. Version 1.1.1.

fireSense_spreadPredict.R

@@ -29,6 +29,10 @@ defineModule(sim, list(
     defineParameter("maxFireSpread", "numeric", default = 0.28,
                     desc = paste("Upper limit on `spreadProb` used when fitting. Here it is only checked",
                                  "to be the same in every module that defines it.")),
+    defineParameter(".rep", "integer", default = 1L, min = 1L, max = NA,
+                    desc = paste("Which replicate this prediction is. It selects parameter set",
+                                 "((rep - 1) %% number of sets) + 1 of each ELF, so replicates cycle through the",
+                                 "fitted sets, each used whole (with its own `yearSpreadSD`), never averaged.")),
     defineParameter(name = ".runInitialTime", class = "numeric", default = start(sim),
                     desc = "Time of the first prediction."),
     defineParameter(name = ".runInterval", class = "numeric", default = 1,
@@ -59,7 +63,8 @@ defineModule(sim, list(
                   desc = "Spread probability of each flammable pixel, this year."),
     createsOutput(objectName = "fireSense_SpreadSD", objectClass = "SpatRaster|numeric",
                   desc = paste("The fitted sd of the per-year random effect on logit spread probability",
-                               "(`yearSpreadSD`; 0 if the fit has none), for `fireSense_burn`. One number with one",
+                               "(`yearSpreadSD` of the parameter set chosen by `.rep`; 0 if the fit has none), for",
+                               "`fireSense_burn`. One number with one",
                                "fitted ELF; with several, a raster blended across ELFs with the weights of",
                                "`fireSense_SpreadPredicted`."))
   ))
@@ -115,11 +120,11 @@ doEvent.fireSense_spreadPredict <- function(sim, eventTime, eventType, debug = F
 #' Predict this year's spread probability
 #'
 #' Rescales `sim$fireSense_SpreadCovariates` with the fit's covariate ranges, computes the spread
-#' probability for each parameter set (row) of the fit, and writes the mean over parameter sets to
-#' `sim$fireSense_SpreadPredicted`.
+#' probability with ONE parameter set (row) of the fit, chosen by the `.rep` parameter (see
+#' `paramSetForRep()`), and writes it to `sim$fireSense_SpreadPredicted`. Parameter sets are not averaged.
 #'
 #' With one fitted ELF (one row of `sim$studyAreaWithSpreadParams`) every pixel uses its parameters and
-#' `sim$covMinMax_spread`. With several, each ELF's model (its parameter sets, its `covMinMax_spread`, and
+#' `sim$covMinMax_spread`. With several, each ELF's model (its parameter set, its `covMinMax_spread`, and
 #' only the covariates it was fitted with) predicts its own pixels, from `sim$rasterToMatchLargeELF`, and
 #' those within `ELFblendWidth` of them; overlapping predictions are averaged with the weights of
 #' `ELFblendWeights()`.
@@ -143,10 +148,11 @@ spreadPredictRun <- function(sim) {
          "module's -- or the shared ledger has no row for this polygon.", call. = FALSE)
 
   if (NROW(sa) == 1L) {
-    pred <- spreadProbOneELF(covs, params = sa$params[[1]], covMinMax = sim$covMinMax_spread,
+    params1 <- paramSetForRep(sa$params[[1]], P(sim)$.rep)
+    pred <- spreadProbOneELF(covs, params = params1, covMinMax = sim$covMinMax_spread,
                              formula = sim$fireSense_spreadFormula, yr = time(sim),
                              maxFireSpread = P(sim)$maxFireSpread, lowerSpreadProb = P(sim)$lowerSpreadProb)
-    sim$fireSense_SpreadSD <- yearSpreadSDOf(sa$params[[1]])
+    sim$fireSense_SpreadSD <- yearSpreadSDOf(params1)
   } else {
     ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]])
     ## each ELF's weight at each pixel: static, so computed once
@@ -162,12 +168,13 @@ spreadPredictRun <- function(sim) {
     for (i in seq_along(ids)) {
       these <- which(w[, i] > 0)
       if (!length(these)) next
-      p <- spreadProbOneELF(covs[these], params = sa$params[[i]], covMinMax = sa$covMinMax_spread[[i]],
+      parsI <- paramSetForRep(sa$params[[i]], P(sim)$.rep)
+      p <- spreadProbOneELF(covs[these], params = parsI, covMinMax = sa$covMinMax_spread[[i]],
                             formula = NULL, yr = time(sim), maxFireSpread = P(sim)$maxFireSpread,
                             lowerSpreadProb = P(sim)$lowerSpreadProb, byParams = TRUE)
       rows <- these[match(p$pixelID, covs$pixelID[these])]
       acc[rows] <- acc[rows] + w[rows, i] * p$spreadProb
-      accSD[rows] <- accSD[rows] + w[rows, i] * yearSpreadSDOf(sa$params[[i]])
+      accSD[rows] <- accSD[rows] + w[rows, i] * yearSpreadSDOf(parsI)
       wsum[rows] <- wsum[rows] + w[rows, i]
     }
     ok <- wsum > 0
@@ -226,14 +233,14 @@ ELFblendWeights <- function(elfRas, template, pixelID, ids, width) {
 #'
 #' @param covs `data.table` of this ELF's pixels: `pixelID` and covariates (as `fireSense_dataPrepPredict`
 #'   makes them; fuel biomass logged).
-#' @param params `data.frame` of the ELF's fitted parameter sets, one per row (ledger `params`).
+#' @param params one-row `data.frame`: the ELF's parameter set for this replicate (`paramSetForRep()`).
 #' @param covMinMax the ELF's `covMinMax_spread`.
 #' @param formula the spread formula, to check the covariates are all there.
 #' @param byParams logical; `TRUE` (several ELFs) ignores `formula`: the covariates are those named in
 #'   `params`, and the table is cut to them, since it holds every ELF's covariates.
 #' @param yr,maxFireSpread,lowerSpreadProb as for `fireSenseUtils::spreadProbFromIntegerCovs()` and
 #'   `fireSenseUtils::logisticAll()`.
-#' @return `data.table` with `pixelID` and `spreadProb`, the mean over parameter sets.
+#' @return `data.table` with `pixelID` and `spreadProb`, from the one parameter set in `params`.
 spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread, lowerSpreadProb,
                              byParams = FALSE) {
   moduleName <- "fireSense_spreadPredict"
@@ -315,36 +322,43 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
 
   mat <- as.matrix(shortAnnDT[, ..colsToUse])
 
-  # for replicate "best" params from DEoptim
-  spreadProbList <- lapply(seq_len(NROW(params)), function(ind) {
-    par <- params[ind, ] |> as.vector() |> unlist()
-    covPars <- intersect(names(par), colsToUse)
-    covPars <- par[covPars]
-    logisticPars <- par[setdiff(names(par), names(covPars))]
-    # Make sure the order is correct in the matrix
-    matching <- intersect(names(covPars), colnames(mat))
-    missingCovs <- setdiff(colnames(mat), names(covPars))
-    if (length(missingCovs))
-      warning("There are covariates in the sim$fireSense_SpreadCovariates: \n",
-              paste0(missingCovs, collapse = ", "),
-              "\n...that are not in the sim$studyAreaWithSpreadParams")
-    logisticAll(logisticPars, mat[, matching, drop = FALSE], covPars, lowerSpreadProb)
-  })
-  spreadProbMat <- do.call(cbind, spreadProbList)
+  par <- params[1L, ] |> as.vector() |> unlist()
+  covPars <- par[intersect(names(par), colsToUse)]
+  logisticPars <- par[setdiff(names(par), names(covPars))]
+  # Make sure the order is correct in the matrix
+  matching <- intersect(names(covPars), colnames(mat))
+  missingCovs <- setdiff(colnames(mat), names(covPars))
+  if (length(missingCovs))
+    warning("There are covariates in the sim$fireSense_SpreadCovariates: \n",
+            paste0(missingCovs, collapse = ", "),
+            "\n...that are not in the sim$studyAreaWithSpreadParams")
 
-  data.table(pixelID = shortAnnDT$pixelID, spreadProb = rowMeans(spreadProbMat))
+  data.table(pixelID = shortAnnDT$pixelID,
+             spreadProb = logisticAll(logisticPars, mat[, matching, drop = FALSE], covPars, lowerSpreadProb))
 }
 
 yearSpreadSDTxt <- "yearSpreadSD"
 
+#' The parameter set a replicate uses
+#'
+#' Replicate `rep` uses row `((rep - 1) %% NROW(params)) + 1` of an ELF's fitted parameter sets, whole: its
+#' coefficients and its `yearSpreadSD` stay together. Each ELF takes the modulo of its own number of rows.
+#'
+#' @param params `data.frame` of an ELF's fitted parameter sets, one per row (ledger `params`).
+#' @param rep integer; the replicate (the module parameter `.rep`).
+#' @return One-row `data.frame`.
+paramSetForRep <- function(params, rep) {
+  params[((rep - 1L) %% NROW(params)) + 1L, , drop = FALSE]
+}
+
 #' The fitted sd of the per-year random effect
 #'
 #' `yearSpreadSD` (fireSense_spreadFit, fireSenseUtils >= 0.2.3.9041) is one eps per year on logit spread
-#' probability. With several retained parameter sets, their mean, as the spread probabilities are averaged.
+#' probability.
 #'
-#' @param params `data.frame` of fitted parameters, one row per retained set.
+#' @param params one-row `data.frame`: the parameter set in use (`paramSetForRep()`).
 #' @return Numeric; 0 when the fit has no `yearSpreadSD`.
 yearSpreadSDOf <- function(params) {
   if (!yearSpreadSDTxt %in% names(params)) return(0)
-  mean(params[[yearSpreadSDTxt]])
+  params[[yearSpreadSDTxt]]
 }

fireSense_spreadPredict.Rmd

@@ -48,8 +48,8 @@ download.file(url = "https://img.shields.io/badge/Made%20with-Markdown-1f425f.pn
 Each year, predicts a raster of fire spread probabilities from the parameters fitted by *fireSense_spreadFit*, for the spread component of fireSense [@Marchal:2017a; @Marchal:2017b; @Marchal:2019].
 
 1. The covariates in `fireSense_SpreadCovariates` are rescaled to [0, 1] using `covMinMax_spread`, the range of the fitting data.
-2. For each parameter set (row) in `studyAreaWithSpreadParams$params[[1]]`, the spread probability is a 2- or 3-parameter logistic of the linear combination of the covariates, with lower asymptote `lowerSpreadProb`.
-3. `fireSense_SpreadPredicted` is the mean over parameter sets, on the `flammableRTM` grid.
+2. With the one parameter set (row of `studyAreaWithSpreadParams$params[[i]]`) chosen by `.rep`, i.e. row `((rep - 1) %% number of sets) + 1` of each ELF, the spread probability is a 2- or 3-parameter logistic of the linear combination of the covariates, with lower asymptote `lowerSpreadProb`.
+3. `fireSense_SpreadPredicted` is that single set's prediction (sets are not averaged), and `fireSense_SpreadSD` its `yearSpreadSD`, on the `flammableRTM` grid.
 
 ### Module inputs and parameters
 

fireSense_spreadPredict.md

@@ -40,8 +40,8 @@ Eliot McIntire <eliot.mcintire@nrcan-rncan.gc.ca> [aut, cre], Tati Micheletti <t
 Each year, predicts a raster of fire spread probabilities from the parameters fitted by *fireSense_spreadFit*, for the spread component of fireSense [@Marchal:2017a; @Marchal:2017b; @Marchal:2019].
 
 1. The covariates in `fireSense_SpreadCovariates` are rescaled to [0, 1] using `covMinMax_spread`, the range of the fitting data.
-2. For each parameter set (row) in `studyAreaWithSpreadParams$params[[1]]`, the spread probability is a 2- or 3-parameter logistic of the linear combination of the covariates, with lower asymptote `lowerSpreadProb`.
-3. `fireSense_SpreadPredicted` is the mean over parameter sets, on the `flammableRTM` grid.
+2. With the one parameter set (row of `studyAreaWithSpreadParams$params[[i]]`) chosen by `.rep`, i.e. row `((rep - 1) %% number of sets) + 1` of each ELF, the spread probability is a 2- or 3-parameter logistic of the linear combination of the covariates, with lower asymptote `lowerSpreadProb`.
+3. `fireSense_SpreadPredicted` is that single set's prediction (sets are not averaged), and `fireSense_SpreadSD` its `yearSpreadSD`, on the `flammableRTM` grid.
 
 ### Module inputs and parameters
 

tests/testthat/setup-toySpread.R

@@ -76,3 +76,24 @@ evOf <- function(dt, type) {
   df <- as.data.frame(dt)
   df[df$moduleName == "fireSense_spreadPredict" & df$eventType == type, , drop = FALSE]
 }
+
+## Two ELFs, "A" and "B", for the multi-ELF tests (see test-multiELF.R for the layout).
+multiInputs <- function() {
+  crs <- "EPSG:3005"
+  flam <- terra::rast(nrows = 1, ncols = 10, xmin = 0, xmax = 50000, ymin = 0, ymax = 5000, crs = crs, vals = 1)
+  elf <- terra::rast(flam); terra::values(elf) <- rep(1:2, each = 5)
+  levels(elf) <- data.frame(id = 1:2, ELFind = c("A", "B"))
+  covs <- data.table::data.table(pixelID = 1:10, MDC = 100, fuelA = 2, fuelB = 3)
+  pars <- function(max, fuel) {
+    d <- data.frame(maxAsymptote = max, hillSlope1 = 2, inflectionPoint1 = 1, MDC = 0, x = 0)
+    names(d)[5] <- fuel
+    d
+  }
+  cmm <- function(fuel) { d <- data.table::data.table(MDC = c(0, 200), f = c(0, 8)); data.table::setnames(d, "f", fuel); d }
+  sa <- data.frame(polygonID = c("A", "B"))
+  sa$params <- list(pars(0.25, "fuelA"), pars(0.27, "fuelB"))
+  sa$covMinMax_spread <- list(cmm("fuelA"), cmm("fuelB"))
+  list(flammableRTM = flam, rasterToMatchLargeELF = elf, fireSense_SpreadCovariates = covs,
+       studyAreaWithSpreadParams = sa, fireSense_spreadFormula = "~ MDC + fuelA - 1",
+       covMinMax_spread = cmm("fuelA"))
+}

tests/testthat/test-metadata.R

@@ -47,10 +47,10 @@ paramTable <- function(md) {
 test_that("parameters have the expected names, classes and defaults", {
   md <- SpaDES.core::moduleMetadata(module = moduleName, path = modulePath)
   expected <- data.frame(
-    class   = c("numeric", "numeric", "numeric", "character", "logical", "numeric", "numeric", "numeric"),
+    class   = c("integer", "numeric", "numeric", "numeric", "character", "logical", "numeric", "numeric", "numeric"),
     ## .runInitialTime defaults to start(sim), which is 0 when only metadata is parsed
-    default = c("0", "1", "NA", "NA", "FALSE", "20000", "0.13", "0.28"),
-    row.names = c(".runInitialTime", ".runInterval", ".saveInitialTime", ".studyAreaName", ".useCache",
+    default = c("1L", "0", "1", "NA", "NA", "FALSE", "20000", "0.13", "0.28"),
+    row.names = c(".rep", ".runInitialTime", ".runInterval", ".saveInitialTime", ".studyAreaName", ".useCache",
                   "ELFblendWidth", "lowerSpreadProb", "maxFireSpread")
   )
   expect_identical(paramTable(md), expected)

tests/testthat/test-multiELF.R

@@ -7,26 +7,6 @@
 ## so each ELF predicts a constant: lower + (maxAsymptote - lower) / 2 = 0.19 for A (max 0.25) and
 ## 0.20 for B (max 0.27). A is fitted on fuelA, B on fuelB; the covariate table holds both.
 
-multiInputs <- function() {
-  crs <- "EPSG:3005"
-  flam <- terra::rast(nrows = 1, ncols = 10, xmin = 0, xmax = 50000, ymin = 0, ymax = 5000, crs = crs, vals = 1)
-  elf <- terra::rast(flam); terra::values(elf) <- rep(1:2, each = 5)
-  levels(elf) <- data.frame(id = 1:2, ELFind = c("A", "B"))
-  covs <- data.table::data.table(pixelID = 1:10, MDC = 100, fuelA = 2, fuelB = 3)
-  pars <- function(max, fuel) {
-    d <- data.frame(maxAsymptote = max, hillSlope1 = 2, inflectionPoint1 = 1, MDC = 0, x = 0)
-    names(d)[5] <- fuel
-    d
-  }
-  cmm <- function(fuel) { d <- data.table::data.table(MDC = c(0, 200), f = c(0, 8)); data.table::setnames(d, "f", fuel); d }
-  sa <- data.frame(polygonID = c("A", "B"))
-  sa$params <- list(pars(0.25, "fuelA"), pars(0.27, "fuelB"))
-  sa$covMinMax_spread <- list(cmm("fuelA"), cmm("fuelB"))
-  list(flammableRTM = flam, rasterToMatchLargeELF = elf, fireSense_SpreadCovariates = covs,
-       studyAreaWithSpreadParams = sa, fireSense_spreadFormula = "~ MDC + fuelA - 1",
-       covMinMax_spread = cmm("fuelA"))
-}
-
 test_that("each ELF predicts its own pixels, and the boundary is a distance-weighted blend", {
   v <- predVals(toyRun(multiInputs(), params = list(ELFblendWidth = 20000)))
   ## raw weights at pixel centres (5 km apart): A = 1 on 1-5, then 0.75, 0.5, 0.25, 0, 0 on 6-10

tests/testthat/test-repParamSet.R

@@ -0,0 +1,51 @@
+## Each prediction replicate uses ONE whole fitted parameter set (ledger row), never an average:
+## row ((rep - 1) %% number of sets) + 1, and fireSense_SpreadSD is that same row's yearSpreadSD.
+
+threeSets <- function() {
+  p <- rbind(toyParams(), toyParams(), toyParams())
+  p$maxAsymptote <- c(0.25, 0.35, 0.20)
+  p$MDC <- c(1.5, 0, 3)
+  p$yearSpreadSD <- c(0.1, 0.2, 0.4)
+  p
+}
+
+test_that("rep selects parameter set ((rep - 1) %% nSets) + 1, and the prediction is that set's, not the mean", {
+  p <- threeSets()
+  for (r in 1:4) {
+    row <- ((r - 1) %% 3) + 1
+    v <- predVals(toyRun(toyInputs(params = p), params = list(.rep = r)))
+    one <- predVals(toyRun(toyInputs(params = p[row, , drop = FALSE])))
+    expect_equal(v, one, tolerance = 1e-12, info = paste("rep", r))
+  }
+  ## the three sets give different maps, so the check above can tell them apart
+  vs <- lapply(1:3, function(r) predVals(toyRun(toyInputs(params = p), params = list(.rep = r))))
+  expect_gt(max(abs(vs[[1]] - vs[[2]]), na.rm = TRUE), 1e-3)
+  expect_gt(max(abs(vs[[1]] - vs[[3]]), na.rm = TRUE), 1e-3)
+})
+
+test_that("fireSense_SpreadSD is the yearSpreadSD of the same row as the coefficients", {
+  p <- threeSets()
+  sds <- vapply(1:4, function(r) toyRun(toyInputs(params = p), params = list(.rep = r))$fireSense_SpreadSD, numeric(1))
+  expect_equal(sds, c(0.1, 0.2, 0.4, 0.1))
+})
+
+test_that("with several ELFs each ELF takes the modulo of its own number of sets", {
+  ins <- multiInputs()
+  sa <- ins$studyAreaWithSpreadParams
+  ## ELF A has 3 sets (max 0.25, 0.30, 0.35), ELF B has 2 (max 0.27, 0.21); sd likewise
+  more <- function(p, max, sd) { p <- p[rep(1L, length(max)), , drop = FALSE]; p$maxAsymptote <- max; p$yearSpreadSD <- sd; p }
+  sa$params <- list(more(sa$params[[1]], c(0.25, 0.30, 0.35), c(0.1, 0.2, 0.3)),
+                    more(sa$params[[2]], c(0.27, 0.21), c(0.5, 0.7)))
+  ins$studyAreaWithSpreadParams <- sa
+  ## all coefficients are 0: each ELF's value is lower + (max - lower) / 2
+  lev <- function(max) 0.13 + (max - 0.13) / 2
+  for (r in 1:4) {
+    sim <- toyRun(ins, params = list(.rep = r, ELFblendWidth = 6000))
+    a <- ((r - 1) %% 3) + 1; b <- ((r - 1) %% 2) + 1
+    v <- predVals(sim); sdv <- terra::values(sim$fireSense_SpreadSD, mat = FALSE)
+    expect_equal(v[1], lev(sa$params[[1]]$maxAsymptote[a]), tolerance = 1e-9, info = paste("rep", r))
+    expect_equal(v[10], lev(sa$params[[2]]$maxAsymptote[b]), tolerance = 1e-9, info = paste("rep", r))
+    expect_equal(sdv[1], sa$params[[1]]$yearSpreadSD[a], tolerance = 1e-9, info = paste("rep", r))
+    expect_equal(sdv[10], sa$params[[2]]$yearSpreadSD[b], tolerance = 1e-9, info = paste("rep", r))
+  }
+})