yearSpreadSD: not a coefficient; passed to fireSense as fireSense_SpreadSD

NEWS.md

@@ -1,5 +1,11 @@
 # fireSense_SpreadPredict (development version)
 
+- The fitted per-year random effect (`yearSpreadSD`, fireSense_SpreadFit / fireSenseUtils >= 0.2.3.9041) is no
+  longer read as a covariate coefficient: with it, a single ELF treated it as a fourth logistic parameter and several
+  ELFs stopped with "'yearSpreadSD' not found". It becomes the new output `fireSense_SpreadSD`, which fireSense
+  scales one draw per year by: one number with one ELF (the mean over retained parameter sets), a raster blended
+  with the spread probabilities' weights with several. 0 when the fit has none.
+
 - Several fitted ELFs in one study area. Every pixel gets a spread probability: each ELF's model (its parameter sets, its `covMinMax_spread` and only the covariates it was fitted with, from its ledger row) predicts its own pixels and those within `ELFblendWidth` (default 20 km) of them, and overlapping predictions are averaged with weights falling linearly from 1 inside an ELF to 0 at `ELFblendWidth` outside it. Each pixel's ELF comes from the new input `rasterToMatchLargeELF` (fireSense_ELFs with a `studyAreaLarge`). One ELF works as before.
 - A fit made with fireSense_SpreadFit's `link = "logistic3pUpper"` stores `upperTail1`; prediction uses the upper-tail link for it, chosen by the parameter's name (`fireSenseUtils::logisticAll()`). Needs fireSenseUtils >= 0.2.3.9038.
 

fireSense_SpreadPredict.R

