From f0ef4e74cfbd0549862899320de4dc14214f99fe Mon Sep 17 00:00:00 2001 From: Eliot McIntire Date: Sun, 18 Oct 2026 17:32:55 -0700 Subject: [PATCH 1/2] Add the "powerLaw" spread link (needs fireSenseUtils PR 167) fireSense_spreadFit: spreadLink = "powerLaw", powerLawGamma, p_c bounds, ledger columns. fireSense_spreadPredict: reads the row's link and gamma and applies powerLawLink via spreadProbFromCovs. Floor fireSenseUtils (>= 0.2.3.9107). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014G18Mik3fH4z5M5CyMSUS6 --- NEWS.md | 2 + fireSense_spreadPredict.R | 60 ++++++++++++++++++++++-- tests/testthat/test-powerLaw.R | 85 ++++++++++++++++++++++++++++++++++ 3 files changed, 142 insertions(+), 5 deletions(-) create mode 100644 tests/testthat/test-powerLaw.R diff --git a/NEWS.md b/NEWS.md index 3841e5e..ee6bce2 100644 --- a/NEWS.md +++ b/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. `fireSense_burn` still adds the year effect (`fireSense_SpreadSD`) on the logit, whereas a power-law fit adds it to the linear predictor, so a power-law row with `yearSpreadSD > 0` gives a warning (once) until `fireSense_burn` follows. 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. diff --git a/fireSense_spreadPredict.R b/fireSense_spreadPredict.R index 3c30bfa..598a19a 100644 --- a/fireSense_spreadPredict.R +++ b/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", @@ -151,10 +151,13 @@ spreadPredictRun <- function(sim) { if (NROW(sa) == 1L) { params1 <- paramSetForRep(sa$params[[1]], P(sim)$.rep) + gamma1 <- ledgerPowerLawGamma(sa, 1L) + warnPowerLawYearEffect(gamma1, params1) 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) } else { ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]]) @@ -172,10 +175,13 @@ spreadPredictRun <- function(sim) { these <- which(w[, i] > 0) if (!length(these)) next parsI <- paramSetForRep(sa$params[[i]], P(sim)$.rep) + gammaI <- ledgerPowerLawGamma(sa, i) + warnPowerLawYearEffect(gammaI, parsI) 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) @@ -251,9 +257,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 +350,31 @@ 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" +#' Warn, once, that `fireSense_burn` adds the year effect on the wrong scale for a power-law fit +#' +#' A fit made with `spreadLink = "powerLaw"` adds its year effect (`yearSpreadSD`) to the linear predictor, not to +#' the logit of `p`; `fireSense_burn` adds `fireSense_SpreadSD` on the logit, so with a power-law fit and a +#' `yearSpreadSD` above 0 the simulated years are not those of the fit. +#' +#' @param powerLawGamma the row's `gamma` (`ledgerPowerLawGamma()`); `NULL` for another link. +#' @param params one-row `data.frame`: the parameter set in use (`paramSetForRep()`). +#' @return `NULL`, invisibly. +warnPowerLawYearEffect <- function(powerLawGamma, params) { + if (!is.null(powerLawGamma) && yearSpreadSDOf(params) > 0 && !isTRUE(mod$warnedPowerLawYearEffect)) { + mod$warnedPowerLawYearEffect <- TRUE + warning("fireSense_spreadPredict: this ledger row is a spreadLink = 'powerLaw' fit with yearSpreadSD > 0; ", + "the fit adds the year effect to the linear predictor, but fireSense_burn adds fireSense_SpreadSD on the ", + "logit of the spread probability, so simulated years will differ from the fit's.", call. = FALSE) + } + invisible(NULL) +} + #' 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 +411,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 diff --git a/tests/testthat/test-powerLaw.R b/tests/testthat/test-powerLaw.R new file mode 100644 index 0000000..e5094cf --- /dev/null +++ b/tests/testthat/test-powerLaw.R @@ -0,0 +1,85 @@ +## 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 with a year effect warns, once, that fireSense_burn adds it on the logit", { + ins <- inputsPL() + ins$studyAreaWithSpreadParams$params <- lapply(ins$studyAreaWithSpreadParams$params, function(p) { p$yearSpreadSD <- 0.3; p }) + expect_warning(sim <- toyRun(ins), "fireSense_burn adds fireSense_SpreadSD on the") + expect_equal(predVals(sim)[cellsPL], pHandPL(1.9), tolerance = 1e-9) + expect_equal(sim$fireSense_SpreadSD, 0.3) + ## no warning without a year effect, or with the default link + expect_no_warning(toyRun(inputsPL())) + ins2 <- toyInputs() + ins2$studyAreaWithSpreadParams$params <- lapply(ins2$studyAreaWithSpreadParams$params, function(p) { p$yearSpreadSD <- 0.3; p }) + expect_no_warning(toyRun(ins2)) +}) From 084b8fac93a1ecbc7af4045fb56b257215bccff5 Mon Sep 17 00:00:00 2001 From: Eliot McIntire Date: Sun, 18 Oct 2026 17:49:03 -0700 Subject: [PATCH 2/2] Output fireSense_SpreadPowerLaw so fireSense_burn adds the year effect to x; drop the warning Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014G18Mik3fH4z5M5CyMSUS6 --- NEWS.md | 2 +- fireSense_spreadPredict.R | 54 ++++++++++++++++++++++------------ tests/testthat/test-metadata.R | 3 +- tests/testthat/test-powerLaw.R | 37 +++++++++++++++++------ 4 files changed, 67 insertions(+), 29 deletions(-) diff --git a/NEWS.md b/NEWS.md index ee6bce2..6bb279b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ # 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. `fireSense_burn` still adds the year effect (`fireSense_SpreadSD`) on the logit, whereas a power-law fit adds it to the linear predictor, so a power-law row with `yearSpreadSD > 0` gives a warning (once) until `fireSense_burn` follows. Needs fireSenseUtils >= 0.2.3.9107. +- 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. diff --git a/fireSense_spreadPredict.R b/fireSense_spreadPredict.R index 598a19a..bac4e2f 100644 --- a/fireSense_spreadPredict.R +++ b/fireSense_spreadPredict.R @@ -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`.")) )) ) @@ -152,13 +158,13 @@ spreadPredictRun <- function(sim) { if (NROW(sa) == 1L) { params1 <- paramSetForRep(sa$params[[1]], P(sim)$.rep) gamma1 <- ledgerPowerLawGamma(sa, 1L) - warnPowerLawYearEffect(gamma1, params1) 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), 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 @@ -171,12 +177,14 @@ 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 <- ledgerPowerLawGamma(sa, i) - warnPowerLawYearEffect(gammaI, parsI) + 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, @@ -185,6 +193,10 @@ spreadPredictRun <- function(sim) { 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 @@ -192,6 +204,13 @@ spreadPredictRun <- function(sim) { ## 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 @@ -356,23 +375,22 @@ spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSp yearSpreadSDTxt <- "yearSpreadSD" -#' Warn, once, that `fireSense_burn` adds the year effect on the wrong scale for a power-law fit +#' 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, not to -#' the logit of `p`; `fireSense_burn` adds `fireSense_SpreadSD` on the logit, so with a power-law fit and a -#' `yearSpreadSD` above 0 the simulated years are not those of the 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 params one-row `data.frame`: the parameter set in use (`paramSetForRep()`). -#' @return `NULL`, invisibly. -warnPowerLawYearEffect <- function(powerLawGamma, params) { - if (!is.null(powerLawGamma) && yearSpreadSDOf(params) > 0 && !isTRUE(mod$warnedPowerLawYearEffect)) { - mod$warnedPowerLawYearEffect <- TRUE - warning("fireSense_spreadPredict: this ledger row is a spreadLink = 'powerLaw' fit with yearSpreadSD > 0; ", - "the fit adds the year effect to the linear predictor, but fireSense_burn adds fireSense_SpreadSD on the ", - "logit of the spread probability, so simulated years will differ from the fit's.", call. = FALSE) - } - invisible(NULL) +#' @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 diff --git a/tests/testthat/test-metadata.R b/tests/testthat/test-metadata.R index 71f7048..919ac87 100644 --- a/tests/testthat/test-metadata.R +++ b/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") ) }) diff --git a/tests/testthat/test-powerLaw.R b/tests/testthat/test-powerLaw.R index e5094cf..e52edfe 100644 --- a/tests/testthat/test-powerLaw.R +++ b/tests/testthat/test-powerLaw.R @@ -71,15 +71,34 @@ test_that("with several ELFs each row applies its own link, blended like the pro expect_equal(v[10], pB, tolerance = 1e-9) }) -test_that("a powerLaw row with a year effect warns, once, that fireSense_burn adds it on the logit", { - ins <- inputsPL() +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_warning(sim <- toyRun(ins), "fireSense_burn adds fireSense_SpreadSD on the") - expect_equal(predVals(sim)[cellsPL], pHandPL(1.9), tolerance = 1e-9) + expect_no_warning(sim <- toyRun(ins)) + expect_equal(predVals(sim)[cellsPL], pHandPL(3), tolerance = 1e-9) expect_equal(sim$fireSense_SpreadSD, 0.3) - ## no warning without a year effect, or with the default link - expect_no_warning(toyRun(inputsPL())) - ins2 <- toyInputs() - ins2$studyAreaWithSpreadParams$params <- lapply(ins2$studyAreaWithSpreadParams$params, function(p) { p$yearSpreadSD <- 0.3; p }) - expect_no_warning(toyRun(ins2)) + 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])) })