4. Networks 2: Shortest Paths and Road Networks¶
ISE 754, Fall 2024
Package Used: Functions from the following package are used in this notebook for the first time:
Graphs: Graphs.jl, "The goal of Graphs.jl is to offer a performant platform for network and graph analysis in Julia, following the example of libraries such as NetworkX in Python."SimpleWeightedGraphs: [SimpleWeightedGraphs.jl (https://juliagraphs.org/SimpleWeightedGraphs.jl/stable/), used to create graphs with undirected (two-way) edges and digraphs with directed (one-way) edges, with all of the edges have and associated weight (e.g., distance, if the edge represents a road segment).DelaunayTriangulation: DelaunayTriangulation.jl, computes Delaunay triangulations of 2-D points.
1. Network Representations¶
The terms network and graph can be used interchangeably. The terms graph/vertex/edge are typically used when the focus is on properties of the mathematical object, and the terms network/node/arc are used when the focus is on the properties of a real-work system. The terms graph and digraph (short for directed graph) are used in mathematics to describe networks in which the arcs are either all undirected or directed, respectively. In many network applications, both directed and undirected arcs are used (e.g., one- and two-way roads), and the arcs have associated weights that represent, for example, the distance or travel time of the road being represented.
| Network | Graph | Digraph |
|---|---|---|
| Node | Vertex | Vertex |
| Arc (undirected/directed) | Edge (undirected) | Edge (directed) |
In Julia, both graphs and digraphs can be created. If a network with a mix of directed and undirected arcs is needed, then a digraph can be used, with two directed arcs in opposite directions used to represent an undirected arc.
There are different network representations that make it easier to operate on the network for specific applications, and it is usually possible to translate from one representation to another. Given a network with $m$ nodes and $n$ arcs:
- Arc list: $n$-element array of the source nodes, n-element array of destination nodes, and an n-element array of the weights of each arc in the network.
- Adjacency matrix: $m \times m$ matrix of arc weights, with weights of $0$ indicating no arc.
- Interlevel matrix: $m_1 \times m_2$ matrix of $m_1 \times m_2$ arc weights, where all nodes can be separated in two groups consisting of $m_1$ arc-source and $m_2$ arc-destination nodes, that can be either be used directly to specify embedded into an adjacency matrix.
- Incidence matrix: $m \times n$ matrix, where values of $-1$ and $1$ in each column indicate the source and destination nodes of the arc represented by the column, that represents of coefficients of the flow-balance equations of the constraints in an LP.
Ex: Digraph (complete bipartite digraph)¶
Interlevel to Adjacency: Embed the $m \times n$ interlevel matrix C into a $(m + n) \times (m + n)$ adjacency matrix:
C = [6. 10.; 0. 13.; 9. 16.] # m x n interlevel matrix from m nodes to n nodes
m,n = size(C)
W = [zeros(m,m) C; zeros(n,m+n)] # (m+n) x (m+n) adjacency matrix
5×5 Matrix{Float64}:
0.0 0.0 0.0 6.0 10.0
0.0 0.0 0.0 0.0 13.0
0.0 0.0 0.0 9.0 16.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
The function SimpleWeightedDiGraph is used to create a directed network using an adjacency matrix, where each element ( i, j ) in the matrix represents an arc directed from a starting node ( i ) to an ending node ( j ):
using Graphs, SimpleWeightedGraphs
g = SimpleWeightedDiGraph(W) # Create digraph from adjacency matrix
{5, 5} directed simple Int64 graph with Float64 weights
adjacency_matrix(g) # Arcs represented by non-zero elements in sparse matrix
5×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 5 stored entries:
⋅ ⋅ ⋅ 6.0 10.0
⋅ ⋅ ⋅ ⋅ 13.0
⋅ ⋅ ⋅ 9.0 16.0
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
Since the adjacency matrix is stored as a sparse matrix, covert it to a regular (dense) matrix for display purposes:
Matrix(adjacency_matrix(g))
5×5 Matrix{Float64}:
0.0 0.0 0.0 6.0 10.0
0.0 0.0 0.0 0.0 13.0
0.0 0.0 0.0 9.0 16.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
An array comprehension can be used to iterate over all of the edges in the graph:
[e for e in edges(g)]
5-element Vector{SimpleWeightedEdge{Int64, Float64}}:
Edge 1 => 4 with weight 6.0
Edge 1 => 5 with weight 10.0
Edge 2 => 5 with weight 13.0
Edge 3 => 4 with weight 9.0
Edge 3 => 5 with weight 16.0
Problem: Arc 2 => 4 of 0 weight not included
replace!(C, 0. => sqrt(eps())) # Use small positive value for all 0 arcs
3×2 Matrix{Float64}:
6.0 10.0
1.49012e-8 13.0
9.0 16.0
W = [zeros(m,m) replace!(C, 0. => sqrt(eps())); zeros(n,m+n)]
g = SimpleWeightedDiGraph(W)
println.(edges(g));
Edge 1 => 4 with weight 6.0 Edge 1 => 5 with weight 10.0 Edge 2 => 4 with weight 1.4901161193847656e-8 Edge 2 => 5 with weight 13.0 Edge 3 => 4 with weight 9.0 Edge 3 => 5 with weight 16.0
src.(edges(g)) # Source or starting nodes of arcs
6-element Vector{Int64}:
1
1
2
2
3
3
dst.(edges(g)) # Destination or ending nodes of arcs
6-element Vector{Int64}:
4
5
4
5
4
5
weights(g) # Weights stored as adjacency matrix
5×5 LinearAlgebra.Adjoint{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}} with 6 stored entries:
⋅ ⋅ ⋅ 6.0 10.0
⋅ ⋅ ⋅ 1.49012e-8 13.0
⋅ ⋅ ⋅ 9.0 16.0
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
Adjacency to Arc list: Iterate over the edge and weight data in the graph to create the arc lists:
i = src.(edges(g))
j = dst.(edges(g))
w = [weights(g)[src(e), dst(e)] for e in edges(g)]
println("i: ", i, "\nj: ", j, "\nw: ", w)
i: [1, 1, 2, 2, 3, 3] j: [4, 5, 4, 5, 4, 5] w: [6.0, 10.0, 1.4901161193847656e-8, 13.0, 9.0, 16.0]
g = SimpleWeightedDiGraph(i, j, w) # Create digraph using arc lists
println.(edges(g));
Edge 1 => 4 with weight 6.0 Edge 1 => 5 with weight 10.0 Edge 2 => 4 with weight 1.4901161193847656e-8 Edge 2 => 5 with weight 13.0 Edge 3 => 4 with weight 9.0 Edge 3 => 5 with weight 16.0
Adjacency/Arc list to Incidence matrix: Output the incidence matrix from the graph, create either from the adjacency matrix or arc lists:
A = Matrix(incidence_matrix(g))
5×6 Matrix{Int64}:
-1 -1 0 0 0 0
0 0 -1 -1 0 0
0 0 0 0 -1 -1
1 0 1 0 1 0
0 1 0 1 0 1
Since the incidence matrix representation does not include arc weight information, the weights need to be determined separately:
w = [weights(g)[src(e), dst(e)] for e in edges(g)]
Matrix([w'; A])
6×6 Matrix{Float64}:
6.0 10.0 1.49012e-8 13.0 9.0 16.0
-1.0 -1.0 0.0 0.0 0.0 0.0
0.0 0.0 -1.0 -1.0 0.0 0.0
0.0 0.0 0.0 0.0 -1.0 -1.0
1.0 0.0 1.0 0.0 1.0 0.0
0.0 1.0 0.0 1.0 0.0 1.0
Ex: (Undirected) Graph:¶
The only change is that each arc in the network is undirected from the starting node ( i ) to the ending node ( j ). The function SimpleWeightedGraph is used to create an undirected network, as opposed to using SimpleWeightedDiGraph to create a directed network:
g = SimpleWeightedGraph(i, j, w) # Create graph using arc lists
println.(edges(g));
Edge 1 => 4 with weight 6.0 Edge 2 => 4 with weight 1.4901161193847656e-8 Edge 3 => 4 with weight 9.0 Edge 1 => 5 with weight 10.0 Edge 2 => 5 with weight 13.0 Edge 3 => 5 with weight 16.0
adjacency_matrix(g) # Both directions represented for each arc
5×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 12 stored entries:
⋅ ⋅ ⋅ 6.0 10.0
⋅ ⋅ ⋅ 1.49012e-8 13.0
⋅ ⋅ ⋅ 9.0 16.0
6.0 1.49012e-8 9.0 ⋅ ⋅
10.0 13.0 16.0 ⋅ ⋅
W
5×5 Matrix{Float64}:
0.0 0.0 0.0 6.0 10.0
0.0 0.0 0.0 1.49012e-8 13.0
0.0 0.0 0.0 9.0 16.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
W + W'
5×5 Matrix{Float64}:
0.0 0.0 0.0 6.0 10.0
0.0 0.0 0.0 1.49012e-8 13.0
0.0 0.0 0.0 9.0 16.0
6.0 1.49012e-8 9.0 0.0 0.0
10.0 13.0 16.0 0.0 0.0
g = SimpleWeightedGraph(W + W') # Create graph using adjacency matrix
println.(edges(g));
Edge 1 => 4 with weight 6.0 Edge 2 => 4 with weight 1.4901161193847656e-8 Edge 3 => 4 with weight 9.0 Edge 1 => 5 with weight 10.0 Edge 2 => 5 with weight 13.0 Edge 3 => 5 with weight 16.0
i, j, w = src.(edges(g)), dst.(edges(g)), [weights(g)[src(e), dst(e)] for e in edges(g)]
println("i: ", i, "\nj: ", j, "\nw: ", w)
i: [1, 2, 3, 1, 2, 3] j: [4, 4, 4, 5, 5, 5] w: [6.0, 1.4901161193847656e-8, 9.0, 10.0, 13.0, 16.0]
A = Matrix(incidence_matrix(g)) # Note: Undirected graph => 1,1 vs digraph => -1,1
5×6 Matrix{Int64}:
1 0 0 1 0 0
0 1 0 0 1 0
0 0 1 0 0 1
1 1 1 0 0 0
0 0 0 1 1 1
2. Dijkstra's Shortest Path Procedure¶
Dijkstra's procedure is used to find the shortest path in a network. It is a greedy procedure that can quickly find the (global) optimal solution.
Ex: Undirected Network¶
Determine the shortest path from node 1 to node 6 in this network:
Pred: $\quad 0 - 3 - 1 - 2 - 4 - 5 $
Path: $\quad 1 \leftarrow 3 \leftarrow 2 \leftarrow 4 \leftarrow 5 \leftarrow 6 $
Cost: $\quad 13$
The network can be created by defining the starting and ending nodes for each arc along with each arc's weight:
using Graphs, SimpleWeightedGraphs
i = [1, 1, 2, 2, 3, 3, 4, 4, 5] # Index of starting node of each arc
j = [2, 3, 3, 4, 4, 5, 5, 6, 6] # Index of ending node of each arc
w = [4, 2, 1, 5, 8, 10, 2, 6, 3] # Weight (e.g., cost, distance, time) of each arc
g = SimpleWeightedGraph(i, j, w) # Create undirected network
{6, 9} undirected simple Int64 graph with Int64 weights
Using Dijkstra's procedure from Graphs (very efficient compared to the code example that follows):
ds = dijkstra_shortest_paths(g, 1) # Shortest paths from node 1 to other all nodes
Graphs.DijkstraState{Int64, Int64}([0, 3, 1, 2, 4, 5], [0, 3, 2, 8, 10, 13], [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], Int64[])
ds.dists # Distance from node 1 to all nodes
6-element Vector{Int64}:
0
3
2
8
10
13
ds.parents # Node predecessors (also called parents)
6-element Vector{Int64}:
0
3
1
2
4
5
enumerate_paths(ds, 6) # Path from node 1 to node 6
6-element Vector{Int64}:
1
3
2
4
5
6
ds.dists[6] # Distance from node 1 to 6
13
To show an implementation of Dijkstra's procedure, we will operate on the adjacency matrix representation of the network. Only inputs are the adjacency matrix and the starting node.
W = Array(float(adjacency_matrix(g)))
6×6 Matrix{Float64}:
0.0 4.0 2.0 0.0 0.0 0.0
4.0 0.0 1.0 5.0 0.0 0.0
2.0 1.0 0.0 8.0 10.0 0.0
0.0 5.0 8.0 0.0 2.0 6.0
0.0 0.0 10.0 2.0 0.0 3.0
0.0 0.0 0.0 6.0 3.0 0.0
replace!(W, 0. => Inf) # Set non-arcs to Inf so never selected
6×6 Matrix{Float64}:
Inf 4.0 2.0 Inf Inf Inf
4.0 Inf 1.0 5.0 Inf Inf
2.0 1.0 Inf 8.0 10.0 Inf
Inf 5.0 8.0 Inf 2.0 6.0
Inf Inf 10.0 2.0 Inf 3.0
Inf Inf Inf 6.0 3.0 Inf
Dijkstra's Shortest Path Procedure:
s = 1 # Starting node
m = size(W, 1) # Number of nodes
d = fill(Inf, m) # Initial distance to all nodes from s = Inf
pred = Int.(zeros(m)) # All predecessors set to 0
S̄ = [1:m;] # Initially, all nodes in set S̄ of unlabeled nodes
d[s] = 0. # Distance from s to s = 0
done = false
while !done
i2 = argmin(d[S̄]) # Find minimum distance unlabeled node (SLOW step!)
i = S̄[i2] # Convert S̄ index (i2) to node index (i)
splice!(S̄,i2) # Label node i (remove from S̄)
pred[d .> d[i] .+ W[i,:]] .= i # Update pred of children of i if dist from i shorter
d = min.(d, d[i] .+ W[i,:]) # Update dist of children of i if dist from i shorter
if isempty(S̄)
done = true
end
end
d, pred
([0.0, 3.0, 2.0, 8.0, 10.0, 13.0], [0, 3, 1, 2, 4, 5])
ds.dists, ds.parents # Using dijkstra_shortest_paths(G, 1)
([0, 3, 2, 8, 10, 13], [0, 3, 1, 2, 4, 5])
Ex: Directed Network¶
The only change is that each arc in the network is directed from the starting node ( i ) to the ending node ( j ). The function SimpleWeightedDiGraph is used to create a directed network, as opposed to using SimpleWeightedGraph to create an undirected network:
i = [1, 1, 2, 2, 3, 3, 4, 4, 5]
j = [2, 3, 3, 4, 4, 5, 5, 6, 6]
w = [4, 2, 1, 5, 8, 10, 2, 6, 3]
g = SimpleWeightedDiGraph(i, j, w)
{6, 9} directed simple Int64 graph with Int64 weights
ds = dijkstra_shortest_paths(g, 1)
println("Min cost path: ", enumerate_paths(ds, 6), ", Cost: ", ds.dists[6])
Min cost path: [1, 2, 4, 5, 6], Cost: 14
Undirected path was [1, 3, 2, 4, 5, 6] at a cost of 13. In the directed graph, there is no edge from 3 to 2, resulting in the different path and higher cost.
Discussion¶
Dijkstra¶
The Dijkstra procedure requires that arcs have all nonnegative lengths on any cycles in the network:
- Negative arc on a cycle would result in a path length $\rightarrow -\infty$.
- It is a “label setting” algorithm since a step to the final solution is made as each node is labeled.
- Can find the longest path (used, e.g., in CPM) by negating all arc lengths
Floyd-Warshall¶
Networks with only some negative arcs on a cycle require slower “label correcting” procedures that repeatedly check for optimality at all nodes or detect a negative cycle.
- Requires $O(m^3)$ via Floyd-Warshall algorithm (cf., $O(m^2)$ Dijkstra), for $m$ nodes.
- Negative arcs used in project scheduling to represent maximum lags between activities.
A* algorithm¶
A* algorithm adds to Dijkstra a heuristic LB estimate of each node's remaining distance to the destination:
- Used in AI search for all types of applications (tic-tac-toe, chess).
- In path planning applications, great circle distance from each node to destination could be used as an LB estimate of the remaining distance.
Time Complexity of Shortest Path Procedures¶
| Procedure | Problem | Time Complexity |
|---|---|---|
| Simplex | LP | $O(2^m)$ |
| Ellipsoid | LP | $O(m^4)$ |
| Hungarian | Transportation | $O(m^3)$ |
| Floyd-Warshall | Shortest Path with Cycles | $O(m^3)$ |
| Dijkstra (linear min) | Shortest Path without Cycles | $O(m^2)$ |
| Dijkstra (Fibonocci heap) | Shortest Path without Cycles | $O(n\,\mbox{log}\,m)$ |
| Number of nodes | $m$ | |
| Number of arcs | $n$ |
The following compares the time complexity of determining the road distance from Raleigh, NC, to Portland, OR, using the U.S. Interstate Highway network:
m = 862 # Number of nodes in Interstate road network
n = 1223 # Number of arcs
@show m^4 # LP for shortest path using ellipsoid
@show m^2 # Dijkstra using linear search for minimum
@show n*log(n) # Dijkstra using Fibonocci heap for minimum
@show (m^4)/(n*log(m)) # Ratio LP to Fibonocci
@show (m^2)/(n*log(m)); # Ratio linear search to Fibonocci
m ^ 4 = 552114385936 m ^ 2 = 743044 n * log(n) = 8694.382991945411 m ^ 4 / (n * log(m)) = 6.6788818050625674e7 m ^ 2 / (n * log(m)) = 89.88541466000086
3. Road Networks¶
National Highway System (NHS) Road Network¶
The files nhs_nodes.csv and nhs_links.csv contain node and link data, respectively, for the U.S. National Highway System (NHS) road network. Each node is defined by an index ( IDX ) and its longitude ( LON ) and latitude ( LAT ), and each link is defined by its starting (source) node index ( SRC ), ending (destination) node index ( DST ), its distance ( DIST ) in miles, and other link attributes like state and county FIPS codes ( STFIP and COFIP ). The network is undirected because all of the arcs can be used for two-way travel. The file roadnetfun.jl contains functions used to work with the road network data.
Ex: Wake County NHS Roads¶
Convert CSV files to data frames and then use STFIP == 37 (NC) and COFIP == 183 (Wake County) to select NHS links, and prune_reindex to create node and link data frames that just contain Wake County roads.
using DataFrames, CSV
include("roadnetfun.jl")
dfnodes = DataFrame(CSV.File("nhs_nodes.csv"))
dflinks = DataFrame(CSV.File("nhs_links.csv"))
dfnodes, dflinks
(369854×3 DataFrame Row │ IDX LON LAT │ Int64 Float64 Float64 ────────┼──────────────────────────── 1 │ 1 -150.289 67.018 2 │ 2 -150.21 67.2448 3 │ 3 -68.0043 46.8593 4 │ 4 -67.9938 46.8208 5 │ 5 -67.9952 46.8308 6 │ 6 -67.9964 46.8142 7 │ 7 -67.9965 46.8125 8 │ 8 -67.9961 46.8112 9 │ 9 -155.111 19.875 10 │ 10 -155.243 19.9873 11 │ 11 -155.813 20.0189 ⋮ │ ⋮ ⋮ ⋮ 369845 │ 369845 -68.01 46.7009 369846 │ 369846 -68.0122 46.694 369847 │ 369847 -68.0161 46.6745 369848 │ 369848 -68.0138 46.686 369849 │ 369849 -68.0158 46.6735 369850 │ 369850 -68.0127 46.689 369851 │ 369851 -68.0087 46.7028 369852 │ 369852 -68.008 46.7052 369853 │ 369853 -68.008 46.7043 369854 │ 369854 -68.0058 46.8537 369833 rows omitted, 387344×10 DataFrame Row │ SRC DST DIST STFIP COFIP URBAN SPEED SIGN N ⋯ │ Int64 Int64 Float64 Int64 Int64 Int64 Float64 String7? S ⋯ ────────┼─────────────────────────────────────────────────────────────────────── 1 │ 1 2 17.9124 2 290 99999 50.0 S11 J ⋯ 2 │ 1 1145 24.2011 2 290 99999 50.0 S11 J 3 │ 2 1522 1.10347 2 290 99999 35.0 S11 J 4 │ 3 19 0.805146 23 3 99998 45.0 U1 m 5 │ 3 369854 0.403536 23 3 99998 35.0 U1 m ⋯ 6 │ 4 5 0.700074 23 3 99999 55.0 U1 m 7 │ 4 6 0.459986 23 3 99999 55.0 U1 m 8 │ 5 13 1.33016 23 3 99999 55.0 U1 m 9 │ 6 7 0.100335 23 3 99999 45.0 U1 m ⋯ 10 │ 7 8 0.090225 23 3 99999 45.0 U1 m 11 │ 8 369842 0.070219 23 3 99999 45.0 U1 m ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 387335 │ 369843 369844 4.8839 23 3 99999 55.0 U1 m 387336 │ 369843 369852 0.180755 23 3 99999 35.0 U1 m ⋯ 387337 │ 369845 369846 0.486638 23 3 99998 35.0 U1 m 387338 │ 369845 369851 0.149565 23 3 99998 35.0 U1 m 387339 │ 369846 369850 0.344989 23 3 99998 25.0 U1 m 387340 │ 369847 369848 0.808124 23 3 99998 25.0 U1 m ⋯ 387341 │ 369847 369849 0.072928 23 3 99999 25.0 U1 m 387342 │ 369848 369850 0.214392 23 3 99998 25.0 U1 m 387343 │ 369851 369853 0.109344 23 3 99999 35.0 U1 m 387344 │ 369852 369853 0.059729 23 3 99999 35.0 U1 m ⋯ 2 columns and 387323 rows omitted)
dfL, dfN = prune_reindex(
filter(r -> (r.STFIP == 37) && (r.COFIP == 183), dflinks), dfnodes)
(446×10 DataFrame Row │ SRC DST DIST STFIP COFIP URBAN SPEED SIGN NAME ⋯ │ Int64 Int64 Float64 Int64 Int64 Int64 Float64 String7? String? ⋯ ─────┼────────────────────────────────────────────────────────────────────────── 1 │ 1 2 0.361991 37 183 73261 70.0 I40 missing ⋯ 2 │ 1 220 3.11067 37 183 73261 70.0 I40 missing 3 │ 2 305 1.09278 37 183 73261 55.0 U70 Us Hwy 4 │ 2 306 0.446956 37 183 73261 70.0 I40 missing 5 │ 3 4 3.39827 37 183 99999 70.0 U1 missing ⋯ 6 │ 4 54 3.85514 37 183 73261 70.0 U1 missing 7 │ 5 11 0.462193 37 183 73261 55.0 U401 missing 8 │ 5 27 0.972324 37 183 73261 55.0 U401 missing 9 │ 6 11 0.798938 37 183 73261 55.0 U401 missing ⋯ 10 │ 6 122 1.25214 37 183 73261 55.0 U401 missing 11 │ 7 8 1.60154 37 183 73261 35.0 S55 missing ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 437 │ 298 299 0.728432 37 183 73261 35.0 U401 missing 438 │ 298 413 0.08392 37 183 99999 35.0 U401 missing ⋯ 439 │ 299 300 0.069977 37 183 73261 35.0 U401 missing 440 │ 300 414 0.82933 37 183 99999 55.0 U401 missing 441 │ 301 415 0.716913 37 183 73261 70.0 U264 missing 442 │ 302 303 0.701596 37 183 73261 70.0 U64 missing ⋯ 443 │ 303 416 2.59257 37 183 73261 70.0 U64 missing 444 │ 303 417 0.507954 37 183 73261 70.0 U264 missing 445 │ 304 415 1.43788 37 183 73261 70.0 U264 missing 446 │ 304 417 0.541973 37 183 73261 70.0 U264 missing ⋯ 2 columns and 425 rows omitted, 417×3 DataFrame Row │ IDX LON LAT │ Int64 Float64 Float64 ─────┼────────────────────────── 1 │ 1 -78.5649 35.6517 2 │ 2 -78.5653 35.6464 3 │ 3 -78.9831 35.6445 4 │ 4 -78.9307 35.6688 5 │ 5 -78.7081 35.6486 6 │ 6 -78.6994 35.6651 7 │ 7 -78.8331 35.6394 8 │ 8 -78.8485 35.6576 9 │ 9 -78.5663 35.6304 10 │ 10 -78.5066 35.7674 11 │ 11 -78.7055 35.6548 ⋮ │ ⋮ ⋮ ⋮ 408 │ 408 -78.6216 35.8017 409 │ 409 -78.5877 35.7878 410 │ 410 -78.5918 35.8207 411 │ 411 -78.5193 36.0089 412 │ 412 -78.5419 35.9653 413 │ 413 -78.4281 35.9375 414 │ 414 -78.4123 35.9578 415 │ 415 -78.2689 35.8204 416 │ 416 -78.275 35.859 417 │ 417 -78.2995 35.834 396 rows omitted)
using Graphs, SimpleWeightedGraphs
g = SimpleWeightedGraph(dfL.SRC, dfL.DST, dfL.DIST)
{417, 446} undirected simple Int64 graph with Float64 weights
using CairoMakie, Logjam.MapTools
fig, ax = makemap(dfN.LON, dfN.LAT)
lines!(ax, x2ln(g, dfN.LON), x2ln(g, dfN.LAT), color=:black, linewidth=.5)
fig
XY = [-78.75864529722959 35.74329140697632;
-78.70352644371913 35.866083946889844]
dfL2, dfN2 = addconnectors(dfL, dfN, XY[:, 1], XY[:, 2])
(452×3 DataFrame Row │ SRC DST DIST │ Int64 Int64 Float64 ─────┼──────────────────────── 1 │ 3 4 0.361991 2 │ 3 222 3.11067 3 │ 4 307 1.09278 4 │ 4 308 0.446956 5 │ 5 6 3.39827 6 │ 6 56 3.85514 7 │ 7 13 0.462193 8 │ 7 29 0.972324 9 │ 8 13 0.798938 10 │ 8 124 1.25214 11 │ 9 10 1.60154 ⋮ │ ⋮ ⋮ ⋮ 443 │ 305 418 2.59257 444 │ 305 419 0.507954 445 │ 306 417 1.43788 446 │ 306 419 0.541973 447 │ 1 73 1.07142 448 │ 1 8 8.24161 449 │ 1 121 4.25563 450 │ 2 190 1.85777 451 │ 2 191 2.36375 452 │ 2 143 0.694065 431 rows omitted, 419×3 DataFrame Row │ IDX LON LAT │ Int64 Float64 Float64 ─────┼────────────────────────── 1 │ 1 -78.7586 35.7433 2 │ 2 -78.7035 35.8661 3 │ 3 -78.5649 35.6517 4 │ 4 -78.5653 35.6464 5 │ 5 -78.9831 35.6445 6 │ 6 -78.9307 35.6688 7 │ 7 -78.7081 35.6486 8 │ 8 -78.6994 35.6651 9 │ 9 -78.8331 35.6394 10 │ 10 -78.8485 35.6576 11 │ 11 -78.5663 35.6304 ⋮ │ ⋮ ⋮ ⋮ 410 │ 410 -78.6216 35.8017 411 │ 411 -78.5877 35.7878 412 │ 412 -78.5918 35.8207 413 │ 413 -78.5193 36.0089 414 │ 414 -78.5419 35.9653 415 │ 415 -78.4281 35.9375 416 │ 416 -78.4123 35.9578 417 │ 417 -78.2689 35.8204 418 │ 418 -78.275 35.859 419 │ 419 -78.2995 35.834 398 rows omitted)
idxN = dfL2.DST[end-5:end]
scatter!(ax, dfN2.LON[idxN], dfN2.LAT[idxN], marker='.', markersize=24, color=:red)
scatter!(ax, XY[:, 1], XY[:, 2])
fig
g = SimpleWeightedGraph(dfL2.SRC, dfL2.DST, dfL2.DIST)
{419, 452} undirected simple Int64 graph with Float64 weights
ds = dijkstra_shortest_paths(g, 1)
println("Min cost path: ", enumerate_paths(ds, 2), ", Cost: ", ds.dists[2])
Min cost path: [1, 73, 340, 120, 123, 338, 68, 69, 175, 168, 320, 39, 321, 170, 169, 146, 145, 192, 364, 144, 143, 2], Cost: 13.689167136202915
dgc(XY[1,:], XY[2,:]) * 1.2
10.834616964839176
Google Maps:
Ex: Road Distances between 100K NC Cities¶
using Logjam.DataTools, DataFrames, CSV
df = filter(r -> (r.STFIP == st2fips(:NC)) && (r.POP > 100_000), usplace())
select!(df, :NAME, :LAT, :LON, :POP)
| Row | NAME | LAT | LON | POP |
|---|---|---|---|---|
| String | Float64 | Float64 | Int64 | |
| 1 | Cary | 35.7814 | -78.819 | 174721 |
| 2 | Charlotte | 35.209 | -80.831 | 874579 |
| 3 | Concord | 35.3921 | -80.6355 | 105240 |
| 4 | Durham | 35.98 | -78.901 | 283506 |
| 5 | Fayetteville | 35.0828 | -78.9735 | 208501 |
| 6 | Greensboro | 36.0951 | -79.827 | 299035 |
| 7 | High Point | 35.9908 | -79.9927 | 114059 |
| 8 | Raleigh | 35.8302 | -78.6415 | 467665 |
| 9 | Wilmington | 34.2092 | -77.8858 | 115451 |
| 10 | Winston-Salem | 36.1029 | -80.2608 | 249545 |
using DataFrames, CSV
using Graphs, SimpleWeightedGraphs
include("roadnetfun.jl")
dfN = DataFrame(CSV.File("nhs_nodes.csv"))
dfL = DataFrame(CSV.File("nhs_links.csv"))
dfL, dfN = prune_reindex(filter(r -> (r.STFIP == st2fips(:NC)), dfL), dfN)
dfL, dfN = addconnectors(dfL, dfN, df.LON, df.LAT)
g = SimpleWeightedGraph(dfL.SRC, dfL.DST, dfL.DIST)
{6893, 7329} undirected simple Int64 graph with Float64 weights
using CairoMakie, Logjam.MapTools
fig, ax = makemap(dfN.LON, dfN.LAT; xexpand=0.05)
lines!(ax, x2ln(g, dfN.LON), x2ln(g, dfN.LAT), color=:black, linewidth=.5)
scatter!(ax, df.LON, df.LAT)
text!(ax, df.LON, df.LAT, text=df.NAME; aligntext(df.LON, df.LAT)..., fontsize=12)
fig
n = nrow(df)
D = hcat([dijkstra_shortest_paths(g, i).dists[1:n] for i in 1:n]...)
10×10 Matrix{Float64}:
0.0 132.816 116.263 16.9174 … 14.4509 139.411 98.21
132.816 0.0 19.5412 134.992 144.194 202.782 81.0455
116.263 19.5412 0.0 118.44 127.641 214.328 62.9545
16.9174 134.992 118.44 0.0 19.36 149.422 83.299
67.853 133.231 133.693 78.0083 73.9081 112.846 116.307
71.8229 92.6482 74.5571 56.9119 … 74.5521 204.614 28.8151
81.9244 79.5903 61.4992 68.2141 85.8543 209.621 20.9838
14.4509 144.194 127.641 19.36 0.0 131.984 100.939
139.411 202.782 214.328 149.422 131.984 0.0 226.573
98.21 81.0455 62.9545 83.299 100.939 226.573 0.0
Distance from Cary (1) to Charlotte (2):
D[1, 2]
132.81556429475137
dgc([df.LON[1], df.LAT[1]],[df.LON[2], df.LAT[2]]) * 1.2
143.86594537973932
Google Maps:
