R/flows-disperse.R

In dodgr: Distances on Directed Graphs

#' @title Aggregate flows dispersed from each point in a network.
#'
#' @description Disperse flows throughout a network based on a input vectors of
#' origin points and associated densities. Dispersal is implemented as an
#' exponential decay, controlled by a parameter, `k`, so that flows decay with
#' `exp(-d / k)`, where `d` is distance. The algorithm allows for efficient
#' fitting of multiple dispersal models for different coefficients to be fitted
#' with a single call. Values of the dispersal coefficients, `k`, may take one
#' of the following forms:
#'
#' \itemize{
#' \item A single numeric value (> 0), with dispersal along all paths
#' calculated with that single value. Return object (see below) will then have
#' a single additional column named "flow".
#' \item A vector of length equal to the number of `from` points, with
#' dispersal from each point then calculated using the corresponding value of
#' `k`. Return object has single additional "flow" column.
#' \item A vector of any other length (that is, >  1 yet different to number of
#' `from` points), in which case different dispersal models will be fitted for
#' each of the `n` specified values, and the resultant return object will have
#' an additional 'n' columns, named 'flow1', 'flow2', ... up to 'n'. These
#' columns must be subsequently matched by the user back on to the
#' corresponding 'k' values.
#' \item A matrix with number of rows equal to the number of `from` points, and
#' any number of columns. Each column will then specify a distinct dispersal
#' model, with different values from each row applied to the corresponding
#' `from` points. The return value will then be the same as the previous
#' version, with an additional `n` columns, "flow1" to "flown".
#' }
#'
#'
#' Flows are calculated by default on contracted graphs, via the `contract =
#' TRUE` parameter. (These are derived by reducing the input graph down to
#' junction vertices only, by joining all intermediate edges between each
#' junction.) If changes to the input graph do not prompt changes to resultant
#' flows, and the default `contract = TRUE` is used, it may be that
#' calculations are using previously cached versions of the contracted graph.
#' If so, please use either \link{clear_dodgr_cache} to remove the cached
#' version, or \link{dodgr_cache_off} prior to initial graph construction to
#' switch the cache off completely.
#'
#' @inheritParams dodgr_flows_aggregate
#' @param graph `data.frame` or equivalent object representing the network
#' graph (see Details)
#' @param from Vector or matrix of points **from** which aggregate dispersed
#' flows are to be calculated (see Details)
#' @param dens Vectors of densities corresponding to the `from` points
#' @param k Width coefficient of exponential diffusion function defined as
#' `exp(-d/k)`, in units of distance column of `graph` (metres by default). Can
#' also be a vector with same length as `from`, giving dispersal coefficients
#' from each point. If value of `k<0` is given, a standard logistic polynomial
#' will be used.
#' @param tol Relative tolerance below which dispersal is considered to have
#' finished. This parameter can generally be ignored; if in doubt, its effect
#' can be removed by setting `tol = 0`.
#' @return Modified version of graph with additional `flow` column added.
#'
#' @family flows
#' @export
#' @examples
#' # This is generally needed to explore different values of `k` on same graph:
#' dodgr_cache_off ()
#'
#' graph <- weight_streetnet (hampi)
#' from <- sample (graph$from_id, size = 10)
#' dens <- rep (1, length (from)) # Uniform densities
#' graph <- dodgr_flows_disperse (graph, from = from, dens = dens)
#' # graph then has an additonal 'flows` column of aggregate flows along all
#' # edges. These flows are directed, and can be aggregated to equivalent
#' # undirected flows on an equivalent undirected graph with:
#' graph_undir <- merge_directed_graph (graph)
#'
#' # Remove `flow` column to avoid warning about over-writing values:
#' graph$flow <- NULL
#' # One dispersal coefficient for each origin point:
#' k <- runif (length (from))
#' graph <- dodgr_flows_disperse (graph, from = from, dens = dens, k = k)
#' grep ("^flow", names (graph), value = TRUE)
#' # single dispersal model; single "flow" column
#'
#' # Multiple models, muliple dispersal coefficients:
#' k <- 1:5
#' graph$flow <- NULL
#' graph <- dodgr_flows_disperse (graph, from = from, dens = dens, k = k)
#' grep ("^flow", names (graph), value = TRUE)
#' # Rm all flow columns:
#' graph [grep ("^flow", names (graph), value = TRUE)] <- NULL
#'
#' # Multiple models with unique coefficient at each origin point:
#' k <- matrix (runif (length (from) * 5), ncol = 5)
#' dim (k)
#' graph <- dodgr_flows_disperse (graph, from = from, dens = dens, k = k)
#' grep ("^flow", names (graph), value = TRUE)
#' # 5 "flow" columns again, but this time different dispersal coefficients each
#' # each origin point.
dodgr_flows_disperse <- function (graph,
                                  from,
                                  dens,
                                  k = 500,
                                  contract = TRUE,
                                  heap = "BHeap",
                                  tol = 1e-12,
                                  quiet = TRUE) {

    if (methods::is (graph, "dodgr_contracted")) {
        contract <- FALSE
    }

    res <- check_k (k, from)
    k <- res$k
    nk <- res$nk

    if (anyNA (dens)) {
        dens [is.na (dens)] <- 0
    }

    check_for_flow_col (graph)

    hps <- get_heap (heap, graph)
    heap <- hps$heap
    graph <- hps$graph

    graph <- preprocess_spatial_cols (graph)
    gr_cols <- dodgr_graph_cols (graph)

    to_from_indices <- to_from_index_with_tp (graph, from, to = NULL)
    if (to_from_indices$compound) {
        graph <- to_from_indices$graph_compound
    }

    if (contract) {
        graph_full <- graph
        graph <- contract_graph_with_pts (
            graph,
            to_from_indices$from$id,
            to = NULL
        )
        hashc <- get_hash (graph, contracted = TRUE)
        fname_c <- fs::path (
            fs::path_temp (),
            paste0 ("dodgr_edge_map_", hashc, ".Rds")
        )
        if (!fs::file_exists (fname_c)) {
            stop ("something went wrong extracting the edge_map ... ")
        } # nocov
        edge_map <- readRDS (fname_c)
    }

    graph2 <- convert_graph (graph, gr_cols)

    if (!is.matrix (dens)) {
        dens <- as.matrix (dens)
    }

    if (!quiet) {
        message ("\nAggregating flows ... ", appendLF = FALSE)
    }

    f <- rcpp_flows_disperse_par (
        graph2,
        to_from_indices$vert_map,
        to_from_indices$from$index,
        k,
        dens,
        tol,
        heap
    )

    if (nk == 1) {
        graph$flow <- f
    } else {
        flowmat <- data.frame (matrix (f, ncol = nk))
        names (flowmat) <- paste0 ("flow", seq_len (nk))
        graph <- cbind (graph, flowmat)
    }

    if (contract) { # map contracted flows back onto full graph
        graph <- uncontract_graph (graph, edge_map, graph_full)
    }

    flow_cols <- grep ("^flow", names (graph), value = TRUE)
    if (to_from_indices$compound) {
        graph <- uncompound_junctions (
            graph,
            flow_cols,
            to_from_indices$compound_junction_map
        )
    }
    graph [, flow_cols] [is.na (graph [, flow_cols])] <- 0

    return (graph)
}