| Title: | Boltzmann–Lotka–Volterra Interaction Model |
|---|---|
| Description: | Estimates Boltzmann–Lotka–Volterra (BLV) interaction model efficiently. Enables programmatic and graphical exploration of the solution space of BLV models when parameters are varied. See Wilson, A. (2008) <dx.doi.org/10.1098/rsif.2007.1288>. |
| Authors: | Fabrice Rossi [aut, cre, cph] (ORCID: <https://orcid.org/0000-0003-4638-1286>), Cyprien Gout [ctb] |
| Maintainer: | Fabrice Rossi <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.1.1.9000 |
| Built: | 2026-10-09 16:47:11 UTC |
| Source: | https://github.com/fabrice-rossi/blvim |
This function creates a data frame with a single column storing its
collection of spatial interaction models. The default name of the column
is "sim" (can be modified using the sim_column parameter). An more
flexible alternative is provided by the sim_df() function.
## S3 method for class 'sim_list' as.data.frame(x, ..., sim_column = "sim")## S3 method for class 'sim_list' as.data.frame(x, ..., sim_column = "sim")
x |
a collection of spatial interaction models, an object of class
|
... |
additional parameters (not used currently) |
sim_column |
the name of the |
a data frame
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows_log <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) as.data.frame(all_flows_log, sim_column = "log flows")distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows_log <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) as.data.frame(all_flows_log, sim_column = "log flows")
Extract the attractivenesses from a spatial interaction model object
attractiveness(sim, ...)attractiveness(sim, ...)
sim |
a spatial interaction model object |
... |
additional parameters |
a vector of attractivenesses at the destination locations
production(), destination_flow()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) attractiveness(model) ## the names of the attractiveness vector are set from the distance matrix ## we remove them for testing equality all.equal(as.numeric(attractiveness(model)), attractiveness)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) attractiveness(model) ## the names of the attractiveness vector are set from the distance matrix ## we remove them for testing equality all.equal(as.numeric(attractiveness(model)), attractiveness)
This function represents graphical the flows of a spatial interaction model, in different direct or aggregated forms.
## S3 method for class 'sim' autoplot( object, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, show_destination = FALSE, show_attractiveness = FALSE, show_production = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, ... )## S3 method for class 'sim' autoplot( object, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, show_destination = FALSE, show_attractiveness = FALSE, show_production = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, ... )
object |
a spatial interaction model object |
flows |
|
with_names |
specifies whether the graphical representation includes
location names ( |
with_positions |
specifies whether the graphical representation is based
on location positions ( |
show_destination |
specifies whether the position based |
show_attractiveness |
specifies whether the position based |
show_production |
specifies whether the position based |
cut_off |
cut off limit for inclusion of a graphical primitive when
|
adjust_limits |
if |
with_labels |
if |
... |
additional parameters, see details |
The graphical representation depends on the values of flows and
with_positions. It is based on the data frame representation produced by
fortify.sim().
If with_position is FALSE (default value), the graphical representations
are "abstract". Depending on flows we have the following representations:
"full": this is the default case for which the full flow matrix is
represented. It is extracted from the spatial interaction model with
flows() and displayed using a matrix representation with origin locations
in rows and destination locations in columns. The colour of a cell
corresponds to the intensity of a flow between the corresponding locations.
To mimic the standard top to bottom reading order of a flow matrix, the top
row of the graphical representation corresponds to the first origin location.
"destination": the function computes the
incoming flows for destination locations (using destination_flow()) and
represents them with a bar plot (each bar is proportional to the incoming
flow);
"attractiveness": the function uses a bar plot to represent the
attractivenesses of the destination locations (as given by
attractiveness()). This is interesting for dynamic models where those
values are updated during the iterations (see blvim() for details). When
the calculation has converged (see sim_converged()), both "destination"
and "attractiveness" graphics should be almost identical.
When the with_names parameter is TRUE, the location names
(location_names()) are used to label the axis of the graphical
representation. If names are not specified, they are replaced by indexes.
When the with_positions parameter is TRUE, the location positions
(location_positions()) are used to produce more "geographically informed"
representations. Notice that if no positions are known for the locations, the
use of with_positions = TRUE is an error. Depending on flows we have the
following representations:
"full": this is the default case for which the full flow matrix is
represented. Positions for both origin and destination locations are needed.
The representation uses arrows from origin location positions to destination
location positions. The thickness of the lines (linewidth aesthetics) is
proportional to the flows. Only segments that carry a flow above the
cut_off value are included. When the spatial interaction model is not
bipartite (see sim_is_bipartite()), zero length segments corresponding to
self exchange are removed. Additional parameters in ... are submitted to
ggplot2::geom_segment(). This can be used to override defaults parameters
used for the arrow shapes, for instance. Those parameters must be named. In
addition to the individual flows, the representation can include location
based information. If show_production is TRUE, the production constraints
(obtained by production()) are shown as disks centred on the origin
locations. If show_destination is TRUE, incoming flows (obtained by
destination_flow()) are shown as disks centred on the destination
locations. If show_attractiveness is TRUE, attractivenesses (obtained by
attractiveness()) are shown as disks centred on the destination locations.
show_destination and show_attractiveness are not compatible (only one can
be TRUE). show_production is compatible with show_destination or
show_attractiveness for bipartite models only (see sim_is_bipartite()).
When disks are used, segments are drawn without arrows, while the default
drawing uses arrows. This can be modified with the additional parameters.
"destination": the function draws a disk at each destination location
using for the size aesthetics the incoming flow at this destination
location (using destination_flow()). Only destinations with an incoming
flow above the cut_off value are included.
"attractiveness": the function draws a disk at each destination location
using for the size aesthetics the attractiveness of the destination. When
the calculation has converged (see sim_converged()), both "destination"
and "attractiveness" graphics should be almost identical. Only destinations
with an attractiveness above the cut_off value are included.
For the position based representations and when with_names is TRUE, the
names of the destinations are added to the graphical representation . If
with_labels is TRUE the names are represented as labels instead of plain
texts (see ggplot2::geom_label()). If the ggrepel package is installed,
its functions are used instead of ggplot2 native functions. When disks are
used to show aggregated flows, the names match the chosen locations: for
destination flow and attractiveness, destination locations are named, while
for production, origin locations are named (they can be both named when the
model is bipartite).
a ggplot object
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) flows <- blvim(distances, production, 1.5, 1 / 150, attractiveness, origin_data = list( names = french_cities$name[1:10], positions = positions ), destination_data = list( names = french_cities$name[1:10], positions = positions ), bipartite = FALSE ) ggplot2::autoplot(flows) ## bar plots should be almost identical if convergence occurred sim_converged(flows) ggplot2::autoplot(flows, "destination") ggplot2::autoplot(flows, "attractiveness") ## names inclusion ggplot2::autoplot(flows, "destination", with_names = TRUE) + ggplot2::coord_flip() ggplot2::autoplot(flows, with_names = TRUE) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ## positions ggplot2::autoplot(flows, "attractiveness", with_positions = TRUE) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, "destination", with_positions = TRUE, with_names = TRUE ) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, "destination", with_positions = TRUE, with_names = TRUE, with_labels = TRUE ) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, with_positions = TRUE) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, with_positions = TRUE, arrow = ggplot2::arrow(length = ggplot2::unit(0.025, "npc")) ) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326") ## individual flows combined with destination flows ## no arrows ggplot2::autoplot(flows, with_positions = TRUE, show_destination = TRUE ) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326") ## readding arrows ggplot2::autoplot(flows, with_positions = TRUE, show_destination = TRUE, arrow = ggplot2::arrow(length = ggplot2::unit(0.025, "npc")) ) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326")positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) flows <- blvim(distances, production, 1.5, 1 / 150, attractiveness, origin_data = list( names = french_cities$name[1:10], positions = positions ), destination_data = list( names = french_cities$name[1:10], positions = positions ), bipartite = FALSE ) ggplot2::autoplot(flows) ## bar plots should be almost identical if convergence occurred sim_converged(flows) ggplot2::autoplot(flows, "destination") ggplot2::autoplot(flows, "attractiveness") ## names inclusion ggplot2::autoplot(flows, "destination", with_names = TRUE) + ggplot2::coord_flip() ggplot2::autoplot(flows, with_names = TRUE) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ## positions ggplot2::autoplot(flows, "attractiveness", with_positions = TRUE) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, "destination", with_positions = TRUE, with_names = TRUE ) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, "destination", with_positions = TRUE, with_names = TRUE, with_labels = TRUE ) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, with_positions = TRUE) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(flows, with_positions = TRUE, arrow = ggplot2::arrow(length = ggplot2::unit(0.025, "npc")) ) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326") ## individual flows combined with destination flows ## no arrows ggplot2::autoplot(flows, with_positions = TRUE, show_destination = TRUE ) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326") ## readding arrows ggplot2::autoplot(flows, with_positions = TRUE, show_destination = TRUE, arrow = ggplot2::arrow(length = ggplot2::unit(0.025, "npc")) ) + ggplot2::scale_linewidth_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326")
This function uses a tile plot from ggplot2 to display a single value for each of the parameter pairs used to produce the collection of spatial interaction models.
## S3 method for class 'sim_df' autoplot(object, value, inverse = TRUE, ...)## S3 method for class 'sim_df' autoplot(object, value, inverse = TRUE, ...)
object |
a data frame of spatial interaction models, an object of class
|
value |
the value to display, default to |
inverse |
whether to use the cost scale parameter (default) |
... |
additional parameters (not used currently) |
The value to display is specified via an expression evaluated in the context
of the data frame. It defaults to the diversity as computed by diversity().
See the below for examples of use.
The horizontal axis is used by default for the cost scale parameter, that is
. This is in general easier to read than using the inverse cost
scale. The inverse parameter can be used to turn off this feature. The
vertical axis is used by default for the return to scale parameter.
a ggplot object
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) ## default display: Shannon diversity ggplot2::autoplot(all_flows_df) ## iterations ggplot2::autoplot(all_flows_df, iterations) ## we leverage non standard evaluation to compute a different diversity ggplot2::autoplot(all_flows_df, diversity(sim, "RW")) ## or to refer to columns of the data frame, either default ones ggplot2::autoplot(all_flows_df, converged) ggplot2::autoplot(all_flows_df, iterations) ggplot2::autoplot(all_flows_df, potential) ## or added ones all_flows_df["Nystuen & Dacey"] <- diversity(sim_column(all_flows_df), "ND") ggplot2::autoplot(all_flows_df, `Nystuen & Dacey`)positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) ## default display: Shannon diversity ggplot2::autoplot(all_flows_df) ## iterations ggplot2::autoplot(all_flows_df, iterations) ## we leverage non standard evaluation to compute a different diversity ggplot2::autoplot(all_flows_df, diversity(sim, "RW")) ## or to refer to columns of the data frame, either default ones ggplot2::autoplot(all_flows_df, converged) ggplot2::autoplot(all_flows_df, iterations) ggplot2::autoplot(all_flows_df, potential) ## or added ones all_flows_df["Nystuen & Dacey"] <- diversity(sim_column(all_flows_df), "ND") ggplot2::autoplot(all_flows_df, `Nystuen & Dacey`)
This function represents graphically the variability of the flows represented
by the spatial interaction models contained in a collection (a sim_list
object).
## S3 method for class 'sim_list' autoplot( object, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, qmin = 0.05, qmax = 0.95, normalisation = c("none", "origin", "full"), ... )## S3 method for class 'sim_list' autoplot( object, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, qmin = 0.05, qmax = 0.95, normalisation = c("none", "origin", "full"), ... )
object |
a collection of spatial interaction models, a |
flows |
|
with_names |
specifies whether the graphical representation includes
location names ( |
with_positions |
specifies whether the graphical representation is based
on location positions ( |
cut_off |
cut off limit for inclusion of a graphical primitive when
|
adjust_limits |
if |
with_labels |
if |
qmin |
lower quantile, see details (default: 0.05) |
qmax |
upper quantile, see details (default: 0.95) |
normalisation |
when |
... |
additional parameters, not used currently |
The graphical representation depends on the values of flows and
with_positions. It is based on the data frame representation produced by
fortify.sim_list(). In all cases, the variations of the flows are
represented via quantiles of their distribution over the collection of models
(computed with quantile.sim_list()). For instance, when flows is
"destination", the function computes the quantiles of the incoming flows
observed in the collection at each destination. We consider three quantiles:
a lower quantile qmin defaulting to 0.05;
the median;
a upper quantile qmax defaulting to 0.95.
If with_position is FALSE (default value), the graphical representations
are "abstract". Depending on flows we have the following representations:
"full": the function displays the quantiles over the collection of
models of the flows using nested squares (flows()). The graph is organised
as matrix with origin locations on rows and destination locations on columns.
At each row and column intersection, three nested squares represent
respectively the lower quantile, the median and the upper quantile of the
distribution of the flows between the corresponding origin and destination
locations over the collection of models. The median square borders are
thicker than the other two squares. The area of each square is proportional
to the represented value.
"destination": the function displays the quantiles over the collection
of models of the incoming flows for each destination location (using
destination_flow()). Quantiles are represented using
ggplot2::geom_crossbar(): each location is represented by a rectangle that
spans from its lower quantile to its upper quantile. An intermediate thicker
bar represents the median quantile.
"attractiveness": the function displays the quantiles over the
collection of models of the attractiveness of each destination location (as
given by attractiveness()). The graphical representation is the same as for
"destination". This is interesting for dynamic models where those values
are updated during the iterations (see blvim() for details). When the
calculation has converged (see sim_converged()), both "destination" and
"attractiveness" graphics should be almost identical.
When the with_names parameter is TRUE, the location names
(location_names()) are used to label the axis of the graphical
representation. If names are not specified, they are replaced by indexes.
When the with_positions parameter is TRUE, the location positions
(location_positions()) are used to produce more "geographically informed"
representations. Notice that if no positions are known for the locations, the
use of with_positions = TRUE is an error. Moreover, flows = "full" is not
supported: the function will issue a warning and revert to the position free
representation if this value is used.
The representations for flows="destination" and flows="attractiveness"
are based on the same principle. Each destination location is represented by
a collection of three nested circles centred on the corresponding location
position, representing respectively the lower quantile, the median and the
upper quantile of the incoming flows or of the attractivenesses. The
diameters of the circles are proportional to the quantities they represent.
The border ot the median circle is thicker than the ones of the other
circles.
When both with_positions and with_names are TRUE, the names of the
destinations are added to the graphical representation. If with_labels is
TRUE the names are represented as labels instead of plain texts (see
ggplot2::geom_label()). If the ggrepel package is installed, its
functions are used instead of ggplot2 native functions.
a ggplot object
fortify.sim_list(), quantile.sim_list()
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, destination_data = list( names = french_cities$name[1:10], positions = positions ), origin_data = list( names = french_cities$name[1:10], positions = positions ) ) ggplot2::autoplot(all_flows, with_names = TRUE) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ggplot2::autoplot(all_flows, with_names = TRUE, normalisation = "none") + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ggplot2::autoplot(all_flows, flow = "destination", with_names = TRUE, qmin = 0, qmax = 1 ) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ggplot2::autoplot(all_flows, flow = "destination", with_positions = TRUE, qmin = 0, qmax = 1 ) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(all_flows, flow = "destination", with_positions = TRUE, qmin = 0, qmax = 1, cut_off = 1.1 ) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(all_flows, flow = "destination", with_positions = TRUE, with_names = TRUE, with_labels = TRUE, qmin = 0, qmax = 1, cut_off = 1.1 ) + ggplot2::coord_sf(crs = "epsg:4326")positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, destination_data = list( names = french_cities$name[1:10], positions = positions ), origin_data = list( names = french_cities$name[1:10], positions = positions ) ) ggplot2::autoplot(all_flows, with_names = TRUE) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ggplot2::autoplot(all_flows, with_names = TRUE, normalisation = "none") + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ggplot2::autoplot(all_flows, flow = "destination", with_names = TRUE, qmin = 0, qmax = 1 ) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) ggplot2::autoplot(all_flows, flow = "destination", with_positions = TRUE, qmin = 0, qmax = 1 ) + ggplot2::scale_size_continuous(range = c(0, 6)) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(all_flows, flow = "destination", with_positions = TRUE, qmin = 0, qmax = 1, cut_off = 1.1 ) + ggplot2::coord_sf(crs = "epsg:4326") ggplot2::autoplot(all_flows, flow = "destination", with_positions = TRUE, with_names = TRUE, with_labels = TRUE, qmin = 0, qmax = 1, cut_off = 1.1 ) + ggplot2::coord_sf(crs = "epsg:4326")
blvim() computes equilibrium flows between origin and destination locations
in the Harris and Wilson spatial interaction model. At equilibrium,
destination attractivenesses satisfy
where is the flow from origin to destination ,
is the attractiveness of destination , and
converts attractiveness into incoming flow. Flows are
computed using the production-constrained entropy-maximising model
implemented by static_blvim().
blvim( costs, X, alpha, beta, Z = NULL, kappa = 1, bipartite = TRUE, origin_data = NULL, destination_data = NULL, epsilon = 0.01, iter_max = 50000, conv_check = 100, precision = 1e-06, algorithm = c("linear", "quadratic", "gradient", "gradient_unit"), control = blvim_gradient_control(), keep_runs = FALSE )blvim( costs, X, alpha, beta, Z = NULL, kappa = 1, bipartite = TRUE, origin_data = NULL, destination_data = NULL, epsilon = 0.01, iter_max = 50000, conv_check = 100, precision = 1e-06, algorithm = c("linear", "quadratic", "gradient", "gradient_unit"), control = blvim_gradient_control(), keep_runs = FALSE )
costs |
a cost matrix |
X |
a vector of production constraints |
alpha |
the return to scale parameter |
beta |
the inverse cost scale parameter |
Z |
Initial destination attractivenesses. This can be |
kappa |
A strictly positive numeric vector of conversion factors
between destination attractivenesses and incoming flows. A scalar is
recycled to all destinations. The default |
bipartite |
when |
origin_data |
|
destination_data |
|
epsilon |
A positive numeric scalar giving the update intensity for
the |
iter_max |
A positive integer giving the maximum number of iterations for each run. |
conv_check |
A positive integer giving the number of iterations between
convergence checks for the |
precision |
A positive numeric scalar giving the convergence tolerance.
It is used in the attractiveness-based convergence test for |
algorithm |
A character string selecting the equilibrium algorithm. It
must be one of |
control |
A list of parameters for the gradient algorithms, normally
produced by |
keep_runs |
A logical scalar controlling whether complete fitted models
are retained when multiple starting points are used. If |
Several iterative algorithms are available. The model can be run from a single initial attractiveness vector or from several starting points. In the latter case, the solution with the largest final potential is returned.
Initial destination attractivenesses are specified through Z:
If Z is NULL, one starting point is constructed automatically.
If Z is a numeric vector, the algorithm is run once from that vector.
If Z is a numeric matrix, each row is used as an independent starting
point.
If Z is an object created by blvim_initialise(), its initialiser is
evaluated using the model context. The resulting vector or matrix is then
handled in the same way as an explicitly supplied Z.
In the bipartite case, automatic initialisation uses the centre of mass of the weighted simplex:
where is the number of destinations. In the non-bipartite case, it
uses
Automatic initialisation therefore produces a single run. Use a matrix or
blvim_initialise() to perform multiple runs.
For multiple runs, all starting points use the same model parameters, algorithm, convergence settings, and algorithm-specific controls. The final potential is evaluated for each run, and the run with the strictly largest potential is selected. If several runs have the same largest potential, the first is selected.
Summaries of all runs are stored in the returned model and can be obtained
with sim_logs(). If keep_runs = TRUE, the complete fitted models are
retained and can be obtained with sim_runs().
At each iteration of a Boltzmann-Lotka-Volterra (BLV) algorithm, flows are
computed with static_blvim(). The total flow received by destination
is
Attractivenesses are then updated to bring closer to
. Two BLV algorithms are available:
"linear" uses the original update rule proposed by Harris and Wilson
(1978):
"quadratic" uses the update rule proposed by Wilson (2008):
For both algorithms, is specified by epsilon and should
generally be smaller than one. Convergence is checked every conv_check
iterations and is attained when
\lVert Z^{t+1} - Z^t \rVert
<
\delta
\left(
\lVert Z^{t+1} \rVert + \delta
\right),
where is specified by precision. Iteration stops at
convergence or after iter_max iterations.
Computing an equilibrium is equivalent to maximising the potential
function; see sim_potential(). The gradient algorithms optimise this
potential directly using projected gradient ascent.
The "gradient" algorithm works directly with the attractivenesses
. The "gradient_unit" algorithm works with rescaled variables
satisfying
The rescaling can accelerate convergence in some settings. The two
algorithms optimise equivalent formulations and should produce the same
solution up to numerical precision. They can nevertheless converge to
solutions different from those obtained with the "linear" and
"quadratic" algorithms.
Both gradient algorithms use projected ascent with a backtracking line
search and an Armijo sufficient-increase condition. Convergence is assessed
using a duality gap with tolerance precision, and the number of iterations
is limited by iter_max. The epsilon and conv_check arguments are not
used. Line-search and other internal parameters can be modified through
control; see blvim_gradient_control().
An object inheriting from sim and sim_blvim. It contains the
equilibrium flow matrix, final destination attractivenesses, initial
attractivenesses of the selected run, and convergence information.
When multiple starting points are used, the returned object represents the
run with the largest final potential. Run summaries are available through
sim_logs() and contain:
run: index of the starting point;
potential: final potential;
iterations: number of reported iterations;
converged: whether the run stopped before iter_max;
status: stopping status for a gradient algorithm;
gap: final duality gap for a gradient algorithm;
evaluations: number of evaluations for a gradient algorithm.
The last three columns are NA for "linear" and "quadratic". If
keep_runs = TRUE, the complete fitted models are additionally available
through sim_runs(), in the same order as the starting points.
While models in this package do not use location data beyond X and Z,
additional data can be stored and used when analysing spatial interaction
models.
Spatial interaction models can store names for origin and destination
locations, using origin_names<-() and destination_names<-(). Names
are taken by default from names of the cost matrix costs. More precisely,
rownames(costs) is used for origin location names and colnames(costs) for
destination location names.
Spatial interaction models can store the positions of the origin and
destination locations, using origin_positions<-() and
destination_positions<-().
In addition to the functions mentioned above, location data can be specified
directly using the origin_data and destination_data parameters. Data are
given by a list whose components are not interpreted excepted the following
ones:
names is used to specify location names and its content has to follow
the restrictions documented in origin_names<-() and
destination_names<-()
positions is used to specify location positions and its content has
to follow the restrictions documented in origin_positions<-() and
destination_positions<-()
Harris, B., & Wilson, A. G. (1978). "Equilibrium Values and Dynamics of Attractiveness Terms in Production-Constrained Spatial-Interaction Models", Environment and Planning A: Economy and Space, 10(4), 371-388. doi:10.1068/a100371
Wilson, A. (2008). "Boltzmann, Lotka and Volterra and spatial structural evolution: an integrated methodology for some dynamical systems", Journal of the Royal Society Interface, 5, 865-871. doi:10.1098/rsif.2007.1288
static_blvim() for the underlying production-constrained interaction
model, blvim_initialise() for generated starting points,
blvim_gradient_control() for gradient and line-search parameters,
sim_potential() for the potential function, and grid_blvim() for
systematic exploration of parameter values. Multiple-run summaries and
fitted models can be obtained with sim_logs() and sim_runs().
distances <- french_cities_distances[1:10, 1:10] / 1000 production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) # Rescale attractivenesses to total production attractiveness <- attractiveness / sum(attractiveness) * sum(production) # Original Harris-Wilson update rule linear_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "linear" ) linear_flows # Wilson's quadratic update rule quadratic_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "quadratic" ) # Gradient ascent on the potential gradient_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "gradient" ) # Gradient ascent using unit-simplex variables unit_gradient_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "gradient_unit" ) # Automatic single-point initialisation automatic_flows <- blvim( distances, production, 1.5, 1 / 250, algorithm = "gradient" ) # Multiple explicitly supplied starting points starts <- rbind( attractiveness, c(rep(0, 5), rep(sum(production) / 5, 5)) ) multistart_flows <- blvim( distances, production, 1.5, 1 / 250, Z = starts, algorithm = "gradient", keep_runs = TRUE ) # Generate starting points when the model context is known generated_flows <- blvim( distances, production, 1.5, 1 / 250, Z = blvim_initialise(n = 20, face_fraction = 0.5), algorithm = "gradient" )distances <- french_cities_distances[1:10, 1:10] / 1000 production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) # Rescale attractivenesses to total production attractiveness <- attractiveness / sum(attractiveness) * sum(production) # Original Harris-Wilson update rule linear_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "linear" ) linear_flows # Wilson's quadratic update rule quadratic_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "quadratic" ) # Gradient ascent on the potential gradient_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "gradient" ) # Gradient ascent using unit-simplex variables unit_gradient_flows <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "gradient_unit" ) # Automatic single-point initialisation automatic_flows <- blvim( distances, production, 1.5, 1 / 250, algorithm = "gradient" ) # Multiple explicitly supplied starting points starts <- rbind( attractiveness, c(rep(0, 5), rep(sum(production) / 5, 5)) ) multistart_flows <- blvim( distances, production, 1.5, 1 / 250, Z = starts, algorithm = "gradient", keep_runs = TRUE ) # Generate starting points when the model context is known generated_flows <- blvim( distances, production, 1.5, 1 / 250, Z = blvim_initialise(n = 20, face_fraction = 0.5), algorithm = "gradient" )
This function returns a list with parameters used to tweak the gradient
algorithms used in blvim(). These parameters control the line search
algorithm which is a standard backtracking method with Armijo stopping
criterion.
blvim_gradient_control( rho_init = 0.5, shrink = 0.5, sufficient_increase = 1e-05, ls_iter_max = 50 )blvim_gradient_control( rho_init = 0.5, shrink = 0.5, sufficient_increase = 1e-05, ls_iter_max = 50 )
rho_init |
Initial value of the gradient step in the line search. |
shrink |
Shrink factor for the gradient step in the line search. |
sufficient_increase |
Relative sufficient increase in the line search. |
ls_iter_max |
Maximum number of iterations of the line search. |
a list containing the values of the parameters
Armijo, Larry (1966), "Minimization of functions having Lipschitz continuous first partial derivatives", Pacific J. Math. 16 (1): 1–3 doi:10.2140/pjm.1966.16.1
blvim_gradient_control(0.5, 0.7, 1e-4, 25)blvim_gradient_control(0.5, 0.7, 1e-4, 25)
blvim_initialise() creates a reusable specification for generating initial
attractivenesses in blvim(). The generator is evaluated only after the
dimensions, production values, and conversion factors of the model are
known.
blvim_initialise(generator = blvim_simplex_initialise, ...)blvim_initialise(generator = blvim_simplex_initialise, ...)
generator |
A function used to generate initial attractivenesses.
Its first argument must be |
... |
Arguments stored in the specification and passed to |
A generator must be a function whose first argument is context. It must
return either:
a numeric vector representing one initial configuration;
a numeric matrix whose rows represent initial configurations; or
a list containing one of these objects in a component named Z.
The context argument is a named list containing:
costs: the cost matrix;
X: the production vector;
alpha and beta: the interaction-model parameters;
kappa: the conversion factors expanded to one value per destination;
C: total production, equal to sum(X);
bipartite: whether the model is bipartite;
n_origins: number of origins;
n_destinations: number of destinations.
Some of the parameters are easily deducible from others (C, n_origins
and n_destinations) but are included to simplify the design of
initialisers.
Additional arguments supplied through ... are stored in the specification
and passed to generator when starting points are generated.
An object of class blvim_initialise.
blvim() for the use of this function, and
blvim_simplex_initialise() for the default initialiser.
# Ten starting points from the default simplex generator initialise <- blvim_initialise(n = 10) # Emphasise proper faces of the simplex boundary_initialise <- blvim_initialise( n = 100, face_fraction = 0.8, points_per_face = 5 ) # A custom generator concentrated_initialise <- function(context, n = 20, concentration = 100) { centre <- rep(1 / context$n_destinations, context$n_destinations) U <- matrix( stats::rgamma(n * context$n_destinations, shape = concentration * centre ), nrow = n, byrow = TRUE ) U <- U / rowSums(U) sweep(U * context$C, 2, context$kappa, "/") } initialise <- blvim_initialise( generator = concentrated_initialise, n = 50, concentration = 20 )# Ten starting points from the default simplex generator initialise <- blvim_initialise(n = 10) # Emphasise proper faces of the simplex boundary_initialise <- blvim_initialise( n = 100, face_fraction = 0.8, points_per_face = 5 ) # A custom generator concentrated_initialise <- function(context, n = 20, concentration = 100) { centre <- rep(1 / context$n_destinations, context$n_destinations) U <- matrix( stats::rgamma(n * context$n_destinations, shape = concentration * centre ), nrow = n, byrow = TRUE ) U <- U / rowSums(U) sweep(U * context$C, 2, context$kappa, "/") } initialise <- blvim_initialise( generator = concentrated_initialise, n = 50, concentration = 20 )
This is the default generator used by blvim_initialise(). It samples initial
attractivenesses from the weighted simplex
blvim_simplex_initialise( context, n = 10L, face_fraction = 0.5, points_per_face = 0L, max_partial_attempts = 1000000L, include_centre = TRUE )blvim_simplex_initialise( context, n = 10L, face_fraction = 0.5, points_per_face = 0L, max_partial_attempts = 1000000L, include_centre = TRUE )
context |
A model context created internally by |
n |
Number of initial configurations to generate. The default is
deliberately modest because |
face_fraction |
Proportion of starting points generated on proper faces of the weighted simplex. |
points_per_face |
Parameter controlling the allocation of points among
simplex faces. See |
max_partial_attempts |
Maximum number of proposals used when selecting
distinct faces in a partial codimension layer. See
|
include_centre |
Whether to include the centre of the weighted simplex
as the first starting point. If |
A list containing:
Z: the matrix of initial attractivenesses;
sampling: diagnostics returned by sample_weighted_simplex(), or NULL
when no random point is requested.
distances <- french_cities_distances[1:10, 1:10] / 1000 production <- rep(1, 10) # Emphasise proper faces of the simplex boundary_initialise <- blvim_initialise( generator = blvim_simplex_initialise, n = 100, face_fraction = 0.8, points_per_face = 5 ) multistart_flows <- blvim( distances, production, 1.5, 1 / 100, Z = boundary_initialise, algorithm = "gradient", keep_runs = TRUE ) head(sim_logs(multistart_flows))distances <- french_cities_distances[1:10, 1:10] / 1000 production <- rep(1, 10) # Emphasise proper faces of the simplex boundary_initialise <- blvim_initialise( generator = blvim_simplex_initialise, n = 100, face_fraction = 0.8, points_per_face = 5 ) multistart_flows <- blvim( distances, production, 1.5, 1 / 100, Z = boundary_initialise, algorithm = "gradient", keep_runs = TRUE ) head(sim_logs(multistart_flows))
This function combines the sim_list and sim objects use as arguments
in a single sim_list, provided they are compatible. Compatibility is
defined as in sim_list(): all spatial interaction models must share
the same costs as well as the same origin and destination data.
## S3 method for class 'sim_list' c(...)## S3 method for class 'sim_list' c(...)
... |
|
A combined object of class sim_list.
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows_log <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) production <- rep(1, 15) attractiveness <- rep(1, 15) all_flows_unit <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) all_flows <- c(all_flows_log, all_flows_unit)distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows_log <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) production <- rep(1, 15) attractiveness <- rep(1, 15) all_flows_unit <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) all_flows <- c(all_flows_log, all_flows_unit)
Extract the cost matrix used to compute this model
costs(sim, ...)costs(sim, ...)
sim |
a spatial interaction model with a cost matrix |
... |
additional parameters |
the cost matrix
positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) costs(model) ## should be equal to distances above all.equal(costs(model), distances)positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) costs(model) ## should be equal to distances above all.equal(costs(model), distances)
Extract the common cost matrix from a collection of spatial interaction models
## S3 method for class 'sim_list' costs(sim, ...)## S3 method for class 'sim_list' costs(sim, ...)
sim |
a collection of spatial interaction models, an object of class
|
... |
additional parameters |
the cost matrix
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness ) ## should be TRUE identical(distances, costs(all_flows))distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness ) ## should be TRUE identical(distances, costs(all_flows))
Compute the flows incoming at each destination location
destination_flow(sim, ...)destination_flow(sim, ...)
sim |
a spatial interaction model object |
... |
additional parameters |
a vector of flows incoming at destination locations
production(), attractiveness()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) destination_flow(model) ## should be different from the attractiveness as the model is static attractiveness(model)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) destination_flow(model) ## should be different from the attractiveness as the model is static attractiveness(model)
Functions to get or set the names of the destination locations in a spatial interaction model (or in a collection of spatial interaction models).
destination_names(sim) destination_names(sim) <- valuedestination_names(sim) destination_names(sim) <- value
sim |
a spatial interaction model object (an object of class |
value |
a character vector of length equal to the number of destination
locations, or |
for destination_names NULL or a character vector with one name
per destination locations in the model. for destination_names<- the
modified sim object or sim_list object.
location_names(), origin_names()
distances <- french_cities_distances[1:10, 1:10] production <- rep(1, 10) attractiveness <- rep(1, 10) ## the row/column names of the cost matrix are used for the location model <- static_blvim(distances, production, 1.5, 1 / 250000, attractiveness) destination_names(model) destination_names(model) <- french_cities$name[1:10] destination_names(model)distances <- french_cities_distances[1:10, 1:10] production <- rep(1, 10) attractiveness <- rep(1, 10) ## the row/column names of the cost matrix are used for the location model <- static_blvim(distances, production, 1.5, 1 / 250000, attractiveness) destination_names(model) destination_names(model) <- french_cities$name[1:10] destination_names(model)
Functions to get or set the positions of the destination locations in a spatial interaction model.
destination_positions(sim) destination_positions(sim) <- valuedestination_positions(sim) destination_positions(sim) <- value
sim |
a spatial interaction model object |
value |
a matrix with as many rows as the number of destination
locations and 2 or 3 columns, or |
for destination_positions NULL or coordinate matrix for the
destination locations. for destination_positions<- the modified sim
object
Location positions are given by numeric matrices with 2 or 3 columns. The
first two columns are assumed to be geographical coordinates while the 3rd
column can be used for instance to store altitude. Coordinates are
interpreted as is in graphical representations (see autoplot.sim()). They
are not matched to the costs as those can be derived from complex movement
models and other non purely geographic considerations.
location_positions(), origin_positions()
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- static_blvim(distances, production, 1.5, 1, attractiveness) destination_positions(model) <- positions destination_positions(model)positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- static_blvim(distances, production, 1.5, 1, attractiveness) destination_positions(model) <- positions destination_positions(model)
This function computes the diversity of the destination locations according
to different definitions that all aim at estimating a number of active
destinations, i.e., the number of destination locations that receive a
"significant fraction" of the total flow or that are attractive enough. The
function applies also to a collection of spatial interaction models as
represented by a sim_list.
diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... ) ## S3 method for class 'sim' diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... ) ## S3 method for class 'sim_list' diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... )diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... ) ## S3 method for class 'sim' diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... ) ## S3 method for class 'sim_list' diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... )
sim |
a spatial interaction model object (an object of class |
definition |
diversity definition |
order |
order of the Rényi entropy, used only when |
activity |
specifies whether the diversity is computed based on the
destination flows (for |
... |
additional parameters |
In general, the activity of a destination location is measured by its
incoming flow a.k.a. its destination flow. If is a flow matrix, the
destination flows are computed as follows
for each destination (see destination_flow()). This is the default
calculation mode in this function (when the parameter activity is set
to "destination").
For dynamic models produced by blvim(), the destination attractivenesses
can be also considered as activity measures. When convergence occurs,
the values are identical, but prior convergence they can be quite different.
When activity is set to "attractiveness", the diversity measures are
computed using the same formula as below but with replaced by
(as given by attractiveness()).
To compute their diversity using entropy based definitions, the activities are first normalised to be interpreted as a probability distribution over the destination locations. For instance for destination flows, we use
The most classic diversity index is given by the exponential of Shannon's
entropy (parameter definition="shannon"). This gives
Rényi generalized entropy can be used to define a collection of other
diversity metrics. The Rényi diversity of order is the
exponential of the Rényi entropy of order of the
distribution, that is
This is defined directly only for , but
extensions to the limit case are straightforward:
is Shannon's entropy/diversity
is the max-entropy, here and thus the
corresponding diversity is the number of locations
is the min-entropy, here and
thhe corresponding diversity is
The definition parameter specifies the diversity used for calculation. The
default value is shannon for Shannon's entropy (in this case the order
parameter is not used). Using renyi gives access to all Rényi diversities
as specified by the order parameter. Large values of order tend to
generate underflows in the calculation that will trigger the use of the
min-entropy instead of the exact Rényi entropy.
In addition to those entropy based definition, terminal based calculations
are also provided. Using any definition supported by the terminals()
function, the diversity is the number of terminals identified. Notice this
applies only to interaction models in which origin and destination locations
are identical, i.e. when the model is not bipartite. In addition, the notion
of terminals is based on destination flows and cannot be used with activities
based on attractivenesses. definition can be:
"ND" for the original Nystuen and Dacey definition
"RW" for the variant by Rihll and Wilson
See terminals() for details.
When applied to a collection of spatial interaction models (an object of
class sim_list) the function uses the same parameters (definition and
order) for all models and returns a vector of diversities. This is
completely equivalent to grid_diversity().
the diversity of destination flows (one value per spatial interaction model)
Jost, L. (2006), "Entropy and diversity", Oikos, 113: 363-375. doi:10.1111/j.2006.0030-1299.14714.x
destination_flow(), attractiveness(), terminals(),
sim_is_bipartite()
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- rep(1, 15) flows <- blvim(distances, production, 1.5, 1 / 100, attractiveness, bipartite = FALSE ) diversity(flows) sim_converged(flows) ## should be identical because of convergence diversity(flows, activity = "attractiveness") diversity(flows, "renyi", 2) diversity(flows, "RW")distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- rep(1, 15) flows <- blvim(distances, production, 1.5, 1 / 100, attractiveness, bipartite = FALSE ) diversity(flows) sim_converged(flows) ## should be identical because of convergence diversity(flows, activity = "attractiveness") diversity(flows, "renyi", 2) diversity(flows, "RW")
Extract the flow matrix from a spatial interaction model object
flows(sim, ...)flows(sim, ...)
sim |
a spatial interaction model object |
... |
additional parameters |
a matrix of flows between origin locations and destination locations
flows_df() for a data frame version of the flows,
destination_flow() for destination flows.
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) ## rescale to production attractiveness <- attractiveness / sum(attractiveness) * sum(production) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness) flows(model)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) ## rescale to production attractiveness <- attractiveness / sum(attractiveness) * sum(production) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness) flows(model)
Extract the flow matrix from a spatial interaction model object in data frame format
flows_df(sim, ...)flows_df(sim, ...)
sim |
a spatial interaction model object |
... |
additional parameters (not used currently) |
This function extracts the flow matrix in a long format. Each row contains the flow between an origin location and a destination location. The resulting data frame has at least three columns:
origin_idx: identifies the origin location by its index from 1 to the number
of origin locations
destination_idx: identifies the destination location by its index from 1
to the number of destination locations
flow: the flow between the corresponding location
In addition, if location information is available, it will be included in the data frame as follows:
location names are included using columns origin_name or destination_name
positions are included using 2 or 3 columns (per location type, origin or
destination) depending on the number of dimensions used for the location.
The names of the columns are by default origin_x, origin_y and
origin_z ( and equivalent names for destination location) unless
coordinate names are specified in the location positions. In this latter
case, the names are prefixed by origin_ or destination_. For instance,
if the destination location position coordinates are named "longitude"
and "latitude", the resulting columns will be destination_longitude and
destination_latitude.
a data frame of flows between origin locations and destination locations with additional content if available (see Details).
location_positions(), location_names(), flows()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) ## rescale to production attractiveness <- attractiveness / sum(attractiveness) * sum(production) ## simple case (no positions and default names) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness) head(flows_df(model)) ## with location data positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness, origin_data = list(positions = positions), destination_data = list(positions = positions) ) head(flows_df(model)) ## with names origin_names(model) <- french_cities$name[1:10] destination_names(model) <- french_cities$name[1:10] head(flows_df(model))distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) ## rescale to production attractiveness <- attractiveness / sum(attractiveness) * sum(production) ## simple case (no positions and default names) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness) head(flows_df(model)) ## with location data positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness, origin_data = list(positions = positions), destination_data = list(positions = positions) ) head(flows_df(model)) ## with names origin_names(model) <- french_cities$name[1:10] destination_names(model) <- french_cities$name[1:10] head(flows_df(model))
This function extracts from a spatial interaction model different types of
data frame that can be used to produce graphical representations.
autoplot.sim() leverages this function to produce its graphical
representations.
## S3 method for class 'sim' fortify( model, data, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, cut_off = 100 * .Machine$double.eps^0.5, ... )## S3 method for class 'sim' fortify( model, data, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, cut_off = 100 * .Machine$double.eps^0.5, ... )
model |
a spatial interaction model object |
data |
not used |
flows |
|
with_names |
specifies whether the extracted data frame includes
location names ( |
with_positions |
specifies whether the extracted data frame is based on
location positions ( |
cut_off |
cut off limit for inclusion of a flow row in the final data frame. |
... |
additional parameters, not used currently |
The data frame produced by the method depends on the values of flows and
with_positions. The general principal is to have one row per flow, either a
single flow from an origin location to a destination location, or an
aggregated flow to a destination location. Flows are stored in one column of
the data frame, while the other columns are used to identify origin and
destination.
If with_position is FALSE (default value), data frames are simple.
Depending on flows, the function extracts different data frames:
"full": this is the default case for which the full flow matrix is extracted.
The data frame has three variables:
origin: identifies the origin location by its index from 1 to the number
of origin locations
destination: identifies the destination location by its index from 1
to the number of destination locations
flow: the flow between the corresponding location
It is recommend to use flows_df() for more control over the extraction
outside of simple graphical representations.
"destination": the data frame has only two or three columns:
destination: identifies the destination location by its index from 1
to the number of destination locations
flow: the incoming flows (see destination_flow())
name: the name of the destination location if with_names is TRUE
"attractiveness": the data frame has also two ot three columns,
destination and name as in the previous case and attractiveness which
contains the attractivenesses of the destinations (see attractiveness()).
When the with_positions parameter is TRUE, the location positions
(location_positions()) are used to produce more "geographically informed"
extractions. Notice that if no positions are known for the locations, the use
of with_positions = TRUE is an error. Depending on flows we have the
following representations:
"full": this is the default case for which the full flow matrix is extracted.
Positions for both origin and destination locations are needed. The data
frame contains five columns:
the first two columns are used for the coordinates of the origin locations (see below for the names of the columns)
xend and yend are used for the coordinates of the destination locations
flow is used for the flows
"destination" and "attractiveness" produce both a data frame with three
or four columns. As when with_positions is FALSE, one column is dedicated either
to the incoming flows (destination_flow()) for flows="destination" (the
name of the column is destination) or to the attractivenesses
(attractiveness()), in which case its name is attractiveness. The other
two columns are used for the positions of the destination locations. Their
names are the names of the columns of the positions
(colnames(destination_location(object))) or "x" and "y", when such
names are not specified. If with_names is TRUE, a name column is included
and contains the names of the destination locations.
In the position based data frames, rows are excluded from the returned data
frames when the flow they represent are small, i.e. when they are smaller
than the cut_off value.
a data frame, see details
positions <- matrix(rnorm(10 * 2), ncol = 2) colnames(positions) <- c("X", "Y") distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) flows <- blvim(distances, production, 1.5, 4, attractiveness, origin_data = list(names = LETTERS[1:10], positions = positions), destination_data = list(names = LETTERS[1:10], positions = positions) ) ggplot2::fortify(flows) ggplot2::fortify(flows, flows = "destination") ggplot2::fortify(flows, flows = "attractiveness") ## positions ggplot2::fortify(flows, flows = "attractiveness", with_positions = TRUE) ## names and positions ggplot2::fortify(flows, flows = "destination", with_positions = TRUE, with_names = TRUE ) ggplot2::fortify(flows, with_positions = TRUE, cut_off = 0.1)positions <- matrix(rnorm(10 * 2), ncol = 2) colnames(positions) <- c("X", "Y") distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) flows <- blvim(distances, production, 1.5, 4, attractiveness, origin_data = list(names = LETTERS[1:10], positions = positions), destination_data = list(names = LETTERS[1:10], positions = positions) ) ggplot2::fortify(flows) ggplot2::fortify(flows, flows = "destination") ggplot2::fortify(flows, flows = "attractiveness") ## positions ggplot2::fortify(flows, flows = "attractiveness", with_positions = TRUE) ## names and positions ggplot2::fortify(flows, flows = "destination", with_positions = TRUE, with_names = TRUE ) ggplot2::fortify(flows, with_positions = TRUE, cut_off = 0.1)
This function extracts from a collection of spatial interaction models
(represented by a sim_list) a data frame in a long format, with one flow
per row. This can be seen a collection oriented version of fortify.sim().
The resulting data frame is used by autoplot.sim_list() to produce summary
graphics.
## S3 method for class 'sim_list' fortify( model, data, flows = c("full", "destination", "attractiveness"), with_names = FALSE, normalisation = c("none", "origin", "full"), ... )## S3 method for class 'sim_list' fortify( model, data, flows = c("full", "destination", "attractiveness"), with_names = FALSE, normalisation = c("none", "origin", "full"), ... )
model |
a collection of spatial interaction models, a |
data |
not used |
flows |
|
with_names |
specifies whether the extracted data frame includes
location names ( |
normalisation |
when |
... |
additional parameters, not used currently |
The data frame produced by the method depends on the values of flows and to
a lesser extent on the value of with_names. In all cases, the data frame
has a configuration column that identify from which spatial interaction
model the other values have been extracted: this is the index of the model in
the original sim_list. Depending on flows we have the following
representations:
if flows="full": this is the default case for which the full flow matrix
of each spatial interaction model is extracted. The data frame contains 4
columns:
origin_idx: identifies the origin location by its index from 1 to the number
of origin locations
destination_idx: identifies the destination location by its index from 1
to the number of destination locations
flow: the flow between the corresponding location. By default, flows
are normalised by origin location (when normalisation="origin"): the total
flows originating from each origin location is normalised to 1. If
normalisation="full", this normalisation is global: the sum of all flows in
each model is normalised to 1. If normalisation="none" flows are not
normalised.
configuration: the spatial interaction model index
if flows="destination" or flows="attractiveness", the data frame contains
3 or 4 columns:
destination: identifies the destination location by its index from 1
to the number of destination locations
flow or attractiveness depending on the value of "flows": this contains
either the destination_flow() or the attractiveness() of the destination
location
configuration: the spatial interaction model index
name: the destination location names if with_names=TRUE (the column is
not present if with_names=FALSE)
The normalisation operated when flows="full" can improve the readability
of the graphical representation proposed in autoplot.sim_list() when the
production constraints differ significantly from one origin location to another.
a data frame, see details
positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) flows_1 <- blvim(distances, production, 1.5, 1, attractiveness) flows_2 <- blvim(distances, production, 1.25, 2, attractiveness) all_flows <- sim_list(list(flows_1, flows_2)) ggplot2::fortify(all_flows) ## somewhat similar to a row bind of sim_df results ggplot2::fortify(all_flows, flows = "destination") destination_names(all_flows) <- letters[1:10] ggplot2::fortify(all_flows, flows = "attractiveness", with_names = TRUE)positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) flows_1 <- blvim(distances, production, 1.5, 1, attractiveness) flows_2 <- blvim(distances, production, 1.25, 2, attractiveness) all_flows <- sim_list(list(flows_1, flows_2)) ggplot2::fortify(all_flows) ## somewhat similar to a row bind of sim_df results ggplot2::fortify(all_flows, flows = "destination") destination_names(all_flows) <- letters[1:10] ggplot2::fortify(all_flows, flows = "attractiveness", with_names = TRUE)
French cities with 50,000 inhabitants or more.
french_citiesfrench_cities
A data.frame with 121 rows and 9 columns
The INSEE code of the city according to the official geographical code of 2025 (OGC)
The name of the city
The code of the department of the city in the OGC (see french_departments)
The code of the region of the city in the OGC (see french_regions)
The population of the city in 2022
Area of the city in squared kilometers
The latitude of the of town hall the city (epsg:4326)
The longitude of the town hall of the city (epsg:4326)
The latitude of the centre of the city (epsg:4326)
The longitude of the centre of the city (epsg:4326)
This data set describes all Metropolitan France cities with 50,000 or more inhabitants in 2022, excluding Corsican cities. It contains 121 cities described by 8 variables. The data frame is sorted in decreasing population order, making it straightforward to select the most populated cities. The same order is used for rows and columns in distance matrices french_cities_distances and french_cities_times.
The population and administrative information was collected from the INSEE open data in November 2025. These data are distributed under the French "Open Licence". Geographical coordinates and areas have been obtained from the Geo API in November 2025 and are also available under the French "Open Licence".
INSEE Population census - Main extraction (2022) https://catalogue-donnees.insee.fr/en/catalogue/recherche/DS_RP_POPULATION_PRINC
INSEE Official Geographical Code (2025) https://www.data.gouv.fr/datasets/code-officiel-geographique-cog/ https://www.data.gouv.fr/api/1/datasets/r/91a95bee-c7c8-45f9-a8aa-f14cc4697545
Geo API (2025) https://geo.api.gouv.fr/
french_departments, french_regions, french_cities_distances and french_cities_times
Distances in meters and in minutes between the French cities with at least 50,000 inhabitants (french_cities).
french_cities_distances french_cities_timesfrench_cities_distances french_cities_times
matrices with 121 rows and 121 columns
An object of class matrix (inherits from array) with 121 rows and 121 columns.
Both data sets are square matrices of dimension (121, 121). If D is one of
the matrix, D[i,j] is the distance from city of id rownames(D)[i] to city
id colnames(D)[j] expressed in meters (french_cities_distances) or in
minutes (french_cities_times). The distance is measured along the fastest
path from i to j on the French road networks as computed using
OSRM engine (see details below). Ids in
column and row names can be used to get information on the cities in the
french_cities data set (column id). Rows and columns are sorted in
decreasing population order, as in french_cities. Note that the matrices
are not symmetric.
The distances and durations have been computed using the OSRM engine in version 6.0.0.
Calculations are based on the car profile included in OSRM sources for
the 6.0.0 version.
France roads are provided by OpenStreetMap under the Open Data Commons Open Database License (ODbL) using the GeoFabrik export dated 2025-11-07T21:20:50Z was used.
Geo API (2025) https://geo.api.gouv.fr/
OpenStreetMap https://www.openstreetmap.org
GeoFabrik https://download.geofabrik.de/europe/france.html
OSRM https://project-osrm.org/
This data set contains the list of all French departments. It can be joined with the french_regions or the french_cities data set. The data set was extracted from the INSEE open data in November 2025. These data are distributed under the French "Open Licence".
french_departmentsfrench_departments
A data.frame with 101 rows and 3 columns
The code of the department in the official geographical code of 2025 (OGC)
The code of the region of the city in the OGC (see french_regions)
The name of the department
INSEE Official Geographical Code (2025) https://www.data.gouv.fr/datasets/code-officiel-geographique-cog/ https://www.data.gouv.fr/api/1/datasets/r/91a95bee-c7c8-45f9-a8aa-f14cc4697545
french_cities and french_regions
This data set contains the list of all French regions It can be joined with the french_departments or the french_cities data set. The data set was extracted from the INSEE open data in November 2025. These data are distributed under the French "Open Licence".
french_regionsfrench_regions
A data.frame with 18 rows and 2 columns
The code of the region in the official geographical code of 2025 (OGC)
The name of the region
INSEE Official Geographical Code (2025) https://www.data.gouv.fr/datasets/code-officiel-geographique-cog/ https://www.data.gouv.fr/api/1/datasets/r/91a95bee-c7c8-45f9-a8aa-f14cc4697545
french_departments and french_cities
The function extract attractivenesses from all the spatial interaction models of the collection and returns them in a matrix in which each row corresponds to a model and each column to a destination location.
grid_attractiveness(sim_list, ...)grid_attractiveness(sim_list, ...)
sim_list |
a collection of spatial interaction models, an object of
class |
... |
additional parameters for the |
a matrix of attractivenesses at the destination locations
attractiveness() and grid_blvim()
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1 ) g_Z <- grid_attractiveness(all_flows) ## should be 12 rows (3 times 4 parameter pairs) and 15 columns (15 ## destination locations) dim(g_Z)distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1 ) g_Z <- grid_attractiveness(all_flows) ## should be 12 rows (3 times 4 parameter pairs) and 15 columns (15 ## destination locations) dim(g_Z)
This function combines spatial interaction model representations similar to
the ones produced by autoplot.sim() into a single ggplot. It provides an
alternative graphical representation to the one produced by
autoplot.sim_df() for collection of spatial interaction models in a
sim_df object.
grid_autoplot( sim_df, key, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, show_destination = FALSE, show_attractiveness = FALSE, show_production = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, max_sims = 25, fw_params = NULL, ... )grid_autoplot( sim_df, key, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, show_destination = FALSE, show_attractiveness = FALSE, show_production = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, max_sims = 25, fw_params = NULL, ... )
sim_df |
a data frame of spatial interaction models, an object of class
|
key |
the wrapping variable which acts as an identifier for spatial interaction models |
flows |
|
with_names |
specifies whether the graphical representation includes
location names ( |
with_positions |
specifies whether the graphical representation is based
on location positions ( |
show_destination |
specifies whether the position based |
show_attractiveness |
specifies whether the position based |
show_production |
specifies whether the position based |
cut_off |
cut off limit for inclusion of a graphical primitive when
|
adjust_limits |
if |
with_labels |
if |
max_sims |
the maximum number of spatial interaction models allowed in
the |
fw_params |
parameters for the ggplot2::facet_wrap call (if non
|
... |
additional (named) parameters passed to |
The rationale of autoplot.sim_df() is to display a single value for each
spatial interaction model (SIM) in the sim_df data frame. On the contrary,
this function produces a full graphical representation of each SIM. It is
therefore limited to small collection of SIMs (as specified by the max_sims
parameter which default to 25).
Under the hood, the function uses fortify.sim() and shares code with
autoplot.sim() to have identical representations. It is simply based on
facet wrapping facility of ggplot2. In particular the key parameter is used
as the wrapping variable in the call to ggplot2::facet_wrap(). If not
specified, the function generates an id variable which ranges from 1 to the
number of SIMs in the sim_df data frame. If specified, it is evaluated in
the context of the data frame and used for wrapping. Notice that if the
expression evaluates to identical values for different SIMs, they will be
drawn on the same panel of the final figure, which may end up with
meaningless representations. Parameters of ggplot2::facet_wrap() can be set
using the fw_params parameter (in a list).
a ggplot object
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.35, by = 0.1), seq(1, 2.5, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, destination_data = list( names = french_cities$name[1:10], positions = positions ), origin_data = list( names = french_cities$name[1:10], positions = positions ) ) all_flows_df <- sim_df(all_flows) ## default display: flows as matrices grid_autoplot(all_flows_df) ## custom wrapping variable grid_autoplot(all_flows_df, paste(alpha, "~", beta)) ## bar plots grid_autoplot(all_flows_df, flows = "attractiveness") ## with positions grid_autoplot(all_flows_df, with_positions = TRUE, show_destination = TRUE) + ggplot2::scale_linewidth_continuous(range = c(0, 1)) + ggplot2::scale_size_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326")positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.35, by = 0.1), seq(1, 2.5, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, destination_data = list( names = french_cities$name[1:10], positions = positions ), origin_data = list( names = french_cities$name[1:10], positions = positions ) ) all_flows_df <- sim_df(all_flows) ## default display: flows as matrices grid_autoplot(all_flows_df) ## custom wrapping variable grid_autoplot(all_flows_df, paste(alpha, "~", beta)) ## bar plots grid_autoplot(all_flows_df, flows = "attractiveness") ## with positions grid_autoplot(all_flows_df, with_positions = TRUE, show_destination = TRUE) + ggplot2::scale_linewidth_continuous(range = c(0, 1)) + ggplot2::scale_size_continuous(range = c(0, 2)) + ggplot2::coord_sf(crs = "epsg:4326")
This function computes a collection of flows between origin locations and
destination locations using blvim() on a grid of parameters. The flows use
the same costs, same production constraints and same attractivenesses. Each
flow is computed using one of all the pairwise combinations between the alpha
values given by alphas and the beta values given by betas. The function
returns an object of class sim_list which contains the resulting flows.
grid_blvim( costs, X, alphas, betas, Z = NULL, kappa = 1, bipartite = TRUE, origin_data = NULL, destination_data = NULL, epsilon = 0.01, iter_max = 50000, conv_check = 100, precision = 1e-06, algorithm = c("linear", "quadratic", "gradient", "gradient_unit"), control = blvim_gradient_control(), progress = FALSE )grid_blvim( costs, X, alphas, betas, Z = NULL, kappa = 1, bipartite = TRUE, origin_data = NULL, destination_data = NULL, epsilon = 0.01, iter_max = 50000, conv_check = 100, precision = 1e-06, algorithm = c("linear", "quadratic", "gradient", "gradient_unit"), control = blvim_gradient_control(), progress = FALSE )
costs |
a cost matrix |
X |
a vector of production constraints |
alphas |
a vector of return to scale parameters |
betas |
a vector of cost inverse scale parameters |
Z |
Initial destination attractivenesses. This can be |
kappa |
A strictly positive numeric vector of conversion factors
between destination attractivenesses and incoming flows. A scalar is
recycled to all destinations. The default |
bipartite |
when |
origin_data |
|
destination_data |
|
epsilon |
A positive numeric scalar giving the update intensity for
the |
iter_max |
A positive integer giving the maximum number of iterations for each run. |
conv_check |
A positive integer giving the number of iterations between
convergence checks for the |
precision |
A positive numeric scalar giving the convergence tolerance.
It is used in the attractiveness-based convergence test for |
algorithm |
A character string selecting the equilibrium algorithm. It
must be one of |
control |
A list of parameters for the gradient algorithms, normally
produced by |
progress |
if TRUE, a progress bar is shown during the calculation (defaults to FALSE) |
an object of class sim_list
While models in this package do not use location data beyond X and Z,
additional data can be stored and used when analysing spatial interaction
models.
Spatial interaction models can store names for origin and destination
locations, using origin_names<-() and destination_names<-(). Names
are taken by default from names of the cost matrix costs. More precisely,
rownames(costs) is used for origin location names and colnames(costs) for
destination location names.
Spatial interaction models can store the positions of the origin and
destination locations, using origin_positions<-() and
destination_positions<-().
In addition to the functions mentioned above, location data can be specified
directly using the origin_data and destination_data parameters. Data are
given by a list whose components are not interpreted excepted the following
ones:
names is used to specify location names and its content has to follow
the restrictions documented in origin_names<-() and
destination_names<-()
positions is used to specify location positions and its content has
to follow the restrictions documented in origin_positions<-() and
destination_positions<-()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness ) all_flows length(all_flows) all_flows[[2]]distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness ) all_flows length(all_flows) all_flows[[2]]
The function extract destination flows from all the spatial interaction models of the collection and returns them in a matrix in which each row corresponds to a model and each column to a destination location.
grid_destination_flow(sim_list, ...)grid_destination_flow(sim_list, ...)
sim_list |
a collection of spatial interaction models, an object of
class |
... |
additional parameters for the |
a matrix of destination flows at the destination locations
destination_flow() and grid_blvim()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1 ) g_df <- grid_destination_flow(all_flows) ## should be 12 rows (3 times 4 parameter pairs) and 10 columns (10 ## destination locations) dim(g_df)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1 ) g_df <- grid_destination_flow(all_flows) ## should be 12 rows (3 times 4 parameter pairs) and 10 columns (10 ## destination locations) dim(g_df)
The function computes for each spatial interaction model of its sim_list
parameter the diversity() of the corresponding destination flows and
returns the values as a vector. The type of diversity and the associated
parameters are identical for all models.
grid_diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... )grid_diversity( sim, definition = c("shannon", "renyi", "ND", "RW"), order = 1L, activity = c("destination", "attractiveness"), ... )
sim |
a collection of spatial interaction models, an object of class
|
definition |
diversity definition |
order |
order of the Rényi entropy, used only when |
activity |
specifies whether the diversity is computed based on the
destination flows (for |
... |
additional parameters |
See diversity() for the definition of the diversities. Notice that
diversity() is generic and can be applied directly to sim_list objects.
The current function is provided to be explicit in R code about what is a
unique model and what is a collection of models (using function names that
start with "grid_")
a vector of diversities, one per spatial interaction model
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE ) diversities <- grid_diversity(all_flows) diversities ## should be a length 12 vector grid_diversity(all_flows, "renyi", 3)distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE ) diversities <- grid_diversity(all_flows) diversities ## should be a length 12 vector grid_diversity(all_flows, "renyi", 3)
The function extract terminal status from all the spatial interaction models
of the collection and returns them in a matrix in which each row corresponds
to a model and each column to a destination location. The value at row i
and column j is TRUE if destination j is a terminal in model i. This
function applies only to non bipartite models.
grid_is_terminal(sim_list, definition = c("ND", "RW"), ...)grid_is_terminal(sim_list, definition = c("ND", "RW"), ...)
sim_list |
a collection of non bipartite spatial interaction models, an
object of class |
definition |
terminal definition, either |
... |
additional parameters for the |
See terminals() for the definition of terminal locations.
a matrix of terminal status at the destination locations
is_terminal() and grid_blvim()
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE ) g_df <- grid_is_terminal(all_flows) ## should be 12 rows (3 times 4 parameter pairs) and 15 columns (15 ## destination locations) dim(g_df)distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE ) g_df <- grid_is_terminal(all_flows) ## should be 12 rows (3 times 4 parameter pairs) and 15 columns (15 ## destination locations) dim(g_df)
The function reports for each spatial interaction model of its sim_list
parameter its convergence status, as defined in sim_converged().
grid_sim_converged(sim, ...)grid_sim_converged(sim, ...)
sim |
a collection of spatial interaction models, an object of class
|
... |
additional parameters |
Notice that sim_converged() is generic and can be applied directly to
sim_list objects. The current function is provided to be explicit in R code
about what is a unique model and what is a collection of models (using
function names that start with "grid_")
a vector of convergence status, one per spatial interaction model
sim_converged(), grid_sim_iterations() and grid_blvim()
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) grid_sim_converged(all_flows)distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) grid_sim_converged(all_flows)
The function reports for each spatial interaction model of its sim_list
parameter the number of iterations used to produce it (see
sim_iterations())
grid_sim_iterations(sim, ...)grid_sim_iterations(sim, ...)
sim |
a collection of spatial interaction models, an object of class
|
... |
additional parameters |
Notice that sim_iterations() is generic and can be applied directly to
sim_list objects. The current function is provided to be explicit in R code
about what is a unique model and what is a collection of models (using
function names that start with "grid_")
a vector of numbers of iteration, one per spatial interaction model
sim_iterations(), grid_sim_converged() and grid_blvim()
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) grid_sim_iterations(all_flows)distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) grid_sim_iterations(all_flows)
The function reports the potentia of each spatial interaction model of its
sim_list parameter (see sim_potential())
grid_sim_potential(sim, ...)grid_sim_potential(sim, ...)
sim |
a collection of spatial interaction models, an object of class
|
... |
additional parameters |
Notice that sim_potential() is generic and can be applied directly to
sim_list objects. The current function is provided to be explicit in R code
about what is a unique model and what is a collection of models (using
function names that start with "grid_")
a vector of numbers of iteration, one per spatial interaction model
sim_potential() and grid_blvim()
distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) grid_sim_potential(all_flows)distances <- french_cities_distances[1:15, 1:15] / 1000 ## convert to km production <- log(french_cities$population[1:15]) attractiveness <- log(french_cities$area[1:15]) all_flows <- grid_blvim( distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1, bipartite = FALSE, iter_max = 750 ) grid_sim_potential(all_flows)
This function combines spatial variability interaction model representations
similar to the ones produced by autoplot.sim_list() into a single ggplot.
It provides an alternative graphical representation to the one produced by
autoplot.sim_df() and by grid_autoplot() for collection of spatial
interaction models in a sim_df object.
grid_var_autoplot( sim_df, key, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, qmin = 0.05, qmax = 0.95, normalisation = c("origin", "full", "none"), fw_params = NULL, ... )grid_var_autoplot( sim_df, key, flows = c("full", "destination", "attractiveness"), with_names = FALSE, with_positions = FALSE, cut_off = 100 * .Machine$double.eps^0.5, adjust_limits = FALSE, with_labels = FALSE, qmin = 0.05, qmax = 0.95, normalisation = c("origin", "full", "none"), fw_params = NULL, ... )
sim_df |
a data frame of spatial interaction models, an object of class
|
key |
the wrapping variable which acts as group identifier for spatial interaction models |
flows |
|
with_names |
specifies whether the graphical representation includes
location names ( |
with_positions |
specifies whether the graphical representation is based
on location positions ( |
cut_off |
cut off limit for inclusion of a graphical primitive when
|
adjust_limits |
if |
with_labels |
if |
qmin |
lower quantile, see details (default: 0.05) |
qmax |
upper quantile, see details (default: 0.95) |
normalisation |
when |
fw_params |
parameters for the ggplot2::facet_wrap call (if non
|
... |
additional parameters passed to |
The rationale of autoplot.sim_df() is to display a single value for each
spatial interaction model (SIM) in the sim_df data frame. On the contrary,
this function produces a graphical representation of the variability of a
partition of the SIMs in the data frame, using autoplot.sim_list() as the
graphical engine.
The key parameter is used to partition the collection of SIMs. It can be
any expression which can be evaluated in the context of the sim_df
parameter. The function uses this parameter as the wrapping variable in a
call to ggplot2::facet_wrap(). It also uses it as a way to specific a
partition of the SIMs: each panel of the final figure is essentially the
variability graph generated by autoplot.sim_list() for the subset of the
SIMs in sim_df that match the corresponding value of key.
Parameters of ggplot2::facet_wrap() can be set using the fw_params
parameter (in a list).
a ggplot object
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, destination_data = list( names = french_cities$name[1:10], positions = positions ), origin_data = list( names = french_cities$name[1:10], positions = positions ) ) all_flows_df <- sim_df(all_flows) ## group models by iteration number grid_var_autoplot(all_flows_df, iterations) ## or by convergence status (showing destination) grid_var_autoplot(all_flows_df, converged, flow = "destination", with_names = TRUE ) + ggplot2::coord_flip() ## using positions grid_var_autoplot(all_flows_df, iterations, flow = "destination", with_positions = TRUE ) + ggplot2::scale_size_continuous(range = c(0, 3)) + ggplot2::coord_sf(crs = "epsg:4326")positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, destination_data = list( names = french_cities$name[1:10], positions = positions ), origin_data = list( names = french_cities$name[1:10], positions = positions ) ) all_flows_df <- sim_df(all_flows) ## group models by iteration number grid_var_autoplot(all_flows_df, iterations) ## or by convergence status (showing destination) grid_var_autoplot(all_flows_df, converged, flow = "destination", with_names = TRUE ) + ggplot2::coord_flip() ## using positions grid_var_autoplot(all_flows_df, iterations, flow = "destination", with_positions = TRUE ) + ggplot2::scale_size_continuous(range = c(0, 3)) + ggplot2::coord_sf(crs = "epsg:4326")
Extracts the destination attractivenesses used to initialise the construction
of a spatial interaction model. For static models, this is identical to their
attractiveness(). For a spatial interaction model produced by an iterative
algorithm, this function returns the attractivenesses used to initialise the
algorithm. attractiveness() returns the final attractivenesses.
initial_attractiveness(sim, ...)initial_attractiveness(sim, ...)
sim |
a spatial interaction model object |
... |
additional parameters |
A vector containing the initial destination attractivenesses.
attractiveness(), static_blvim()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) initial_attractiveness(model) identical( initial_attractiveness(model), attractiveness(model) )distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) initial_attractiveness(model) identical( initial_attractiveness(model), attractiveness(model) )
Extract the inverse cost scale parameter used to compute this model
inverse_cost(sim, ...)inverse_cost(sim, ...)
sim |
a spatial interaction model with a inverse cost scale parameter |
... |
additional parameters |
the inverse cost scale parameter
positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) inverse_cost(model) ## should be 1positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) inverse_cost(model) ## should be 1
This function returns a logical vector whose length equals the number of
locations. The value in position i is TRUE if location number i is
a terminal and FALSE if it is not. For the definition of terminals,
see terminals().
is_terminal(sim, definition = c("ND", "RW"), ...)is_terminal(sim, definition = c("ND", "RW"), ...)
sim |
a spatial interaction model object |
definition |
terminal definition, either |
... |
additional parameters |
a logical vector with TRUE at the positions of locations that are
terminals and FALSE for other locations.
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- blvim(distances, production, 1.3, 1 / 500, attractiveness, bipartite = FALSE ) destination_names(model) <- french_cities$name[1:10] is_terminal(model) dist_times <- french_cities_times[1:10, 1:10] tmodel <- blvim(dist_times, production, 1.3, 1 / 10000, attractiveness, bipartite = FALSE ) destination_names(tmodel) <- french_cities$name[1:10] is_terminal(tmodel)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- blvim(distances, production, 1.3, 1 / 500, attractiveness, bipartite = FALSE ) destination_names(model) <- french_cities$name[1:10] is_terminal(model) dist_times <- french_cities_times[1:10, 1:10] tmodel <- blvim(dist_times, production, 1.3, 1 / 10000, attractiveness, bipartite = FALSE ) destination_names(tmodel) <- french_cities$name[1:10] is_terminal(tmodel)
Those functions provide low level access to origin and destination local
names. It is recommended to use origin_names() and destination_names()
instead of location_names and location_names<-.
location_names(sim) location_names(sim) <- valuelocation_names(sim) location_names(sim) <- value
sim |
a spatial interaction model object (an object of class |
value |
a list with two components (see the returned value) or |
for location_names NULL or a list with two components: origin
for the origin location names and destination for the destination
location names. For location_names<-() the modified sim object or
sim_list object.
origin_names(), destination_names()
distances <- french_cities_distances[1:10, 1:10] production <- rep(1, 10) attractiveness <- rep(1, 10) ## the row/column names of the cost matrix are used for the location model <- static_blvim(distances, production, 1.5, 1 / 250000, attractiveness) location_names(model) location_names(model) <- NULL location_names(model) location_names(model) <- list( origin = french_cities$name[1:10], destination = LETTERS[1:10] ) destination_names(model) origin_names(model)distances <- french_cities_distances[1:10, 1:10] production <- rep(1, 10) attractiveness <- rep(1, 10) ## the row/column names of the cost matrix are used for the location model <- static_blvim(distances, production, 1.5, 1 / 250000, attractiveness) location_names(model) location_names(model) <- NULL location_names(model) location_names(model) <- list( origin = french_cities$name[1:10], destination = LETTERS[1:10] ) destination_names(model) origin_names(model)
These functions provide low level access to origin and destination local
positions. It is recommended to use origin_positions() and
destination_positions() instead of location_positions and
location_positions<-.
location_positions(sim) location_positions(sim) <- valuelocation_positions(sim) location_positions(sim) <- value
sim |
a spatial interaction model object |
value |
a list with two components (see the returned value) or |
for location_positions NULL or a list with two components:
origin for the origin location positions and destination for the
destination location positions. For location_positions<-() the modified
sim object.
Location positions are given by numeric matrices with 2 or 3 columns. The
first two columns are assumed to be geographical coordinates while the 3rd
column can be used for instance to store altitude. Coordinates are
interpreted as is in graphical representations (see autoplot.sim()). They
are not matched to the costs as those can be derived from complex movement
models and other non purely geographic considerations.
origin_positions(), destination_positions()
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) ## No positions location_positions(model) <- list( origin = positions, destination = positions ) destination_positions(model) origin_positions(model)positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) ## No positions location_positions(model) <- list( origin = positions, destination = positions ) destination_positions(model) origin_positions(model)
This function computes all pairwise distances between spatial interaction
models (SIMs) of its x parameter, using sim_distance() with the specified
distance parameters. Then it returns the "median" of the collection defined
as the SIM that is in average the closest to all the other SIMs. Tie breaking
uses the order of the SIMs in the collection.
## S3 method for class 'sim_list' median( x, na.rm = FALSE, flows = c("full", "destination", "attractiveness"), method = c("euclidean"), return_distances = FALSE, ... )## S3 method for class 'sim_list' median( x, na.rm = FALSE, flows = c("full", "destination", "attractiveness"), method = c("euclidean"), return_distances = FALSE, ... )
x |
a collection of spatial interaction models, an object of class
|
na.rm |
not used |
flows |
|
method |
the distance measure to be used. Currently only |
return_distances |
should the distances computed to find the median be
returned as a |
... |
additional parameters (not used currently) |
As distance calculation can be slow in a large collection of SIMs, the
distance object can be returned as a distances attribute of the median SIM
by setting the return_distances parameter to TRUE. In addition, the
returned SIM has always two attributes:
index gives the index of the mode in the original sim_list
distortion gives the mean of the distances from the median SIM to all
the other SIMs
a spatial interaction model, an object of class sim with
additional attributes
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_median <- median(all_flows) attr(all_flows_median, "index") attr(all_flows_median, "distortion") median(all_flows, flows = "destination") median(all_flows, flows = "attractiveness")distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_median <- median(all_flows) attr(all_flows_median, "index") attr(all_flows_median, "distortion") median(all_flows, flows = "destination") median(all_flows, flows = "attractiveness")
Set the column names of a SIM data frame. Renaming the sim_list column is
supported and tracked, but renaming any core column turns the sim_df
into a standard data.frame.
## S3 replacement method for class 'sim_df' names(x) <- value## S3 replacement method for class 'sim_df' names(x) <- value
x |
data frame of spatial interaction models, an object of class
|
value |
unique names for the columns of the data frame |
a sim_df data frame if possible, a standard data.frame when
this is not possible.
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) names(all_flows_df) names(all_flows_df)[7] <- "my_sim" names(all_flows_df) ## still a sim_df class(all_flows_df) names(all_flows_df)[1] <- "return to scale" names(all_flows_df) ## not a sim_df class(all_flows_df)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) names(all_flows_df) names(all_flows_df)[7] <- "my_sim" names(all_flows_df) ## still a sim_df class(all_flows_df) names(all_flows_df)[1] <- "return to scale" names(all_flows_df) ## not a sim_df class(all_flows_df)
This function computes the most important flows in a spatial interaction
model according to the approach outlined by J. D. Nystuen and M. F. Dacey
(Nystuen & Dacey, 1961. In this work, a nodal flow is the largest flow sent
from a non terminal location (based on the definition of terminals recalled
in terminals()). The nodal structure is the collection of those flows.
They form an oriented graph that has interesting properties. In particular
each weakly connected component contains a single terminal location which can
be seen as the dominant location of the component. Notice that because nodal
flows are based on terminals, this function applies only to the non bipartite
setting.
nd_graph(sim, definition = c("ND", "RW"), ...) ## S3 method for class 'sim' nd_graph(sim, definition = c("ND", "RW"), ...)nd_graph(sim, definition = c("ND", "RW"), ...) ## S3 method for class 'sim' nd_graph(sim, definition = c("ND", "RW"), ...)
sim |
a spatial interaction model object |
definition |
terminal definition, either |
... |
additional parameters |
In practice, the function computes first the terminals and non terminals according to either Nystuen & Dacey (1961) or Rihll and Wilson (1991). Then it extracts the nodal flows. The result of the analysis is returned as a data frame with three columns:
from: the index of the non terminal origin location
to: the index of destination location of the nodal flow of from
flow: the value of the nodal flow
An important aspect of the node structure is that is does not contain isolated terminals. If a location is a terminal but is never the receiver of a nodal flow it will not appear in the collection of nodal flows. It constitutes a a trivial connected component in itself.
a data frame describing the Nystuen and Dacey graph a.k.a. the nodal structure of a spatial interaction model
Nystuen, J.D. and Dacey, M.F. (1961), "A graph theory interpretation of nodal regions", Papers and Proceedings of the Regional Science Association 7: 29-42. doi:10.1007/bf01969070
Rihll, T., and Wilson, A. (1991), "Modelling settlement structures in ancient Greece: new approaches to the polis", In City and Country in the Ancient World, Vol. 3, Edited by J. Rich and A. Wallace-Hadrill, 58-95. London: Routledge.
sim_is_bipartite(), is_terminal(), terminals()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- blvim(distances, production, 1.3, 1 / 250, attractiveness, bipartite = FALSE ) destination_names(model) <- french_cities$name[1:10] nd_graph(model) dist_times <- french_cities_times[1:15, 1:15] tmodel <- blvim(dist_times, rep(1, 15), 1.3, 1 / 5000, rep(1, 15), bipartite = FALSE ) destination_names(tmodel) <- french_cities$name[1:15] terminals(tmodel, definition = "RW") nd_graph(tmodel, "RW")distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- blvim(distances, production, 1.3, 1 / 250, attractiveness, bipartite = FALSE ) destination_names(model) <- french_cities$name[1:10] nd_graph(model) dist_times <- french_cities_times[1:15, 1:15] tmodel <- blvim(dist_times, rep(1, 15), 1.3, 1 / 5000, rep(1, 15), bipartite = FALSE ) destination_names(tmodel) <- french_cities$name[1:15] terminals(tmodel, definition = "RW") nd_graph(tmodel, "RW")
Functions to get or set the names of the origin locations in a spatial interaction model (or in a collection of spatial interaction models).
origin_names(sim) origin_names(sim) <- valueorigin_names(sim) origin_names(sim) <- value
sim |
a spatial interaction model object (an object of class |
value |
a character vector of length equal to the number of origin
locations, or |
for origin_names NULL or a character vector with one name per
origin locations in the model. for origin_names<- the modified sim
object or sim_list object.
location_names(), destination_names()
distances <- french_cities_distances[1:10, 1:10] production <- rep(1, 10) attractiveness <- rep(1, 10) ## the row/column names of the cost matrix are used for the location model <- static_blvim(distances, production, 1.5, 1 / 250000, attractiveness) origin_names(model) origin_names(model) <- french_cities$name[1:10] origin_names(model)distances <- french_cities_distances[1:10, 1:10] production <- rep(1, 10) attractiveness <- rep(1, 10) ## the row/column names of the cost matrix are used for the location model <- static_blvim(distances, production, 1.5, 1 / 250000, attractiveness) origin_names(model) origin_names(model) <- french_cities$name[1:10] origin_names(model)
Functions to get or set the positions of the origin locations in a spatial interaction model.
origin_positions(sim) origin_positions(sim) <- valueorigin_positions(sim) origin_positions(sim) <- value
sim |
a spatial interaction model object |
value |
a matrix with as many rows as the number of origin locations
and 2 or 3 columns, or |
for origin_positions NULL or the coordinate matrix for the
origin locations. for origin_positions<- the modified sim object
Location positions are given by numeric matrices with 2 or 3 columns. The
first two columns are assumed to be geographical coordinates while the 3rd
column can be used for instance to store altitude. Coordinates are
interpreted as is in graphical representations (see autoplot.sim()). They
are not matched to the costs as those can be derived from complex movement
models and other non purely geographic considerations.
location_positions(), destination_positions()
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) origin_positions(model) <- positions origin_positions(model)positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) origin_positions(model) <- positions origin_positions(model)
Extract the production constraints from a spatial interaction model object
production(sim, ...)production(sim, ...)
sim |
a spatial interaction model object |
... |
additional parameters |
a vector of production constraints at the origin locations
attractiveness(), destination_flow()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) production(model) ## the names of the production vector are set from the distance matrix ## we remove them for testing equality all.equal(as.numeric(production(model)), production)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) production(model) ## the names of the production vector are set from the distance matrix ## we remove them for testing equality all.equal(as.numeric(production(model)), production)
The function computes the specified quantiles for individual or aggregated flows in a collection of spatial interaction models.
## S3 method for class 'sim_list' quantile( x, flows = c("full", "destination", "attractiveness"), probs = seq(0, 1, 0.25), normalisation = c("none", "origin", "full"), ... )## S3 method for class 'sim_list' quantile( x, flows = c("full", "destination", "attractiveness"), probs = seq(0, 1, 0.25), normalisation = c("none", "origin", "full"), ... )
x |
a collection of spatial interaction models, a |
flows |
|
probs |
numeric vector of probabilities with values in |
normalisation |
when |
... |
additional parameters, not used currently |
The exact behaviour depends on the value of flows. In all cases, the
function returns a data frame. The data frame has one column per quantile,
named after the quantile using a percentage based named, as in the default
stats::quantile(). For a graphical representation of those quantiles, see
autoplot.sim_list().
if flows="full": this is the default case in which the function
computes for each pair of origin and destination locations
the quantiles of the distribution of the flow from to in the
collection of models (the flows maybe normalised before quantile
calculations, depending on the value of normalisation). In addition to the
quantiles, the data frame has two columns:
origin_idx: identifies the origin location by its index from 1 to the number
of origin locations;
destination_idx: identifies the destination location by its index from 1
to the number of destination locations.
if flows="destination", the function computes quantiles for each
destination location of the distribution of its incoming flow
(destination_flow()) in the collection of models. In addition to the
quantiles, the data frame has one column destination that identifies the
destination location by its index from 1 to the number of destination
locations.
if flows="attractiveness", the function computes quantiles for each
destination location of the distribution of its attractiveness
(attractiveness()) in the collection of models. In addition to the
quantiles, the data frame has one column destination that identifies the
destination location by its index from 1 to the number of destination
locations.
a date frame, see details
fortify.sim_list() autoplot.sim_list()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000 ) head(quantile(all_flows)) head(quantile(all_flows, flows = "destination")) head(quantile(all_flows, flows = "attractiveness", probs = c(0.1, 0.3, 0.6, 0.9) ))distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.1), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000 ) head(quantile(all_flows)) head(quantile(all_flows, flows = "destination")) head(quantile(all_flows, flows = "attractiveness", probs = c(0.1, 0.3, 0.6, 0.9) ))
Extract the return to scale parameter used to compute this model
return_to_scale(sim, ...)return_to_scale(sim, ...)
sim |
a spatial interaction model with a return to scale parameter |
... |
additional parameters |
the return to scale parameter
positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) return_to_scale(model) ## should be 1.5positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) return_to_scale(model) ## should be 1.5
Generates random points in the weighted simplex
The sample is divided between points in the relative interior of the simplex and points lying on proper faces, that is, points having at least one coordinate equal to zero.
Interior points are sampled uniformly after mapping the weighted simplex to the standard simplex. Points on faces are allocated preferentially to faces of increasing codimension: codimension-one faces are selected first, followed by codimension-two faces, and so on.
sample_weighted_simplex( kappa, C = 1, n = 1000L, face_fraction = 0.5, points_per_face = 0L, max_partial_attempts = 1000000L )sample_weighted_simplex( kappa, C = 1, n = 1000L, face_fraction = 0.5, points_per_face = 0L, max_partial_attempts = 1000000L )
kappa |
A numeric vector of length at least two containing finite,
strictly positive weights. These weights define the simplex constraint
|
C |
A finite, strictly positive numeric scalar giving the right-hand side of the weighted-simplex constraint. |
n |
A finite positive integer giving the total number of sample points. |
face_fraction |
A finite numeric scalar in |
points_per_face |
A nonnegative integer controlling the allocation of face points. If zero, the function maximizes the number of distinct selected faces. If positive, it is the requested maximum number of points on each selected nonvertex face. The function may increase this value, with a warning, when the requested value is incompatible with the number and capacities of the available faces. |
max_partial_attempts |
A positive integer giving the maximum number of random proposals allowed when selecting distinct faces in the final, partially covered codimension layer. Increasing this value may be necessary when most faces in that layer must be selected. |
The number of points allocated to proper faces is
where is face_fraction. The remaining points are sampled in the
relative interior.
The points_per_face argument controls the tradeoff between covering many
distinct faces and sampling each selected face more thoroughly:
If points_per_face = 0, the function selects as many different faces as
possible. If the number of requested face points does not exceed the
number of available faces, one point is generated on each selected face.
Otherwise, all available faces are selected and the remaining points are
distributed as evenly as possible among the nonvertex faces.
If points_per_face > 0, the function attempts to generate at most
points_per_face points on each selected nonvertex face. It selects the
smallest number of faces that provides sufficient capacity, prioritizing
low-codimension faces.
A vertex has support size one and consists of a single point. Consequently,
a selected vertex can receive only one sample point, irrespective of the
value of points_per_face.
If the requested positive value of points_per_face is too small to
accommodate all face points, even after selecting all available faces, the
function increases it to the smallest feasible value and issues a warning.
The adjusted value is returned in the points_per_face component.
For each selected nonvertex face, points are sampled uniformly with respect
to the intrinsic Lebesgue measure on that face. If is the support
of the selected face, independent standard exponential variables
, , are generated and transformed according to
with for .
A list with the following components:
Z: An n by length(kappa) numeric matrix. Each row is nonnegative and
satisfies sum(kappa * Z[i, ]) = C, up to floating-point error.
codim: An integer vector of length n. The value for a point is the
number of coordinates constrained to zero. Interior points have
codimension zero.
zero: A list of length n. Each element contains the one-based indices
of the zero coordinates of the corresponding row of Z.
face_id: An integer vector of length n. Interior points have identifier
zero. Face points have a positive identifier referring to an entry in
faces.
faces: A list describing the selected faces, with components:
codim: codimension of each selected face;
zero: one-based zero-coordinate indices for each selected face;
count: number of sample points generated on each selected face.
n_interior: Number of points sampled in the relative interior.
n_face_points: Number of points allocated to proper faces.
selected_face_count: Number of distinct selected faces.
points_per_face: Effective maximum number of points per selected
nonvertex face. This may be larger than the value requested by the user.
layer_count: Integer vector whose kth element is the number of selected
codimension-k faces.
partial_codim: Codimension of the final partially selected face layer,
or NA_integer_ if all selected codimension layers are complete.
max_constraint_error: Maximum absolute residual
over the generated points.
In dimension , there are
proper nonempty faces, including the vertices. There are therefore
proper nonvertex faces.
The number of codimension- faces is
The function exhausts these layers in increasing order of . Only the
final selected codimension layer can be partial.
Sampling uses R's random-number generator. Use set.seed() before calling
the function to obtain reproducible results.
Coordinates outside the selected support are set exactly to zero. Positive coordinates are generated using normalized exponential random variables. The weighted equality constraint is satisfied up to floating-point roundoff.
# Mix interior points and face points. set.seed(123) sample <- sample_weighted_simplex( kappa = c(1, 2, 3, 4, 5), C = 2, n = 1000, face_fraction = 0.4, points_per_face = 16 ) dim(sample$Z) sample$n_interior sample$n_face_points sample$selected_face_count # Verify the weighted equality constraint. max(abs(drop(sample$Z %*% c(1, 2, 3, 4, 5)) - 2)) # Maximize coverage across distinct faces. set.seed(123) broad_sample <- sample_weighted_simplex( kappa = rep(1, 5), n = 100, face_fraction = 1, points_per_face = 0 ) broad_sample$selected_face_count broad_sample$faces$count # Inspect the selected codimension layers. broad_sample$layer_count# Mix interior points and face points. set.seed(123) sample <- sample_weighted_simplex( kappa = c(1, 2, 3, 4, 5), C = 2, n = 1000, face_fraction = 0.4, points_per_face = 16 ) dim(sample$Z) sample$n_interior sample$n_face_points sample$selected_face_count # Verify the weighted equality constraint. max(abs(drop(sample$Z %*% c(1, 2, 3, 4, 5)) - 2)) # Maximize coverage across distinct faces. set.seed(123) broad_sample <- sample_weighted_simplex( kappa = rep(1, 5), n = 100, face_fraction = 1, points_per_face = 0 ) broad_sample$selected_face_count broad_sample$faces$count # Inspect the selected codimension layers. broad_sample$layer_count
Get the collection of spatial interaction models from a SIM data frame
sim_column(sim_df)sim_column(sim_df)
sim_df |
a data frame of spatial interaction models, an object of class
|
the collection of spatial interaction models in the sim_df object,
as a sim_list object
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows, sim_colum = "my_col") names(all_flows_df) sim_column(all_flows_df)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows, sim_colum = "my_col") names(all_flows_df) sim_column(all_flows_df)
Some spatial interaction models are the result of an iterative calculation,
see for instance blvim(). This calculation may have been interrupted before
convergence. The present function returns TRUE if the calculation converged,
FALSE if this was not the case and NA if the spatial interaction model
is not the result of an iterative calculation. The function applies also to a
collection of spatial interaction models as represented by a sim_list.
sim_converged(sim, ...)sim_converged(sim, ...)
sim |
a spatial interaction model object (an object of class |
... |
additional parameters |
TRUE, FALSE or NA, as described above. In the case of a sim_list
the function returns a logical vector with one value per model.
sim_iterations(), blvim(), grid_blvim()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) destination_flow(model) sim_converged(model) ## must be NAdistances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) destination_flow(model) sim_converged(model) ## must be NA
Some spatial interaction models use conversion factors between attractivenesses
and incoming flows to check whether they are at equilibrium (see for instance
blvim()). This function returns the conversion factors used to build this
sim is applicable, and NA when not applicable.
sim_conversion(sim, ...)sim_conversion(sim, ...)
sim |
a spatial interaction model object (an object of class |
... |
additional parameters |
NA or a vector of conversion factors, see above.
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_conversion(model) ## must NAdistances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_conversion(model) ## must NA
This function build a data frame from a collection of spatial interaction
models. The data frame has a list column sim of type sim_list which
stores the collection of models and classical columns that contain
characteristics of the models. The name of the list column can be set to
something else than sim (but not a name used by other default columns). See
details for the default columns.
sim_df(x, sim_column = "sim")sim_df(x, sim_column = "sim")
x |
a collection of spatial interaction models, an object of class
|
sim_column |
the name of the |
The data frame has one row per spatial interaction model and the following columns:
sim (default name): the last column that contains the models
alpha: the return to scale parameter used to build the model
beta: the cost inverse scale parameter used to build the model
diversity: model default diversity() (Shannon's diversity)
potential: model potential sim_potential()
iterations: the number of iterations used to produce the model (1 for a
static model)
converged: TRUE is the iterative calculation of the model converged
(for models produced by blvim() and related approaches), FALSE for no
convergence and NA for static models
The resulting object behaves mostly like a data.frame and support standard
extraction and replacement operators. The object tries to keep its sim_df
class during modifications. In particular, names<-.sim_df() tracks
name change for the sim_list column. If a modification or an extraction
operation changes the type of the sim_list column or drops it, the
resulting object is a standard data.frame. See [.sim_df
and names<-.sim_df() for details.
a data frame representation of the spatial interaction model
collection with classes sim_df and data.frame
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) all_flows_df$converged ## change the name of the sim column names(all_flows_df)[7] <- "models" ## still a sim_df class(all_flows_df) ## get the models sim_column(all_flows_df)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) all_flows_df$converged ## change the name of the sim column names(all_flows_df)[7] <- "models" ## still a sim_df class(all_flows_df) ## get the models sim_column(all_flows_df)
Extract or replace subsets of SIM data frames. The behaviour of the functions
is very close to the one of the corresponding data.frame functions, see
[.data.frame. However, modifications of the SIM columns
or suppression of core columns will turn the object into a standard
data.frame to void issues in e.g. graphical representations, see below for
details.
## S3 replacement method for class 'sim_df' x$name <- value ## S3 replacement method for class 'sim_df' x[[i, j]] <- value ## S3 method for class 'sim_df' x[i, j, ..., drop] ## S3 replacement method for class 'sim_df' x[i, j] <- value## S3 replacement method for class 'sim_df' x$name <- value ## S3 replacement method for class 'sim_df' x[[i, j]] <- value ## S3 method for class 'sim_df' x[i, j, ..., drop] ## S3 replacement method for class 'sim_df' x[i, j] <- value
x |
data frame of spatial interaction models, an object of class
|
name |
a literal character string |
value |
a suitable replacement value: it will be repeated a whole number
of times if necessary and it may be coerced: see the Coercion section in
|
i, j, ...
|
elements to extract or replace. For |
drop |
If |
In a sim_df, the core columns are derived from the sim_list column.
Replacement functions maintain this property by updating the columns after
any modification of the sim_list column. Modifications of the core columns
are rejected (removing a core column is accepted but this turns the sim_df
into a standard data.frame).
In addition, the sim_list column obeys to restriction on sim_list
modifications (i.e, a sim_list contains a homogeneous collection of spatial
interaction models).
Extraction functions keep the sim_df class only if the result is a data
frame with a sim_list column, the core columns and potentially additional
columns.
For [ a sim_df, a data.frame or a single column depending on the
values of i and j.
For [[ a column of the sim_df (or NULL) or an element of a column when
two indices are used.
For $ a column of the sim_df (or NULL).
For [<-, [[<-, and $<- a sim_df or a data frame (see
details).
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) ## the models as a sim_list all_flows_df[, "sim"] ## sub data frame, a sim_df all_flows_df[1:5, ] ## sub data frame, not a sim_df (alpha is missing) all_flows_df[6:10, 2:7] all_flows_2 <- grid_blvim(distances, log(french_cities$population[1:10]), seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) ## replace the sim_list column by the new models ## before all_flows_df$diversity all_flows_df$sim <- all_flows_2 ## after (all core columns have been updated) all_flows_df$diversitydistances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim(distances, production, seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) all_flows_df <- sim_df(all_flows) ## the models as a sim_list all_flows_df[, "sim"] ## sub data frame, a sim_df all_flows_df[1:5, ] ## sub data frame, not a sim_df (alpha is missing) all_flows_df[6:10, 2:7] all_flows_2 <- grid_blvim(distances, log(french_cities$population[1:10]), seq(1.05, 1.45, by = 0.2), seq(1, 3, by = 0.5) / 400, attractiveness, bipartite = FALSE, epsilon = 0.1, iter_max = 1000, ) ## replace the sim_list column by the new models ## before all_flows_df$diversity all_flows_df$sim <- all_flows_2 ## after (all core columns have been updated) all_flows_df$diversity
Some spatial interaction models are the result of an iterative calculation,
see for instance blvim(). This function returns a diagnostic object that
gives details about the behaviour of the iterative algorithm.
sim_diagnostic(sim, ...)sim_diagnostic(sim, ...)
sim |
a spatial interaction model object (an object of class |
... |
additional parameters |
a diagnostic object of class deriving from sim_diagnostic
sim_iterations(), sim_converged(), blvim(), grid_blvim()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_diagnostic(model)distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_diagnostic(model)
This function extracts from each spatial interaction model of the collection a vector representation derived from its flow matrix (see details). This vector is then used to compute distances between the models.
sim_distance( sim_list, flows = c("full", "destination", "attractiveness"), method = c("euclidean"), ... )sim_distance( sim_list, flows = c("full", "destination", "attractiveness"), method = c("euclidean"), ... )
sim_list |
a collection of spatial interaction models, an object of
class |
flows |
|
method |
the distance measure to be used. Currently only |
... |
additional parameters (not used currently) |
The vector representation is selected using the flows parameters. Possible
values are
"full" (default value): the representation is obtained by considering
the matrix of flows() as a vector (with the standard as.vector()
function);
"destination": the representation is the destination_flow() vector
associated to each spatial interaction model;
"attractiveness": the representation is the attractiveness() vector
associated to each spatial interaction model.
an object of class "dist"
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1 ) flows_distances <- sim_distance(all_flows) inflows_distances <- sim_distance(all_flows, "destination")distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) all_flows <- grid_blvim( distances, production, c(1.25, 1.5), c(1, 2, 3, 4) / 500, attractiveness, epsilon = 0.1 ) flows_distances <- sim_distance(all_flows) inflows_distances <- sim_distance(all_flows, "destination")
The function computes the Jacobian matrix of the fixed-point map with respect
to the attractivenesses. The function returns NA when applied to spatial
interaction models for which the map is not defined.
sim_fp_jacobian(sim, ...)sim_fp_jacobian(sim, ...)
sim |
a spatial interaction model, an object of class |
... |
additional parameters |
The fixed point map is given by
where is the destination flow (destination_flow()),
the attractiveness (attractiveness()) and the
conversion factor (sim_conversion()), all attached to destination .
The derivatives of w.r.t. can be used to
analyse the stability of the solutions, among other applications.
a p x p matrix or NA
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_fp_jacobian(model) ## must NAdistances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_fp_jacobian(model) ## must NA
The function returns TRUE is the spatial interaction model (SIM) is
bipartite, that is if the origin locations are distinct from the destination
locations (at least from the analysis point of view). The function return
FALSE when the SIM uses the same locations for origin and destination.
sim_is_bipartite(sim)sim_is_bipartite(sim)
sim |
a spatial interaction model object |
TRUE if the spatial interaction model is bipartite, FALSE if
not.
positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) ## returns TRUE despite the use of a single set of positions sim_is_bipartite(model) ## now we are clear about the non bipartite nature of the model model <- static_blvim(distances, production, 1.5, 1, attractiveness, bipartite = FALSE ) sim_is_bipartite(model)positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) model <- static_blvim(distances, production, 1.5, 1, attractiveness) ## returns TRUE despite the use of a single set of positions sim_is_bipartite(model) ## now we are clear about the non bipartite nature of the model model <- static_blvim(distances, production, 1.5, 1, attractiveness, bipartite = FALSE ) sim_is_bipartite(model)
Returns the number of iterations used to produce this spatial interaction model
sim_iterations(sim, ...)sim_iterations(sim, ...)
sim |
a spatial interaction model object (an object of class |
... |
additional parameters |
a number of iterations that may be one if the spatial interaction model
has been obtained using a static model (see static_blvim()). In the case of a sim_list
the function returns a vector with iteration number per model.
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) destination_flow(model) sim_iterations(model) ## must be onedistances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) destination_flow(model) sim_iterations(model) ## must be one
The collection of sim objects represented by a sim_list object is assumed
to be homogeneous, that is to correspond to a fix set of origin and
destination locations, associated to a fixed cost matrix.
sim_list(sims, validate = TRUE)sim_list(sims, validate = TRUE)
sims |
a list of homogeneous spatial interaction objects |
validate |
should the function validate the homogeneity of the list of
spatial interaction objects (defaults to |
a sim_list object
positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) flows_1 <- blvim(distances, production, 1.5, 1, attractiveness) flows_2 <- blvim(distances, production, 1.25, 2, attractiveness) all_flows <- sim_list(list(flows_1, flows_2))positions <- matrix(rnorm(10 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 10) attractiveness <- c(2, rep(1, 9)) flows_1 <- blvim(distances, production, 1.5, 1, attractiveness) flows_2 <- blvim(distances, production, 1.25, 2, attractiveness) all_flows <- sim_list(list(flows_1, flows_2))
Extracts the summary information recorded for the individual algorithm runs
performed by blvim().
sim_logs(sim)sim_logs(sim)
sim |
A spatial interaction model object produced by |
The rows are in the same order as the starting points supplied to
blvim(). Run summaries are retained independently of the value of
keep_runs.
See blvim() for details on how the nature of Z determines the number of
runs. Use sim_runs() to extract the complete spatial interaction models
produced by the runs when they are available.
A data frame with one row per starting point and the following columns:
run: index of the run from 1 to number of initial point;
potential: final value of the potential;
iterations: number of iterations;
converged: whether the run converged;
status: stopping status reported by a gradient algorithm;
gap: final duality gap reported by a gradient algorithm;
evaluations: number of evaluations reported by a gradient algorithm.
The gradient-specific columns are NA for the "linear" and
"quadratic" algorithms.
blvim(), sim_runs(), sim_diagnostic()
The potential is a measure of stability introduced by Osawa, Akamatsu, and
Kogure. It is only defined for spatial interaction models produced by BLV
framework. The function returns NA when applied to spatial interaction
models for which a potential cannot be computed. The function applies also to
a collection of spatial interaction models as represented by a sim_list.
sim_potential(sim, ...)sim_potential(sim, ...)
sim |
a spatial interaction model object (an object of class |
... |
additional parameters |
The potential of a spatial interaction model is given by
where:
is the cost matrix costs()
is the production constraint vector (production())
is the return to scale parameter (return_to_scale())
is the inverse of a cost scale parameter (inverse_cost())
is the attractiveness vector (attractiveness())
is the conversion factor vector (sim_conversion())
the scalar value of the potential. In the case of a sim_list
the function returns a numerical vector with one value per model.
Osawa, M., Akamatsu, T., & Kogure, Y. (2025). "Most likely retail agglomeration patterns: Potential maximization and stochastic stability of spatial equilibria". https://arxiv.org/abs/2011.06778v2
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_conversion(model) ## must NAdistances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- log(french_cities$population[1:10]) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 250, attractiveness) sim_conversion(model) ## must NA
Extracts the spatial interaction model objects produced by the individual
algorithm runs performed by blvim().
sim_runs(sim)sim_runs(sim)
sim |
A spatial interaction model object produced by |
See blvim() for details on how the Z and keep_runs arguments
determine the number of runs and whether the resulting spatial interaction
models are retained.
Use sim_logs() to extract the summary information recorded for multiple
runs.
A sim_list containing the available spatial interaction models
produced by the algorithm runs.
For a single run, a one-element sim_list containing x is returned. For
multiple runs retained with keep_runs = TRUE, the models are returned in
the same order as the starting points supplied to blvim().
The function returns NULL when multiple runs were performed with
keep_runs = FALSE, because only the model associated with the largest
final potential was retained.
blvim() for computing equilibrium spatial interaction models,
sim_logs() for summaries of the individual runs, and sim_list() for
collections of spatial interaction models.
distances <- french_cities_distances[1:10, 1:10] / 1000 production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) attractiveness <- attractiveness / sum(attractiveness) * sum(production) # A single run produces a one-element sim_list single <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "gradient" ) length(sim_runs(single)) # Multiple runs, with all spatial interaction models retained starts <- rbind( attractiveness, rev(attractiveness) ) multiple <- blvim( distances, production, 1.5, 1 / 250, starts, algorithm = "gradient", keep_runs = TRUE ) runs <- sim_runs(multiple) length(runs)distances <- french_cities_distances[1:10, 1:10] / 1000 production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) attractiveness <- attractiveness / sum(attractiveness) * sum(production) # A single run produces a one-element sim_list single <- blvim( distances, production, 1.5, 1 / 250, attractiveness, algorithm = "gradient" ) length(sim_runs(single)) # Multiple runs, with all spatial interaction models retained starts <- rbind( attractiveness, rev(attractiveness) ) multiple <- blvim( distances, production, 1.5, 1 / 250, starts, algorithm = "gradient", keep_runs = TRUE ) runs <- sim_runs(multiple) length(runs)
This function computes flows between origin locations and destination locations according to the production constrained entropy maximising model proposed by A. Wilson.
static_blvim( costs, X, alpha, beta, Z, bipartite = TRUE, origin_data = NULL, destination_data = NULL )static_blvim( costs, X, alpha, beta, Z, bipartite = TRUE, origin_data = NULL, destination_data = NULL )
costs |
a cost matrix |
X |
a vector of production constraints |
alpha |
the return to scale parameter |
beta |
the inverse cost scale parameter |
Z |
a vector of destination attractivenesses |
bipartite |
when |
origin_data |
|
destination_data |
|
The model computes flows using the following parameters:
costs () is a matrix whose entry is
the cost of having a "unitary" flow from origin location to
destination location
X () is a vector of size containing non negative
production constraints for the origin locations
alpha () is a return to scale parameter that enhance (or
reduce if smaller that 1) the attractivenesses of destination locations when
they are larger than 1
beta () is the inverse of a cost scale parameter, i.e., costs
are multiplied by beta in the model
Z () is a vector of size containing the positive
attractivenesses of the destination locations
According to Wilson's model, the flow from origin location to
destination location , , is given by
The model is production constrained because
that is the origin location sends a total flow of exactly .
an object of class sim (and sim_wpc) for spatial interaction
model that contains the matrix of flows from the origin locations to the
destination locations (see
above) and the attractivenesses of the destination locations.
While models in this package do not use location data beyond X and Z,
additional data can be stored and used when analysing spatial interaction
models.
Spatial interaction models can store names for origin and destination
locations, using origin_names<-() and destination_names<-(). Names
are taken by default from names of the cost matrix costs. More precisely,
rownames(costs) is used for origin location names and colnames(costs) for
destination location names.
Spatial interaction models can store the positions of the origin and
destination locations, using origin_positions<-() and
destination_positions<-().
In addition to the functions mentioned above, location data can be specified
directly using the origin_data and destination_data parameters. Data are
given by a list whose components are not interpreted excepted the following
ones:
names is used to specify location names and its content has to follow
the restrictions documented in origin_names<-() and
destination_names<-()
positions is used to specify location positions and its content has
to follow the restrictions documented in origin_positions<-() and
destination_positions<-()
Wilson, A. (1971), "A family of spatial interaction models, and associated developments", Environment and Planning A: Economy and Space, 3(1), 1-32 doi:10.1068/a030001
origin_names(), destination_names(), origin_positions(),
destination_positions()
positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness, origin_data = list( names = french_cities$name[1:10], positions = positions ), destination_data = list( names = french_cities$name[1:10], positions = positions ) ) model location_names(model) location_positions(model)positions <- as.matrix(french_cities[1:10, c("th_longitude", "th_latitude")]) distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- log(french_cities$area[1:10]) model <- static_blvim(distances, production, 1.5, 1 / 500, attractiveness, origin_data = list( names = french_cities$name[1:10], positions = positions ), destination_data = list( names = french_cities$name[1:10], positions = positions ) ) model location_names(model) location_positions(model)
This function computes summary statistics on a collection of spatial
interaction models (in a sim_list).
## S3 method for class 'sim_list' summary(object, ...) ## S3 method for class 'summary_sim_list' print(x, ...)## S3 method for class 'sim_list' summary(object, ...) ## S3 method for class 'summary_sim_list' print(x, ...)
object |
a collection of spatial interaction models, an object of class
|
... |
additional parameters (not used currently) |
x |
an object of class |
The list returned by the function contains the following elements:
median: the median of the collection, as return by the median.sim_list()
function
distortion: the average distance of all elements of the collection to
the median model
withinss: the sum of all pairwise distances between the elements of
the collection
nb_sims: the size of the collection
In addition, if the collection contains non bipartite models, the result
has another element, nb_configurations which gives the number of distinct
terminal sets in the collection, where the terminals are computed by
terminals(), using the "RW" definition.
an object of class summary_sim_list and list with a set of
summary statistics computed on the collection of spatial interaction models
median.sim_list(), terminals()
positions <- matrix(rnorm(15 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 15) attractiveness <- rep(1, 15) all_flows <- grid_blvim(distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3), attractiveness, epsilon = 0.1, bipartite = FALSE, ) summary(all_flows)positions <- matrix(rnorm(15 * 2), ncol = 2) distances <- as.matrix(dist(positions)) production <- rep(1, 15) attractiveness <- rep(1, 15) all_flows <- grid_blvim(distances, production, c(1.1, 1.25, 1.5), c(1, 2, 3), attractiveness, epsilon = 0.1, bipartite = FALSE, ) summary(all_flows)
This function identifies terminals in the locations underlying the given
spatial interaction model. Terminals are locally dominating locations that
essentially send less to other locations than they receive (see details for
formal definitions). As we compare incoming flows to outgoing flows, terminal
computation is restricted to interaction models in which origin and
destination locations are identical, i.e. models that are not bipartite.
terminals(sim, definition = c("ND", "RW"), ...)terminals(sim, definition = c("ND", "RW"), ...)
sim |
a spatial interaction model object |
definition |
terminal definition, either |
... |
additional parameters |
The notion of terminal used in this function is based on seminal work by J.
D. Nystuen and M. F. Dacey (Nystuen & Dacey, 1961), as well as on the follow
up variation from Rihll & Wislon (1987 and 1991). We assume given a square
flow matrix . The incoming flow
at location is given by
and is used as a measure of importance of this location. Then in Nystuen &
Dacey (1961), location is a "terminal point" (or a "central city") if
where is such that
In words, is a terminal if the location to which it sends
its largest flow is less important than itself, in terms of incoming
flows. This is the definition used by the function when definition is
"ND".
Rihll & Wilson (1987) use a modified version of this definition described in
details in Rihll and Wilson (1991). With this relaxed version, location
is a terminal if
In words, is a terminal if it receives more flows than it is sending
to each other location. It is easy to see that each Nystuen & Dacey terminal
is a Rihll & Wilson terminal, but the reverse is false in general. The
function use the Rihll & Wilson definition when definition is "RW"
a vector containing the indexes of the terminals identified from the flow matrix of the interaction model.
Nystuen, J.D. and Dacey, M.F. (1961), "A graph theory interpretation of nodal regions", Papers and Proceedings of the Regional Science Association 7: 29-42. doi:10.1007/bf01969070
Rihll, T.E., and Wilson, A.G. (1987). "Spatial interaction and structural models in historical analysis: some possibilities and an example", Histoire & Mesure 2: 5-32. doi:10.3406/hism.1987.1300
Rihll, T., and Wilson, A. (1991), "Modelling settlement structures in ancient Greece: new approaches to the polis", In City and Country in the Ancient World, Vol. 3, Edited by J. Rich and A. Wallace-Hadrill, 58-95. London: Routledge.
sim_is_bipartite(), is_terminal(), grid_is_terminal()
distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- blvim(distances, production, 1.3, 1 / 250, attractiveness, bipartite = FALSE ) destination_names(model) <- french_cities$name[1:10] terminals(model) dist_times <- french_cities_times[1:10, 1:10] tmodel <- blvim(dist_times, production, 1.3, 1 / 5000, attractiveness, bipartite = FALSE ) destination_names(tmodel) <- french_cities$name[1:10] terminals(tmodel) terminals(tmodel, definition = "RW")distances <- french_cities_distances[1:10, 1:10] / 1000 ## convert to km production <- rep(1, 10) attractiveness <- rep(1, 10) model <- blvim(distances, production, 1.3, 1 / 250, attractiveness, bipartite = FALSE ) destination_names(model) <- french_cities$name[1:10] terminals(model) dist_times <- french_cities_times[1:10, 1:10] tmodel <- blvim(dist_times, production, 1.3, 1 / 5000, attractiveness, bipartite = FALSE ) destination_names(tmodel) <- french_cities$name[1:10] terminals(tmodel) terminals(tmodel, definition = "RW")