Add the "powerLaw" spread link

NEWS.md

@@ -1,5 +1,7 @@
 # fireSense_spreadPredict (development version)
 
+- A ledger row from a fit made with `fireSense_spreadFit`'s new `spreadLink = "powerLaw"` is predicted with that link: `spreadProbOneELF()` hands the row's `gamma` (`fireSenseUtils::spreadFitPowerLawGammaTxt`, `ledgerPowerLawGamma()`) to `fireSenseUtils::spreadProbFromCovs()` with `link = "powerLaw"`, so `p = pc - (pc - lowerSpreadProb) * exp(-x / gamma)` with `pc` the row's `maxAsymptote`, exactly as in the fit; with several ELFs each row uses its own. A `"powerLaw"` row without a valid `gamma` stops. Every other row, such as all rows written so far, predicts exactly as before. The year effect of such a fit is added to the linear predictor, not to the logit, so the module also outputs `fireSense_SpreadPowerLaw` (`pc`, `gamma`, `floor`; blended like `fireSense_SpreadSD` with several ELFs) for `fireSense_burn` (its follow-up PR), which recovers `x` from `p`; it is `NULL` for any other link. With mixed links it is `NA` where only a non-power-law ELF predicts. Needs fireSenseUtils >= 0.2.3.9107.
+
 - A ledger row from a fit made with climate outside the `maxAsymptote` ceiling (`fireSense_spreadFit`'s new `spreadLink = "climateOutside"`) is predicted with that link: `spreadProbOneELF()` hands the row's stored climate reference (`fireSenseUtils::spreadFitClimateRefTxt`, `ledgerClimateRef()`) to `fireSenseUtils::spreadProbFromCovs()`, the one function the fit computes its spread probability with, so the climate term is added on the logit after the fuel ceiling exactly as in the fit; with several ELFs each row uses its own. `fireSense_burn` adds the year effect on the same logit, so nothing changes there. A row without a reference, such as every row written so far, predicts exactly as before (`spreadProbFromCovs()` with no reference is `logisticAll()`). Needs fireSenseUtils >= 0.2.3.9106.
 
 - A single fitted ELF now rescales its covariates with the `covMinMax_spread` of its ledger row, as several ELFs already did, not with `sim$covMinMax_spread`, which `fireSense_spreadFit` rebuilds in every run from that run's tables (in a hindcast, from a formula rebuilt for another year, which can lack a fitted term). A row without ranges (an older ledger) still uses `sim$covMinMax_spread`. A fitted covariate with no range now stops, naming it; before, it silently stayed on its raw scale.

fireSense_spreadPredict.R

@@ -17,7 +17,7 @@ defineModule(sim, list(
   documentation = list("README.txt", "fireSense_spreadPredict.Rmd"),
   reqdPkgs = list("data.table", "magrittr", "Matrix", "methods", "terra", "SpaDES.core (>=3.0.4)", "stats",
                   "ggplot2", "viridis",
-                  "PredictiveEcology/fireSenseUtils@development (>= 0.2.3.9106)"),
+                  "PredictiveEcology/fireSenseUtils@development (>= 0.2.3.9107)"),
   parameters = bindrows(
     defineParameter(name = "lowerSpreadProb", class = "numeric", default = fireSenseUtils::spreadProbFloor,
                     desc = paste("Lower asymptote of the 2- and 3-parameter logistic. Default",
@@ -68,7 +68,13 @@ defineModule(sim, list(
                                "(`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`."))
+                               "`fireSense_SpreadPredicted`.")),
+    createsOutput(objectName = "fireSense_SpreadPowerLaw", objectClass = "SpatRaster|numeric",
+                  desc = paste("Only for a fit made with `spreadLink = \"powerLaw\"` (else NULL): `pc` (the row's",
+                               "`maxAsymptote`), `gamma` and `floor` (`lowerSpreadProb`) of",
+                               "`fireSenseUtils::powerLawLink()`, for `fireSense_burn`, which adds the year effect to the",
+                               "linear predictor, as the fit does. A numeric vector with one fitted ELF; with several, a",
+                               "3-layer raster blended like `fireSense_SpreadSD`."))
   ))
 )
 
