diff --git a/DESCRIPTION b/DESCRIPTION
index c4f68347..9e1d05ed 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: dodgr
 Title: Distances on Directed Graphs
-Version: 0.4.3.030
+Version: 0.4.3.033
 Authors@R: c(
     person("Mark", "Padgham", , "mark.padgham@email.com", role = c("aut", "cre")),
     person("Andreas", "Petutschnig", role = "aut"),
diff --git a/R/flows-aggregate.R b/R/flows-aggregate.R
new file mode 100644
index 00000000..51affb83
--- /dev/null
+++ b/R/flows-aggregate.R
@@ -0,0 +1,237 @@
+#' @title Aggregate flows throughout a network.
+#'
+#' @description Aggregate flows throughout a network based on an input matrix
+#' of flows between all pairs of `from` and `to` points. 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.
+#'
+#' @param graph `data.frame` or equivalent object representing the network
+#' graph (see Details)
+#' @param flows Matrix of flows with `nrow(flows)==length(from)` and
+#' `ncol(flows)==length(to)`.
+#' @param pairwise If `TRUE`, aggregate flows only only paths connecting the
+#' ordered pairs of `from` and `to`. In this case, both `from` and `to` must be
+#' of the same length, and `flows` must be either a vector of the same length,
+#' or a matrix with only one column and same number of rows. `flows` then
+#' quantifies the flows between each pair of `from` and `to` points.
+#' @param contract If `TRUE` (default), calculate flows on contracted graph
+#' before mapping them back on to the original full graph (recommended as this
+#' will generally be much faster). `FALSE` should only be used if the `graph`
+#' has already been contracted.
+#' @param heap Type of heap to use in priority queue. Options include
+#' Fibonacci Heap (default; `FHeap`), Binary Heap (`BHeap`),
+#' Trinomial Heap (`TriHeap`), Extended Trinomial Heap
+#' (`TriHeapExt`, and 2-3 Heap (`Heap23`).
+#' @param tol Relative tolerance below which flows towards `to` vertices are not
+#' considered. This will generally have no effect, but can provide speed gains
+#' when flow matrices represent spatial interaction models, in which case this
+#' parameter effectively reduces the radius from each `from` point over which
+#' flows are aggregated. To remove any such effect, set `tol = 0`.
+#' @param norm_sums Standardise sums from all origin points, so sum of flows
+#' throughout entire network equals sum of densities from all origins (see
+#' Note).
+#' @inheritParams dodgr_dists
+#' @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.
+#'
+#' Flows are calculated by default using parallel computation with the maximal
+#' number of available cores or threads. This number can be reduced by
+#' specifying a value via
+#' `RcppParallel::setThreadOptions (numThreads = <desired_number>)`.
+#'
+#' @family flows
+#' @export
+#' @examples
+#' graph <- weight_streetnet (hampi)
+#' from <- sample (graph$from_id, size = 10)
+#' to <- sample (graph$to_id, size = 5)
+#' to <- to [!to %in% from]
+#' flows <- matrix (10 * runif (length (from) * length (to)),
+#'     nrow = length (from)
+#' )
+#' graph <- dodgr_flows_aggregate (graph, from = from, to = to, flows = flows)
+#' # 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
+#'
+#' # The following code can be used to convert the resultant graph to an `sf`
+#' # object suitable for plotting
+#' \dontrun{
+#' gsf <- dodgr_to_sf (graph_undir)
+#'
+#' # example of plotting with the 'mapview' package
+#' library (mapview)
+#' flow <- gsf$flow / max (gsf$flow)
+#' ncols <- 30
+#' cols <- c ("lawngreen", "red")
+#' colranmp <- colorRampPalette (cols) (ncols) [ceiling (ncols * flow)]
+#' mapview (gsf, color = colranmp, lwd = 10 * flow)
+#' }
+#'
+#' # An example of flow aggregation across a generic (non-OSM) highway,
+#' # represented as the `routes_fast` object of the \pkg{stplanr} package,
+#' # which is a SpatialLinesDataFrame containing commuter densities along
+#' # components of a street network.
+#' \dontrun{
+#' library (stplanr)
+#' # merge all of the 'routes_fast' lines into a single network
+#' r <- overline (routes_fast, attrib = "length", buff_dist = 1)
+#' r <- sf::st_as_sf (r)
+#' # then extract the start and end points of each of the original 'routes_fast'
+#' # lines and use these for routing with `dodgr`
+#' l <- lapply (routes_fast@lines, function (i) {
+#'     c (
+#'         sp::coordinates (i) [[1]] [1, ],
+#'         tail (sp::coordinates (i) [[1]], 1)
+#'     )
+#' })
+#' l <- do.call (rbind, l)
+#' xy_start <- l [, 1:2]
+#' xy_end <- l [, 3:4]
+#' # Then just specify a generic OD matrix with uniform values of 1:
+#' flows <- matrix (1, nrow = nrow (l), ncol = nrow (l))
+#' # We need to specify both a `type` and `id` column for the
+#' # \link{weight_streetnet} function.
+#' r$type <- 1
+#' r$id <- seq (nrow (r))
+#' graph <- weight_streetnet (
+#'     r,
+#'     type_col = "type",
+#'     id_col = "id",
+#'     wt_profile = 1
+#' )
+#' f <- dodgr_flows_aggregate (
+#'     graph,
+#'     from = xy_start,
+#'     to = xy_end,
+#'     flows = flows
+#' )
+#' # Then merge directed flows and convert to \pkg{sf} for plotting as before:
+#' f <- merge_directed_graph (f)
+#' geoms <- dodgr_to_sfc (f)
+#' gc <- dodgr_contract_graph (f)
+#' gsf <- sf::st_sf (geoms)
+#' gsf$flow <- gc$flow
+#' # sf plot:
+#' plot (gsf ["flow"])
+#' }
+dodgr_flows_aggregate <- function (graph,
+                                   from,
+                                   to,
+                                   flows,
+                                   pairwise = FALSE,
+                                   contract = TRUE,
+                                   heap = "BHeap",
+                                   tol = 1e-12,
+                                   norm_sums = TRUE,
+                                   quiet = TRUE) {
+
+    if (methods::is (graph, "dodgr_contracted")) {
+        contract <- FALSE
+    }
+
+    if (anyNA (flows)) {
+        flows [is.na (flows)] <- 0
+    }
+    hps <- get_heap (heap, graph)
+    heap <- hps$heap
+    graph <- hps$graph
+
+    if (!identical (class (from), class (to))) {
+        stop ("from and to must be the same class of object.")
+    }
+    check_for_flow_col (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 (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 (!is.matrix (flows)) {
+        flows <- matrix (flows, nrow = length (to_from_indices$from$index))
+    } else if (!(nrow (flows) == length (to_from_indices$from$index) &&
+        ncol (flows) == length (to_from_indices$to$index))) {
+        stop ("flows matrix is not compatible with 'from'/'to' arguments")
+    }
+    if (pairwise) {
+        check_pairwise_from_to (from, to, flows)
+    }
+
+    if (!quiet) {
+        message ("\nAggregating flows ... ", appendLF = FALSE)
+    }
+
+    if (pairwise) {
+        graph$flow <- rcpp_flows_aggregate_pairwise (
+            graph2, to_from_indices$vert_map,
+            to_from_indices$from$index, to_from_indices$to$index,
+            flows, norm_sums, tol, heap
+        )
+    } else {
+        graph$flow <- rcpp_flows_aggregate_par (
+            graph2, to_from_indices$vert_map,
+            to_from_indices$from$index, to_from_indices$to$index,
+            flows, norm_sums, tol, heap
+        )
+    }
+
+    if (contract) { # map contracted flows back onto full graph
+        graph <- uncontract_graph (graph, edge_map, graph_full)
+        graph$flow [is.na (graph$flow)] <- 0
+    }
+    if (to_from_indices$compound) {
+        graph <- uncompound_junctions (
+            graph,
+            "flow",
+            to_from_indices$compound_junction_map
+        )
+        graph$flow [is.na (graph$flow)] <- 0
+    }
+
+    return (graph)
+}
diff --git a/R/flows-disperse.R b/R/flows-disperse.R
new file mode 100644
index 00000000..edfe61e0
--- /dev/null
+++ b/R/flows-disperse.R
@@ -0,0 +1,192 @@
+#' @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)
+}
diff --git a/R/flows-si.R b/R/flows-si.R
new file mode 100644
index 00000000..29256dae
--- /dev/null
+++ b/R/flows-si.R
@@ -0,0 +1,248 @@
+#' @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)
+}
diff --git a/R/flows-utils.R b/R/flows-utils.R
new file mode 100644
index 00000000..5887934d
--- /dev/null
+++ b/R/flows-utils.R
@@ -0,0 +1,99 @@
+check_for_flow_col <- function (graph) {
+    if ("flow" %in% names (graph)) {
+        warning (
+            "graph already has a 'flow' column; ",
+            "this will be overwritten"
+        )
+    }
+}
+
+check_k <- function (k, from) {
+
+    nk <- 1
+
+    if (is.data.frame (k)) {
+        k <- as.matrix (k)
+    }
+
+    if (is.matrix (k)) {
+        if (nrow (k) != length (from)) {
+            stop ("nrow(k) must equal length of 'from' points")
+        }
+        nk <- ncol (k)
+    } else if (is.numeric (k)) {
+        if (length (k) == 1) {
+            k <- rep (k, length (from))
+        } else if (length (k) != length (from)) {
+            # convert to matrix
+            nk <- length (k)
+            k <- array (rep (k, each = length (from)),
+                dim = c (length (from), nk)
+            )
+        }
+    } else {
+        stop ("'k' must be either a single value, a vector, or a matrix")
+    }
+
+    list (k = k, nk = nk)
+}
+
+check_pairwise_from_to <- function (from, to, flows) {
+
+    nf <- ifelse (is.vector (from), length (from), nrow (from))
+    nt <- ifelse (is.vector (to), length (to), nrow (to))
+
+    if (nf != nt) {
+        stop (
+            "'from' and 'to' must be the same length or dimensions",
+            "when using 'pairwise = TRUE'.",
+            call. = FALSE
+        )
+    }
+    if (nrow (flows) != nf) {
+        stop (
+            "'flows' must be the same length or dimensions as 'from' and 'to'",
+            call. = FALSE
+        )
+    }
+}
+
+get_random_prefix <- function (prefix = "flow",
+                               n = 5) {
+
+    charvec <- c (letters, LETTERS, 0:9)
+    rand_prefix <- paste (sample (charvec, n, replace = TRUE), collapse = "")
+    fs::path (fs::path_temp (), paste0 (prefix, "_", rand_prefix))
+}
+
+nodes_arg_to_pts <- function (nodes,
+                              graph) {
+
+    if (!is.matrix (nodes) && is.numeric (nodes)) {
+        nodes <- matrix (nodes, ncol = 2)
+    }
+    if (is.vector (nodes)) {
+        # non-numeric, so must be vector of node IDs
+        return (nodes)
+    }
+    if (ncol (nodes) == 2) {
+        verts <- dodgr_vertices (graph)
+        nodes <- verts$id [match_pts_to_verts (verts, nodes)]
+    }
+    return (nodes)
+}
+
+
+# keep from and to routing points in contracted graph
+contract_graph_with_pts <- function (graph,
+                                     from,
+                                     to) {
+
+    pts <- NULL
+    if (!missing (from)) {
+        pts <- c (pts, from)
+    }
+    if (!missing (to)) {
+        pts <- c (pts, to)
+    }
+    dodgr_contract_graph (graph, unique (pts))
+}
diff --git a/R/flows.R b/R/flows.R
deleted file mode 100644
index b75aa3fe..00000000
--- a/R/flows.R
+++ /dev/null
@@ -1,779 +0,0 @@
-#' @title Aggregate flows throughout a network.
-#'
-#' @description Aggregate flows throughout a network based on an input matrix
-#' of flows between all pairs of `from` and `to` points. 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.
-#'
-#' @param graph `data.frame` or equivalent object representing the network
-#' graph (see Details)
-#' @param flows Matrix of flows with `nrow(flows)==length(from)` and
-#' `ncol(flows)==length(to)`.
-#' @param pairwise If `TRUE`, aggregate flows only only paths connecting the
-#' ordered pairs of `from` and `to`. In this case, both `from` and `to` must be
-#' of the same length, and `flows` must be either a vector of the same length,
-#' or a matrix with only one column and same number of rows. `flows` then
-#' quantifies the flows between each pair of `from` and `to` points.
-#' @param contract If `TRUE` (default), calculate flows on contracted graph
-#' before mapping them back on to the original full graph (recommended as this
-#' will generally be much faster). `FALSE` should only be used if the `graph`
-#' has already been contracted.
-#' @param heap Type of heap to use in priority queue. Options include
-#' Fibonacci Heap (default; `FHeap`), Binary Heap (`BHeap`),
-#' Trinomial Heap (`TriHeap`), Extended Trinomial Heap
-#' (`TriHeapExt`, and 2-3 Heap (`Heap23`).
-#' @param tol Relative tolerance below which flows towards `to` vertices are not
-#' considered. This will generally have no effect, but can provide speed gains
-#' when flow matrices represent spatial interaction models, in which case this
-#' parameter effectively reduces the radius from each `from` point over which
-#' flows are aggregated. To remove any such effect, set `tol = 0`.
-#' @param norm_sums Standardise sums from all origin points, so sum of flows
-#' throughout entire network equals sum of densities from all origins (see
-#' Note).
-#' @inheritParams dodgr_dists
-#' @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.
-#'
-#' Flows are calculated by default using parallel computation with the maximal
-#' number of available cores or threads. This number can be reduced by
-#' specifying a value via
-#' `RcppParallel::setThreadOptions (numThreads = <desired_number>)`.
-#'
-#' @family flows
-#' @export
-#' @examples
-#' graph <- weight_streetnet (hampi)
-#' from <- sample (graph$from_id, size = 10)
-#' to <- sample (graph$to_id, size = 5)
-#' to <- to [!to %in% from]
-#' flows <- matrix (10 * runif (length (from) * length (to)),
-#'     nrow = length (from)
-#' )
-#' graph <- dodgr_flows_aggregate (graph, from = from, to = to, flows = flows)
-#' # 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
-#'
-#' # The following code can be used to convert the resultant graph to an `sf`
-#' # object suitable for plotting
-#' \dontrun{
-#' gsf <- dodgr_to_sf (graph_undir)
-#'
-#' # example of plotting with the 'mapview' package
-#' library (mapview)
-#' flow <- gsf$flow / max (gsf$flow)
-#' ncols <- 30
-#' cols <- c ("lawngreen", "red")
-#' colranmp <- colorRampPalette (cols) (ncols) [ceiling (ncols * flow)]
-#' mapview (gsf, color = colranmp, lwd = 10 * flow)
-#' }
-#'
-#' # An example of flow aggregation across a generic (non-OSM) highway,
-#' # represented as the `routes_fast` object of the \pkg{stplanr} package,
-#' # which is a SpatialLinesDataFrame containing commuter densities along
-#' # components of a street network.
-#' \dontrun{
-#' library (stplanr)
-#' # merge all of the 'routes_fast' lines into a single network
-#' r <- overline (routes_fast, attrib = "length", buff_dist = 1)
-#' r <- sf::st_as_sf (r)
-#' # then extract the start and end points of each of the original 'routes_fast'
-#' # lines and use these for routing with `dodgr`
-#' l <- lapply (routes_fast@lines, function (i) {
-#'     c (
-#'         sp::coordinates (i) [[1]] [1, ],
-#'         tail (sp::coordinates (i) [[1]], 1)
-#'     )
-#' })
-#' l <- do.call (rbind, l)
-#' xy_start <- l [, 1:2]
-#' xy_end <- l [, 3:4]
-#' # Then just specify a generic OD matrix with uniform values of 1:
-#' flows <- matrix (1, nrow = nrow (l), ncol = nrow (l))
-#' # We need to specify both a `type` and `id` column for the
-#' # \link{weight_streetnet} function.
-#' r$type <- 1
-#' r$id <- seq (nrow (r))
-#' graph <- weight_streetnet (
-#'     r,
-#'     type_col = "type",
-#'     id_col = "id",
-#'     wt_profile = 1
-#' )
-#' f <- dodgr_flows_aggregate (
-#'     graph,
-#'     from = xy_start,
-#'     to = xy_end,
-#'     flows = flows
-#' )
-#' # Then merge directed flows and convert to \pkg{sf} for plotting as before:
-#' f <- merge_directed_graph (f)
-#' geoms <- dodgr_to_sfc (f)
-#' gc <- dodgr_contract_graph (f)
-#' gsf <- sf::st_sf (geoms)
-#' gsf$flow <- gc$flow
-#' # sf plot:
-#' plot (gsf ["flow"])
-#' }
-dodgr_flows_aggregate <- function (graph,
-                                   from,
-                                   to,
-                                   flows,
-                                   pairwise = FALSE,
-                                   contract = TRUE,
-                                   heap = "BHeap",
-                                   tol = 1e-12,
-                                   norm_sums = TRUE,
-                                   quiet = TRUE) {
-
-    if (methods::is (graph, "dodgr_contracted")) {
-        contract <- FALSE
-    }
-
-    if (anyNA (flows)) {
-        flows [is.na (flows)] <- 0
-    }
-    hps <- get_heap (heap, graph)
-    heap <- hps$heap
-    graph <- hps$graph
-
-    if (!identical (class (from), class (to))) {
-        stop ("from and to must be the same class of object.")
-    }
-    check_for_flow_col (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 (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 (!is.matrix (flows)) {
-        flows <- matrix (flows, nrow = length (to_from_indices$from$index))
-    } else if (!(nrow (flows) == length (to_from_indices$from$index) &&
-        ncol (flows) == length (to_from_indices$to$index))) {
-        stop ("flows matrix is not compatible with 'from'/'to' arguments")
-    }
-    if (pairwise) {
-        check_pairwise_from_to (from, to, flows)
-    }
-
-    if (!quiet) {
-        message ("\nAggregating flows ... ", appendLF = FALSE)
-    }
-
-    if (pairwise) {
-        graph$flow <- rcpp_flows_aggregate_pairwise (
-            graph2, to_from_indices$vert_map,
-            to_from_indices$from$index, to_from_indices$to$index,
-            flows, norm_sums, tol, heap
-        )
-    } else {
-        graph$flow <- rcpp_flows_aggregate_par (
-            graph2, to_from_indices$vert_map,
-            to_from_indices$from$index, to_from_indices$to$index,
-            flows, norm_sums, tol, heap
-        )
-    }
-
-    if (contract) { # map contracted flows back onto full graph
-        graph <- uncontract_graph (graph, edge_map, graph_full)
-        graph$flow [is.na (graph$flow)] <- 0
-    }
-    if (to_from_indices$compound) {
-        graph <- uncompound_junctions (
-            graph,
-            "flow",
-            to_from_indices$compound_junction_map
-        )
-        graph$flow [is.na (graph$flow)] <- 0
-    }
-
-    return (graph)
-}
-
-#' @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)
-}
-
-#' @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)
-}
-
-check_for_flow_col <- function (graph) {
-    if ("flow" %in% names (graph)) {
-        warning (
-            "graph already has a 'flow' column; ",
-            "this will be overwritten"
-        )
-    }
-}
-
-check_k <- function (k, from) {
-
-    nk <- 1
-
-    if (is.data.frame (k)) {
-        k <- as.matrix (k)
-    }
-
-    if (is.matrix (k)) {
-        if (nrow (k) != length (from)) {
-            stop ("nrow(k) must equal length of 'from' points")
-        }
-        nk <- ncol (k)
-    } else if (is.numeric (k)) {
-        if (length (k) == 1) {
-            k <- rep (k, length (from))
-        } else if (length (k) != length (from)) {
-            # convert to matrix
-            nk <- length (k)
-            k <- array (rep (k, each = length (from)),
-                dim = c (length (from), nk)
-            )
-        }
-    } else {
-        stop ("'k' must be either a single value, a vector, or a matrix")
-    }
-
-    list (k = k, nk = nk)
-}
-
-check_pairwise_from_to <- function (from, to, flows) {
-
-    nf <- ifelse (is.vector (from), length (from), nrow (from))
-    nt <- ifelse (is.vector (to), length (to), nrow (to))
-
-    if (nf != nt) {
-        stop (
-            "'from' and 'to' must be the same length or dimensions",
-            "when using 'pairwise = TRUE'.",
-            call. = FALSE
-        )
-    }
-    if (nrow (flows) != nf) {
-        stop (
-            "'flows' must be the same length or dimensions as 'from' and 'to'",
-            call. = FALSE
-        )
-    }
-}
-
-get_random_prefix <- function (prefix = "flow",
-                               n = 5) {
-
-    charvec <- c (letters, LETTERS, 0:9)
-    rand_prefix <- paste (sample (charvec, n, replace = TRUE), collapse = "")
-    fs::path (fs::path_temp (), paste0 (prefix, "_", rand_prefix))
-}
-
-nodes_arg_to_pts <- function (nodes,
-                              graph) {
-
-    if (!is.matrix (nodes) && is.numeric (nodes)) {
-        nodes <- matrix (nodes, ncol = 2)
-    }
-    if (is.vector (nodes)) {
-        # non-numeric, so must be vector of node IDs
-        return (nodes)
-    }
-    if (ncol (nodes) == 2) {
-        verts <- dodgr_vertices (graph)
-        nodes <- verts$id [match_pts_to_verts (verts, nodes)]
-    }
-    return (nodes)
-}
-
-
-# keep from and to routing points in contracted graph
-contract_graph_with_pts <- function (graph,
-                                     from,
-                                     to) {
-
-    pts <- NULL
-    if (!missing (from)) {
-        pts <- c (pts, from)
-    }
-    if (!missing (to)) {
-        pts <- c (pts, to)
-    }
-    dodgr_contract_graph (graph, unique (pts))
-}
diff --git a/codemeta.json b/codemeta.json
index d07e29d6..ba6e983e 100644
--- a/codemeta.json
+++ b/codemeta.json
@@ -11,7 +11,7 @@
   "codeRepository": "https://github.com/UrbanAnalyst/dodgr",
   "issueTracker": "https://github.com/UrbanAnalyst/dodgr/issues",
   "license": "https://spdx.org/licenses/GPL-3.0",
-  "version": "0.4.3.030",
+  "version": "0.4.3.033",
   "programmingLanguage": {
     "@type": "ComputerLanguage",
     "name": "R",
diff --git a/man/dodgr_flows_aggregate.Rd b/man/dodgr_flows_aggregate.Rd
index ddf7fbcc..26a26195 100644
--- a/man/dodgr_flows_aggregate.Rd
+++ b/man/dodgr_flows_aggregate.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/flows.R
+% Please edit documentation in R/flows-aggregate.R
 \name{dodgr_flows_aggregate}
 \alias{dodgr_flows_aggregate}
 \title{Aggregate flows throughout a network.}
diff --git a/man/dodgr_flows_disperse.Rd b/man/dodgr_flows_disperse.Rd
index 3d4e50b3..9df6ed9f 100644
--- a/man/dodgr_flows_disperse.Rd
+++ b/man/dodgr_flows_disperse.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/flows.R
+% Please edit documentation in R/flows-disperse.R
 \name{dodgr_flows_disperse}
 \alias{dodgr_flows_disperse}
 \title{Aggregate flows dispersed from each point in a network.}
diff --git a/man/dodgr_flows_si.Rd b/man/dodgr_flows_si.Rd
index d3f76bab..9f5426f2 100644
--- a/man/dodgr_flows_si.Rd
+++ b/man/dodgr_flows_si.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/flows.R
+% Please edit documentation in R/flows-si.R
 \name{dodgr_flows_si}
 \alias{dodgr_flows_si}
 \title{Aggregate flows throughout a network using a spatial interaction model.}
diff --git a/src/pathfinders.cpp b/src/pathfinders.cpp
index 81354529..439a1215 100644
--- a/src/pathfinders.cpp
+++ b/src/pathfinders.cpp
@@ -91,7 +91,7 @@ void PF::PathFinder::scan_edges_heur (const DGraphEdge *edge,
         bool *m_open_vec,
         const bool *m_closed_vec,
         const size_t &v0,
-        const std::vector<double> &heur)    // heuristic for A*
+        const PF::AStarHeuristic &heur)    // heuristic for A*
 {
     while (edge) {
         size_t et = edge->target;
@@ -243,7 +243,7 @@ void PF::PathFinder::DijkstraLimit (
 void PF::PathFinder::AStar (std::vector<double>& d,
         std::vector<double>& w,
         std::vector<long int>& prev,
-        const std::vector<double>& heur,
+        const PF::AStarHeuristic& heur,
         const size_t v0,
         const std::vector <size_t> &to_index)
 {
diff --git a/src/pathfinders.h b/src/pathfinders.h
index 5de38baa..238af77e 100644
--- a/src/pathfinders.h
+++ b/src/pathfinders.h
@@ -22,6 +22,33 @@ class DGraph;
 
 namespace PF{
 
+// Lazily-evaluated Euclidean A* heuristic: straight-line distances from the
+// search origin are computed on demand as each vertex is first encountered,
+// rather than pre-computed for every vertex in the graph up front. A* with
+// early termination on reaching all targets often visits only a small
+// fraction of vertices, so eager pre-computation of the full array is wasted
+// work on large graphs.
+struct AStarHeuristic
+{
+    const std::vector <double> &vx;
+    const std::vector <double> &vy;
+    const double x0;
+    const double y0;
+
+    AStarHeuristic (const std::vector <double> &vx_in,
+            const std::vector <double> &vy_in,
+            const size_t v0) :
+        vx (vx_in), vy (vy_in), x0 (vx_in [v0]), y0 (vy_in [v0])
+    {
+    }
+
+    double operator[] (const size_t v) const
+    {
+        const double dx = vx [v] - x0, dy = vy [v] - y0;
+        return sqrt (dx * dx + dy * dy);
+    }
+};
+
 /* DijkstraEdge is used for the std::set implementation, everything else is for
  * Shane Saunders's heap sort versions */
 struct DijkstraEdge
@@ -85,7 +112,7 @@ class PathFinder {
                 bool *m_open_vec,
                 const bool *m_closed_vec,
                 const size_t &v0,
-                const std::vector<double> &heur);
+                const AStarHeuristic &heur);
         // with A* heuristic for dists-categorical
         void scan_edge_types_heur (
                 const DGraphEdge *edge,
@@ -95,7 +122,7 @@ class PathFinder {
                 bool *m_open_vec,
                 const bool *m_closed_vec,
                 const size_t &v0,
-                const std::vector<double> &heur);
+                const AStarHeuristic &heur);
         // run_sp_categorical for threshold dists
         void scan_edge_types (
                 const DGraphEdge *edge,
@@ -133,13 +160,13 @@ class PathFinder {
         void AStar (std::vector<double>& d,
                 std::vector<double>& w,
                 std::vector<long int>& prev,
-                const std::vector<double>& heur,
+                const AStarHeuristic& heur,
                 const size_t v0,
                 const std::vector <size_t> &to_index);
         void AStarEdgeType (std::vector<double>& d,
                 std::vector<double>& w,
                 std::vector<long int>& prev,
-                const std::vector<double>& heur,
+                const AStarHeuristic& heur,
                 const size_t v0,
                 const std::vector <size_t> &to_index);
         void Dijkstra_set (std::vector <double>& d,
diff --git a/src/run_sp.cpp b/src/run_sp.cpp
index 094b21dd..17220e7a 100644
--- a/src/run_sp.cpp
+++ b/src/run_sp.cpp
@@ -100,7 +100,6 @@ struct OneDist : public RcppParallel::Worker
         std::vector <double> w (nverts);
         std::vector <double> d (nverts);
         std::vector <long int> prev (nverts);
-        std::vector <double> heuristic (nverts, 0.0);
 
         for (std::size_t i = begin; i < end; i++)
         {
@@ -108,12 +107,10 @@ struct OneDist : public RcppParallel::Worker
 
             if (is_spatial)
             {
-                for (size_t j = 0; j < nverts; j++)
-                {
-                    const double dx = vx [j] - vx [from_i],
-                        dy = vy [j] - vy [from_i];
-                    heuristic [j] = sqrt (dx * dx + dy * dy);
-                }
+                // Heuristic values are computed lazily, only for vertices
+                // actually visited by AStar, rather than pre-computed for
+                // every vertex in the graph.
+                const PF::AStarHeuristic heuristic (vx, vy, from_i);
                 pathfinder->AStar (d, w, prev, heuristic, from_i, toi);
             } else if (heap_type.find ("set") == std::string::npos)
                 pathfinder->Dijkstra (d, w, prev, from_i, toi);
@@ -228,7 +225,6 @@ struct OneDistPaired : public RcppParallel::Worker
         std::vector <double> w (nverts);
         std::vector <double> d (nverts);
         std::vector <long int> prev (nverts);
-        std::vector <double> heuristic (nverts, 0.0);
 
         for (std::size_t i = begin; i < end; i++)
         {
@@ -237,18 +233,21 @@ struct OneDistPaired : public RcppParallel::Worker
 
             if (is_spatial)
             {
+                const PF::AStarHeuristic heuristic (vx, vy, from_i);
+
                 // need to set an additional target vertex that is somewhat
                 // beyond the single actual target vertex. Default here is max
-                // heuristic, but reduced in following loop.
+                // heuristic, but reduced in following loop. This selection
+                // needs the heuristic value of every vertex, so gains
+                // nothing from lazy evaluation, unlike the AStar search
+                // itself below.
                 long int max_h_index = -1;
                 double max_h_value = -1.0;
                 for (size_t j = 0; j < nverts; j++)
                 {
-                    const double dx = vx [j] - vx [from_i],
-                        dy = vy [j] - vy [from_i];
-                    heuristic [j] = sqrt (dx * dx + dy * dy);
-                    if (heuristic [j] > max_h_value) {
-                        max_h_value = heuristic [j];
+                    const double h = heuristic [j];
+                    if (h > max_h_value) {
+                        max_h_value = h;
                         max_h_index = static_cast <long int> (j);
                     }
                 }
@@ -260,8 +259,9 @@ struct OneDistPaired : public RcppParallel::Worker
                 // be adjusted?
                 const double thr = 0.1;
                 for (size_t j = 0; j < nverts; j++) {
-                    if ((heuristic [j] < (thr * htemp)) && (heuristic [j] > min_h_value)) {
-                        min_h_value = heuristic [j];
+                    const double h = heuristic [j];
+                    if ((h < (thr * htemp)) && (h > min_h_value)) {
+                        min_h_value = h;
                         min_h_index = static_cast <long int> (j);
                     }
                 }
diff --git a/src/run_sp_categorical.cpp b/src/run_sp_categorical.cpp
index a868b639..e6655386 100644
--- a/src/run_sp_categorical.cpp
+++ b/src/run_sp_categorical.cpp
@@ -72,17 +72,11 @@ struct OneCategoricalDist : public RcppParallel::Worker
             std::vector <double> d (nverts * (num_edge_types + 1));
             std::vector <long int> prev (nverts);
 
-            std::vector <double> heuristic (nverts, 0.0);
-
             size_t from_i = static_cast <size_t> (dp_fromi [i]);
 
-            // only implemented for spatial graphs
-            for (size_t j = 0; j < nverts; j++)
-            {
-                const double dx = vx [j] - vx [from_i],
-                    dy = vy [j] - vy [from_i];
-                heuristic [j] = sqrt (dx * dx + dy * dy);
-            }
+            // only implemented for spatial graphs; heuristic values are
+            // computed lazily, only for vertices actually visited by AStar.
+            const PF::AStarHeuristic heuristic (vx, vy, from_i);
             pathfinder->AStarEdgeType (d, w, prev, heuristic, from_i, toi);
 
             for (size_t j = 0; j < toi.size (); j++)
@@ -97,7 +91,7 @@ struct OneCategoricalDist : public RcppParallel::Worker
             }
         }
     }
-                                   
+
 };
 
 struct OnePairedCategoricalDist : public RcppParallel::Worker
@@ -146,18 +140,12 @@ struct OnePairedCategoricalDist : public RcppParallel::Worker
             std::vector <double> d (nverts * (num_edge_types + 1));
             std::vector <long int> prev (nverts);
 
-            std::vector <double> heuristic (nverts, 0.0);
-
             const size_t from_i = static_cast <size_t> (dp_fromi [i]);
             const std::vector <size_t> to_i = { toi [i] };
 
-            // only implemented for spatial graphs
-            for (size_t j = 0; j < nverts; j++)
-            {
-                const double dx = vx [j] - vx [from_i],
-                    dy = vy [j] - vy [from_i];
-                heuristic [j] = sqrt (dx * dx + dy * dy);
-            }
+            // only implemented for spatial graphs; heuristic values are
+            // computed lazily, only for vertices actually visited by AStar.
+            const PF::AStarHeuristic heuristic (vx, vy, from_i);
             pathfinder->AStarEdgeType (d, w, prev, heuristic, from_i, to_i);
 
             for (size_t k = 0; k <= num_edge_types; k++)
@@ -218,17 +206,11 @@ struct OneCategory : public RcppParallel::Worker
             std::vector <double> d (nverts * (num_edge_types + 1));
             std::vector <long int> prev (nverts);
 
-            std::vector <double> heuristic (nverts, 0.0);
-
             size_t from_i = static_cast <size_t> (dp_fromi [i]);
 
-            // only implemented for spatial graphs
-            for (size_t j = 0; j < nverts; j++)
-            {
-                const double dx = vx [j] - vx [from_i],
-                    dy = vy [j] - vy [from_i];
-                heuristic [j] = sqrt (dx * dx + dy * dy);
-            }
+            // only implemented for spatial graphs; heuristic values are
+            // computed lazily, only for vertices actually visited by AStar.
+            const PF::AStarHeuristic heuristic (vx, vy, from_i);
             pathfinder->AStarEdgeType (d, w, prev, heuristic, from_i, toi);
 
             for (size_t j = 0; j < toi.size (); j++)
@@ -350,7 +332,7 @@ size_t categorical::get_num_edge_types (const std::vector <size_t> &edge_type)
 void PF::PathFinder::AStarEdgeType (std::vector<double>& d,
         std::vector<double>& w,
         std::vector<long int>& prev,
-        const std::vector<double>& heur,
+        const PF::AStarHeuristic& heur,
         const size_t v0,
         const std::vector <size_t> &to_index)
 {
@@ -400,7 +382,7 @@ void PF::PathFinder::scan_edge_types_heur (const DGraphEdge *edge,
         bool *m_open_vec,
         const bool *m_closed_vec,
         const size_t &v0,
-        const std::vector<double> &heur)    // heuristic for A*
+        const PF::AStarHeuristic &heur)    // heuristic for A*
 {
     const size_t nverts = w.size ();
     const size_t num_edge_types = d.size () / nverts - 1L;