Software 8.6.1 (Shortest paths). Dijkstra’s algorithm with a heap on the network of Example 8.3.2, the path recovered from the predecessor labels, and the same problem as a minimum-cost flow LP with the node–arc incidence matrix:

>>> import heapq
>>> import numpy as np
>>> from scipy.optimize import linprog
>>> arcs = {(1, 2): 7, (1, 3): 9, (1, 6): 14, (2, 3): 10, (2, 4): 15,
...         (3, 4): 11, (3, 6): 2, (4, 5): 6, (6, 5): 9}  # length of each arc
>>> def dijkstra(arcs, s):
...     dist, prev, done = {s: 0}, {}, set()
...     heap = [(0, s)]
...     while heap:
...         d, i = heapq.heappop(heap)
...         if i in done: continue
...         done.add(i)  # d is now the shortest distance to i
...         for (a, b), w in arcs.items():
...             if a == i and d + w < dist.get(b, np.inf):
...                 dist[b] = d + w; prev[b] = i
...                 heapq.heappush(heap, (d + w, b))
...     return dist, prev
...
>>> dist, prev = dijkstra(arcs, 1)
>>> dict(sorted(dist.items()))
{1: 0, 2: 7, 3: 9, 4: 20, 5: 20, 6: 11}
>>> path, k = [5], 5
>>> while k != 1: k = prev[k]; path.append(k)
>>> path[::-1]                             # the shortest route from 1 to 5
[1, 3, 6, 5]
>>> nodes = sorted({v for a in arcs for v in a}); E = list(arcs)
>>> N = np.zeros((len(nodes), len(E)))     # node-arc incidence matrix
>>> for e, (a, b) in enumerate(E): N[a - 1, e], N[b - 1, e] = 1, -1
>>> supply = np.zeros(len(nodes)); supply[0], supply[4] = 1, -1  # 1 -> 5
>>> r = linprog(list(arcs.values()), A_eq=N, b_eq=supply, method='highs')
>>> r.fun, [E[e] for e in range(len(E)) if r.x[e] > 0.5]
(20.0, [(1, 3), (3, 6), (6, 5)])

Reading the output. The distances from node \(1\) are \(0, 7, 9, 20, 20, 11\), and the shortest route to \(5\) is \(1, 3, 6, 5\). The LP, which knows nothing about paths, returns a flow of one unit on exactly the arcs \((1, 3)\), \((3, 6)\), \((6, 5)\), with cost \(20\): integral, as total unimodularity guarantees.

MATLAB/Octave: shortestpath(digraph(s, t, w), 1, 5).

Software 8.6.2 (Augmenting paths). The Edmonds–Karp algorithm: breadth-first search in the residual network, augmenting along the shortest path found, until none remains; it returns the flow value and the set \(S\) of nodes still reachable, and the capacity of that cut:

>>> cap = {(1, 2): 10, (1, 3): 5, (1, 6): 15, (2, 3): 4, (2, 4): 9,
...        (3, 4): 15, (3, 6): 4, (4, 5): 10, (6, 5): 15}
>>> def max_flow(cap, s, t):  # Edmonds-Karp: shortest augmenting paths
...     R = dict(cap)                # residual capacities, with reverses
...     R.update({(b, a): 0 for (a, b) in cap if (b, a) not in cap})
...     total = 0
...     while True:
...         prev, queue = {s: None}, [s]
...         for i in queue:  # breadth-first search in the residual
...             for (a, b), c in R.items():
...                 if a == i and c > 0 and b not in prev:
...                     prev[b] = a; queue.append(b)
...         if t not in prev:              # no augmenting path is left
...             return total, set(prev)    # (prev: the nodes reachable)
...         path, k = [], t
...         while prev[k] is not None:
...             path.append((prev[k], k)); k = prev[k]
...         delta = min(R[e] for e in path)
...         for (a, b) in path: R[(a, b)] -= delta; R[(b, a)] += delta
...         total += delta
...         print('augment', delta, 'along', [a for a, _ in path[::-1]] + [t])
...
>>> value, S = max_flow(cap, 1, 5)
augment 15 along [1, 6, 5]
augment 9 along [1, 2, 4, 5]
augment 1 along [1, 3, 4, 5]
>>> value, S, sum(c for (a, b), c in cap.items() if a in S and b not in S)
(25, {1, 2, 3, 4, 6}, 25)

