R/spread.R

Defines functions spread

Documented in spread

utils::globalVariables(c(
  ".", ".I", "dists", "dup", "id", "indices", "initialLocus"
))

###############################################################################
#' Simulate a spread process on a landscape.
#'
#' This can be used to simulate fires, seed dispersal, calculation of iterative,
#' concentric landscape values (symmetric or asymmetric) and many other things.
#' Essentially, it starts from a collection of cells (`loci`) and spreads
#' to neighbours, according to the `directions` and `spreadProb` arguments.
#' This can become quite general, if `spreadProb` is 1 as it will expand
#' from every loci until all cells in the landscape have been covered.
#' With `id` set to `TRUE`, the resulting map will be classified
#' by the index of the cell where that event propagated from.
#' This can be used to examine things like fire size distributions.
#' **NOTE:** See also [spread2()], which is more robust and can be
#' used to build custom functions.
#' However, under some conditions, this `spread` function is faster.
#' The two functions can accomplish many of the same things, and key differences are internal.
#'
#' For large rasters, a combination of `lowMemory = TRUE` and
#' `returnIndices = TRUE` or `returnIndices = 2`
#' will be fastest and use the least amount of memory.
#' **2022-07-25: `lowMemory = TRUE` is deprecated due to removal of package `ffbase` from CRAN.**
#'
#' This function can be interrupted before all active cells are exhausted if
#' the `iterations` value is reached before there are no more active
#' cells to spread into. If this is desired, `returnIndices` should be
#' `TRUE` and the output of this call can be passed subsequently as an input
#' to this same function. This is intended to be used for situations where external
#' events happen during a spread event, or where one or more arguments to the spread
#' function change before a spread event is completed. For example, if it is
#' desired that the `spreadProb` change before a spread event is completed because,
#' for example, a fire is spreading, and a new set of conditions arise due to
#' a change in weather.
#'
#' `asymmetry` is currently used to modify the `spreadProb` in the following way.
#' First for each active cell, spreadProb is converted into a length 2 numeric of Low and High
#' spread probabilities for that cell:
#' `spreadProbsLH <- (spreadProb*2) // (asymmetry+1)*c(1,asymmetry)`,
#' whose ratio is equal to
#' `asymmetry`.
#' Then, using `asymmetryAngle`, the angle between the
#' initial starting point of the event and all potential
#' cells is found. These are converted into a proportion of the angle from
#' `-asymmetryAngle`
#' to
#' `asymmetryAngle`
#' using:
#' `angleQuality <- (cos(angles - rad2(asymmetryAngle))+1)/2`
#' where `rad2 <- function (degree) (degree * pi)/180`
#'
#' These are then converted to multiple `spreadProbs` by
#' `spreadProbs <- lowSpreadProb + (angleQuality * diff(spreadProbsLH))`
#' To maintain an expected `spreadProb` that is the same as the asymmetric
#' `spreadProbs`, these are then rescaled so that the mean of the
#' asymmetric spreadProbs is always equal to spreadProb at every iteration:
#' `spreadProbs <- spreadProbs - diff(c(spreadProb, mean(spreadProbs)))`
#'
#' @section Breaking out of spread events:
#'
#' There are 4 ways for the spread to "stop" spreading. Here, each "event" is defined as
#' all cells that are spawned from a single starting loci. So, one spread call can have
#' multiple spreading "events". The ways outlines below are all acting at all times,
#' i.e., they are not mutually exclusive. Therefore, it is the user's
#' responsibility to make sure the different rules are interacting with
#' each other correctly. Using `spreadProb` or `maxSize` are computationally
#' fastest, sometimes dramatically so.
#'
#' \tabular{ll}{
#'   `spreadProb` \tab Probabilistically, if spreadProb is low enough,
#'                          active spreading events will stop. In practice,
#'                          active spreading events will stop. In practice,
#'                          this number generally should be below 0.3 to actually
#'                          see an event stop\cr
#'   `maxSize` \tab This is the number of cells that are "successfully" turned
#'                       on during a spreading event. This can be vectorized, one value
#'                       for each event   \cr
#'   `circleMaxRadius` \tab If `circle` is TRUE, then this will be the maximum
#'                       radius reached, and then the event will stop. This is
#'                       vectorized, and if length is >1, it will be matched
#'                       in the order of `loci`\cr
#'   `stopRule` \tab This is a function that can use "landscape", "id", "cells",
#'                       or any named vector passed into `spread` in the `...`.
#'                       This can take on relatively complex functions.
#'                       Passing in, say, a `RasterLayer` to `spread`
#'                       can access the individual values on that arbitrary
#'                       `RasterLayer` using "cells".
#'                       These will be calculated within all the cells of the individual
#'                       event (equivalent to a "group_by(event)" in `dplyr`.
#'                       So, `sum(arbitraryRaster[cells])` would sum up all
#'                       the raster values on the `arbitraryRaster` raster
#'                       that are overlaid by the individual event.
#'                       This can then be used in a logical statement. See examples.
#'                       To confirm the cause of stopping, the user can assess the values
#'                       after the function has finished.\cr
#' }
#'
#' The spread function does not return the result of this `stopRule`.
#' If, say, an event has both `circleMaxRadius` and `stopRule`,
#' and it is the `circleMaxRadius` that caused the event spreading to stop,
#' there will be no indicator returned from this function that indicates
#' which rule caused the stop.
#'
#' `stopRule` has many use cases. One common use case is evaluating
#' a neighbourhood around a focal set of points. This provides,
#' therefore, an alternative to the [terra::buffer()] function or
#' [terra::focal()] function.
#' In both of those cases, the window/buffer size must be an input to the function. Here,
#' the resulting size can be emergent based on the incremental growing and calculating
#' of the `landscape` values underlying the spreading event.
#'
#' @section `stopRuleBehavior`:
#' This determines how the `stopRule` should be implemented. Because
#' spreading occurs outwards in concentric circles or shapes, one cell width at a time, there
#' are 4 possible ways to interpret the logical inequality defined in `stopRule`.
#' In order of number of cells included in resulting events, from most cells to fewest cells:
#'
#' \tabular{ll}{
#'   `"includeRing"` \tab Will include the entire ring of cells that, as a group,
#'                             caused `stopRule` to be `TRUE`.\cr
#'   `"includePixel"` \tab Working backwards from the entire ring that caused the
#'                              `stopRule` to be `TRUE`, this will iteratively
#'                              random cells in the final ring
#'                              until the `stopRule` is `FALSE`. This will add back
#'                              the last removed cell and include it in the return result
#'                              for that event.\cr
#'   `"excludePixel"` \tab Like `"includePixel"`, but it will not add back the cell
#'                        that causes `stopRule` to be `TRUE`\cr
#'   `"excludeRing"` \tab Analogous to `"excludePixel"`, but for the entire final
#'                             ring of cells added. This will exclude the entire ring of cells
#'                             that caused the `stopRule` to be `TRUE`\cr
#' }
#'
#' @seealso [spread2()] for a different implementation of the same algorithm.
#' It is more robust, meaning, there will be fewer unexplainable errors, and the behaviour
#' has been better tested, so it is more likely to be exactly as described under all
#' argument combinations.
#' Also, [rings()] which uses `spread` but with specific argument
#' values selected for a specific purpose.
#' [terra::distance()].
#' [cir()] to create "circles"; it is fast for many small problems.
#'
#' @param landscape     A `RasterLayer` or `SpatRaster` object. This defines the possible
#'                      locations for spreading events to start and spread into.
#'                      This can also be used as part of `stopRule`.
#'
#' @param loci          A vector of locations in `landscape`.
#'                      These should be cell indices.
#'                      If user has x and y coordinates, these can be converted
#'                      with [`cellFromXY()`][terra::cellFromXY].
#'
#' @param spreadProb    Numeric, `RasterLayer`, or `SpatRaster`.
#'                      If numeric of length 1, then this is the global probability
#'                      of spreading into each cell from a neighbour.
#'                      If a raster (or a vector of length `terra::ncell(landscape)`,
#'                      resolution and extent of `landscape`), then this will
#'                      be the cell-specific probability. Default is `0.23`.
#'                      If a `spreadProbLater` is provided, then this is
#'                      only used for the first iteration. Also called "escape
#'                      probability". See section on "Breaking out of spread events".
#'
#' @param persistence   A length 1 probability that an active cell will continue
#'                      to burn, per time step.
#'
#' @param mask          `RasterLayer` or `SpatRaster` object congruent with `landscape`,
#'                      whose elements are `0,1`, where `1` indicates "cannot spread to".
#'                      Currently not implemented, but identical behaviour can be
#'                      achieved if `spreadProb` has zeros in all unspreadable
#'                      locations.
#'
#' @param maxSize       Numeric. Maximum number of cells for a single or all events to be spread.
#'                      Recycled to match `loci` length, if it is not as long as `loci`.
#'                      See section on `Breaking out of spread events`.
#'
#' @param directions    The number of adjacent cells in which to look;
#'                      default is 8 (Queen case). Can only be 4 or 8.
#'
#' @param iterations    Number of iterations to spread.
#'                      Leaving this `NULL` allows the spread to continue
#'                      until stops spreading itself (i.e., exhausts itself).
#'
#' @param lowMemory     Deprecated.
#'
#' @param returnIndices Logical or numeric. If `1` or `TRUE`, will return a `data.table`
#'                      with indices and values of successful spread events.
#'                      If `2`, it will simply return a vector of pixel indices of
#'                      all cells that were touched. This will be the fastest option.
#'                      If `FALSE`, then it will return a raster with values.
#'                      See Details.
#'
#' @param returnDistances Logical. Should the function include a column with the
#'                      individual cell distances from the locus where that event
#'                      started. Default is `FALSE`.
#'                      See Details.
#'
#' @param spreadProbLater Numeric, or `RasterLayer`. If provided, then this
#'                      will become the `spreadProb` after the first iteration.
#'                      See Details.
#'
#' @param spreadState   `data.table`. This should be the output of a previous call
#'                      to `spread`, where `returnIndices` was `TRUE`.
#'                      Default `NA`, meaning the spread is starting from `loci`.
#'                      See Details.
#'
#' @param circle        Logical. If `TRUE`, then outward spread will be by
#'                      equidistant rings, rather than solely by adjacent cells
#'                      (via `directions` arg.). Default is `FALSE`.
#'                      Using `circle = TRUE` can be dramatically slower for large problems.
#'                      Note, this should usually be used with `spreadProb = 1`.
#'
#' @param circleMaxRadius Numeric. A further way to stop the outward spread of events.
#'                      If `circle` is `TRUE`, then it will grow to this maximum radius.
#'                      See section on `Breaking out of spread events`.
#'                      Default is `NA`.
#'
#' @param stopRule      A function which will be used to assess whether each
#'                      individual cluster should stop growing.
#'                      This function can be an argument of `"landscape"`,
#'                      `"id"`, `"cells"`, and any other variables passed to
#'                      `spread` in the `...`. `cells` and `landscape` will both
#'                      be numeric vectors of length of active cells. `cells` will
#'                      be the raster index, so can be used to extract values
#'                      from another raster passed via ... .
#'                      Default `NA`, meaning that spreading will not stop
#'                      as a function of the landscape.
#'                      See section on "Breaking out of spread events" and examples.
#'
#' @param stopRuleBehavior Character. Can be one of `"includePixel"`,
#'                      `"excludePixel"`, `"includeRing"`, or
#'                      `"excludeRing"`.
#'                      If `stopRule` contains a function, this argument is
#'                      used determine what to do with the cell(s) that caused
#'                      the rule to be `TRUE`. See details.
#'                      Default is `"includeRing"` which means to accept the
#'                      entire ring of cells that caused the rule to be `TRUE`.
#'
#' @param allowOverlap  Logical. If `TRUE`, then individual events can overlap
#'                      with one another, i.e., they do not interact (this is slower
#'                      than if `allowOverlap = FALSE`).
#'                      Default is `FALSE`.
#'
#' @param asymmetry     A numeric indicating the ratio of the asymmetry to be used.
#'                      Default is `NA`, indicating no asymmetry.
#'                      See details. This is still experimental.
#'                      **Use with caution.**
#'
#' @param asymmetryAngle A numeric indicating the angle in degrees (0 is "up",
#'                       as in North on a map), that describes which way the `asymmetry` is.
#'
#' @param quick  Logical. If `TRUE`, then several potentially time consuming
#'               checking (such as `inRange`) will be skipped.
#'               This should only be used if there is no concern about checking
#'               to ensure that inputs are legal.
#'
#' @param neighProbs A numeric vector, whose sum is 1.
#'                   It indicates the probabilities an individual spread iteration
#'                   spreading to `1:length(neighProbs)` neighbours.
#'
#' @param exactSizes Logical. If `TRUE`, then the `maxSize` will be
#'                   treated as exact sizes, i.e., the spread events will continue
#'                   until they are `floor(maxSize)`.
#'                   This is overridden by `iterations`, but if `iterations`
#'                   is run, and individual events haven't reached `maxSize`,
#'                   then the returned `data.table` will still have at least
#'                   one active cell per event that did not achieve `maxSize`,
#'                   so that the events can continue if passed into `spread`
#'                   with `spreadState`.
#'
#' @param relativeSpreadProb Logical. If `TRUE`, then `spreadProb` will
#'                      be rescaled *within* the `directions` neighbours, such that
#'                      the sum of the probabilities of all neighbours will be 1. Default
#'                      `FALSE`, unless `spreadProb` values are not contained
#'                      between 0 and 1, which will force `relativeSpreadProb`
#'                      to be `TRUE`.
#'
#' @param ...           Additional named vectors or named list of named vectors
#'                      required for `stopRule`. These
#'                      vectors should be as long as required e.g., length
#'                      `loci` if there is one value per event.
#' @param plot.it  If `TRUE`, then plot the raster at every iteration,
#'                 so one can watch the spread event grow.
#'
#' @param mapID    Deprecated. Use `id`.
#'
#' @param id    Logical. If `TRUE`, returns a raster of events ids.
#'              If `FALSE`, returns a raster of iteration numbers,
#'              i.e., the spread history of one or more events.
#'              NOTE: this is overridden if `returnIndices` is `TRUE`
#'              or `1` or `2`.
#'
#' @return Either a `RasterLayer` indicating the spread of the process in
#' the landscape or a `data.table` if `returnIndices` is `TRUE`.
#' If a `RasterLayer`, then it represents
#' every cell in which a successful spread event occurred. For the case of, say, a fire
#' this would represent every cell that burned. If `allowOverlap` is `TRUE`,
#' This `RasterLayer` will represent the sum of the individual event ids
#' (which are numerics `seq_along(loci)`.
#' This will generally be of minimal use because it won't be possible to distinguish
#' if event 2 overlapped with event 5 or if it was just event 7.
#'
#' If `returnIndices` is `TRUE`,
#' then this function returns a `data.table` with columns:
#'
#' \tabular{ll}{
#'   `id` \tab an arbitrary ID `1:length(loci)` identifying
#'                      unique clusters of spread events, i.e., all cells
#'                      that have been spread into that have a
#'                      common initial cell.\cr
#'   `initialLocus` \tab the initial cell number of that particular
#'                            spread event.\cr
#'   `indices` \tab The cell indices of cells that have
#'                        been touched by the spread algorithm.\cr
#'   `active` \tab a logical indicating whether the cell is active (i.e.,
#'                        could still be a source for spreading) or not (no
#'                        spreading will occur from these cells).\cr
#' }
#'
#' This will generally be more useful when `allowOverlap` is `TRUE`.
#'
#' @example inst/examples/example_spread.R
#'
#' @author Eliot McIntire and Steve Cumming
#' @export
#' @importFrom data.table := data.table setcolorder set
#' @importFrom fpCompare %<=%
#' @importFrom reproducible maxFn minFn
#' @importFrom stats runif
#' @importFrom terra ext ncell rast res setValues
#' @importFrom utils assignInMyNamespace
#' @rdname spread
spread <- function(landscape, loci = NA_real_, spreadProb = 0.23, persistence = 0,
                   mask = NA, maxSize = 1e8L, directions = 8L, iterations = 1e6L,
                   lowMemory = NULL, # getOption("spades.lowMemory"),
                   returnIndices = FALSE,
                   returnDistances = FALSE, mapID = NULL, id = FALSE, plot.it = FALSE,
                   spreadProbLater = NA_real_, spreadState = NA,
                   circle = FALSE, circleMaxRadius = NA_real_,
                   stopRule = NA, stopRuleBehavior = "includeRing", allowOverlap = FALSE,
                   asymmetry = NA_real_, asymmetryAngle = NA_real_, quick = FALSE,
                   neighProbs = NULL, exactSizes = FALSE, relativeSpreadProb = FALSE, ...) {

    if (!is.null(neighProbs)) {
      if (isTRUE(allowOverlap))
        stop("Can't use neighProbs and allowOverlap = TRUE together")
    }
    if (requireNamespace("dqrng", quietly = TRUE)) {
      samInt <- dqrng::dqsample.int
      # set dqrng seed from base state
      dqrng::dqset.seed(sample.int(1e9, 2))
    } else {
      samInt <- sample.int
    }

    if (!is.null(mapID)) {
      warning("mapID is deprecated, use id")
      id <- mapID
    }
    if (!quick) {
      allowedRules <- c("includePixel", "excludePixel", "includeRing", "excludeRing")
      if (!any(stopRuleBehavior %in% allowedRules))
        stop("stopRuleBehaviour must be one of '", paste(allowedRules, collapse = "', '"), "'.")
    }
    if (isTRUE(lowMemory)) {
      stop("lowMemory is no longer supported due to removal of ffbase from CRAN.")
      # .requireNamespace("ff")
      # .requireNamespace("ffbase")
    }

    spreadStateExists <- is(spreadState, "data.table")
    spreadProbLaterExists <- TRUE

    if (!is(spreadProbLater, "Raster")) {
      if (anyNA(spreadProbLater)) {
        spreadProbLaterExists <- FALSE
        spreadProbLater <- spreadProb
      }
    }

    ### should sanity check map extents
    if (any(is.na(loci)))  {
      # start it in the centre cell, if there is no spreadState
      if (!spreadStateExists)
        loci <- middlePixel(landscape)
    }
    if (!is.integer(loci)) {
      loci <- as.integer(loci)
    }
    if (!quick) {
      dupLoci <- duplicated(loci)
      if (any(duplicated(loci))) {
        message("duplicate initial loci are provided")
        # loci <- loci[dupLoci]
      }
    }

    if (length(loci) == 0) stop("No loci. Nothing to do")

    if (any(!is.na(maxSize))) {
      msEqZero <- maxSize < 1
      if (any(msEqZero)) {
        loci <- loci[!msEqZero]
        maxSize <- maxSize[!msEqZero]
      }
    }

    if (spreadStateExists) {
      keepers <- spreadState$active == TRUE
      loci <- initialActiveCells <- spreadState[keepers, indices]
      initialLoci <- unique(spreadState$initialLocus)
    } else {
      initialLoci <- loci
    }
    lenInitialLoci <- length(initialLoci)
    sequenceInitialLoci <- seq(lenInitialLoci)

    # Check for probabilities
    if (!quick) {
      if (inherits(spreadProbLater, "RasterLayer") || inherits(spreadProb, "Rasterlayer") ||
          inherits(spreadProbLater, "SpatRaster") || inherits(spreadProb, "SpatRaster")) {
        if ((minFn(spreadProb) > 1L) || (maxFn(spreadProb) < 0L) ||
            (maxFn(spreadProb) > 1L) || (minFn(spreadProb) < 0L)) {
          relativeSpreadProb <- TRUE
        }
        if (spreadProbLaterExists)
          if ((minFn(spreadProbLater) > 1L) || (maxFn(spreadProbLater) < 0L) ||
              (maxFn(spreadProbLater) > 1L) || (minFn(spreadProbLater) < 0L)) {
            relativeSpreadProb <- TRUE
          }
      } else {
        if (!all(inRange(na.omit(spreadProb)))) {
          relativeSpreadProb <- TRUE
          stop("spreadProb is not a probability")
        }
        if (spreadProbLaterExists) {
          relativeSpreadProb <- TRUE
          if (!all(inRange(na.omit(spreadProbLater)))) stop("spreadProbLater is not a probability")
        }
      }
    }

    ncells <- as.integer(terra::ncell(landscape))

    #browser(expr = exists("aaaaa"))
    allowOverlapOrReturnDistances <- allowOverlap | returnDistances
    useMatrixVersionSpreads <- allowOverlapOrReturnDistances | spreadStateExists
    if (useMatrixVersionSpreads) {
      if (spreadStateExists) {
        spreads <- as.matrix(spreadState[, list(initialLocus, indices, id, active)])
      } else {
        spreads <- cbind(initialLocus = initialLoci, indices = initialLoci,
                         id = seq_along(loci), active = 1)
      }
    } else {
      if (!is.null(lowMemory)) {
        message("lowMemory argument is now deprecated; using standard spread.")
      }
      # The experimental new spread function has some changes for speed. 1) The
      # bottleneck amazingly, was the creation of a new empty vector of length
      # ncell(landscape) ... it took >50% of the time of the spread function
      # when called 100,000s of times on a variety of spreadProb situations. 2) I
      # found that the only way to stop instantiating this was to have a
      # data.table object that uses reference semantics. 3) Put a simple, 1 column
      # data.table object into the SpaDES.tools namespace. It will contain the
      # former spreads object which was 0 everywhere the events hadn't spread
      # to, and a non-zero integer otherwise. 4) The function has to make sure that
      # it is "correct" on leaving the function. Two different cases: A) it
      # exits improperly --> action is delete this object; B) it exits correctly
      # --> action is to change all the values that were non-zero back to zero,
      # rather than delete the object. The whole point is to keep the object
      # intact after it has exited spread, so that it is available again
      # immediately for reuse.
      needEmptySpreads <- TRUE
      stNamespace <- asNamespace("SpaDES.tools")
      if (exists("spreadsDTInNamespace", envir = stNamespace)) {
        spreadsDT <- getFromNamespace("spreadsDTInNamespace", "SpaDES.tools")
        # set(spreadsDT, NULL, "spreads", 0L)
        # spreads <- spreadsDT$spreads
        if (identical(NROW(spreadsDT), ncells)) {
          needEmptySpreads <- FALSE
        }
      }
      if (needEmptySpreads) {
        spreadsDT <- data.table(spreads = vector("integer", ncells))
        set(spreadsDT, NULL, "spreads", 0L)
        # put the empty data.table into the SpaDES.tools namespace
        assignInMyNamespace("spreadsDTInNamespace", spreadsDT)
        on.exit(assignInMyNamespace("spreadsDTInNamespace", integer()), add = TRUE)
      }
    }

    n <- 1L

    # circle needs directions to be 8
    if (circle | !is.na(asymmetry)) {
      if (circle) directions <- 8L # only required for circle
      initialLociXY <- cbind(id = seq_along(initialLoci), xyFromCell(landscape, initialLoci))
      id <- TRUE
      if (allowOverlapOrReturnDistances) {
        spreads <- cbind(spreads, dists = 0)
      }
    }

    # determine ... variables
    otherVars <- list(...)
    anyList <- unlist(lapply(otherVars, is.list))

    if (any(anyList)) {
      otherVarsLists <- unlist(unname(otherVars), recursive = FALSE)
      otherVars[anyList] <- NULL
      otherVars <- append(otherVars, otherVarsLists)
    }

    # check validity of stopRule
    if (is.function(stopRule)) {
      id <- TRUE
      stopRuleObjs <- names(formals(stopRule))
      if (!quick) {
        if (any(is.na(match(stopRuleObjs,
                            c("id", "landscape", "cells", names(otherVars)))))) {
          stop("Arguments in stopRule not valid.\n",
               "The function definition must be a function of built-in options,",
               " (id, landscape, or cells) or user supplied variables.",
               " If user supplied, the variables",
               " must be passed as named vectors, or lists or data.frames.",
               " See examples.")
        }
      }
      landRasNeeded <- any(stopRuleObjs == "landscape")
      colNamesPotentials <- c("id", "landscape"[landRasNeeded], "cells", "prev")
      argNames <- c(colNamesPotentials, names(otherVars))
      whArgs <- match(names(formals(stopRule)), argNames)

      # Raster indexing is slow. If there is are Rasters submitted with the stopRule
      #  then this will convert them to vectors first. Clearly, this will have
      #  memory consequences if the Rasters are on disk, but spread is optimized for speed

      # rasters <- unlist(lapply(otherVars[names(otherVars)], function(x) is(x, "Raster")))
      # if (any(rasters)) {
      #   for (i in 1:which(rasters)) {
      #     otherVars[[names(rasters[i])]] <- otherVars[[names(rasters[i])]][]
      #   }
      # }
      landRas <- landscape[] # For speed
    }

    if (!allowOverlap && !returnDistances) {
      if (id | returnIndices > 0 | relativeSpreadProb) {
        if (!spreadStateExists) {
          set(spreadsDT, loci, "spreads", seq_along(loci))
          ##DT spreads[loci] <- seq_along(loci)
          # give values to spreads vector at initialLoci
        }
      } else {
        spreadsDT$spreads[loci] <- n
      }
      spreadsIndices <- unname(loci)
      #browser(expr = exists("aaaaa"))
      length(spreadsIndices) <- length(loci) * 100
      prevSpreadIndicesActiveLen <- length(loci)
      prevSpreadIndicesFullLen <- length(spreadsIndices)
    }

    # Convert mask and NAs to 0 on the spreadProb Raster
    if (is(spreadProb, "Raster")) {
      # convert NA to 0s
      #isNASpreadProb <- is.na(spreadProb[])
      # if (anyNA(spreadProb[])) {
      #   isNASpreadProb <- is.na(spreadProb[])
      #   spreadProb[isNASpreadProb] <- 0L
      # }
    } else if (is.numeric(spreadProb)) {
      # Translate numeric spreadProb into a Raster, if there is a mask
      if (is(mask, "Raster")) {
        spreadProb <- terra::rast(terra::ext(landscape), res = res(landscape), vals = spreadProb)
      }
    }

    # Convert mask and NAs to 0 on the spreadProbLater Raster
    if (is(spreadProbLater, "Raster")) {
    } else if (is.numeric(spreadProbLater)) {
      # Translate numeric spreadProbLater into a Raster, if there is a mask
      if (is(mask, "Raster")) {
        spreadProbLater <- terra::rast(terra::ext(landscape), res = res(landscape),
                                       vals = spreadProbLater)
      }
    }

    # Mask spreadProbLater and spreadProb
    if (is(mask, "Raster")) {
      spreadProbLater[mask[] == 1L] <- 0L
      spreadProb[mask[] == 1L] <- 0L
    }

    if (spreadStateExists) {
      if (allowOverlapOrReturnDistances) {
        stop("Using spreadState with either allowOverlap = TRUE",
             " or returnDistances = TRUE is not implemented")
      } else {
        if (sum(colnames(spreadState) %in% c("indices", "id", "active", "initialLocus")) != 4) {
          stop("spreadState must have at least columns: ",
               "indices, id, active, and initialLocus.")
        }
      }
    }

    if (!quick) {
      if (any(loci > ncells)) stop("loci indices are not on landscape")
    }

    # Recycling maxSize as needed
    if (any(!is.na(maxSize))) {
      if (!is.integer(maxSize)) maxSize <- floor(maxSize)
      if (spreadStateExists) {
        sizeAll <- spreadState[, list(len = .N), by = id]
        size <- c(sizeAll[, len])
      } else {
        maxSize <- rep_len(maxSize, length(loci))
        size <- rep_len(1L, length(loci))
      }
    } else {
      maxSize <- ncells
      size <- length(loci)
    }

    #browser(expr = exists("aaaaa"))
    noMaxSize <- all(maxSize >= ncells) # will be used to omit testing for maxSize
    if (is.null(neighProbs)) {
      numNeighs <- NULL
    }

    if (!exists("numRetries", envir = .pkgEnv)) {
      assign("numRetries", rep(0, lenInitialLoci), envir = .pkgEnv)
    }

    toColumn <- c("to", "indices")

    #browser(expr = exists("aaaaa"))
    # while there are active cells
    while (length(loci) & (n <= iterations)) {
      if (!is.null(neighProbs)) {
        numNeighs <- if (is.list(neighProbs)) {
          unlist(lapply(neighProbs, function(x) {
            sample.int(length(x), size = 1, replace = TRUE, prob = x)
          }))
        } else {
          sample.int(length(neighProbs), size = length(loci), replace = TRUE, prob = neighProbs)
        }
      }

      # identify neighbours
      if (useMatrixVersionSpreads) {
        whActive <- spreads[, "active"] == 1 # spreads carries over
        potentials <- adj(landscape, loci, directions, pairs = TRUE,
                          id = spreads[whActive, "id"])#, numNeighs = numNeighs)
        spreads[whActive, "active"] <- 0
        potentials <- cbind(potentials, active = 1)
      } else {
        if (id | returnIndices > 0 | circle | relativeSpreadProb | !is.null(neighProbs)) {
          potentials <- adj(landscape, loci, directions, pairs = TRUE)
        } else {
          # must pad the first column of potentials
          newAdj <- adj(landscape, loci, directions, pairs = FALSE)
          potentials <- cbind(NA_integer_, newAdj)
        }
      }

      if (circle) {
        potentials <- cbind(potentials, dists = 0)
      }

      # keep only neighbours that have not been spread to yet
      if (useMatrixVersionSpreads) {
        # data.table version is faster for potentials > 2000 or so
        if (NROW(potentials) > 2000) {
          spreadsDT <- as.data.table(spreads)
          potentialsDT <- as.data.table(potentials)
          potentialsDT[, initialLocus := initialLoci[potentialsDT$id]]
          colnamesPot <- colnames(potentialsDT)
          whIL <- which(colnamesPot == "initialLocus")
          whFrom <- which(colnamesPot == "from")
          setcolorder(potentialsDT,
                      c(colnamesPot[whIL], colnamesPot[-c(whIL, whFrom)], colnamesPot[whFrom]))
          setnames(potentialsDT, old = "to", new = "indices")
          newPot <- potentialsDT[!spreadsDT, on = c("id", "indices")]
          potentials <- as.matrix(newPot)
        } else {
          potentials <- cbind(initialLocus = initialLoci[potentials[, "id"]], potentials)
          colnames(potentials)[which(colnames(potentials) == "to")] <- "indices"
          colnamesPot <- colnames(potentials)
          whIL <- which(colnamesPot == "initialLocus")
          whFrom <- which(colnamesPot == "from")
          potentials <- potentials[,c(colnamesPot[whIL],
                                      colnamesPot[-c(whIL, whFrom)],
                                      colnamesPot[whFrom])]

          # These next lines including the lapply are the rate limiting
          #   step and it has been heavily worked to speed it up March 31, 2020
          seq2 <- sequenceInitialLoci[sequenceInitialLoci %in% potentials[,"id"]]
          out <- lapply(seq2, function(ind) {
            hasID <- potentials[, "id"] == ind
            po <- potentials[hasID, ]
            hasID2 <- spreads[, "id"] == ind
            inds <- spreads[hasID2, "indices"]
            vals <- po[, 2L] %in% inds
            po[!vals,]
          })
          potentials <- do.call(rbind, out)
        }
      } else {
        # Keep only the ones where it hasn't been spread to yet
        keep <- spreadsDT$spreads[potentials[, 2L]] == 0L
        # keep <- spreads[potentials[, 2L]] == 0L
        potentials <- potentials[keep, , drop = FALSE]
      }

      if (n == 2) {
        spreadProb <- spreadProbLater
      }

      # extract spreadProb values from spreadProb argument
      if (is.numeric(spreadProb)) {
        if (!(length(spreadProb) == 1 || length(spreadProb) == terra::ncell(landscape)))
          stop("spreadProb must be length 1 or length terra::ncell(landscape), or a raster")
        if (n == 1 & spreadProbLaterExists) {
          # need cell specific values
          spreadProbs <- rep(spreadProb, NROW(potentials))
          spreadProb <- spreadProbLater
        } else {
          if (length(spreadProb) > 1) {
            spreadProbs <- spreadProb[potentials[, 2L]]
          } else {
            spreadProbs <- rep(spreadProb, NROW(potentials))
          }
        }
      } else {
        # here for raster spreadProb
        if (n == 1 & spreadProbLaterExists) {
          # need cell specific values
          spreadProbs <- spreadProb[][potentials[, 2L]]
          spreadProb <- spreadProbLater
        } else {
          spreadProbs <- spreadProb[][potentials[, 2L]]
        }
      }

      if (anyNA(spreadProbs)) spreadProbs[is.na(spreadProbs)] <- 0

      if (!is.na(asymmetry)) {
        if (allowOverlapOrReturnDistances) {
          a <- cbind(id = potentials[, 3L], to = potentials[, 2L],
                     xyFromCell(landscape, potentials[, 2L]))
        } else {
          if (useMatrixVersionSpreads) {
            a <- cbind(id = spreads[potentials[, 1L]], to = potentials[, 2L],
                       xyFromCell(landscape, potentials[, 2L]))
          } else {
            a <- cbind(id = spreadsDT$spreads[potentials[, 1L]], to = potentials[, 2L],
                       xyFromCell(landscape, potentials[, 2L]))
          }
        }
        d <- directionFromEachPoint(from = initialLociXY, to = a)
        newSpreadProbExtremes <- (spreadProb[] * 2) / (asymmetry + 1) * c(1, asymmetry)
        angleQuality <- (cos(d[, "angles"] - rad2(asymmetryAngle)) + 1) / 2
        spreadProbs <- newSpreadProbExtremes[1] + (angleQuality * diff(newSpreadProbExtremes))
        spreadProbs <- spreadProbs - diff(c(spreadProb[], mean(spreadProbs)))
      }

      if (!is.null(neighProbs) | relativeSpreadProb) {
        aaa <- split(seq_along(potentials[, toColumn[spreadStateExists + 1]]),
                     potentials[, "from"])
        if (length(aaa) != length(numNeighs)) {
          activeCellContinue <- loci %in% unique(potentials[, "from"])
          numNeighs <- numNeighs[activeCellContinue]
        }

        tmpA <- unlist(lapply(aaa, length))
        tmpB <- which(tmpA < numNeighs)
        if (length(tmpB) > 0)
          numNeighs[tmpB] <- unname(tmpA[tmpB])

        if (relativeSpreadProb) {
          rescaledProbs <- tapply(spreadProbs, potentials[, "from"], function(x) {
            x / sum(x, na.rm = TRUE)
          }, simplify = FALSE)
          neighIndexToKeep <- unlist(lapply(seq_along(aaa), function(x) {
            resample(aaa[[x]], size = numNeighs[x], prob = rescaledProbs[[x]])
          }))
        } else {
          neighIndexToKeep <- unlist(lapply(seq_along(aaa), function(x)
            resample(aaa[[x]], size = numNeighs[x])))
        }
        potentials <- potentials[neighIndexToKeep, , drop = FALSE]
        spreadProbs <- spreadProbs[neighIndexToKeep]
        spreadProbs[spreadProbs > 0] <- 1
      }

      randomSuccesses <- runifC(NROW(potentials)) <= spreadProbs
      potentials <- potentials[randomSuccesses, , drop = FALSE]

      # random ordering so not always same:
      lenPot <- NROW(potentials)
      if (lenPot) {
        reorderVals <- samInt(lenPot) ## TODO: uses sample.int(..., replace = FALSE)
        potentials <- potentials[reorderVals, , drop = FALSE]
      }
      if (!allowOverlap) {
        # here is where allowOverlap and returnDistances are different ##### NOW OBSOLETE, I BELIEVE ELIOT March 2020
        potentials <- potentials[!duplicated(potentials[, 2L]), , drop = FALSE]
      } else {
        pots <- potentials[, c("id", "indices"), drop = FALSE]
        potentials <- potentials[!duplicated(pots), , drop = FALSE]
      }

      # increment iteration
      n <- n + 1L

      # potentials can become zero because all active cells are edge cells
      if (length(potentials) > 0) {
        # implement circle
        if (!missing(circle)) {
          if (circle) {
            if (allowOverlapOrReturnDistances) {
              a <- cbind(potentials, xyFromCell(landscape, potentials[, 2L]))
            } else {
              #browser(expr = exists("aaaaa"))
              if (useMatrixVersionSpreads) {
                a <- cbind(potentials, id = spreads[potentials[, "from"]],
                           xyFromCell(landscape, potentials[, "to"]))
              } else {
                a <- cbind(potentials, id = spreadsDT$spreads[potentials[, "from"]],
                           xyFromCell(landscape, potentials[, "to"]))
              }
            }
            # need to remove dists column because distanceFromEachPoint, adds one back
            a <- a[, !(colnames(a) %in% c("dists")), drop = FALSE]
            # need 3 columns, id, x, y in both initialLociXY and a
            d <- distanceFromEachPoint(initialLociXY, a, angles = asymmetry) # d is sorted
            cMR <- (n - 1) * res(landscape)[1]
            if (!any(is.na(circleMaxRadius))) {
              # don't bother proceeding if circleMaxRadius is larger than current iteration
              if (any(circleMaxRadius <= ((n - 1) * res(landscape)[1]))) {
                if (length(circleMaxRadius) > 1) {
                  # if it is a vector of values
                  cMR <- circleMaxRadius[d[, "id"]]
                } else {
                  cMR <- circleMaxRadius
                }
              }
            }
            potentials <- d[, !(colnames(d) %in% xycolNames), drop = FALSE]
            potentials <- potentials[(d[, "dists"] %<=% cMR), , drop = FALSE]
          }
        }

        # If potentials has distances in it, it will be a numeric matrix; events should be integer
        events <- if (!is.integer(potentials)) as.integer(potentials[, 2L]) else potentials[, 2L]

        if (!noMaxSize) {
          if (useMatrixVersionSpreads) {
            len <- tabulate(potentials[, 3L], length(maxSize))
          } else {
            # actually interested in potential[,2L], but they don't have values yet..
            #  can use their source
            len <- tabulate(spreadsDT$spreads[potentials[, 1L]], length(maxSize))
          }
          if (any((size + len) > maxSize & size <= maxSize)) {
            whichID <- which(size + len > maxSize)

            # remove some active cells, if more than maxSize
            toRm <- (size + len)[whichID] - maxSize[whichID]

            for (i in seq_len(length(whichID))) {
              if (useMatrixVersionSpreads) {
                thisID <- which(potentials[, 3L] == whichID[i])
              } else {
                thisID <- which(spreadsDT$spreads[potentials[, 1L]] == whichID[i])
              }

              # some unusual cases where there are none on the spreads. Unsure how this occurs
              if (length(thisID))
                potentials <- potentials[-resample(thisID, toRm[i]), , drop = FALSE]
            }
            events <- as.integer(potentials[, 2L])
          }
          size <- pmin(size + len, maxSize) ## Quick? and dirty. fast but loose (too flexible)
        }

        # Implement stopRule section
        if (is.function(stopRule) & length(events) > 0) {
          if (allowOverlapOrReturnDistances) {
            prevCells <- cbind(
              id = spreads[, "id"],
              landscape = if (landRasNeeded) landRas[spreads[, "indices"]] else NULL,
              cells = spreads[, "indices"], prev = 1)
            eventCells <- cbind(id = potentials[, "id"],
                                landscape = if (landRasNeeded) landRas[events] else NULL,
                                cells = events, prev = 0)
          } else {
            whgtZero <- spreadsIndices
            #browser(expr = exists("aaaaa"))
            if (useMatrixVersionSpreads) {
              prevCells <- cbind(id = spreads[whgtZero],
                                 landscape = if (landRasNeeded) landRas[whgtZero] else NULL,
                                 cells = whgtZero, prev = 1)
              eventCells <- cbind(id = spreads[potentials[, 1L]],
                                  landscape = if (landRasNeeded) landRas[potentials[, 2L]] else NULL,
                                  cells = potentials[, 2L], prev = 0)
            } else {
              prevCells <- cbind(id = spreadsDT$spreads[whgtZero],
                                 landscape = if (landRasNeeded) landRas[whgtZero] else NULL,
                                 cells = whgtZero, prev = 1)
              eventCells <- cbind(id = spreadsDT$spreads[potentials[, 1L]],
                                  landscape = if (landRasNeeded) landRas[potentials[, 2L]] else NULL,
                                  cells = potentials[, 2L], prev = 0)
            }
          }
          if (circle) {
            prevCells <- cbind(prevCells, dist = NA)
            eventCells <- cbind(eventCells, dist = potentials[, "dists"])
          }
          # don't need to continue doing ids that are not active
          tmp <- rbind(prevCells[prevCells[, "id"] %in% unique(eventCells[, "id"]), ], eventCells)

          ids <- unique(tmp[, "id"])

          shouldStopList <- lapply(ids, function(id) {
            shortTmp <- tmp[tmp[, "id"] == id, ]
            args <- append(list(id = id),
                           lapply(colNamesPotentials[-1], function(j) shortTmp[, j]))
            names(args) <- colNamesPotentials
            args <- append(args, otherVars)
            do.call(stopRule, args[whArgs])
          })

          if (any(lapply(shouldStopList, length) > 1))
            stop("stopRule does not return a length-one logical.",
                 " Perhaps stopRule need indexing by cells or id?")

          shouldStop <- unlist(shouldStopList)

          names(shouldStop) <- ids

          if (any(shouldStop)) {
            if (stopRuleBehavior != "includeRing") {
              if (stopRuleBehavior != "excludeRing") {
                whStop <- as.numeric(names(shouldStop)[shouldStop])
                whStopAll <- tmp[, "id"] %in% whStop
                tmp2 <- tmp[whStopAll, ]

                whStopEvents <- eventCells[, "id"] %in% whStop

                # If an event needs to stop, then must identify which cells are included
                out <- lapply(whStop, function(id) {
                  tmp3 <- tmp2[tmp2[, "id"] == id, ]
                  newOnes <- tmp3[, "prev"] == 0
                  ord <- seq_along(newOnes)

                  # because of undesired behaviour of sample when length(x) == 1
                  if (sum(newOnes) > 1) {
                    ord[newOnes] <- sample(ord[newOnes])
                    if (circle) ord[newOnes] <- ord[newOnes][order(tmp3[ord[newOnes], "dist"])]
                    tmp3 <- tmp3[ord, ]
                  }
                  startLen <- sum(!newOnes)
                  addIncr <- 1
                  done <- FALSE
                  args <- append(list(id = id),
                                 lapply(colNamesPotentials[-1], function(j) {
                                   tmp3[1:startLen, j]
                                 })) # instead of as.data.frame
                  names(args) <- colNamesPotentials
                  args <- append(args, otherVars)
                  argsSeq <- seq_along(colNamesPotentials[-1]) + 1

                  while (!done) {
                    args[argsSeq] <- lapply(colNamesPotentials[-1], function(j) {
                      unname(c(args[[j]], tmp3[(startLen + addIncr), j]))
                    }) # instead of as.data.frame
                    done <- do.call(stopRule, args[whArgs])
                    addIncr <- addIncr + 1
                  }
                  if (stopRuleBehavior == "excludePixel") addIncr <- addIncr - 1
                  firstInd <- startLen + addIncr
                  lastInd <- NROW(tmp3)
                  sequ <- if (firstInd > lastInd) 0 else firstInd:lastInd
                  tmp3[sequ, , drop = FALSE]
                })
                eventRm <- do.call(rbind, out)[, "cells"]
                cellsKeep <- !(potentials[, 2L] %in% eventRm)
              } else {
                cellsKeep <- rep(FALSE, NROW(potentials))
              }
              potentials <- potentials[cellsKeep, , drop = FALSE]
              events <- as.integer(potentials[, 2L])
              eventCells <- eventCells[cellsKeep, , drop = FALSE]
            }
            toKeepSR <- !(eventCells[, "id"] %in% as.numeric(names(which((shouldStop)))))
          }
        }

        if (length(events) > 0) {
          # place new value at new cells that became active
          if (useMatrixVersionSpreads) {
            fromCol <- colnames(potentials) == "from"

            spreads <- rbind(spreads, potentials[, !fromCol])
            if ((returnDistances | spreadStateExists) & !allowOverlap) {
              # 2nd place where allowOverlap and returnDistances differ
              notDups <- !duplicated(spreads[, "indices"])
              nrSpreads <- NROW(spreads)
              nrPotentials <- NROW(potentials)
              notDupsEvents <- notDups[-(1:(nrSpreads - nrPotentials))]
              spreads <- spreads[notDups, , drop = FALSE]
              events <- events[notDupsEvents]
            }
          } else {
            if (id | returnIndices > 0 | relativeSpreadProb) {
              # give new cells, the id of the source cell
              set(spreadsDT, events, "spreads", spreadsDT$spreads[potentials[, 1L]])
            } else {
              set(spreadsDT, events, "spreads", n)
            }
            curEventsLen <- length(events)
            addedIndices <- prevSpreadIndicesActiveLen + 1:curEventsLen

            if (sum(curEventsLen, prevSpreadIndicesActiveLen) > prevSpreadIndicesFullLen) {
              length(spreadsIndices) <- (prevSpreadIndicesActiveLen + curEventsLen) * 2
              prevSpreadIndicesFullLen <- length(spreadsIndices)
            }
            spreadsIndices[addedIndices] <- events
            prevSpreadIndicesActiveLen <- prevSpreadIndicesActiveLen + curEventsLen
            # spreadsIndices <- c(spreadsIndices, events)
          }
        }

        # remove the cells from "events" that push it over maxSize
        #  There are some edge cases that aren't captured above ... identify them...
        if (length(maxSize) > 1L) {
          if (exists("whichID", inherits = FALSE)) {
            # must update toKeepSR in case that is a second reason to stop event
            if (exists("toKeepSR", inherits = FALSE)) {
              if (allowOverlapOrReturnDistances) {
                maxSizeKeep <- !(spreads[spreads[, "active"] == 1, "id"] %in% whichID)
                spreads <- spreads[c(rep(TRUE, sum(spreads[, "active"] == 0)), maxSizeKeep), ]
              } else {
                if (useMatrixVersionSpreads) {
                  maxSizeKeep <- !spreads[events] %in% whichID
                } else {
                  maxSizeKeep <- !spreadsDT$spreads[events] %in% whichID
                }
              }
              events <- events[maxSizeKeep]
              toKeepSR <- toKeepSR[maxSizeKeep]
            }
            rm(whichID)
          }
        } else {
          if (all(size >= maxSize)) {
            potentials <- potentials[0L, ] ## remove any potential cells, as size is met
            events <- NULL
          }
        }

        # Remove cells that were stopped by stopRule
        if (is.function(stopRule)) {
          if (exists("toKeepSR", inherits = FALSE)) {
            events <- events[toKeepSR]
            if (allowOverlapOrReturnDistances) {
              spreads[c(rep(TRUE, sum(spreads[, "active"] == 0)), !toKeepSR), "active"] <- 0
            }
            rm(toKeepSR)
          }
        }
      } else {
        # there are no potentials -- possibly from failed runif, or spreadProbs all 0
        events <- NULL
      }

      if (exactSizes) {
        if (all(get("numRetries", inherits = FALSE, envir = .pkgEnv) < 10)) {
          if (spreadStateExists) {
            tooSmall <- tabulate(spreads[, "id"], length(maxSize)) < maxSize
            inactive <- tabulate(spreads[spreads[, "active"] == 1, "id"], length(maxSize)) == 0
          } else {
            if (useMatrixVersionSpreads) {
              tooSmall <- tabulate(spreads, length(maxSize)) < maxSize
              inactive <- tabulate(spreads[events], length(maxSize)) == 0
            } else {
              tooSmall <- tabulate(spreadsDT$spreads, length(maxSize)) < maxSize
              inactive <- tabulate(spreadsDT$spreads[events], length(maxSize)) == 0
            }
          }

          # these are ones that are stuck ... i.e., too small, and inactive
          needPersist <- tooSmall & inactive
          needPersistJump <- TRUE
          if (any(needPersist)) {
            assign("numRetries", envir = .pkgEnv,
                   get("numRetries", inherits = FALSE, envir = .pkgEnv) + needPersist)

            if (spreadStateExists) {
              whSmallInactive <- which(tooSmall & inactive)
              spreadsSmallInactive <- spreads[spreads[, "id"] %in% whSmallInactive, , drop = FALSE]
              if (needPersistJump) {
                message("Jumping to new active location, up to 1000 m away")
                mmm <- rings(landscape, loci = spreadsSmallInactive[, "indices"],
                             maxRadius = 1000, minRadius = 1, returnIndices = TRUE)
                wh <- mmm[, list(whKeepLoci = resample(.I, 1)), by = id]$whKeepLoci
              } else {
                for (whSI in whSmallInactive) {
                  wh <- which(spreads[, "id"] == whSI)
                  wh <- tail(wh, 2) # pick last two ones from all inactive cells
                  keepLoci <- spreads[wh, "indices"]
                  events <- c(keepLoci, events)
                  spreads[wh, "active"] <- 1
                }
              }
            } else {
              keepLoci <- spreads[loci] %in% which(tooSmall & inactive)
              events <- c(loci[keepLoci], events)
            }
          }
        }
      }

      # drop or keep loci
      if (is.na(persistence) | persistence == 0L) {
        loci <- NULL
      } else {
        if (inRange(persistence)) {
          loci <- loci[runif(length(loci)) <= persistence]
        } else {
          # here is were we would handle methods for raster* or functions
          stop("Unsupported type: persistence")
        }
      }

      if (plot.it) {
        # if (requireNamespace("quickPlot", quietly = TRUE)) {
        #   if (n == 2 & !spreadStateExists) quickPlot::clearPlot()
        #   if (allowOverlapOrReturnDistances) {
        #     spreadsDT <- data.table(spreads)
        #     hab2 <- landscape
        #     hab2[] <- 0
        #     pixVal <- spreadsDT[, sum(id), by = indices]
        #     hab2[pixVal$indices] <- pixVal$V1
        #     quickPlot::Plot(hab2, legendRange = c(0, sum(seq_along(initialLoci))))
        #   } else {
        #     plotCur <- terra::rast(landscape)
        #     plotCur <- setValues(plotCur, spreads)
        #     quickPlot::Plot(plotCur)
        #   }
        # } else {
             plotCur <- terra::rast(landscape)
             plotCur <- setValues(plotCur, spreadsDT$spreads)
             terra::plot(plotCur)
        # }
      }

      # new loci list for next while loop, concat of persistent and new events
      loci <- c(loci, events)
    } # end of while loop

    # Reset the base R seed so it is deterministic
    if (requireNamespace("dqrng", quietly = TRUE)) {
      set.seed(dqrng::dqsample.int(1e9, 1) + sample.int(1e9, 1))
    }

    if (!allowOverlap & !returnDistances) {
      spreadsIndices <- spreadsIndices[1:prevSpreadIndicesActiveLen]
    }

    # Convert the data back to raster
    if (!allowOverlap & !returnDistances & !spreadStateExists) {
      # if (lowMemory) {
      #   wh <- as.ram(ffwhich(spreads, spreads > 0))
      #   if (returnIndices > 0) {
      #     wh <- wh[!(wh %in% potentials[,2L])]
      #     completed <- data.table(indices = wh, id = spreads[wh], active = FALSE)
      #     if (NROW(potentials) > 0) {
      #       active <- data.table(indices = potentials[, 2L],
      #                            id = spreads[potentials[, 1L]],
      #                            active = TRUE)
      #     } else {
      #       active <- data.table(indices = numeric(0), id = numeric(0),
      #                            active = logical(0))
      #     }
      #   }
      # } else {
      wh <- if (spreadStateExists) {
        c(spreadState[!keepers]$indices, spreadsIndices)
      } else {
        spreadsIndices
      }
      if (returnIndices > 0) {
        # wh already contains the potentials for next iteration -- these should be not duplicated
        #   inside "completed"
        wh <- wh[!(wh %in% potentials[, 2L])]
        completed <- data.table(indices = wh, id = spreadsDT$spreads[wh], active = FALSE)
        if (NROW(potentials) > 0) {
          active <- data.table(indices = as.integer(potentials[, 2L]),
                               id = spreadsDT$spreads[potentials[, 1L]],
                               active = TRUE)
        } else {
          active <- data.table(indices = integer(0), id = integer(0), active = logical(0))
        }
      }
    }

    if (returnIndices == 1) {
      if (useMatrixVersionSpreads) {
        keepCols <- c(3, 1, 2, 4)
        if (circle) keepCols <- c(keepCols, 5)

        # change column order to match non allowOverlap
        allCells <- as.data.table(spreads[, keepCols, drop = FALSE])
        set(allCells, NULL, j = "active", as.logical(allCells$active))
        # setkeyv(allCells, "id")
      } else {
        #browser(expr = exists("aaaaa"))
        allCells <- rbindlist(list(completed, active)) # active first; next line will keep active
        if (spreadStateExists) {
          initEventID <- unique(spreadState$id)
        } else {
          initEventID <- allCells[indices %in% initialLoci, id]
        }
        if (!all(is.na(initialLoci))) {
          attr(initialLoci, ".match.hash") <- NULL # something in data.table put this
          dtToJoin <- data.table(id = sort(initEventID), initialLocus = initialLoci)
        } else {
          dtToJoin <- data.table(id = numeric(0), initialLocus = numeric(0))
        }
        setkeyv(dtToJoin, "id")
        setkeyv(allCells, "id")

        # tack on initialLoci
        allCells <- dtToJoin[allCells]
      }
      allCells[]
      if (exactSizes)
        if (exists("numRetries", envir = .pkgEnv)) {
          if (sum(allCells$active) == 0) rm("numRetries", envir = .pkgEnv)
        }
      if  (!(useMatrixVersionSpreads)) {
        #browser(expr = exists("aaaaa"))
        set(spreadsDT, allCells$indices, "spreads", 0L)
        # remove the previous on.exit which had the effect of deleting the contents
        #   completely on a failed `spread`. Here, we want to delete the previous
        #   on.exit --> allowing the object to stay intact, but with only zeros.
        on.exit()
        # spreadsIndices <- spreadsIndices[1:prevSpreadIndicesActiveLen]
      }

      return(allCells)
    }
    if (returnIndices == 2) {
      return(wh)
    }

    landscape[] <- 0

    ## remove colour table
    if (inherits(landscape, "RasterLayer")) {
      .requireNamespace("raster")
      raster::colortable(landscape) <- logical(0)
    } else if (inherits(landscape, "SpatRaster")) {
      terra::coltab(landscape) <- NULL
    }

    if (allowOverlapOrReturnDistances) {
      if (returnDistances & !allowOverlap) {
        landscape[spreads[, "indices"]] <- spreads[, "dists"]
      } else {
        spreadsDTFinal <- data.table(spreads)
        if (returnDistances & allowOverlap) {
          pixVal <- spreadsDTFinal[, min(dists), by = indices]
          message("returnDistances is TRUE, allowOverlap is TRUE, but returnIndices is FALSE; ",
                  "returning minimum distance raster.")
        } else {
          pixVal <- spreadsDTFinal[, sum(id), by = indices]
        }
        landscape[pixVal$indices] <- pixVal$V1
      }
    } else {
      landscape[wh] <- spreadsDT$spreads[wh]
      if (exists("potentials")) {
        if (NROW(potentials) > 0) {
          landscape[potentials[, 1L]] <- spreadsDT$spreads[potentials[, 2L]]
          set(spreadsDT, as.integer(potentials[, 1L]), "spreads", 0L)
        }
      }
      set(spreadsDT, wh, "spreads", 0L)

    }
    return(landscape)
}

spreadsDTInNamespace <- integer()

Try the SpaDES.tools package in your browser

Any scripts or data that you put into this service are public.

SpaDES.tools documentation built on Nov. 11, 2023, 1:07 a.m.