diff --git a/NEWS.md b/NEWS.md index 654b4db..3841e5e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # fireSense_spreadPredict (development version) +- 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. - A single fitted ELF now predicts with the covariates its fitted parameters name, as several ELFs already did, not the terms of `sim$fireSense_spreadFormula`. A hindcast whose formula had a term the fit did not (e.g. `treedWetland_agb`) stopped with "not found in data objects". Predictions where the formula matched the fit do not change. A fitted coefficient with no covariate still stops, naming it. diff --git a/fireSense_spreadPredict.R b/fireSense_spreadPredict.R index 4652677..f5ffcd4 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.9084)"), + "PredictiveEcology/fireSenseUtils@development (>= 0.2.3.9106)"), parameters = bindrows( defineParameter(name = "lowerSpreadProb", class = "numeric", default = fireSenseUtils::spreadProbFloor, desc = paste("Lower asymptote of the 2- and 3-parameter logistic. Default", @@ -154,7 +154,7 @@ spreadPredictRun <- function(sim) { 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)) + covCentre = ledgerCovCentre(sa, 1L), climateRef = ledgerClimateRef(sa, 1L)) sim$fireSense_SpreadSD <- yearSpreadSDOf(params1) } else { ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]]) @@ -175,7 +175,7 @@ spreadPredictRun <- function(sim) { 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)) + covCentre = ledgerCovCentre(sa, i), climateRef = ledgerClimateRef(sa, i)) 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) @@ -247,9 +247,13 @@ ELFblendWeights <- function(elfRas, template, pixelID, ids, width) { #' from the rescaled covariates; `NULL` (a fit without an intercept, or an older ledger row) centres nothing. #' A parameter set with `fireSenseUtils::spreadInterceptTxt` has an intercept, added as a column of 1s after #' the covariates are rescaled and centred, as in the fit. +#' @param climateRef the ELF's `climateRef_spread` (`ledgerClimateRef()`): with a fit made with climate outside the +#' `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. #' @return `data.table` with `pixelID` and `spreadProb`, from the one parameter set in `params`. spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSpreadProb, - covCentre = NULL) { + covCentre = NULL, climateRef = NULL) { moduleName <- "fireSense_spreadPredict" covs <- copy(covs) ## the per-year random effect is not a covariate coefficient: fireSense_burn applies it (fireSense_SpreadSD) @@ -333,8 +337,10 @@ spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSp # Make sure the order is correct in the matrix matching <- intersect(names(covPars), colnames(mat)) + ## the one function the fit computes its spread probability with, so the prediction cannot differ from it data.table(pixelID = shortAnnDT$pixelID, - spreadProb = logisticAll(logisticPars, mat[, matching, drop = FALSE], covPars, lowerSpreadProb)) + spreadProb = spreadProbFromCovs(mat[, matching, drop = FALSE], logisticPars, covPars, lowerSpreadProb, + climateRef = climateRef)) } yearSpreadSDTxt <- "yearSpreadSD" @@ -357,6 +363,24 @@ ledgerCovMinMax <- function(sa, i, default = NULL) { cmm } +#' The climate reference of one ledger row +#' +#' A fit made with climate outside the `maxAsymptote` ceiling (`fireSense_spreadFit`'s `spreadLink = "climateOutside"`) +#' stores the reference value of each climate covariate in the ledger row (`fireSenseUtils::spreadFitClimateRefTxt`), +#' with `spreadLink` beside it. Prediction applies the same link with the same reference +#' (`fireSenseUtils::spreadProbFromCovs()`). +#' +#' @param sa `sim$studyAreaWithSpreadParams`, one row per ELF. +#' @param i integer; the row. +#' @return A named list, or `NULL` for a row without a reference: a fit with climate under the ceiling, or a row +#' written before the link existed. Such a row predicts exactly as it always did. +ledgerClimateRef <- function(sa, i) { + col <- fireSenseUtils::spreadFitClimateRefTxt + if (!col %in% names(sa)) return(NULL) + ref <- sa[[col]][[i]] + if (is.list(ref) && length(ref)) ref else NULL +} + #' 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-climateOutside.R b/tests/testthat/test-climateOutside.R new file mode 100644 index 0000000..ab47a35 --- /dev/null +++ b/tests/testthat/test-climateOutside.R @@ -0,0 +1,76 @@ +## A ledger row from a fit made with climate outside the maxAsymptote ceiling (fireSense_spreadFit's +## spreadLink = "climateOutside") carries, in fireSenseUtils::spreadFitClimateRefTxt, the fit's reference value of +## each climate covariate (on the rescaled scale), and "climateOutside" in fireSenseUtils::spreadLinkTxt. Prediction +## applies the same link through fireSenseUtils::spreadProbFromCovs(), the function the fit used: +## p0 = logistic of the other covariates, under the ceiling; p = plogis(qlogis(p0) + beta_clim * (climate - ref)), +## uncapped. A row without a reference predicts as it always did. +## +## The toy covariates are those of test-spreadProb-values.R: MDC / 200 is the climate covariate, youngAge and +## fuelA / 8 the rest; the toy parameters have hillSlope1 = 2, so p0 = 0.13 + 0.12 * plogis(2 * x_fuel), and the +## climate shift is beta_clim * (MDC / 200 - ref) with no slope factor (the fit fixes hillSlope1 at 1). + +refCol <- fireSenseUtils::spreadFitClimateRefTxt +linkCol <- fireSenseUtils::spreadLinkTxt +cellsCO <- c(1, 2, 3, 4, 6, 7, 8) + +inputsCO <- function(ref = list(MDC = 0.4), link = "climateOutside") { + ins <- toyInputs() + if (!is.null(link)) ins$studyAreaWithSpreadParams[[linkCol]] <- link + ins$studyAreaWithSpreadParams[[refCol]] <- list(ref) + ins +} + +## by hand, in cell order +pHandCO <- function(ref = 0.4) { + cv <- toyCovariates(); o <- order(cv$pixelID) + mdc <- cv$MDC[o] / 200; ya <- cv$youngAge[o]; fuel <- cv$fuelA[o] / 8 + p0 <- handLogistic3(-0.5 * ya + 1 * fuel) + plogis(qlogis(p0) + 1.5 * (mdc - ref)) +} + +test_that("a row with a climate reference predicts climate on the logit after the fuel ceiling, uncapped", { + v <- predVals(toyRun(inputsCO())) + expect_equal(v[cellsCO], pHandCO(0.4), tolerance = 1e-7) + expect_true(all(is.na(v[c(5, 9)]))) + ## cell 4 (MDC 200, fuelA 8): p0 = 0.13 + 0.12 * plogis(2) = 0.2358, then + 1.5 * (1 - 0.4) on the logit + expect_equal(v[4], plogis(qlogis(handLogistic3(1)) + 0.9), tolerance = 1e-7) + expect_gt(v[4], 0.25) # above maxAsymptote: the ceiling no longer binds for climate + ## cell 1 (MDC 0, no fuel): p0 = 0.19, then 1.5 * (0 - 0.4) = -0.6 on the logit: below p0 + expect_equal(v[1], plogis(qlogis(0.19) - 0.6), tolerance = 1e-7) + ## at climate = the reference the new link is p0: cell 3 has MDC 100 -> 0.5, so with ref 0.5 it is the + ## fuel-only logistic (fuelA 4 -> 0.5, x_fuel = 0.5) + v2 <- predVals(toyRun(inputsCO(ref = list(MDC = 0.5)))) + expect_equal(v2[3], handLogistic3(0.5), tolerance = 1e-7) + expect_equal(v2[cellsCO], pHandCO(0.5), tolerance = 1e-7) +}) + +test_that("a row without a reference predicts exactly as before", { + base <- predVals(toyRun()) + ## the column is there but empty: a fit with climate under the ceiling, written beside one without + expect_identical(predVals(toyRun(inputsCO(ref = list(), link = "ceiling"))), base) + expect_identical(predVals(toyRun(inputsCO(ref = NULL, link = NULL))), base) + ## and the reference changes the prediction only where climate differs from it + withRef <- predVals(toyRun(inputsCO(ref = list(MDC = 0.5)))) + expect_false(isTRUE(all.equal(withRef, base))) + sa <- toyInputs()$studyAreaWithSpreadParams + expect_null(ledgerClimateRef(sa, 1L)) + sa[[refCol]] <- list(list()) + expect_null(ledgerClimateRef(sa, 1L)) + sa[[refCol]] <- list(list(MDC = 0.4)) + expect_identical(ledgerClimateRef(sa, 1L), list(MDC = 0.4)) +}) + +test_that("with several ELFs each row applies its own reference, blended like the probabilities", { + ins <- multiInputs() + ## a climate effect for both ELFs; MDC is 100 everywhere, i.e. 0.5 rescaled + ins$studyAreaWithSpreadParams$params <- lapply(ins$studyAreaWithSpreadParams$params, function(p) { p$MDC <- 1; p }) + ins$studyAreaWithSpreadParams[[linkCol]] <- "climateOutside" + ins$studyAreaWithSpreadParams[[refCol]] <- list(list(MDC = 0.5), list(MDC = 0.25)) + v <- predVals(toyRun(ins, params = list(ELFblendWidth = 20000))) + pA <- 0.19 # at its reference: the fuel-only value (all other coefficients 0) + pB <- plogis(qlogis(0.20) + 1 * (0.5 - 0.25)) # 0.25 above its reference on the logit + 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) +})