@@ -10,7 +10,7 @@ defineModule(sim, list(
     person("Alex M.", "Chubaty", email = "achubaty@for-cast.ca", role = "ctb")
   ),
   childModules = character(),
-  version = list(fireSense_SpreadPredict = "1.0.0.9003", SpaDES.core = "0.1.0"),
+  version = list(fireSense_SpreadPredict = "1.0.0.9004", SpaDES.core = "0.1.0"),
   timeframe = as.POSIXlt(c(NA, NA)),
   timeunit = "year",
   citation = list("citation.bib"),
@@ -54,7 +54,12 @@ defineModule(sim, list(
   ),
   outputObjects = bindrows(
     createsOutput(objectName = "fireSense_SpreadPredicted", objectClass = "SpatRaster",
-                  desc = "Spread probability of each flammable pixel, this year.")
+                  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`. One number with one",
+                               "fitted ELF; with several, a raster blended across ELFs with the weights of",
+                               "`fireSense_SpreadPredicted`."))
   ))
 )
 
@@ -139,6 +144,7 @@ spreadPredictRun <- function(sim) {
     pred <- spreadProbOneELF(covs, params = sa$params[[1]], 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]])
   } else {
     ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]])
     ## each ELF's weight at each pixel: static, so computed once
@@ -150,7 +156,7 @@ spreadPredictRun <- function(sim) {
     if (any(none))
       warning("fireSense_SpreadPredict: ", sum(none), " flammable pixels are more than ",
               P(sim)$ELFblendWidth, " m from every fitted ELF; they get no spread probability", call. = FALSE)
-    acc <- numeric(NROW(covs)); wsum <- numeric(NROW(covs))
+    acc <- numeric(NROW(covs)); wsum <- numeric(NROW(covs)); accSD <- numeric(NROW(covs))
     for (i in seq_along(ids)) {
       these <- which(w[, i] > 0)
       if (!length(these)) next
@@ -159,10 +165,14 @@ spreadPredictRun <- function(sim) {
                             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]])
       wsum[rows] <- wsum[rows] + w[rows, i]
     }
     ok <- wsum > 0
     pred <- data.table(pixelID = covs$pixelID[ok], spreadProb = acc[ok] / wsum[ok])
+    ## each ELF's year effect sd, blended like the probabilities: fireSense scales one z per year by it
+    sim$fireSense_SpreadSD <- rast(sim$flammableRTM)
+    sim$fireSense_SpreadSD[covs$pixelID[ok]] <- accSD[ok] / wsum[ok]
   }
 
   # Return to raster format
@@ -226,6 +236,11 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
                              byParams = FALSE) {
   moduleName <- "fireSense_SpreadPredict"
   covs <- copy(covs)
+  ## the per-year random effect is not a covariate coefficient: fireSense applies it (fireSense_SpreadSD)
+  if (yearSpreadSDTxt %in% names(params)) {
+    params <- as.data.frame(params)
+    params[[yearSpreadSDTxt]] <- NULL
+  }
 
   ## Fuel biomass arrives logged (fireSenseUtils::logMinB()). A fit made on LINEAR fuel biomass has
   ## fireSenseUtils::fuelLinearRange, c(0, 1e4), as that covariate's covMinMax_spread, and its
@@ -290,3 +305,17 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
 
   data.table(pixelID = shortAnnDT$pixelID, spreadProb = rowMeans(spreadProbMat))
 }
+
+yearSpreadSDTxt <- "yearSpreadSD"
+
+#' 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.
+#'
+#' @param params `data.frame` of fitted parameters, one row per retained set.
+#' @return Numeric; 0 when the fit has no `yearSpreadSD`.
+yearSpreadSDOf <- function(params) {
+  if (!yearSpreadSDTxt %in% names(params)) return(0)
+  mean(params[[yearSpreadSDTxt]])
+}

tests/testthat/test-metadata.R

@@ -29,7 +29,7 @@ test_that("outputs are the expected names and classes", {
   outputs <- stats::setNames(md$outputObjects$objectClass, md$outputObjects$objectName)
   expect_identical(
     outputs[order(names(outputs))],
-    c(fireSense_SpreadPredicted = "SpatRaster")
+    c(fireSense_SpreadPredicted = "SpatRaster", fireSense_SpreadSD = "SpatRaster|numeric")
   )
 })
 

tests/testthat/test-multiELF.R

@@ -49,3 +49,25 @@ test_that("with several ELFs, a missing ELF raster stops with a message that say
   ins <- multiInputs(); ins$rasterToMatchLargeELF <- NULL
   expect_error(toyRun(ins), "rasterToMatchLargeELF")
 })
+
+## ---- the per-year random effect (yearSpreadSD) ----
+withSD <- function(inputs, sds) {
+  inputs$studyAreaWithSpreadParams$params <- Map(function(p, s) { p$yearSpreadSD <- s; p },
+                                                 inputs$studyAreaWithSpreadParams$params, sds)
+  inputs
+}
+
+test_that("yearSpreadSD does not change the spread probability, and is blended into fireSense_SpreadSD", {
+  base <- toyRun(multiInputs(), params = list(ELFblendWidth = 20000))
+  sim <- toyRun(withSD(multiInputs(), c(0.2, 0.6)), params = list(ELFblendWidth = 20000))
+  expect_equal(predVals(sim), predVals(base), tolerance = 1e-12)
+  sdv <- terra::values(sim$fireSense_SpreadSD, mat = FALSE)
+  rA <- c(1, 1, 1, 1, 1, 0.75, 0.5, 0.25, 0, 0); rB <- rev(rA)
+  expect_equal(sdv, (rA * 0.2 + rB * 0.6) / (rA + rB), tolerance = 1e-9)   # the probabilities' weights
+  expect_equal(sdv[c(1, 10)], c(0.2, 0.6))
+})
+
+test_that("without yearSpreadSD in the fits the sd is 0 everywhere", {
+  sim <- toyRun(multiInputs(), params = list(ELFblendWidth = 20000))
+  expect_true(all(terra::values(sim$fireSense_SpreadSD, mat = FALSE) == 0))
+})

tests/testthat/test-spreadProb-values.R

@@ -215,3 +215,13 @@ test_that("a fit with upperTail1 is predicted with the upper-tail link, from the
   v0 <- predVals(toyRun(toyInputs(params = toyParams(upperTail1 = 0))))
   expect_equal(v0, predVals(toyRun()), tolerance = 1e-12)
 })
+
+test_that("one ELF: yearSpreadSD is not a coefficient; fireSense_SpreadSD is its mean over parameter sets", {
+  base <- toyRun()
+  p2 <- rbind(toyParams(), toyParams())
+  p2$yearSpreadSD <- c(0.3, 0.5)
+  sim <- toyRun(toyInputs(params = p2))
+  expect_equal(predVals(sim), predVals(base), tolerance = 1e-12)
+  expect_equal(sim$fireSense_SpreadSD, 0.4)
+  expect_identical(toyRun()$fireSense_SpreadSD, 0)
+})