diff --git a/NEWS.md b/NEWS.md index c57f0d4..a760895 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# fireSense_SpreadPredict (development version) + +- 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 1.0.0 First release from `development` since `master` was last updated (2021-01-27). Full history: https://github.com/PredictiveEcology/fireSense_SpreadPredict/compare/a5b41f9...v1.0.0 diff --git a/fireSense_SpreadPredict.R b/fireSense_SpreadPredict.R index ef9e09f..3a0c6d3 100644 --- a/fireSense_SpreadPredict.R +++ b/fireSense_SpreadPredict.R @@ -10,17 +10,22 @@ defineModule(sim, list( person("Alex M.", "Chubaty", email = "achubaty@for-cast.ca", role = "ctb") ), childModules = character(), - version = list(fireSense_SpreadPredict = "1.0.0.9001", SpaDES.core = "0.1.0"), + version = list(fireSense_SpreadPredict = "1.0.0.9003", SpaDES.core = "0.1.0"), timeframe = as.POSIXlt(c(NA, NA)), timeunit = "year", citation = list("citation.bib"), documentation = list("README.txt", "fireSense_SpreadPredict.Rmd"), reqdPkgs = list("magrittr", "Matrix", "methods", "terra", "SpaDES.core (>=3.0.4)", "stats", "ggplot2", "viridis", - "PredictiveEcology/fireSenseUtils@development (>= 0.2.3.9029)"), + "PredictiveEcology/fireSenseUtils@development (>= 0.2.3.9038)"), parameters = bindrows( defineParameter(name = "lowerSpreadProb", class = "numeric", default = 0.13, desc = "Lower asymptote of the 2- and 3-parameter logistic."), + defineParameter("ELFblendWidth", "numeric", default = 20000, + desc = paste("With several ELFs: each ELF's model also predicts this far (m) outside its own", + "pixels, and where predictions overlap they are averaged with weights that fall", + "linearly from 1 inside the ELF to 0 at this distance outside it. 50/50 at a", + "boundary. The default is the buffer fireSenseUtils::makeELFs() puts around ELFs.")), 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.")), @@ -41,6 +46,9 @@ defineModule(sim, list( expectsInput(objectName = "fireSense_SpreadCovariates", objectClass = "data.table", desc = paste("This year's covariates, from `fireSense_dataPrepPredict`.", "`pixelID` is the cell index of `flammableRTM`.")), + expectsInput(objectName = "rasterToMatchLargeELF", objectClass = "SpatRaster", sourceURL = NA, + desc = paste("Only with several fitted ELFs: each pixel's ELF (`ELFind`), on the grid of", + "`flammableRTM`, from `fireSense_ELFs` with a `studyAreaLarge`.")), expectsInput(objectName = "flammableRTM", objectClass = "SpatRaster", sourceURL = NA, desc = "Binary raster, 1 where the pixel is flammable. Template for `fireSense_SpreadPredicted`.") ), @@ -99,105 +107,186 @@ doEvent.fireSense_SpreadPredict <- function(sim, eventTime, eventType, debug = F #' Predict this year's spread probability #' -#' Rescales `sim$fireSense_SpreadCovariates` with `sim$covMinMax_spread`, computes the spread -#' probability for each parameter set (row) in `sim$studyAreaWithSpreadParams$params[[1]]`, -#' and writes the mean over parameter sets to `sim$fireSense_SpreadPredicted`. +#' 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`. +#' +#' 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 +#' 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()`. #' #' @param sim A `simList`. #' #' @return The `simList`, invisibly. spreadPredictRun <- function(sim) { - moduleName <- current(sim)$moduleName + covs <- copy(sim$fireSense_SpreadCovariates) + sa <- sim$studyAreaWithSpreadParams - fireSense_SpreadCovariates <- copy(sim$fireSense_SpreadCovariates) + # Without fitted parameters there is nothing to predict from; say so instead of + # dying in rowMeans() on an empty matrix (which is what an unfitted ELF produced + # when fireSense_SpreadFit had not run first). This must come before anything + # indexes `params[[1]]`: with zero rows that fails first, "subscript out of bounds". + nPar <- tryCatch(NROW(sa$params[[1]]), error = function(e) 0L) + if (NROW(sa) == 0L || is.null(nPar) || nPar == 0L) + stop("fireSense_SpreadPredict: sim$studyAreaWithSpreadParams holds no fitted spread ", + "parameters for this run (", if (!is.null(sim$.runName)) sim$.runName else "unknown", + "). Either fireSense_SpreadFit has not run yet -- its `run` event must precede this ", + "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, + formula = sim$fireSense_spreadFormula, yr = time(sim), + maxFireSpread = P(sim)$maxFireSpread, lowerSpreadProb = P(sim)$lowerSpreadProb) + } else { + ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]]) + ## each ELF's weight at each pixel: static, so computed once + if (is.null(mod$ELFweights) || !identical(attr(mod$ELFweights, "key"), list(ids, covs$pixelID))) + mod$ELFweights <- ELFblendWeights(sim$rasterToMatchLargeELF, sim$flammableRTM, covs$pixelID, ids, + width = P(sim)$ELFblendWidth) + w <- mod$ELFweights + none <- rowSums(w) == 0 + 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)) + 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]], + 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 + wsum[rows] <- wsum[rows] + w[rows, i] + } + ok <- wsum > 0 + pred <- data.table(pixelID = covs$pixelID[ok], spreadProb = acc[ok] / wsum[ok]) + } + + # Return to raster format + sim$fireSense_SpreadPredicted <- rast(sim$flammableRTM) ## use flammableRTM as template + sim$fireSense_SpreadPredicted[pred$pixelID] <- pred$spreadProb + + invisible(sim) +} + +#' Each ELF's weight at each pixel +#' +#' An ELF's raw weight is 1 on its own pixels and falls linearly to 0 at `width` metres from them (distance +#' to the nearest pixel of the ELF). Weights are then normalised to sum to 1 at each pixel, so a pixel on a +#' boundary between two ELFs is 50/50 and one `width` inside an ELF uses that ELF only. +#' +#' @param elfRas `sim$rasterToMatchLargeELF`: `ELFind` per pixel, as a categorical raster (or the id). +#' @param template `sim$flammableRTM`; `pixelID` indexes its cells, so `elfRas` must share its grid. +#' @param pixelID integer; the cells to weight. +#' @param ids character; the ELFs, in the order of the rows of `studyAreaWithSpreadParams`. +#' @param width numeric; metres. +#' @return matrix, one row per `pixelID`, one column per ELF (`ids`), rows summing to 1 (or 0 where no +#' ELF is within `width`). `attr(, "key")` records `ids` and `pixelID`. +ELFblendWeights <- function(elfRas, template, pixelID, ids, width) { + if (is.null(elfRas)) + stop("fireSense_SpreadPredict: sim$studyAreaWithSpreadParams has several ELFs, so each pixel's ELF ", + "must come from sim$rasterToMatchLargeELF (fireSense_ELFs, with a studyAreaLarge); it is missing", + call. = FALSE) + if (!isTRUE(terra::compareGeom(elfRas, template, stopOnError = FALSE))) + stop("fireSense_SpreadPredict: sim$rasterToMatchLargeELF is not on the grid of sim$flammableRTM", + call. = FALSE) + r <- elfRas[[1]] + v <- terra::values(r, mat = FALSE) + lv <- terra::levels(r)[[1]] + lab <- if (is.data.frame(lv) && NCOL(lv) >= 2) as.character(lv[[2]][match(v, lv[[1]])]) else as.character(v) + w <- vapply(ids, function(id) { + core <- r; terra::values(core) <- ifelse(lab %in% id, 1, NA) + if (all(is.na(terra::values(core, mat = FALSE)))) return(numeric(length(pixelID))) + d <- terra::values(terra::distance(core), mat = FALSE)[pixelID] + pmax(0, 1 - d / width) + }, numeric(length(pixelID))) + w <- matrix(w, ncol = length(ids), dimnames = list(NULL, ids)) + tot <- rowSums(w) + w[tot > 0, ] <- w[tot > 0, , drop = FALSE] / tot[tot > 0] + attr(w, "key") <- list(ids, pixelID) + w +} + +#' Spread probability for the pixels of one ELF +#' +#' @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 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. +spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread, lowerSpreadProb, + byParams = FALSE) { + moduleName <- "fireSense_SpreadPredict" + covs <- copy(covs) ## 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 ## coefficients only mean anything for biomass / 1e4: undo the log with the function the fit used. ## A fit made on the log scale has the log range there, and its covariates are left as they are. - for (cn in intersect(names(sim$covMinMax_spread), names(fireSense_SpreadCovariates))) { - if (fireSenseUtils::isLinearFuelRange(sim$covMinMax_spread[[cn]])) - fireSense_SpreadCovariates[[cn]] <- fireSenseUtils::fuelLogToLinear(fireSense_SpreadCovariates[[cn]]) + for (cn in intersect(names(covMinMax), names(covs))) { + if (fireSenseUtils::isLinearFuelRange(covMinMax[[cn]])) + covs[[cn]] <- fireSenseUtils::fuelLogToLinear(covs[[cn]]) } - # Load inputs in the data container - mod_env <- new.env(parent = globalenv()) - list2env(fireSense_SpreadCovariates, envir = mod_env) - ## In case there is a response in the formula remove it - - terms <- as.formula(sim$fireSense_spreadFormula) %>% - terms.formula() %>% - delete.response() - - formula <- reformulate(attr(terms, "term.labels"), intercept = attr(terms, "intercept")) - allxy <- all.vars(formula) - - missing <- !allxy %in% ls(mod_env, all.names = TRUE) + ## the covariates this ELF was fitted with + needed <- if (isTRUE(byParams)) { + setdiff(names(params), unlist(fireSenseUtils::logisticParamNames)) + } else { + terms <- delete.response(terms.formula(as.formula(formula))) + all.vars(reformulate(attr(terms, "term.labels"), intercept = attr(terms, "intercept"))) + } + missing <- !needed %in% names(covs) if (s <- sum(missing)) { stop( - moduleName, "> '", allxy[missing][1L], "'", + moduleName, "> '", needed[missing][1L], "'", if (s > 1) paste0(" (and ", s - 1L, " other", if (s > 2) "s", ")"), " not found in data objects." ) } + ## with several ELFs the table holds every ELF's covariates; use this ELF's only + if (isTRUE(byParams)) covs <- covs[, c("pixelID", needed), with = FALSE] # integers x 1000, the form `spreadProbFromIntegerCovs` expects - shortAnnDTx1000 <- toX1000(list(fireSense_SpreadCovariates))[[1]] |> setDT() - colsToUse <- setdiff(names(fireSense_SpreadCovariates), "pixelID") - - # Without fitted parameters there is nothing to predict from; say so instead of - # dying in rowMeans() on an empty matrix (which is what an unfitted ELF produced - # when fireSense_SpreadFit had not run first). This must come before anything - # indexes `params[[1]]`: with zero rows that fails first, "subscript out of bounds". - nPar <- tryCatch(NROW(sim$studyAreaWithSpreadParams$params[[1]]), error = function(e) 0L) - if (NROW(sim$studyAreaWithSpreadParams) == 0L || is.null(nPar) || nPar == 0L) - stop("fireSense_SpreadPredict: sim$studyAreaWithSpreadParams holds no fitted spread ", - "parameters for this run (", if (!is.null(sim$.runName)) sim$.runName else "unknown", - "). Either fireSense_SpreadFit has not run yet -- its `run` event must precede this ", - "module's -- or the shared ledger has no row for this polygon.", call. = FALSE) + shortAnnDTx1000 <- toX1000(list(covs))[[1]] |> setDT() + colsToUse <- setdiff(names(covs), "pixelID") - logisticPars <- sim$studyAreaWithSpreadParams$params[[1]] - shortAnnDT <- spreadProbFromIntegerCovs(shortAnnDTx1000 = shortAnnDTx1000, - yr = time(sim), - covMinMax = sim$covMinMax_spread, + yr = yr, + covMinMax = covMinMax, mutuallyExclusive = NULL, # alraedy done in dataPrepPredict colsToUse = colsToUse, doAssertions = FALSE, - logisticPars = logisticPars, - maxFireSpread = Par$maxFireSpread + logisticPars = params, + maxFireSpread = maxFireSpread ) - parsModel <- length(colsToUse) mat <- as.matrix(shortAnnDT[, ..colsToUse]) # for replicate "best" params from DEoptim - spreadProbList <- purrr::pmap(.l = list(ind = seq(NROW(sim$studyAreaWithSpreadParams$params[[1]]))), - sa = sim$studyAreaWithSpreadParams, function(ind, sa) { - par <- sa$params[[1]][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") - mat <- mat[, matching] - - logisticAll(logisticPars, - mat, covPars, P(sim)$lowerSpreadProb) - }) + 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) - - set(shortAnnDT, NULL, "spreadProb", rowMeans(spreadProbMat)) - # Return to raster format - sim$fireSense_SpreadPredicted <- rast(sim$flammableRTM) ## use flammableRTM as template - sim$fireSense_SpreadPredicted[shortAnnDT$pixelID] <- shortAnnDT$spreadProb - - invisible(sim) + data.table(pixelID = shortAnnDT$pixelID, spreadProb = rowMeans(spreadProbMat)) } diff --git a/tests/testthat/test-metadata.R b/tests/testthat/test-metadata.R index 90f3616..9ba88d7 100644 --- a/tests/testthat/test-metadata.R +++ b/tests/testthat/test-metadata.R @@ -19,7 +19,8 @@ test_that("inputs are the expected names and classes", { inputs[order(names(inputs))], c(covMinMax_spread = "data.table", fireSense_SpreadCovariates = "data.table", - flammableRTM = "SpatRaster") + flammableRTM = "SpatRaster", + rasterToMatchLargeELF = "SpatRaster") ) }) @@ -46,11 +47,11 @@ 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", "logical", "numeric", "numeric"), + class = c("numeric", "numeric", "numeric", "logical", "numeric", "numeric", "numeric"), ## .runInitialTime defaults to start(sim), which is 0 when only metadata is parsed - default = c("0", "1", "NA", "FALSE", "0.13", "0.28"), + default = c("0", "1", "NA", "FALSE", "20000", "0.13", "0.28"), row.names = c(".runInitialTime", ".runInterval", ".saveInitialTime", ".useCache", - "lowerSpreadProb", "maxFireSpread") + "ELFblendWidth", "lowerSpreadProb", "maxFireSpread") ) expect_identical(paramTable(md), expected) }) diff --git a/tests/testthat/test-multiELF.R b/tests/testthat/test-multiELF.R new file mode 100644 index 0000000..15cf7e5 --- /dev/null +++ b/tests/testthat/test-multiELF.R @@ -0,0 +1,51 @@ +## Several fitted ELFs in one study area (the 2-ELF Mackenzie forecast, 2026-09). Every pixel gets a spread +## probability. Each ELF's model predicts its own pixels and those within ELFblendWidth of them; where two +## overlap they are averaged with weights falling linearly from 1 inside an ELF to 0 at ELFblendWidth +## outside it, normalised to sum to 1 (Eliot, 2026-09-22). +## +## Toy: one row of 10 pixels, 5 km wide; ELF "A" is columns 1-5, "B" columns 6-10. All coefficients are 0, +## 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 + rA <- c(1, 1, 1, 1, 1, 0.75, 0.5, 0.25, 0, 0); rB <- rev(rA) + expect_equal(v, (rA * 0.19 + rB * 0.20) / (rA + rB), tolerance = 1e-9) + expect_equal(v[1], 0.19, tolerance = 1e-9) # 25 km from B: A only + expect_equal(v[10], 0.20, tolerance = 1e-9) + expect_equal(v[5], (0.19 + 0.75 * 0.20) / 1.75, tolerance = 1e-9) + expect_true(all(is.finite(v))) # every pixel has a probability +}) + +test_that("a narrow blend width leaves only the boundary pixels mixed", { + v <- predVals(toyRun(multiInputs(), params = list(ELFblendWidth = 6000))) + ## 5 km from the other ELF: raw weight 1 - 5/6 + expect_equal(v[c(1:4, 7:10)], rep(c(0.19, 0.20), each = 4), tolerance = 1e-9) + expect_equal(v[5], (0.19 + (1/6) * 0.20) / (1 + 1/6), tolerance = 1e-9) +}) + +test_that("with several ELFs, a missing ELF raster stops with a message that says what is needed", { + ins <- multiInputs(); ins$rasterToMatchLargeELF <- NULL + expect_error(toyRun(ins), "rasterToMatchLargeELF") +}) diff --git a/tests/testthat/test-spreadProb-values.R b/tests/testthat/test-spreadProb-values.R index 5e84705..573de75 100644 --- a/tests/testthat/test-spreadProb-values.R +++ b/tests/testthat/test-spreadProb-values.R @@ -195,3 +195,23 @@ test_that("a fit made on the log scale is predicted on the log scale, as before" ## (x = 1.25 -> 0.2409) but (8.517 - 3.605) / 6.395 = 0.768 on the log scale (x = 1.518 -> 0.2445) expect_gt(abs(v[3] - handLogistic3(1.25)), 0.003) }) + +test_that("a fit with upperTail1 is predicted with the upper-tail link, from the parameter's name", { + ## fireSense_SpreadFit's link "logistic3pUpper" stores one more logistic parameter, upperTail1 + ## (fireSenseUtils::logistic3pUpper(), Stukel 1988). It changes the curve only where + ## hillSlope1 * x > 0: there u = hillSlope1 * x becomes -log(1 - a * u) / a for a < 0. + v <- predVals(toyRun(toyInputs(params = toyParams(upperTail1 = -0.5)))) + handUpper <- function(x, a = -0.5) { + u <- 2 * x + u[u > 0] <- -log(1 - a * u[u > 0]) / a + 0.13 + 0.12 / (1 + exp(-u)) + } + ## cell 4: x = 2.5, u = 5 -> 2 * log(3.5) = 2.505526 -> 0.13 + 0.12 * 0.9245283 = 0.2409434 + expect_equal(v[4], 0.2409434, tolerance = 1e-6) + expect_equal(v[c(1, 3, 4, 6, 7, 8)], handUpper(c(0, 1.25, 2.5, 1, 1.125, 1.25)), tolerance = 1e-7) + ## below the inflection nothing changes: cell 2 (x = -0.125) is the logistic3p value + expect_equal(v[2], 0.1825388, tolerance = 1e-6) + ## and upperTail1 = 0 is logistic3p everywhere + v0 <- predVals(toyRun(toyInputs(params = toyParams(upperTail1 = 0)))) + expect_equal(v0, predVals(toyRun()), tolerance = 1e-12) +})