@@ -151,11 +157,14 @@ spreadPredictRun <- function(sim) {
 
   if (NROW(sa) == 1L) {
     params1 <- paramSetForRep(sa$params[[1]], P(sim)$.rep)
+    gamma1 <- ledgerPowerLawGamma(sa, 1L)
     pred <- spreadProbOneELF(covs, params = params1, covMinMax = ledgerCovMinMax(sa, 1L, sim$covMinMax_spread),
                              yr = time(sim),
                              maxFireSpread = P(sim)$maxFireSpread, lowerSpreadProb = P(sim)$lowerSpreadProb,
-                             covCentre = ledgerCovCentre(sa, 1L), climateRef = ledgerClimateRef(sa, 1L))
+                             covCentre = ledgerCovCentre(sa, 1L), climateRef = ledgerClimateRef(sa, 1L),
+                             powerLawGamma = gamma1)
     sim$fireSense_SpreadSD <- yearSpreadSDOf(params1)
+    sim$fireSense_SpreadPowerLaw <- powerLawParsOf(params1, gamma1, P(sim)$lowerSpreadProb)
   } else {
     ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]])
     ## each ELF's weight at each pixel: static, so computed once
@@ -168,24 +177,40 @@ spreadPredictRun <- function(sim) {
       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)); accSD <- numeric(NROW(covs))
+    accPL <- matrix(0, NROW(covs), length(powerLawParNames), dimnames = list(NULL, powerLawParNames))
+    gammas <- lapply(seq_along(ids), function(i) ledgerPowerLawGamma(sa, i))
+    wsumPL <- numeric(NROW(covs))   ## weight of the power-law ELFs only
     for (i in seq_along(ids)) {
       these <- which(w[, i] > 0)
       if (!length(these)) next
       parsI <- paramSetForRep(sa$params[[i]], P(sim)$.rep)
+      gammaI <- gammas[[i]]
       p <- spreadProbOneELF(covs[these], params = parsI, covMinMax = sa$covMinMax_spread[[i]],
                             yr = time(sim), maxFireSpread = P(sim)$maxFireSpread,
                             lowerSpreadProb = P(sim)$lowerSpreadProb,
-                            covCentre = ledgerCovCentre(sa, i), climateRef = ledgerClimateRef(sa, i))
+                            covCentre = ledgerCovCentre(sa, i), climateRef = ledgerClimateRef(sa, i),
+                            powerLawGamma = gammaI)
       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(parsI)
+      if (!is.null(gammas[[i]])) {
+        accPL[rows, ] <- accPL[rows, ] + w[rows, i] * rep(powerLawParsOf(parsI, gammas[[i]], P(sim)$lowerSpreadProb), each = length(rows))
+        wsumPL[rows] <- wsumPL[rows] + w[rows, 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_burn scales one z per year by it
     sim$fireSense_SpreadSD <- rast(sim$flammableRTM)
     sim$fireSense_SpreadSD[covs$pixelID[ok]] <- accSD[ok] / wsum[ok]
+    if (any(wsumPL > 0)) {   ## a power-law fit: pc, gamma and floor blended over the power-law ELFs; NA where there is none
+      sim$fireSense_SpreadPowerLaw <- rast(sim$flammableRTM, nlyrs = length(powerLawParNames))
+      names(sim$fireSense_SpreadPowerLaw) <- powerLawParNames
+      okPL <- ok & wsumPL > 0
+      for (nm in powerLawParNames)
+        sim$fireSense_SpreadPowerLaw[[nm]][covs$pixelID[okPL]] <- accPL[okPL, nm] / wsumPL[okPL]
+    }
   }
 
   # Return to raster format
@@ -251,9 +276,13 @@ ELFblendWeights <- function(elfRas, template, pixelID, ids, width) {
 #'   `maxAsymptote` ceiling (`fireSense_spreadFit`'s `spreadLink = "climateOutside"`), the reference value of each
 #'   climate covariate; `fireSenseUtils::spreadProbFromCovs()` then applies that link, exactly as the fit did.
 #'   `NULL` (a fit with climate under the ceiling, or an older ledger row) is the link as it always was.
+#' @param powerLawGamma the ELF's `powerLawGamma_spread` (`ledgerPowerLawGamma()`): with a fit made with
+#'   `fireSense_spreadFit`'s `spreadLink = "powerLaw"`, its `gamma`; `fireSenseUtils::spreadProbFromCovs()` then
+#'   applies that link (`fireSenseUtils::powerLawLink()`, with `maxAsymptote` as `p_c`), exactly as the fit did.
+#'   `NULL` (any other link) is the link as it always was.
 #' @return `data.table` with `pixelID` and `spreadProb`, from the one parameter set in `params`.
 spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSpreadProb,
-                             covCentre = NULL, climateRef = NULL) {
+                             covCentre = NULL, climateRef = NULL, powerLawGamma = NULL) {
   moduleName <- "fireSense_spreadPredict"
   covs <- copy(covs)
   ## the per-year random effect is not a covariate coefficient: fireSense_burn applies it (fireSense_SpreadSD)
@@ -340,11 +369,30 @@ spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSp
   ## the one function the fit computes its spread probability with, so the prediction cannot differ from it
   data.table(pixelID = shortAnnDT$pixelID,
              spreadProb = spreadProbFromCovs(mat[, matching, drop = FALSE], logisticPars, covPars, lowerSpreadProb,
-                                             climateRef = climateRef))
+                                             link = if (!is.null(powerLawGamma)) "powerLaw",
+                                             climateRef = climateRef, powerLawGamma = powerLawGamma))
 }
 
 yearSpreadSDTxt <- "yearSpreadSD"
 