Reading the output. Three augmentations, of \(15\), \(9\) and \(1\), give a flow of \(25\). At the end nodes \(1, 2, 3, 4, 6\) are reachable and \(5\) is not; the arcs leaving that set have capacity \(25\), equal to the flow, which proves it maximal.

MATLAB/Octave: maxflow(digraph(s, t, cap), 1, 5).

Software 8.6.3 (CPM and PERT). Expected times and variances, the forward and backward passes, the critical path and floats, the normal approximation to \(P(\text {finish by } 21)\), and a simulation of \(200\,000\) projects with beta-PERT durations:

>>> from scipy.stats import norm
>>> acts = {'A': ([], 2, 4, 6), 'B': ([], 3, 5, 9), 'C': (['A'], 4, 5, 7),
...         'D': (['A'], 4, 6, 10), 'E': (['B', 'C'], 4, 5, 7),
...         'F': (['D'], 3, 4, 8), 'G': (['E'], 3, 5, 8)}  # (preds, a, m, b)
>>> te = {k: (a + 4*m + b) / 6 for k, (_, a, m, b) in acts.items()}
>>> var = {k: ((b - a) / 6)**2 for k, (_, a, m, b) in acts.items()}
>>> ES, EF = {}, {}
>>> for k in acts:                         # forward pass (already in order)
...     ES[k] = max([EF[p] for p in acts[k][0]], default=0)
...     EF[k] = ES[k] + te[k]
...
>>> T = max(EF.values()); LF, LS = {}, {}
>>> for k in reversed(list(acts)):         # backward pass
...     succ = [j for j in acts if k in acts[j][0]]
...     LF[k] = min([LS[j] for j in succ], default=T)
...     LS[k] = LF[k] - te[k]
...
>>> crit = [k for k in acts if abs(LS[k] - ES[k]) < 1e-9]
>>> round(T, 3), crit
(19.5, ['A', 'C', 'E', 'G'])
>>> {k: round(LS[k] - ES[k], 3) for k in acts}   # total floats
{'A': 0.0, 'B': 3.833, 'C': 0.0, 'D': 4.667, 'E': 0.0, 'F': 4.667, 'G': 0.0}
>>> sd = np.sqrt(sum(var[k] for k in crit)); round(sd, 4)
1.2802
>>> round(norm.cdf((21 - T) / sd), 4)      # P(project done by time 21)
0.8793
>>> rng = np.random.default_rng(0)         # simulate with beta-PERT durations
>>> def beta_pert(a, m, b, size):          # the beta law with mean (a+4m+b)/6
...     al, be = 1 + 4*(m - a)/(b - a), 1 + 4*(b - m)/(b - a)
...     return a + (b - a) * rng.beta(al, be, size)
...
>>> sims = {k: beta_pert(a, m, b, 200000) for k, (_, a, m, b) in acts.items()}
>>> fin = {}
>>> for k in acts:                         # finish time of each activity
...     pre = [fin[p] for p in acts[k][0]]
...     fin[k] = (np.maximum.reduce(pre) if pre else 0) + sims[k]
...
>>> done = np.maximum.reduce(list(fin.values()))
>>> round(done.mean(), 3), round((done <= 21).mean(), 4)
(19.507, 0.8477)

Reading the output. \(T = 19.5\) along \(A, C, E, G\), with floats \(3.83\) for \(B\) and \(4.67\) for \(D\) and \(F\); \(\sigma = 1.28\), and the normal approximation gives \(0.879\). The simulation’s mean duration, \(19.51\), is close to \(T\) because the other paths have large floats, but the probability of finishing by \(21\) is \(0.848\): the beta durations are skewed, and the path \(B \to E \to G\) is occasionally the longest, so the PERT formula overstates the chance by \(3\) percentage points.

MATLAB/Octave: the same passes over a sorted activity list; normcdf for \(\Phi \).