Net Mod 2: Shortest Paths and CPM

Packages Used: Functions from the following packages are used in this notebook for the first time:

• GraphPlot
• LightGraphs
• SimpleWeightedGraphs

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 1: 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 LightGraphs, 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 adjoint(::SparseArrays.SparseMatrixCSC{Float64, Int64}) with eltype 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

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 2: (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  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

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 3: 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 LightGraphs, 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

Since the location of each node in the network is not specified, the function gplot can be used to try to find a "good" set of node locations. The procedure used starts with a random set of locations for the nodes and then tries to improve it. As a result, it is good to try a variety of different random seed values to find a good layout of the network.

using GraphPlot, Random
Random.seed!(234212)
gplot(G, nodelabel=1:nv(G), edgelabel=[weights(G)[src(e),dst(e)] for e in edges(G)])

Using Dijkstra's procedure from LightGraphs (very efficient compared to the code example that follows):

ds = dijkstra_shortest_paths(G, 1)    # Shortest paths from node 1 to other all nodes

LightGraphs.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 4: 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)
Random.seed!(234212)
gplot(G, nodelabel=1:nv(G), edgelabel=Int.([weights(G)[src(e),dst(e)] for e in edges(G)]))

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

Ex 5: Road Distance from Raleigh to Portland

The files _NetMod-2-Data-node.csv and _NetMod-2-Data-arc.csv contain 862-node and 1223-arc data, respectively, for the U.S. Interstate road network. Each node is defined by its longitude ( x ) and latitude ( y ), and each arc is defined by its starting node index ( i ), ending node index ( j ), and its distance ( d ) in miles. The network is undirected because all of the arcs can be used for two-way travel.

The network can be used to find the shortest distance form Raleigh, NC (node 565) to Portland, OR (node 629):

using DataFrames, CSV
df = DataFrame(CSV.File("Net_Mod-2-Data-arc.csv"))
i, j, d = df.i, df.j, df.d
idxRaleigh = 565
idxPortland = 629
G = SimpleWeightedGraph(i, j, d)
ds = dijkstra_shortest_paths(G, idxRaleigh)
dist = ds.dists[idxPortland]

2937.2999999999997

Using Google Maps, the distance of the shortest travel time between Raleigh and Portland was found to be 2,870 miles, as compared to 2,937. (One reason for the difference could be new Interstate segments that are in Google and are not in the older road data in _NetMod-2-Data.xlsx.)

Since the location of each node in the network is specified, the function plot can be used to display the network instead of using gplot:

using Plots
df = DataFrame(CSV.File("Net_Mod-2-Data-node.csv"))
x, y = df.x, df.y
plot([x[i]'; x[j]'], [y[i]'; y[j]'], label=false)
scatter!([x[idxRaleigh]], [y[idxRaleigh]], label="Raleigh")
scatter!([x[idxPortland]], [y[idxPortland]], label="Portland")

Ex 6: Scheduling with Alternate Routings

In this example, eight different machines are available for the production of a product. There is some flexibility in the set of machines that can be used to perform all of the operations needed to produce the product, and so the product can be produced with several alternate machine routings. Machine 1 is required to be used for the first set of operations and will require two hours to complete (due to processing time and estimated waiting time at the machine). After Machine 1, either Machine 2 or Machine 3 can be used for the next operations, each machine requiring different times because both the number and type of operations are different and the waiting times are different. For each machine, the following table lists its predecessor(s) and time. In the case of more than a single predecessor, either machine listed is a feasible predecessor in a routing.

Machine │ Predecessor   Time
────────┼────────────────────
      1 │ 0                2
      2 │ 1                3
      3 │ 1                5
      4 │ 2                4
      5 │ 2                1
      6 │ [3, 5]           4
      7 │ 4                3
      8 │ [6, 7]           2

A directed network can be used to indicate the predecessor relationships, where the weight of each arc is the time associated with the arc's starting (or outgoing) node:

