-
Notifications
You must be signed in to change notification settings - Fork 16
/
fund-cycles.R
386 lines (357 loc) · 14 KB
/
fund-cycles.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
# flatten lists of lists to single list
flatten_list <- function (x) {
x2 <- list ()
for (i in x) {
x2 <- c (x2, i)
}
return (x2)
}
#' Calculate fundamental cycles in a graph.
#'
#' @param graph `data.frame` or equivalent object representing the contracted
#' network graph (see Details).
#' @param vertices `data.frame` returned from \link{dodgr_vertices}`(graph)`.
#' Will be calculated if not provided, but it's quicker to pass this if it has
#' already been calculated.
#' @param graph_max_size Maximum size submitted to the internal C++ routines as
#' a single chunk. Warning: Increasing this may lead to computer meltdown!
#' @param expand For large graphs which must be broken into chunks, this factor
#' determines the relative overlap between chunks to ensure all cycles are
#' captured. (This value should only need to be modified in special cases.)
#' @return List of cycle paths, in terms of vertex IDs in `graph` and, for
#' spatial graphs, the corresponding coordinates.
#'
#' @note Calculation of fundamental cycles is VERY computationally demanding,
#' and this function should only be executed on CONTRACTED graphs (that is,
#' graphs returned from \link{dodgr_contract_graph}), and even than may take a
#' long time to execute. Results for full graphs can be obtained with the
#' function \link{dodgr_full_cycles}. The computational complexity can also not
#' be calculated in advance, and so the parameter `graph_max_size` will lead to
#' graphs larger than that (measured in numbers of edges) being cut into smaller
#' parts. (Note that that is only possible for spatial graphs, meaning that it
#' is not at all possible to apply this function to large, non-spatial graphs.)
#' Each of these smaller parts will be expanded by the specified amount
#' (`expand`), and cycles found within. The final result is obtained by
#' aggregating all of these cycles and removing any repeated ones arising due to
#' overlap in the expanded portions. Finally, note that this procedure of
#' cutting graphs into smaller, computationally manageable sub-graphs provides
#' only an approximation and may not yield all fundamental cycles.
#'
#' @family misc
#' @export
#' @examples
#' net <- weight_streetnet (hampi)
#' graph <- dodgr_contract_graph (net)
#' verts <- dodgr_vertices (graph)
#' cyc <- dodgr_fundamental_cycles (graph, verts)
dodgr_fundamental_cycles <- function (graph,
vertices = NULL,
graph_max_size = 10000,
expand = 0.05) {
if (missing (graph)) {
stop ("graph must be provided")
}
if (!inherits (graph, "data.frame")) {
stop ("graph must be a data.frame object")
}
if (is.null (vertices)) {
vertices <- dodgr_vertices (graph)
}
if (!"flow" %in% names (graph)) { # makes no difference
graph$flow <- 1
}
graph <- merge_directed_graph (graph) # uses fast C++ routines
graph$flow <- NULL
bb <- get_graph_bb (graph)
if (nrow (graph) <= graph_max_size) {
bb_indices <- list (bb)
} else {
ndivs <- get_ndivs (graph, graph_max_size)
bb_list <- get_bb_list (bb, ndivs, expand = expand)
bb_data <- subdivide_bb (graph, bb_list, graph_max_size, expand)
bb_indices <- bb_data$bb_indices
}
graphc <- convert_graph (graph, dodgr_graph_cols (graph))
if (length (bb_indices) == 1) {
res <- rcpp_fundamental_cycles (graphc, verts = vertices)
} else {
message (
"Now computing fundamental cycles by breaking graph with ",
nrow (graphc),
" edges into ",
length (bb_indices),
" components ..."
)
pb <- utils::txtProgressBar (style = 3)
res <- list ()
for (i in seq (bb_indices)) {
graphi <- graphc [bb_indices [[i]], ]
verti <- dodgr_vertices (graphi)
res [[i]] <- rcpp_fundamental_cycles (graphi, verts = verti)
utils::setTxtProgressBar (pb, i / length (bb_indices))
}
close (pb)
# each element of res is a list, so flatten these:
res <- flatten_list (res)
# These hash each and remove any duplicated ones:
dig <- unlist (lapply (res, digest::digest))
res <- res [which (!duplicated (dig))]
}
if (is_graph_spatial (graph)) {
if (length (bb_indices) > 1) {
message (
"Generating spatial coordinates of polygons ",
"(this should be fairly quick ...)"
)
}
res <- lapply (res, function (i) {
index <- match (i, vertices$id)
data.frame (
id = i,
x = vertices$x [index],
y = vertices$y [index],
stringsAsFactors = FALSE
)
})
}
return (res)
}
# ********** FUNCTIONS TO BREAK SPATIAL GRAPHS INTO SUB-COMPONENTS **********
# Initial estimate of how many divisions needed
get_ndivs <- function (graph, graph_max_size) {
ndivs <- ceiling (nrow (graph) / graph_max_size)
ceiling (sqrt (ndivs)) # num of grid rows and cols
}
# get boundingn box of graph
get_graph_bb <- function (graph) {
gr_cols <- dodgr_graph_cols (graph)
from_lon <- graph [, gr_cols$xfr]
from_lat <- graph [, gr_cols$yfr]
to_lon <- graph [, gr_cols$xto]
to_lat <- graph [, gr_cols$yto]
apply (cbind (c (from_lon, to_lon), c (from_lat, to_lat)), 2, range)
}
# divide a bb into rectangular grid of sub-boxes and return list of
# corresponding bboxes
get_bb_list <- function (bb, ndivs, expand = 0.05) {
# divide one column of bb: either lons or lats
divide_bb_vec <- function (bb, ndivs, colnum = 2, expand) {
bb <- c (bb [1, colnum], vapply (
seq (ndivs), function (i) {
bb [1, colnum] + i / ndivs *
diff (bb [, colnum])
},
numeric (1)
))
bb <- cbind (bb [1:ndivs], bb [2:(ndivs + 1)])
t (apply (bb, 1, function (i) {
mean (i) + c (-0.5 - expand, 0.5 + expand) * diff (i)
}))
}
bb_lons <- divide_bb_vec (bb, ndivs, colnum = 1, expand = expand)
bb_lats <- divide_bb_vec (bb, ndivs, colnum = 2, expand = expand)
bb_list <- list ()
for (i in seq (ndivs)) {
for (j in seq (ndivs)) {
bb_list [[length (bb_list) + 1]] <-
cbind (bb_lons [i, ], bb_lats [j, ])
}
}
return (bb_list)
}
# get indices into graph of edges lying within each bit of a bb_list
get_bb_indices <- function (graph, bb_list) {
res <- list ()
for (i in seq (bb_list)) {
res [[i]] <- which (graph$from_lon > bb_list [[i]] [1, 1] &
graph$from_lon < bb_list [[i]] [2, 1] &
graph$from_lat > bb_list [[i]] [1, 2] &
graph$from_lat < bb_list [[i]] [2, 2] &
graph$to_lon > bb_list [[i]] [1, 1] &
graph$to_lon < bb_list [[i]] [2, 1] &
graph$to_lat > bb_list [[i]] [1, 2] &
graph$to_lat < bb_list [[i]] [2, 2])
}
return (res)
}
# Rectangularly subdivide any components of bb_list that are > graph_max_size
# into 4 sub-components.
subdivide_bb <- function (graph, bb_list, graph_max_size, expand) {
bb_indices <- get_bb_indices (graph, bb_list)
lens <- unlist (lapply (bb_indices, length))
while (any (lens > graph_max_size)) {
indx <- which (lens > graph_max_size)
bbs <- bb_list [indx]
bb_list [indx] <- NULL
bb_indices [indx] <- NULL
for (i in bbs) {
bb_list <- c (bb_list, get_bb_list (i, ndivs = 2, expand = expand))
bb_indices <- get_bb_indices (graph, bb_list)
lens <- unlist (lapply (bb_indices, length))
}
bb_list [which (lens == 0)] <- NULL
bb_indices [which (lens == 0)] <- NULL
}
list (bb_list = bb_list, bb_indices = bb_indices)
}
#' Calculate fundamental cycles on a FULL (that is, non-contracted) graph.
#'
#' @inheritParams dodgr_fundamental_cycles
#' @note This function converts the `graph` to its contracted form, calculates
#' the fundamental cycles on that version, and then expands these cycles back
#' onto the original graph. This is far more computationally efficient than
#' calculating fundamental cycles on a full (non-contracted) graph.
#'
#' @family misc
#' @export
#' @examples
#' \dontrun{
#' net <- weight_streetnet (hampi)
#' graph <- dodgr_contract_graph (net)
#' cyc1 <- dodgr_fundamental_cycles (graph)
#' cyc2 <- dodgr_full_cycles (net)
#' }
#' # cyc2 has same number of cycles, but each one is generally longer, through
#' # including all points intermediate to junctions; cyc1 has cycles composed of
#' # junction points only.
dodgr_full_cycles <- function (graph,
graph_max_size = 10000,
expand = 0.05) {
graph$flow <- 1
graph <- merge_directed_graph (graph)
graph$flow <- NULL
# graph <- graph [graph$component == 1, ]
graphc <- dodgr_contract_graph (graph)
v <- dodgr_vertices (graphc)
# edge_map <- get_edge_map (graphc) # TODO: Implement this
hashc <- get_hash (graphc, contracted = TRUE)
fname_c <- fs::path (
fs::path_temp (),
paste0 ("dodgr_edge_map_", hashc, ".Rds")
)
if (!fs::file_exists (fname_c)) {
stop ("something unexpected went wrong extracting the edge map")
} # nocov
edge_map <- readRDS (fname_c)
x <- dodgr_fundamental_cycles (graphc,
vertices = v,
graph_max_size = graph_max_size,
expand = expand
)
gr_cols <- dodgr_graph_cols (graphc)
from_to <- paste0 (graphc [[gr_cols$from]], "-", graphc [[gr_cols$to]])
to_from <- paste0 (graphc [[gr_cols$to]], "-", graphc [[gr_cols$from]])
ids <- lapply (x, function (i) {
idpairs <- paste0 (i$id [-length (i$id)], "-", i$id [-1])
# Get the edge pairs that match the idpairs, whether as from->to or
# to->from
edges_c1 <- graphc [[gr_cols$edge_id]] [match (idpairs, from_to)]
edges_c2 <- graphc [[gr_cols$edge_id]] [match (idpairs, to_from)]
edges_c <- apply (rbind (edges_c1, edges_c2), 2, function (i) {
i [which (!is.na (i))] [1]
})
edges_new <- lapply (as.list (edges_c), function (j) {
if (j %in% edge_map$edge_new) {
index <- which (edge_map$edge_new %in% j)
edges <- edge_map$edge_old [index]
j <- match (edges, graph [[gr_cols$edge_id]])
if ((graph [[gr_cols$from]] [j [1]] == # nolint
graph [[gr_cols$to]] [utils::tail (j, 1)]) ||
(graph [[gr_cols$from]] [j [1]] == # nolint
graph [[gr_cols$to]] [j [2]])) { # nolint
j <- rev (j)
}
} else {
j <- match (j, graph [[gr_cols$edge_id]])
}
return (as.numeric (j))
}) # end edges_new lapply
unlist (edges_new)
}) # end ids lapply
# ids at that point is a sequences of indices into graph. This is then
# converted to a sequence of vertex IDs, through just adding the last vertex
# of the sequence on to close the polygon
res <- lapply (ids, function (i) {
graph [c (i, i [1]), gr_cols$from, drop = TRUE]
})
if (is_graph_spatial (graph)) {
vertices <- dodgr_vertices (graph)
res <- lapply (res, function (i) {
index <- match (i, vertices$id)
data.frame (
id = i,
x = vertices$x [index],
y = vertices$y [index],
stringsAsFactors = FALSE
)
})
}
return (res)
}
#' Convert \pkg{sf} `LINESTRING` objects to `POLYGON` objects representing all
#' fundamental cycles within the `LINESTRING` objects.
#'
#' @inheritParams dodgr_fundamental_cycles
#' @param sflines An \pkg{sf} `LINESTRING` object representing a network.
#' @return An `sf::sfc` collection of `POLYGON` objects.
#' @family misc
#' @export
dodgr_sflines_to_poly <- function (sflines,
graph_max_size = 10000,
expand = 0.05) {
if (!(methods::is (sflines, "sf") || methods::is (sflines, "sf"))) {
stop ("lines must be an object of class 'sf' or 'sfc'")
}
if (!methods::is (
sflines [[attr (sflines, "sf_column")]],
"sfc_LINESTRING"
)) {
stop ("lines must be an 'sfc_LINESTRING' object")
}
graph <- weight_streetnet (sflines, wt_profile = 1)
# Different graph components need to be analysed seperately, and an
# arbitrary decision is made here to only consider components with > 100
# edges - smaller ones are unlikely to have any cycles.
comps <- table (graph$component)
comps <- as.numeric (names (comps) [which (comps > 100)])
x <- list ()
for (i in seq (comps)) {
graphi <- graph [graph$component == comps [i], ]
graphi$edge_id <- seq (nrow (graphi))
attr (graphi, "hash") <- NULL
x [[i]] <- dodgr_full_cycles (graphi,
graph_max_size = graph_max_size,
expand = expand
)
}
x <- flatten_list (x)
polys_to_sfc (x, sflines)
}
# convert list of polygon coordinates to `sf::sfc` collection
# code adapted from osmdata/tests/testthat/test-sf-construction.R
polys_to_sfc <- function (x, sflines) {
g <- sflines [[attr (sflines, "sf_column")]]
crs <- attr (g, "crs")
xy <- do.call (rbind, x)
xvals <- xy [, 2]
yvals <- xy [, 3]
bb <- structure (rep (NA_real_, 4),
names = c ("xmin", "ymin", "xmax", "ymax")
)
bb [1:4] <- c (min (xvals), min (yvals), max (xvals), max (yvals))
class (bb) <- "bbox"
attr (bb, "crs") <- crs
x <- lapply (x, function (i) {
res <- as.matrix (i [, 2:3])
colnames (res) <- NULL
structure (list (res),
class = c ("XY", "POLYGON", "sfg")
)
})
attr (x, "n_empty") <- 0
attr (x, "precision") <- 0.0
class (x) <- c ("sfc_POLYGON", "sfc")
attr (x, "bbox") <- bb
attr (x, "crs") <- crs
x
}