+#' The names of `fireSense_SpreadPowerLaw`
+powerLawParNames <- c("pc", "gamma", "floor")
+
+#' What `fireSense_burn` needs to add the year effect on the linear predictor of a power-law fit
+#'
+#' A fit made with `spreadLink = "powerLaw"` adds its year effect (`yearSpreadSD`) to the linear predictor `x`, not to
+#' the logit of `p`. `fireSense_burn` recovers `x` from `p` with the link's `pc`, `gamma` and `floor`, adds the effect,
+#' and maps back with `fireSenseUtils::powerLawLink()`.
+#'
+#' @param params one-row `data.frame`: the parameter set in use (`paramSetForRep()`); `maxAsymptote` is `pc`.
+#' @param powerLawGamma the row's `gamma` (`ledgerPowerLawGamma()`); `NULL` for another link.
+#' @param lowerSpreadProb the `floor`.
+#' @return A named numeric vector (`powerLawParNames`), or `NULL` for another link.
+powerLawParsOf <- function(params, powerLawGamma, lowerSpreadProb) {
+  if (is.null(powerLawGamma)) return(NULL)
+  c(pc = params[["maxAsymptote"]][1L], gamma = powerLawGamma, floor = lowerSpreadProb)[powerLawParNames]
+}
+
 #' The covariate ranges of one ledger row
 #'
 #' A fit rescales its covariates with fixed ranges and stores them in its ledger row (`covMinMax_spread`).
@@ -381,6 +429,26 @@ ledgerClimateRef <- function(sa, i) {
   if (is.list(ref) && length(ref)) ref else NULL
 }
 
+#' The power-law exponent of one ledger row
+#'
+#' A fit made with `fireSense_spreadFit`'s `spreadLink = "powerLaw"` stores its `gamma` in the ledger row
+#' (`fireSenseUtils::spreadFitPowerLawGammaTxt`), with `"powerLaw"` in `fireSenseUtils::spreadLinkTxt`. Prediction
+#' applies the same link with the same `gamma` (`fireSenseUtils::spreadProbFromCovs()`).
+#'
+#' @param sa `sim$studyAreaWithSpreadParams`, one row per ELF.
+#' @param i integer; the row.
+#' @return `gamma`, or `NULL` for a row of any other link, such as every row written before the link existed. Such
+#'   a row predicts exactly as it always did. A `"powerLaw"` row without a valid `gamma` stops.
+ledgerPowerLawGamma <- function(sa, i) {
+  linkCol <- fireSenseUtils::spreadLinkTxt
+  if (!linkCol %in% names(sa) || !identical(as.character(sa[[linkCol]][[i]]), "powerLaw")) return(NULL)
+  gamma <- sa[[fireSenseUtils::spreadFitPowerLawGammaTxt]][[i]]
+  if (!(is.numeric(gamma) && length(gamma) == 1L && is.finite(gamma) && gamma > 0))
+    stop("fireSense_spreadPredict: ledger row ", i, " has spreadLink 'powerLaw' but no valid '",
+         fireSenseUtils::spreadFitPowerLawGammaTxt, "' (", format(gamma), ")", call. = FALSE)
+  gamma
+}
+
 #' The covariate centres of one ledger row
 #'
 #' A fit with an intercept centres its covariates and stores the centres in the ledger row