i = [1,1,2,2,3,5,4,6,7]
j = [2,3,4,5,6,6,7,8,8]
time = [2,3,5,4,1,4,3,2]
t = time[i]
[i'; j'; t']

3×9 Matrix{Int64}:
 1  1  2  2  3  5  4  6  7
 2  3  4  5  6  6  7  8  8
 2  2  3  3  5  1  4  4  3

Machine 8 is the last machine in the network, and as a result, its time has not been used as the weight on an outgoing arc from node 8. An additional node (9) can be added to the network to serve as the ending node for an arc from node 8 with a weight equal to Machine 8's time:

push!(i,8)
push!(j,9)
push!(t,time[end])
[i'; j'; t']

3×10 Matrix{Int64}:
 1  1  2  2  3  5  4  6  7  8
 2  3  4  5  6  6  7  8  8  9
 2  2  3  3  5  1  4  4  3  2

Displaying and then determining the shortest path from node 1 to 9 results in a routing (or path) that visits machines 1, 2, 5, 6, and 8 at a cost of 12 time units:

Random.seed!(378933)
g = SimpleWeightedDiGraph(i, j, t)
gplot(g, nodelabel=1:nv(g), edgelabel=Int.([weights(g)[src(e),dst(e)] for e in edges(g)]))

ds = dijkstra_shortest_paths(g, 1)
path = enumerate_paths(ds, 9)
println("Min time routing: ", enumerate_paths(ds, 9), ", Time: ", ds.dists[9])

Min time routing: [1, 2, 5, 6, 8, 9], Time: 12

3. Critical Path Method (CPM)

The critical path method (CPM) is a procedure for scheduling the activities of a project. Unlike scheduling with alternate routings (Ex 4), all of the predecessors of a node must be complete before the activity represented by the node can occur. As a result, the minimum time to complete all of the activities represented in a network corresponds to the longest path through the network. The activities along this path are the critical activities that, if delayed, would delay the entire project. The longest path can be determined by finding the shortest path through the network using the negative of all arc times.

Ex 7: House Construction Time

The table created below shows the activities associated with the construction of a house. Listed for each activity are the activities that must immediately precede the start of the activity and the time required to complete the activity (in working days). Estimate how many days will be required to construct the house.

using DataFrames
name = ["Excavate","Foundation","Rough wall","Rough exterior plumbing",
    "Roof","Rough electrical work","Rough interior plumbing","Exterior siding",
    "Wall board","Exterior painting","Flooring","Interior painting",
    "Exterior fixtures","Interior fixtures"]
act = [1:14;]   # Note: Assuming all activites numbered consecutivly starting at one
pred = [0,1,2,3,3,3,4,5,[6,7],[4,8],9,9,10,[11,12]]
time = [2,4,10,4,6,7,5,7,8,9,4,5,2,6]
df = DataFrame(Activity = act, Description = name, Predecessor = pred, Time = time)
show(df, show_row_number=false)

14×4 DataFrame
 Activity  Description              Predecessor  Time  
 Int64     String                   Any          Int64 
───────────────────────────────────────────────────────
        1  Excavate                 0                2
        2  Foundation               1                4
        3  Rough wall               2               10
        4  Rough exterior plumbing  3                4
        5  Roof                     3                6
        6  Rough electrical work    3                7
        7  Rough interior plumbing  4                5
        8  Exterior siding          5                7
        9  Wall board               [6, 7]           8
       10  Exterior painting        [4, 8]           9
       11  Flooring                 9                4
       12  Interior painting        9                5
       13  Exterior fixtures        10               2
       14  Interior fixtures        [11, 12]         6

Activity 1 is the single initial activity because it is the only activity without any predecessors. Activities 13 and 14 both represent final activities because they are not the predecessors of any other activities. A directed network can be used to indicate the predecessor relationships, where the weight of each arc is the time associated with the arc's starting (or outgoing) node:

i = [1,2,3,3,3,4,5,6,7, 4, 8, 9, 9,10,11,12]
j = [2,3,4,5,6,7,8,9,9,10,10,11,12,13,14,14]
t = time[i]
[i'; j'; t']

3×16 Matrix{Int64}:
 1  2   3   3   3  4  5  6  7   4   8   9   9  10  11  12
 2  3   4   5   6  7  8  9  9  10  10  11  12  13  14  14
 2  4  10  10  10  4  6  7  5   4   7   8   8   9   4   5

Nodes 13 and 14 represent final activities, and, as a result, their times have not been used as the weight on any outgoing arcs. An additional node (15) can be added to the network to serve as the ending node for arcs from nodes 13 and 14 with weights equal to their activity times:

push!(act,15)
push!(i,13,14)
push!(j,15,15)
push!(t,time[13],time[14])
[i'; j'; t']

3×18 Matrix{Int64}:
 1  2   3   3   3  4  5  6  7   4   8   9   9  10  11  12  13  14
 2  3   4   5   6  7  8  9  9  10  10  11  12  13  14  14  15  15
 2  4  10  10  10  4  6  7  5   4   7   8   8   9   4   5   2   6

Creating the directed graph, displaying it, and then negativing the arc times so that the longest path from node 1 to 15 can be determined results in a critical path of activities in the construction of the house that will require a minimum of 44 working days to complete:

using LightGraphs, SimpleWeightedGraphs
using GraphPlot, Random
G = SimpleWeightedDiGraph(i, j, t)
Random.seed!(3784933)
gplot(G, nodelabel=act, edgelabel=[weights(G)[src(e),dst(e)] for e in edges(G)])

G = SimpleWeightedDiGraph(i, j, -t)
ds = dijkstra_shortest_paths(G, 1)
criticalpath = enumerate_paths(ds, 15)

9-element Vector{Int64}:
  1
  2
  3
  4
  7
  9
 12
 14
 15

mintime = -ds.dists[15]

44

Issue: Multiple Initial Activities

When there is a single initial activity, it can be used as the first node of the critical path. If there are multiple initial activities, a new single dummy initial activity can be used as the predecessor for each, and the arc connecting the dummy node to the real initial nodes can have a time of 0. For example, if an additional initial activity 15 (Curb cut) has no predecessor, takes one day, and is the predecessor of activity 2, dummy activity 16 can be added as its predecessor and the predecessor of activity 1:

push!(name,"Curb cut")
act = [1:15;]   # Note: Assuming all activites numbered consecutivly starting at one
pred = [0,[1,15],2,3,3,3,4,5,[6,7],[4,8],9,9,10,[11,12],0]
time = [2,4,10,4,6,7,5,7,8,9,4,5,2,6,1]
df = DataFrame(Act = act, Description = name, Pred = pred, Time = time)
show(df, show_row_number=false)

15×4 DataFrame
 Act    Description              Pred      Time  
 Int64  String                   Any       Int64 
─────────────────────────────────────────────────
     1  Excavate                 0             2
     2  Foundation               [1, 15]       4
     3  Rough wall               2            10
     4  Rough exterior plumbing  3             4
     5  Roof                     3             6
     6  Rough electrical work    3             7
     7  Rough interior plumbing  4             5
     8  Exterior siding          5             7
     9  Wall board               [6, 7]        8
    10  Exterior painting        [4, 8]        9
    11  Flooring                 9             4
    12  Interior painting        9             5
    13  Exterior fixtures        10            2
    14  Interior fixtures        [11, 12]      6
    15  Curb cut                 0             1

initact = act[pred .== 0]          # Determine initial activities

2-element Vector{Int64}:
  1
 15

if length(initact) > 1             # Add dummy start activity if multiple initial activities
    push!(act, maximum(act) + 1)
    push!(pred, 0)
    push!(time, 0)
    push!(name, "(Start)")
    pred[initact] .= maximum(act)  # Make dummy predecessor of original initial activities
    initact = maximum(act)         # Make dummy the new initital activity
end
i = vcat(pred[pred .!= 0]...)      # Flatten pred nested array to single-level array
j = vcat([repeat([x],length(pred[x])) for x in act[pred .!= 0]]...)   # Duplicate to match
t = time[i]
[i'; j'; t']

3×19 Matrix{Int64}:
 16  1  15  2   3   3   3  4  5  6  7   4   8   9   9  10  11  12  16
  1  2   2  3   4   5   6  7  8  9  9  10  10  11  12  13  14  14  15
  0  2   1  4  10  10  10  4  6  7  5   4   7   8   8   9   4   5   0

finalact = setdiff(j,i)            # Determine final activities

2-element Vector{Int64}:
 13
 14

push!(i, finalact...)              # Always need to add dummy final activity so that an
push!(act, maximum(act) + 1)       # arc can extend from original final activities that
push!(time, 0)                     # includes their times
push!(name, "(Finish)")
push!(pred, finalact)
push!(j, repeat([act[end]], length(finalact))...)
push!(t, time[finalact]...)
finalact = maximum(act)
[i'; j'; t']

3×21 Matrix{Int64}:
 16  1  15  2   3   3   3  4  5  6  7   4   8   9   9  10  11  12  16  13  14
  1  2   2  3   4   5   6  7  8  9  9  10  10  11  12  13  14  14  15  17  17
  0  2   1  4  10  10  10  4  6  7  5   4   7   8   8   9   4   5   0   2   6

df = DataFrame(Act = act, Description = name, Pred = pred, Time = time)
show(df, show_row_number=false)

17×4 DataFrame
 Act    Description              Pred      Time  
 Int64  String                   Any       Int64 
─────────────────────────────────────────────────
     1  Excavate                 16            2
     2  Foundation               [1, 15]       4
     3  Rough wall               2            10
     4  Rough exterior plumbing  3             4
     5  Roof                     3             6
     6  Rough electrical work    3             7
     7  Rough interior plumbing  4             5
     8  Exterior siding          5             7
     9  Wall board               [6, 7]        8
    10  Exterior painting        [4, 8]        9
    11  Flooring                 9             4
    12  Interior painting        9             5
    13  Exterior fixtures        10            2
    14  Interior fixtures        [11, 12]      6
    15  Curb cut                 16            1
    16  (Start)                  0             0
    17  (Finish)                 [13, 14]      0

G = SimpleWeightedDiGraph(i, j, t)
Random.seed!(24111)
gplot(G, nodelabel=act, edgelabel=[weights(G)[src(e),dst(e)] for e in edges(G)])

Determine Critical Path

G = SimpleWeightedDiGraph(i, j, -t)            # New graph with negative arc weights
ds = dijkstra_shortest_paths(G, initact)
criticalpath = enumerate_paths(ds, finalact)

10-element Vector{Int64}:
 16
  1
  2
  3
  4
  7
  9
 12
 14
 17

mintime = -ds.dists[finalact]

44

Determine Earliest Start-Finish and Latest Start-Finish for Activities

A forward breadth-first search starting with the initial node can be used to determine the earliest start time (es) for each activity, which, after adding the time of each activity, gives the earliest finish time (ef). Similarly, a backward breadth-first search starting with the final node can be used to determine the latest finish time (lf) for each activity, which, after subtracting the time of each activity, gives the latest start time (ls).

es = zeros(length(time))           # Initalize earliest start time to 0 for all activities
q = [initact]                      # Initialize queue with intial node
done = false
while !done
    x = popfirst!(q)               # Dequeue: Remove node from front of queue
    y = outneighbors(G, x)         # Get nodes connected by arcs starting from x
    display((x,y))
    if !isempty(y)                 # Final node has no out neighbors
        for yi ∈ y
            if yi ∉ q
                push!(q,yi)        # Enqueue: Add node to end of queue
            end
            es[yi] = max(es[yi], es[x]+time[x])   # If increasing, update earliest start
        end
    end
    if isempty(q)
        done = true
    end
end
ef = es .+ time;

(16, [1, 15])

(1, [2])

(15, [2])

(2, [3])

(3, [4, 5, 6])

(4, [7, 10])

(5, [8])

(6, [9])

(7, [9])

(10, [13])

(8, [10])

(9, [11, 12])

(13, [17])

(10, [13])

(11, [14])

(12, [14])

(17, Int64[])

(13, [17])

(14, [17])

(17, Int64[])

lf = fill(ef[finalact],length(time))  # Initialize latest finish time to final activity ef
q = [finalact]                        # Initialize queue with final node
done = false
while !done
    x = popfirst!(q)                  # Dequeue: Remove node from front of queue
    y = inneighbors(G, x)             # Get nodes connected by arcs ending at x
    display((x,y))
    if !isempty(y)                    # Initial node has no in neighbors
        for yi ∈ y
            if yi ∉ q
                push!(q,yi)           # Enqueue: Add node to end of queue
            end
            lf[yi] = min(lf[yi], lf[x]-time[x]) # If decreasing, update latest finish
        end
    end
    if isempty(q)
        done = true
    end
end
ls = lf .- time;

(17, [13, 14])

(13, [10])

(14, [11, 12])

(10, [4, 8])

(11, [9])

(12, [9])

(4, [3])

(8, [5])

(9, [6, 7])

(3, [2])

(5, [3])

(6, [3])

(7, [4])

(2, [1, 15])

(3, [2])

(4, [3])

(1, [16])

(15, [16])

(2, [1, 15])

(3, [2])

(16, Int64[])

(1, [16])

(15, [16])

(2, [1, 15])

(16, Int64[])

(1, [16])

(15, [16])

(16, Int64[])

Activities that are on the critical path can be identified by having identical earliest and latest start times. Any positive difference between the latest and earliest start time of an activity represents the slack available in starting the activity.

df.EStart = es
df.EFinish = ef
df.LStart = ls
df.LFinish = lf
df.Slack = ls .- es
df.CP = isapprox.(ls,es)
show(df[!,Not(:Description)], show_row_number=false)

17×9 DataFrame
 Act    Pred      Time   EStart   EFinish  LStart   LFinish  Slack    CP    
 Int64  Any       Int64  Float64  Float64  Float64  Float64  Float64  Bool  
────────────────────────────────────────────────────────────────────────────
     1  16            2      0.0      2.0      0.0      2.0      0.0   true
     2  [1, 15]       4      2.0      6.0      2.0      6.0      0.0   true
     3  2            10      6.0     16.0      6.0     16.0      0.0   true
     4  3             4     16.0     20.0     16.0     20.0      0.0   true
     5  3             6     16.0     22.0     20.0     26.0      4.0  false
     6  3             7     16.0     23.0     18.0     25.0      2.0  false
     7  4             5     20.0     25.0     20.0     25.0      0.0   true
     8  5             7     22.0     29.0     26.0     33.0      4.0  false
     9  [6, 7]        8     25.0     33.0     25.0     33.0      0.0   true
    10  [4, 8]        9     29.0     38.0     33.0     42.0      4.0  false
    11  9             4     33.0     37.0     34.0     38.0      1.0  false
    12  9             5     33.0     38.0     33.0     38.0      0.0   true
    13  10            2     38.0     40.0     42.0     44.0      4.0  false
    14  [11, 12]      6     38.0     44.0     38.0     44.0      0.0   true
    15  16            1      0.0      1.0      1.0      2.0      1.0  false
    16  0             0      0.0      0.0      0.0      0.0      0.0   true
    17  [13, 14]      0     44.0     44.0     44.0     44.0      0.0   true

Ex 8: Another CPM Example

Activity 1 requires three days to complete and can start after activities 2, 3, and 4 are complete. Activity 2 requires four days and can begin after activity 5, and activity 3 takes nine days and can begin after both activities 5 and 6 are complete. Activities 4, 5, 6, and 7 each take seven, 12, six, and four days and each can begin after activities 7, 8, 9, and 9 are complete, respectively, and activities 8 and 9 each take ten and eight days to complete, respectively, and can begin at any time. Determine the minimum number of days required to complete all of the activities.

Data Input

using DataFrames
act = [1:9;]
pred = [[2,3,4],5,[5,6],7,8,9,9,0,0]
time = [3,4,9,7,12,6,4,10,8]
df = DataFrame(Act = act, Pred = pred, Time = time)
show(df, show_row_number=false)

9×3 DataFrame
 Act    Pred       Time  
 Int64  Any        Int64 
─────────────────────────
     1  [2, 3, 4]      3
     2  5              4
     3  [5, 6]         9
     4  7              7
     5  8             12
     6  9              6
     7  9              4
     8  0             10
     9  0              8

Add Inital and Final Activities

initact = act[pred .== 0]          # Determine initial activities
if length(initact) > 1             # Add dummy start activity if multiple initial activities
    push!(act, maximum(act) + 1)
    push!(pred, 0)
    push!(time, 0)
    push!(name, "(Start)")
    pred[initact] .= maximum(act)  # Make dummy predecessor of original initial activities
    initact = maximum(act)         # Make dummy the new initital activity
end
i = vcat(pred[pred .!= 0]...)      # Flatten pred nested array to single-level array
j = vcat([repeat([x],length(pred[x])) for x in act[pred .!= 0]]...)   # Duplicate to match
t = time[i]

finalact = setdiff(j,i)            # Determine final activities
push!(i, finalact...)              # Always need to add dummy final activity so that an
push!(act, maximum(act) + 1)       # arc can extend from original final activities that
push!(time, 0)                     # includes their times
push!(name, "(Finish)")
push!(pred, finalact)
push!(j, repeat([act[end]], length(finalact))...)
push!(t, time[finalact]...)
finalact = maximum(act)

df = DataFrame(Act = act, Pred = pred, Time = time)
show(df, show_row_number=false)

G = SimpleWeightedDiGraph(i, j, t)
Random.seed!(24111)
gplot(G, nodelabel=act, edgelabel=[weights(G)[src(e),dst(e)] for e in edges(G)])

11×3 DataFrame
 Act    Pred       Time  
 Int64  Any        Int64 
─────────────────────────
     1  [2, 3, 4]      3
     2  5              4
     3  [5, 6]         9
     4  7              7
     5  8             12
     6  9              6
     7  9              4
     8  10            10
     9  10             8
    10  0              0
    11  [1]            0

Determine Critical Path

G = SimpleWeightedDiGraph(i, j, -t)            # New graph with negative arc weights
ds = dijkstra_shortest_paths(G, initact)
println("Minimum days required to complete all activities: ", -ds.dists[finalact])
println("Critical path activities: ", enumerate_paths(ds, finalact))

Minimum days required to complete all activities: 34
Critical path activities: [10, 8, 5, 3, 1, 11]

Report Activity Slack

es = zeros(length(time))           # Initalize earliest start time to 0 for all activities
q = [initact]                      # Initialize queue with intial node
done = false
while !done
    x = popfirst!(q)               # Dequeue: Remove node from front of queue
    y = outneighbors(G, x)         # Get nodes connected by arcs starting from x
    if !isempty(y)                 # Final node has no out neighbors
        for yi ∈ y
            if yi ∉ q
                push!(q,yi)        # Enqueue: Add node to end of queue
            end
            es[yi] = max(es[yi], es[x]+time[x])   # If increasing, update earliest start
        end
    end
    if isempty(q)
        done = true
    end
end
ef = es .+ time;

lf = fill(ef[finalact],length(time))  # Initialize latest finish time to final activity ef
q = [finalact]                        # Initialize queue with final node
done = false
while !done
    x = popfirst!(q)                  # Dequeue: Remove node from front of queue
    y = inneighbors(G, x)             # Get nodes connected by arcs ending at x
    if !isempty(y)                    # Initial node has no in neighbors
        for yi ∈ y
            if yi ∉ q
                push!(q,yi)           # Enqueue: Add node to end of queue
            end
            lf[yi] = min(lf[yi], lf[x]-time[x]) # If decreasing, update latest finish
        end
    end
    if isempty(q)
        done = true
    end
end
ls = lf .- time;

df.EStart = es
df.EFinish = ef
df.LStart = ls
df.LFinish = lf
df.Slack = ls .- es
df.CP = isapprox.(ls,es)
show(df, show_row_number=false)

11×9 DataFrame
 Act    Pred       Time   EStart   EFinish  LStart   LFinish  Slack    CP    
 Int64  Any        Int64  Float64  Float64  Float64  Float64  Float64  Bool  
─────────────────────────────────────────────────────────────────────────────
     1  [2, 3, 4]      3     31.0     34.0     31.0     34.0      0.0   true
     2  5              4     22.0     26.0     27.0     31.0      5.0  false
     3  [5, 6]         9     22.0     31.0     22.0     31.0      0.0   true
     4  7              7     12.0     19.0     24.0     31.0     12.0  false
     5  8             12     10.0     22.0     10.0     22.0      0.0   true
     6  9              6      8.0     14.0     16.0     22.0      8.0  false
     7  9              4      8.0     12.0     20.0     24.0     12.0  false
     8  10            10      0.0     10.0      0.0     10.0      0.0   true
     9  10             8      0.0      8.0      8.0     16.0      8.0  false
    10  0              0      0.0      0.0      0.0      0.0      0.0   true
    11  [1]            0     34.0     34.0     34.0     34.0      0.0   true

4. 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:

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