R/flows-si.R

In dodgr: Distances on Directed Graphs

#' @title Aggregate flows throughout a network using a spatial interaction model.
#'
#' @description Aggregate flows throughout a network using an exponential
#' Spatial Interaction (SI) model between a specified set of origin and
#' destination points, and associated vectors of densities. Spatial
#' interactions are implemented using an exponential decay, controlled by a
#' parameter, `k`, so that interactions decay with `exp(-d / k)`, where `d` is
#' distance. The algorithm allows for efficient fitting of multiple interaction
#' models for different coefficients to be fitted with a single call. Values of
#' the interaction coefficients, `k`, may take one of the following forms:
#'
#' \itemize{
#' \item A single numeric value (> 0), with interactions 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
#' interactions 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 interaction 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 interaction
#' 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 k Width of exponential spatial interaction function (exp (-d / k)),
#' in units of 'd', specified in one of 3 forms: (i) a single value; (ii) a
#' vector of independent values for each origin point (with same length as
#' 'from' points); or (iii) an equivalent matrix with each column holding values
#' for each 'from' point, so 'nrow(k)==length(from)'. See Note.
#' @param dens_from Vector of densities at origin ('from') points
#' @param dens_to Vector of densities at destination ('to') points
#' @return Modified version of graph with additional `flow` column added.
#'
#' @note The `norm_sums` parameter should be used whenever densities at origins
#' and destinations are absolute values, and ensures that the sum of resultant
#' flow values throughout the entire network equals the sum of densities at all
#' origins. For example, with `norm_sums = TRUE` (the default), a flow from a
#' single origin with density one to a single destination along two edges will
#' allocate flows of one half to each of those edges, such that the sum of flows
#' across the network will equal one, or the sum of densities from all origins.
#' The `norm_sums = TRUE` option is appropriate where densities are relative
#' values, and ensures that each edge maintains relative proportions. In the
#' above example, flows along each of two edges would equal one, for a network
#' sum of two, or greater than the sum of densities.
#'
#' With `norm_sums = TRUE`, the sum of network flows (`sum(output$flow)`) should
#' equal the sum of origin densities (`sum(dens_from)`). This may nevertheless
#' not always be the case, because origin points may simply be too far from any
#' destination (`to`) points for an exponential model to yield non-zero values
#' anywhere in a network within machine tolerance. Such cases may result in sums
#' of output flows being less than sums of input densities.
#'
#' @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)
#' to <- sample (graph$from_id, size = 20)
#' dens_from <- runif (length (from))
#' dens_to <- runif (length (to))
#' graph <- dodgr_flows_si (
#'     graph,
#'     from = from,
#'     to = to,
#'     dens_from = dens_from,
#'     dens_to = dens_to
#' )
#' # 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)
#' # This graph will only include those edges having non-zero flows, and so:
#' nrow (graph)
#' nrow (graph_undir) # the latter is much smaller
#'
#' # ----- One dispersal coefficient for each origin point:
#' # Remove `flow` column to avoid warning about over-writing values:
#' graph$flow <- NULL
#' k <- runif (length (from))
#' graph <- dodgr_flows_si (
#'     graph,
#'     from = from,
#'     to = to,
#'     dens_from = dens_from,
#'     dens_to = dens_to,
#'     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_si (
#'     graph,
#'     from = from,
#'     to = to,
#'     dens_from = dens_from,
#'     dens_to = dens_to,
#'     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_si (
#'     graph,
#'     from = from,
#'     to = to,
#'     dens_from = dens_from,
#'     dens_to = dens_to,
#'     k = k
#' )
#' grep ("^flow", names (graph), value = TRUE)
#' # 5 "flow" columns again, but this time different dispersal coefficients each
#' # each origin point.
dodgr_flows_si <- function (graph,
                            from,
                            to,
                            k = 500,
                            dens_from = NULL,
                            dens_to = NULL,
                            contract = TRUE,
                            norm_sums = TRUE,
                            heap = "BHeap",
                            tol = 1e-12,
                            quiet = TRUE) {

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

    if (missing (from)) {
        stop ("'from' must be provided for spatial interaction models.")
    }

    check_for_flow_col (graph)

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

    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)
    if (to_from_indices$compound) {
        graph <- to_from_indices$graph_compound
    }

    if (is.null (dens_from)) {
        dens_from <- rep (1, length (from))
    }
    if (is.null (dens_to)) {
        dens_to <- rep (1, length (to))
    }

    if (contract) {
        graph_full <- graph
        graph <- contract_graph_with_pts (
            graph,
            to_from_indices$from$id,
            to_from_indices$to$id
        )
        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 (!quiet) {
        message ("\nAggregating flows ... ", appendLF = FALSE)
    }

    f <- rcpp_flows_si (
        graph2,
        to_from_indices$vert_map,
        to_from_indices$from$index,
        to_from_indices$to$index,
        k,
        dens_from,
        dens_to,
        norm_sums,
        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) {
        flow_cols <- grep ("^flow", names (graph), value = TRUE)
        graph <- uncompound_junctions (
            graph,
            flow_cols,
            to_from_indices$compound_junction_map
        )
    }
    graph [, flow_cols] [is.na (graph [, flow_cols])] <- 0

    return (graph)
}