tests/testthat/test-metadata.R

@@ -29,7 +29,8 @@ 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", fireSense_SpreadSD = "SpatRaster|numeric")
+    c(fireSense_SpreadPowerLaw = "SpatRaster|numeric", fireSense_SpreadPredicted = "SpatRaster",
+      fireSense_SpreadSD = "SpatRaster|numeric")
   )
 })
 

tests/testthat/test-powerLaw.R

@@ -0,0 +1,104 @@
+## A ledger row from a fit made with fireSense_spreadFit's spreadLink = "powerLaw" carries "powerLaw" in
+## fireSenseUtils::spreadLinkTxt and the fit's gamma in fireSenseUtils::spreadFitPowerLawGammaTxt (the columns
+## fireSense_spreadFit's spreadFitLedgerRow() writes). Prediction applies the same link through
+## fireSenseUtils::spreadProbFromCovs(): p = pc - (pc - 0.13) * exp(-x / gamma), x = covariates %*% coefficients,
+## pc the row's maxAsymptote. A row of any other link predicts as it always did.
+##
+## The toy covariates are those of test-spreadProb-values.R: MDC / 200, youngAge and fuelA / 8; the toy parameters
+## have maxAsymptote = 0.25, hillSlope1 = 2, and coefficients 1.5, -0.5, 1 (the power law ignores hillSlope1:
+## the fit fixes it at 1, and it is the linear predictor that p depends on).
+
+linkColPL <- fireSenseUtils::spreadLinkTxt
+gammaColPL <- fireSenseUtils::spreadFitPowerLawGammaTxt
+cellsPL <- c(1, 2, 3, 4, 6, 7, 8)
+
+inputsPL <- function(gamma = 1.9, link = "powerLaw") {
+  ins <- toyInputs()
+  if (!is.null(link)) ins$studyAreaWithSpreadParams[[linkColPL]] <- link
+  if (!is.null(gamma)) ins$studyAreaWithSpreadParams[[gammaColPL]] <- gamma
+  ins
+}
+
+## by hand, in cell order
+pHandPL <- function(gamma = 1.9, pc = 0.25, floor = 0.13) {
+  cv <- toyCovariates(); o <- order(cv$pixelID)
+  x <- 1.5 * cv$MDC[o] / 200 - 0.5 * cv$youngAge[o] + 1 * cv$fuelA[o] / 8
+  pmax(pc - (pc - floor) * exp(-x / gamma), 0)
+}
+
+test_that("a powerLaw row predicts p = pc - (pc - floor) * exp(-x / gamma), with its own gamma", {
+  v <- predVals(toyRun(inputsPL()))
+  expect_equal(v[cellsPL], pHandPL(1.9), tolerance = 1e-9)
+  expect_true(all(is.na(v[c(5, 9)])))
+  ## cell 1: x = 0, so the floor; cell 2 (youngAge 1, x < 0) is below it; cell 4 (x = 2.5) is near pc and below it
+  expect_equal(v[1], 0.13, tolerance = 1e-9)
+  expect_lt(v[2], 0.13)
+  expect_lt(v[4], 0.25)
+  expect_equal(predVals(toyRun(inputsPL(gamma = 4)))[cellsPL], pHandPL(4), tolerance = 1e-9)
+  ## the logit link of the same row would give something else
+  expect_false(isTRUE(all.equal(v, predVals(toyRun()))))
+})
+
+test_that("a row of any other link predicts exactly as before", {
+  base <- predVals(toyRun())
+  expect_identical(predVals(toyRun(inputsPL(gamma = NULL, link = NULL))), base)
+  expect_identical(predVals(toyRun(inputsPL(gamma = 4, link = "ceiling"))), base)   # gamma ignored
+  sa <- toyInputs()$studyAreaWithSpreadParams
+  expect_null(ledgerPowerLawGamma(sa, 1L))
+  sa[[linkColPL]] <- "climateOutside"
+  expect_null(ledgerPowerLawGamma(sa, 1L))
+  sa[[linkColPL]] <- "powerLaw"; sa[[gammaColPL]] <- 2.2
+  expect_identical(ledgerPowerLawGamma(sa, 1L), 2.2)
+  ## a powerLaw row without a usable gamma stops, instead of predicting with the wrong link
+  sa[[gammaColPL]] <- NA_real_
+  expect_error(ledgerPowerLawGamma(sa, 1L), "no valid")
+  sa[[gammaColPL]] <- NULL
+  expect_error(ledgerPowerLawGamma(sa, 1L), "no valid")
+})
+
+test_that("with several ELFs each row applies its own link, blended like the probabilities", {
+  ins <- multiInputs()
+  ## ELF A is a power-law fit (gamma 3), ELF B the default link; every other coefficient is 0 (as in test-climateOutside.R)
+  ins$studyAreaWithSpreadParams$params <- lapply(ins$studyAreaWithSpreadParams$params, function(p) { p$MDC <- 1; p })
+  ins$studyAreaWithSpreadParams[[linkColPL]] <- c("powerLaw", "ceiling")
+  ins$studyAreaWithSpreadParams[[gammaColPL]] <- c(3, NA)
+  v <- predVals(toyRun(ins, params = list(ELFblendWidth = 20000)))
+  pA <- pmax(0.25 - (0.25 - 0.13) * exp(-0.5 / 3), 0)      # MDC is 100 everywhere: x = 1 * 0.5; pc = ELF A's maxAsymptote
+  pB <- handLogistic3(0.5, upper = 0.27)                    # ELF B keeps the logistic (its hillSlope1 is 2)
+  rA <- c(1, 1, 1, 1, 1, 0.75, 0.5, 0.25, 0, 0); rB <- rev(rA)
+  expect_equal(v, (rA * pA + rB * pB) / (rA + rB), tolerance = 1e-9)
+  expect_equal(v[1], pA, tolerance = 1e-9)
+  expect_equal(v[10], pB, tolerance = 1e-9)
+})
+
+test_that("a powerLaw row hands fireSense_burn pc, gamma and floor, so it can add the year effect to x", {
+  ins <- inputsPL(gamma = 3)
+  ins$studyAreaWithSpreadParams$params <- lapply(ins$studyAreaWithSpreadParams$params, function(p) { p$yearSpreadSD <- 0.3; p })
+  expect_no_warning(sim <- toyRun(ins))
+  expect_equal(predVals(sim)[cellsPL], pHandPL(3), tolerance = 1e-9)
+  expect_equal(sim$fireSense_SpreadSD, 0.3)
+  expect_identical(sim$fireSense_SpreadPowerLaw, c(pc = 0.25, gamma = 3, floor = 0.13))
+  ## any other link: nothing for fireSense_burn, which keeps the year effect on the logit
+  expect_null(toyRun(inputsPL(gamma = NULL, link = NULL))$fireSense_SpreadPowerLaw)
+  expect_null(toyRun(inputsPL(gamma = 4, link = "ceiling"))$fireSense_SpreadPowerLaw)
+})
+
+test_that("with several ELFs fireSense_SpreadPowerLaw is blended like the probabilities", {
+  ins <- multiInputs()
+  ins$studyAreaWithSpreadParams$params <- lapply(ins$studyAreaWithSpreadParams$params, function(p) { p$MDC <- 1; p })
+  ins$studyAreaWithSpreadParams[[linkColPL]] <- c("powerLaw", "powerLaw")
+  ins$studyAreaWithSpreadParams[[gammaColPL]] <- c(2, 4)
+  pl <- toyRun(ins, params = list(ELFblendWidth = 20000))$fireSense_SpreadPowerLaw
+  expect_s4_class(pl, "SpatRaster")
+  expect_identical(names(pl), c("pc", "gamma", "floor"))
+  g <- terra::values(pl[["gamma"]], mat = FALSE)
+  expect_equal(g[1], 2); expect_equal(g[10], 4)                 # each ELF's own gamma, in its own interior
+  expect_true(all(g >= 2 & g <= 4, na.rm = TRUE))
+  expect_equal(terra::values(pl[["floor"]], mat = FALSE), rep(0.13, 10))
+  ins$studyAreaWithSpreadParams[[linkColPL]] <- c("powerLaw", "ceiling")
+  ins$studyAreaWithSpreadParams[[gammaColPL]] <- c(2, NA)
+  ## one power-law ELF and one other: the parameters are NA where only the other ELF predicts (burn then uses the logit)
+  pl2 <- toyRun(ins, params = list(ELFblendWidth = 20000))$fireSense_SpreadPowerLaw
+  g2 <- terra::values(pl2[["gamma"]], mat = FALSE)
+  expect_equal(g2[1], 2); expect_true(is.na(g2[10]))
+})