restructured networksimplex.py and added test_networksimplex.py (#4685)

* restructured networksimplex.py and added test_networksimplex.py

* formatting changed

* formatting changed

* formatting changed

* removed extra underscore

* added fixture in test_networksimplex.py

* formatting changed

* initilized -> initialized.

* fix typo in comment

Co-authored-by: Ross Barnowski <rossbar@berkeley.edu>
Co-authored-by: Dan Schult <dschult@colgate.edu>

2 files changed, 768 insertions, 305 deletions

diff --git a/networkx/algorithms/flow/networksimplex.py b/networkx/algorithms/flow/networksimplex.py
index 561dbb8a..0de9be80 100644
--- a/networkx/algorithms/flow/networksimplex.py
+++ b/networkx/algorithms/flow/networksimplex.py
@@ -10,6 +10,321 @@ import networkx as nx
 from networkx.utils import not_implemented_for
 
 
+class _DataEssentialsAndFunctions:
+    def __init__(
+        self, G, multigraph, demand="demand", capacity="capacity", weight="weight"
+    ):
+
+        # Number all nodes and edges and hereafter reference them using ONLY their numbers
+        self.node_list = list(G)  # nodes
+        self.node_indices = {u: i for i, u in enumerate(self.node_list)}  # node indices
+        self.node_demands = [
+            G.nodes[u].get(demand, 0) for u in self.node_list
+        ]  # node demands
+
+        self.edge_sources = []  # edge sources
+        self.edge_targets = []  # edge targets
+        if multigraph:
+            self.edge_keys = []  # edge keys
+        self.edge_indices = {}  # edge indices
+        self.edge_capacities = []  # edge capacities
+        self.edge_weights = []  # edge weights
+
+        if not multigraph:
+            edges = G.edges(data=True)
+        else:
+            edges = G.edges(data=True, keys=True)
+
+        inf = float("inf")
+        edges = (e for e in edges if e[0] != e[1] and e[-1].get(capacity, inf) != 0)
+        for i, e in enumerate(edges):
+            self.edge_sources.append(self.node_indices[e[0]])
+            self.edge_targets.append(self.node_indices[e[1]])
+            if multigraph:
+                self.edge_keys.append(e[2])
+            self.edge_indices[e[:-1]] = i
+            self.edge_capacities.append(e[-1].get(capacity, inf))
+            self.edge_weights.append(e[-1].get(weight, 0))
+
+        # spanning tree specific data to be initialized
+
+        self.edge_count = None  # number of edges
+        self.edge_flow = None  # edge flows
+        self.node_potentials = None  # node potentials
+        self.parent = None  # parent nodes
+        self.parent_edge = None  # edges to parents
+        self.subtree_size = None  # subtree sizes
+        self.next_node_dft = None  # next nodes in depth-first thread
+        self.prev_node_dft = None  # previous nodes in depth-first thread
+        self.last_descendent_dft = None  # last descendants in depth-first thread
+        self._spanning_tree_initialized = (
+            False  # False until initialize_spanning_tree() is called
+        )
+
+    def initialize_spanning_tree(self, n, faux_inf):
+        self.edge_count = len(self.edge_indices)  # number of edges
+        self.edge_flow = list(
+            chain(repeat(0, self.edge_count), (abs(d) for d in self.node_demands))
+        )  # edge flows
+        self.node_potentials = [
+            faux_inf if d <= 0 else -faux_inf for d in self.node_demands
+        ]  # node potentials
+        self.parent = list(chain(repeat(-1, n), [None]))  # parent nodes
+        self.parent_edge = list(
+            range(self.edge_count, self.edge_count + n)
+        )  # edges to parents
+        self.subtree_size = list(chain(repeat(1, n), [n + 1]))  # subtree sizes
+        self.next_node_dft = list(
+            chain(range(1, n), [-1, 0])
+        )  # next nodes in depth-first thread
+        self.prev_node_dft = list(range(-1, n))  # previous nodes in depth-first thread
+        self.last_descendent_dft = list(
+            chain(range(n), [n - 1])
+        )  # last descendants in depth-first thread
+        self._spanning_tree_initialized = True  # True only if all the assignments pass
+
+    def find_apex(self, p, q):
+        """
+        Find the lowest common ancestor of nodes p and q in the spanning tree.
+        """
+        size_p = self.subtree_size[p]
+        size_q = self.subtree_size[q]
+        while True:
+            while size_p < size_q:
+                p = self.parent[p]
+                size_p = self.subtree_size[p]
+            while size_p > size_q:
+                q = self.parent[q]
+                size_q = self.subtree_size[q]
+            if size_p == size_q:
+                if p != q:
+                    p = self.parent[p]
+                    size_p = self.subtree_size[p]
+                    q = self.parent[q]
+                    size_q = self.subtree_size[q]
+                else:
+                    return p
+
+    def trace_path(self, p, w):
+        """
+        Returns the nodes and edges on the path from node p to its ancestor w.
+        """
+        Wn = [p]
+        We = []
+        while p != w:
+            We.append(self.parent_edge[p])
+            p = self.parent[p]
+            Wn.append(p)
+        return Wn, We
+
+    def find_cycle(self, i, p, q):
+        """
+        Returns the nodes and edges on the cycle containing edge i == (p, q)
+        when the latter is added to the spanning tree.
+
+        The cycle is oriented in the direction from p to q.
+        """
+        w = self.find_apex(p, q)
+        Wn, We = self.trace_path(p, w)
+        Wn.reverse()
+        We.reverse()
+        if We != [i]:
+            We.append(i)
+        WnR, WeR = self.trace_path(q, w)
+        del WnR[-1]
+        Wn += WnR
+        We += WeR
+        return Wn, We
+
+    def augment_flow(self, Wn, We, f):
+        """
+        Augment f units of flow along a cycle represented by Wn and We.
+        """
+        for i, p in zip(We, Wn):
+            if self.edge_sources[i] == p:
+                self.edge_flow[i] += f
+            else:
+                self.edge_flow[i] -= f
+
+    def trace_subtree(self, p):
+        """
+        Yield the nodes in the subtree rooted at a node p.
+        """
+        yield p
+        l = self.last_descendent_dft[p]
+        while p != l:
+            p = self.next_node_dft[p]
+            yield p
+
+    def remove_edge(self, s, t):
+        """
+        Remove an edge (s, t) where parent[t] == s from the spanning tree.
+        """
+        size_t = self.subtree_size[t]
+        prev_t = self.prev_node_dft[t]
+        last_t = self.last_descendent_dft[t]
+        next_last_t = self.next_node_dft[last_t]
+        # Remove (s, t).
+        self.parent[t] = None
+        self.parent_edge[t] = None
+        # Remove the subtree rooted at t from the depth-first thread.
+        self.next_node_dft[prev_t] = next_last_t
+        self.prev_node_dft[next_last_t] = prev_t
+        self.next_node_dft[last_t] = t
+        self.prev_node_dft[t] = last_t
+        # Update the subtree sizes and last descendants of the (old) acenstors
+        # of t.
+        while s is not None:
+            self.subtree_size[s] -= size_t
+            if self.last_descendent_dft[s] == last_t:
+                self.last_descendent_dft[s] = prev_t
+            s = self.parent[s]
+
+    def make_root(self, q):
+        """
+        Make a node q the root of its containing subtree.
+        """
+        ancestors = []
+        while q is not None:
+            ancestors.append(q)
+            q = self.parent[q]
+        ancestors.reverse()
+        for p, q in zip(ancestors, islice(ancestors, 1, None)):
+            size_p = self.subtree_size[p]
+            last_p = self.last_descendent_dft[p]
+            prev_q = self.prev_node_dft[q]
+            last_q = self.last_descendent_dft[q]
+            next_last_q = self.next_node_dft[last_q]
+            # Make p a child of q.
+            self.parent[p] = q
+            self.parent[q] = None
+            self.parent_edge[p] = self.parent_edge[q]
+            self.parent_edge[q] = None
+            self.subtree_size[p] = size_p - self.subtree_size[q]
+            self.subtree_size[q] = size_p
+            # Remove the subtree rooted at q from the depth-first thread.
+            self.next_node_dft[prev_q] = next_last_q
+            self.prev_node_dft[next_last_q] = prev_q
+            self.next_node_dft[last_q] = q
+            self.prev_node_dft[q] = last_q
+            if last_p == last_q:
+                self.last_descendent_dft[p] = prev_q
+                last_p = prev_q
+            # Add the remaining parts of the subtree rooted at p as a subtree
+            # of q in the depth-first thread.
+            self.prev_node_dft[p] = last_q
+            self.next_node_dft[last_q] = p
+            self.next_node_dft[last_p] = q
+            self.prev_node_dft[q] = last_p
+            self.last_descendent_dft[q] = last_p
+
+    def add_edge(self, i, p, q):
+        """
+        Add an edge (p, q) to the spanning tree where q is the root of a subtree.
+        """
+        last_p = self.last_descendent_dft[p]
+        next_last_p = self.next_node_dft[last_p]
+        size_q = self.subtree_size[q]
+        last_q = self.last_descendent_dft[q]
+        # Make q a child of p.
+        self.parent[q] = p
+        self.parent_edge[q] = i
+        # Insert the subtree rooted at q into the depth-first thread.
+        self.next_node_dft[last_p] = q
+        self.prev_node_dft[q] = last_p
+        self.prev_node_dft[next_last_p] = last_q
+        self.next_node_dft[last_q] = next_last_p
+        # Update the subtree sizes and last descendants of the (new) ancestors
+        # of q.
+        while p is not None:
+            self.subtree_size[p] += size_q
+            if self.last_descendent_dft[p] == last_p:
+                self.last_descendent_dft[p] = last_q
+            p = self.parent[p]
+
+    def update_potentials(self, i, p, q):
+        """
+        Update the potentials of the nodes in the subtree rooted at a node
+        q connected to its parent p by an edge i.
+        """
+        if q == self.edge_targets[i]:
+            d = self.node_potentials[p] - self.edge_weights[i] - self.node_potentials[q]
+        else:
+            d = self.node_potentials[p] + self.edge_weights[i] - self.node_potentials[q]
+        for q in self.trace_subtree(q):
+            self.node_potentials[q] += d
+
+    def reduced_cost(self, i):
+        """Returns the reduced cost of an edge i."""
+        c = (
+            self.edge_weights[i]
+            - self.node_potentials[self.edge_sources[i]]
+            + self.node_potentials[self.edge_targets[i]]
+        )
+        return c if self.edge_flow[i] == 0 else -c
+
+    def find_entering_edges(self):
+        """Yield entering edges until none can be found."""
+        if self.edge_count == 0:
+            return
+
+        # Entering edges are found by combining Dantzig's rule and Bland's
+        # rule. The edges are cyclically grouped into blocks of size B. Within
+        # each block, Dantzig's rule is applied to find an entering edge. The
+        # blocks to search is determined following Bland's rule.
+        B = int(ceil(sqrt(self.edge_count)))  # pivot block size
+        M = (self.edge_count + B - 1) // B  # number of blocks needed to cover all edges
+        m = 0  # number of consecutive blocks without eligible
+        # entering edges
+        f = 0  # first edge in block
+        while m < M:
+            # Determine the next block of edges.
+            l = f + B
+            if l <= self.edge_count:
+                edges = range(f, l)
+            else:
+                l -= self.edge_count
+                edges = chain(range(f, self.edge_count), range(l))
+            f = l
+            # Find the first edge with the lowest reduced cost.
+            i = min(edges, key=self.reduced_cost)
+            c = self.reduced_cost(i)
+            if c >= 0:
+                # No entering edge found in the current block.
+                m += 1
+            else:
+                # Entering edge found.
+                if self.edge_flow[i] == 0:
+                    p = self.edge_sources[i]
+                    q = self.edge_targets[i]
+                else:
+                    p = self.edge_targets[i]
+                    q = self.edge_sources[i]
+                yield i, p, q
+                m = 0
+        # All edges have nonnegative reduced costs. The current flow is
+        # optimal.
+
+    def residual_capacity(self, i, p):
+        """Returns the residual capacity of an edge i in the direction away
+        from its endpoint p.
+        """
+        return (
+            self.edge_capacities[i] - self.edge_flow[i]
+            if self.edge_sources[i] == p
+            else self.edge_flow[i]
+        )
+
+    def find_leaving_edge(self, Wn, We):
+        """Returns the leaving edge in a cycle represented by Wn and We."""
+        j, s = min(
+            zip(reversed(We), reversed(Wn)),
+            key=lambda i_p: self.residual_capacity(*i_p),
+        )
+        t = self.edge_targets[j] if self.edge_sources[j] == s else self.edge_sources[j]
+        return j, s, t
+
+
 @not_implemented_for("undirected")
 def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
     r"""Find a minimum cost flow satisfying all demands in digraph G.
@@ -180,43 +495,23 @@ def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
     if len(G) == 0:
         raise nx.NetworkXError("graph has no nodes")
 
-    # Number all nodes and edges and hereafter reference them using ONLY their
-    # numbers
-
-    N = list(G)  # nodes
-    I = {u: i for i, u in enumerate(N)}  # node indices
-    D = [G.nodes[u].get(demand, 0) for u in N]  # node demands
-
-    inf = float("inf")
-    for p, b in zip(N, D):
-        if abs(b) == inf:
-            raise nx.NetworkXError(f"node {p!r} has infinite demand")
-
     multigraph = G.is_multigraph()
-    S = []  # edge sources
-    T = []  # edge targets
-    if multigraph:
-        K = []  # edge keys
-    E = {}  # edge indices
-    U = []  # edge capacities
-    C = []  # edge weights
 
-    if not multigraph:
-        edges = G.edges(data=True)
-    else:
-        edges = G.edges(data=True, keys=True)
-    edges = (e for e in edges if e[0] != e[1] and e[-1].get(capacity, inf) != 0)
-    for i, e in enumerate(edges):
-        S.append(I[e[0]])
-        T.append(I[e[1]])
-        if multigraph:
-            K.append(e[2])
-        E[e[:-1]] = i
-        U.append(e[-1].get(capacity, inf))
-        C.append(e[-1].get(weight, 0))
+    # extracting data essential to problem
+    DEAF = _DataEssentialsAndFunctions(
+        G, multigraph, demand=demand, capacity=capacity, weight=weight
+    )
 
-    for e, c in zip(E, C):
-        if abs(c) == inf:
+    ###########################################################################
+    # Quick Error Detection
+    ###########################################################################
+
+    inf = float("inf")
+    for u, d in zip(DEAF.node_list, DEAF.node_demands):
+        if abs(d) == inf:
+            raise nx.NetworkXError(f"node {u!r} has infinite demand")
+    for e, w in zip(DEAF.edge_indices, DEAF.edge_weights):
+        if abs(w) == inf:
             raise nx.NetworkXError(f"edge {e!r} has infinite weight")
     if not multigraph:
         edges = nx.selfloop_edges(G, data=True)
@@ -227,13 +522,13 @@ def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
             raise nx.NetworkXError(f"edge {e[:-1]!r} has infinite weight")
 
     ###########################################################################
-    # Quick infeasibility detection
+    # Quick Infeasibility Detection
     ###########################################################################
 
-    if sum(D) != 0:
+    if sum(DEAF.node_demands) != 0:
         raise nx.NetworkXUnfeasible("total node demand is not zero")
-    for e, u in zip(E, U):
-        if u < 0:
+    for e, c in zip(DEAF.edge_indices, DEAF.edge_capacities):
+        if c < 0:
             raise nx.NetworkXUnfeasible(f"edge {e!r} has negative capacity")
     if not multigraph:
         edges = nx.selfloop_edges(G, data=True)
@@ -252,293 +547,65 @@ def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
     # spanning tree of the network simplex method. The new edges will used to
     # trivially satisfy the node demands and create an initial strongly
     # feasible spanning tree.
-    n = len(N)  # number of nodes
-    for p, d in enumerate(D):
+    for i, d in enumerate(DEAF.node_demands):
         # Must be greater-than here. Zero-demand nodes must have
-        # edges pointing towards the root to ensure strong
-        # feasibility.
+        # edges pointing towards the root to ensure strong feasibility.
         if d > 0:
-            S.append(-1)
-            T.append(p)
+            DEAF.edge_sources.append(-1)
+            DEAF.edge_targets.append(i)
         else:
-            S.append(p)
-            T.append(-1)
+            DEAF.edge_sources.append(i)
+            DEAF.edge_targets.append(-1)
     faux_inf = (
         3
         * max(
             chain(
-                [sum(u for u in U if u < inf), sum(abs(c) for c in C)],
-                (abs(d) for d in D),
+                [
+                    sum(c for c in DEAF.edge_capacities if c < inf),
+                    sum(abs(w) for w in DEAF.edge_weights),
+                ],
+                (abs(d) for d in DEAF.node_demands),
             )
         )
         or 1
     )
-    C.extend(repeat(faux_inf, n))
-    U.extend(repeat(faux_inf, n))
+
+    n = len(DEAF.node_list)  # number of nodes
+    DEAF.edge_weights.extend(repeat(faux_inf, n))
+    DEAF.edge_capacities.extend(repeat(faux_inf, n))
 
     # Construct the initial spanning tree.
-    e = len(E)  # number of edges
-    x = list(chain(repeat(0, e), (abs(d) for d in D)))  # edge flows
-    pi = [faux_inf if d <= 0 else -faux_inf for d in D]  # node potentials
-    parent = list(chain(repeat(-1, n), [None]))  # parent nodes
-    edge = list(range(e, e + n))  # edges to parents
-    size = list(chain(repeat(1, n), [n + 1]))  # subtree sizes
-    next = list(chain(range(1, n), [-1, 0]))  # next nodes in depth-first thread
-    prev = list(range(-1, n))  # previous nodes in depth-first thread
-    last = list(chain(range(n), [n - 1]))  # last descendants in depth-first thread
+    DEAF.initialize_spanning_tree(n, faux_inf)
 
     ###########################################################################
     # Pivot loop
     ###########################################################################
 
-    def reduced_cost(i):
-        """Returns the reduced cost of an edge i."""
-        c = C[i] - pi[S[i]] + pi[T[i]]
-        return c if x[i] == 0 else -c
-
-    def find_entering_edges():
-        """Yield entering edges until none can be found."""
-        if e == 0:
-            return
-
-        # Entering edges are found by combining Dantzig's rule and Bland's
-        # rule. The edges are cyclically grouped into blocks of size B. Within
-        # each block, Dantzig's rule is applied to find an entering edge. The
-        # blocks to search is determined following Bland's rule.
-        B = int(ceil(sqrt(e)))  # pivot block size
-        M = (e + B - 1) // B  # number of blocks needed to cover all edges
-        m = 0  # number of consecutive blocks without eligible
-        # entering edges
-        f = 0  # first edge in block
-        while m < M:
-            # Determine the next block of edges.
-            l = f + B
-            if l <= e:
-                edges = range(f, l)
-            else:
-                l -= e
-                edges = chain(range(f, e), range(l))
-            f = l
-            # Find the first edge with the lowest reduced cost.
-            i = min(edges, key=reduced_cost)
-            c = reduced_cost(i)
-            if c >= 0:
-                # No entering edge found in the current block.
-                m += 1
-            else:
-                # Entering edge found.
-                if x[i] == 0:
-                    p = S[i]
-                    q = T[i]
-                else:
-                    p = T[i]
-                    q = S[i]
-                yield i, p, q
-                m = 0
-        # All edges have nonnegative reduced costs. The current flow is
-        # optimal.
-
-    def find_apex(p, q):
-        """Find the lowest common ancestor of nodes p and q in the spanning
-        tree.
-        """
-        size_p = size[p]
-        size_q = size[q]
-        while True:
-            while size_p < size_q:
-                p = parent[p]
-                size_p = size[p]
-            while size_p > size_q:
-                q = parent[q]
-                size_q = size[q]
-            if size_p == size_q:
-                if p != q:
-                    p = parent[p]
-                    size_p = size[p]
-                    q = parent[q]
-                    size_q = size[q]
-                else:
-                    return p
-
-    def trace_path(p, w):
-        """Returns the nodes and edges on the path from node p to its ancestor
-        w.
-        """
-        Wn = [p]
-        We = []
-        while p != w:
-            We.append(edge[p])
-            p = parent[p]
-            Wn.append(p)
-        return Wn, We
-
-    def find_cycle(i, p, q):
-        """Returns the nodes and edges on the cycle containing edge i == (p, q)
-        when the latter is added to the spanning tree.
-
-        The cycle is oriented in the direction from p to q.
-        """
-        w = find_apex(p, q)
-        Wn, We = trace_path(p, w)
-        Wn.reverse()
-        We.reverse()
-        if We != [i]:
-            We.append(i)
-        WnR, WeR = trace_path(q, w)
-        del WnR[-1]
-        Wn += WnR
-        We += WeR
-        return Wn, We
-
-    def residual_capacity(i, p):
-        """Returns the residual capacity of an edge i in the direction away
-        from its endpoint p.
-        """
-        return U[i] - x[i] if S[i] == p else x[i]
-
-    def find_leaving_edge(Wn, We):
-        """Returns the leaving edge in a cycle represented by Wn and We."""
-        j, s = min(
-            zip(reversed(We), reversed(Wn)), key=lambda i_p: residual_capacity(*i_p)
-        )
-        t = T[j] if S[j] == s else S[j]
-        return j, s, t
-
-    def augment_flow(Wn, We, f):
-        """Augment f units of flow along a cycle represented by Wn and We."""
-        for i, p in zip(We, Wn):
-            if S[i] == p:
-                x[i] += f
-            else:
-                x[i] -= f
-
-    def trace_subtree(p):
-        """Yield the nodes in the subtree rooted at a node p."""
-        yield p
-        l = last[p]
-        while p != l:
-            p = next[p]
-            yield p
-
-    def remove_edge(s, t):
-        """Remove an edge (s, t) where parent[t] == s from the spanning tree."""
-        size_t = size[t]
-        prev_t = prev[t]
-        last_t = last[t]
-        next_last_t = next[last_t]
-        # Remove (s, t).
-        parent[t] = None
-        edge[t] = None
-        # Remove the subtree rooted at t from the depth-first thread.
-        next[prev_t] = next_last_t
-        prev[next_last_t] = prev_t
-        next[last_t] = t
-        prev[t] = last_t
-        # Update the subtree sizes and last descendants of the (old) acenstors
-        # of t.
-        while s is not None:
-            size[s] -= size_t
-            if last[s] == last_t:
-                last[s] = prev_t
-            s = parent[s]
-
-    def make_root(q):
-        """Make a node q the root of its containing subtree."""
-        ancestors = []
-        while q is not None:
-            ancestors.append(q)
-            q = parent[q]
-        ancestors.reverse()
-        for p, q in zip(ancestors, islice(ancestors, 1, None)):
-            size_p = size[p]
-            last_p = last[p]
-            prev_q = prev[q]
-            last_q = last[q]
-            next_last_q = next[last_q]
-            # Make p a child of q.
-            parent[p] = q
-            parent[q] = None
-            edge[p] = edge[q]
-            edge[q] = None
-            size[p] = size_p - size[q]
-            size[q] = size_p
-            # Remove the subtree rooted at q from the depth-first thread.
-            next[prev_q] = next_last_q
-            prev[next_last_q] = prev_q
-            next[last_q] = q
-            prev[q] = last_q
-            if last_p == last_q:
-                last[p] = prev_q
-                last_p = prev_q
-            # Add the remaining parts of the subtree rooted at p as a subtree
-            # of q in the depth-first thread.
-            prev[p] = last_q
-            next[last_q] = p
-            next[last_p] = q
-            prev[q] = last_p
-            last[q] = last_p
-
-    def add_edge(i, p, q):
-        """Add an edge (p, q) to the spanning tree where q is the root of a
-        subtree.
-        """
-        last_p = last[p]
-        next_last_p = next[last_p]
-        size_q = size[q]
-        last_q = last[q]
-        # Make q a child of p.
-        parent[q] = p
-        edge[q] = i
-        # Insert the subtree rooted at q into the depth-first thread.
-        next[last_p] = q
-        prev[q] = last_p
-        prev[next_last_p] = last_q
-        next[last_q] = next_last_p
-        # Update the subtree sizes and last descendants of the (new) ancestors
-        # of q.
-        while p is not None:
-            size[p] += size_q
-            if last[p] == last_p:
-                last[p] = last_q
-            p = parent[p]
-
-    def update_potentials(i, p, q):
-        """Update the potentials of the nodes in the subtree rooted at a node
-        q connected to its parent p by an edge i.
-        """
-        if q == T[i]:
-            d = pi[p] - C[i] - pi[q]
-        else:
-            d = pi[p] + C[i] - pi[q]
-        for q in trace_subtree(q):
-            pi[q] += d
-
-    # Pivot loop
-    for i, p, q in find_entering_edges():
-        Wn, We = find_cycle(i, p, q)
-        j, s, t = find_leaving_edge(Wn, We)
-        augment_flow(Wn, We, residual_capacity(j, s))
+    for i, p, q in DEAF.find_entering_edges():
+        Wn, We = DEAF.find_cycle(i, p, q)
+        j, s, t = DEAF.find_leaving_edge(Wn, We)
+        DEAF.augment_flow(Wn, We, DEAF.residual_capacity(j, s))
         # Do nothing more if the entering edge is the same as the leaving edge.
         if i != j:
-            if parent[t] != s:
+            if DEAF.parent[t] != s:
                 # Ensure that s is the parent of t.
                 s, t = t, s
             if We.index(i) > We.index(j):
                 # Ensure that q is in the subtree rooted at t.
                 p, q = q, p
-            remove_edge(s, t)
-            make_root(q)
-            add_edge(i, p, q)
-            update_potentials(i, p, q)
+            DEAF.remove_edge(s, t)
+            DEAF.make_root(q)
+            DEAF.add_edge(i, p, q)
+            DEAF.update_potentials(i, p, q)
 
     ###########################################################################
     # Infeasibility and unboundedness detection
     ###########################################################################
 
-    if any(x[i] != 0 for i in range(-n, 0)):
+    if any(DEAF.edge_flow[i] != 0 for i in range(-n, 0)):
         raise nx.NetworkXUnfeasible("no flow satisfies all node demands")
 
-    if any(x[i] * 2 >= faux_inf for i in range(e)) or any(
+    if any(DEAF.edge_flow[i] * 2 >= faux_inf for i in range(DEAF.edge_count)) or any(
         e[-1].get(capacity, inf) == inf and e[-1].get(weight, 0) < 0
         for e in nx.selfloop_edges(G, data=True)
     ):
@@ -548,9 +615,9 @@ def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
     # Flow cost calculation and flow dict construction
     ###########################################################################
 
-    del x[e:]
-    flow_cost = sum(c * x for c, x in zip(C, x))
-    flow_dict = {n: {} for n in N}
+    del DEAF.edge_flow[DEAF.edge_count :]
+    flow_cost = sum(w * x for w, x in zip(DEAF.edge_weights, DEAF.edge_flow))
+    flow_dict = {n: {} for n in DEAF.node_list}
 
     def add_entry(e):
         """Add a flow dict entry."""
@@ -564,14 +631,27 @@ def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
                 d = t
         d[e[-2]] = e[-1]
 
-    S = (N[s] for s in S)  # Use original nodes.
-    T = (N[t] for t in T)  # Use original nodes.
+    DEAF.edge_sources = (
+        DEAF.node_list[s] for s in DEAF.edge_sources
+    )  # Use original nodes.
+    DEAF.edge_targets = (
+        DEAF.node_list[t] for t in DEAF.edge_targets
+    )  # Use original nodes.
     if not multigraph:
-        for e in zip(S, T, x):
+        for e in zip(
+            DEAF.edge_sources,
+            DEAF.edge_targets,
+            DEAF.edge_flow,
+        ):
             add_entry(e)
         edges = G.edges(data=True)
     else:
-        for e in zip(S, T, K, x):
+        for e in zip(
+            DEAF.edge_sources,
+            DEAF.edge_targets,
+            DEAF.edge_keys,
+            DEAF.edge_flow,
+        ):
             add_entry(e)
         edges = G.edges(data=True, keys=True)
     for e in edges:
@@ -579,12 +659,12 @@ def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
             if e[-1].get(capacity, inf) == 0:
                 add_entry(e[:-1] + (0,))
         else:
-            c = e[-1].get(weight, 0)
-            if c >= 0:
+            w = e[-1].get(weight, 0)
+            if w >= 0:
                 add_entry(e[:-1] + (0,))
             else:
-                u = e[-1][capacity]
-                flow_cost += c * u
-                add_entry(e[:-1] + (u,))
+                c = e[-1][capacity]
+                flow_cost += w * c
+                add_entry(e[:-1] + (c,))
 
     return flow_cost, flow_dict
diff --git a/networkx/algorithms/flow/tests/test_networksimplex.py b/networkx/algorithms/flow/tests/test_networksimplex.py
new file mode 100644
index 00000000..9a9287bb
--- /dev/null
+++ b/networkx/algorithms/flow/tests/test_networksimplex.py
@@ -0,0 +1,383 @@
+import pytest

+import networkx as nx

+import os

+

+

+@pytest.fixture

+def simple_flow_graph():

+    G = nx.DiGraph()

+    G.add_node("a", demand=0)

+    G.add_node("b", demand=-5)

+    G.add_node("c", demand=50000000)

+    G.add_node("d", demand=-49999995)

+    G.add_edge("a", "b", weight=3, capacity=4)

+    G.add_edge("a", "c", weight=6, capacity=10)

+    G.add_edge("b", "d", weight=1, capacity=9)

+    G.add_edge("c", "d", weight=2, capacity=5)

+    return G

+

+

+@pytest.fixture

+def simple_no_flow_graph():

+    G = nx.DiGraph()

+    G.add_node("s", demand=-5)

+    G.add_node("t", demand=5)

+    G.add_edge("s", "a", weight=1, capacity=3)

+    G.add_edge("a", "b", weight=3)

+    G.add_edge("a", "c", weight=-6)

+    G.add_edge("b", "d", weight=1)

+    G.add_edge("c", "d", weight=-2)

+    G.add_edge("d", "t", weight=1, capacity=3)

+    return G

+

+

+def get_flowcost_from_flowdict(G, flowDict):

+    """Returns flow cost calculated from flow dictionary"""

+    flowCost = 0

+    for u in flowDict.keys():

+        for v in flowDict[u].keys():

+            flowCost += flowDict[u][v] * G[u][v]["weight"]

+    return flowCost

+

+

+def test_infinite_demand_raise(simple_flow_graph):

+    G = simple_flow_graph

+    inf = float("inf")

+    nx.set_node_attributes(G, {"a": {"demand": inf}})

+    pytest.raises(nx.NetworkXError, nx.network_simplex, G)

+

+

+def test_neg_infinite_demand_raise(simple_flow_graph):

+    G = simple_flow_graph

+    inf = float("inf")

+    nx.set_node_attributes(G, {"a": {"demand": -inf}})

+    pytest.raises(nx.NetworkXError, nx.network_simplex, G)

+

+

+def test_infinite_weight_raise(simple_flow_graph):

+    G = simple_flow_graph

+    inf = float("inf")

+    nx.set_edge_attributes(

+        G, {("a", "b"): {"weight": inf}, ("b", "d"): {"weight": inf}}

+    )

+    pytest.raises(nx.NetworkXError, nx.network_simplex, G)

+

+

+def test_nonzero_net_demand_raise(simple_flow_graph):

+    G = simple_flow_graph

+    nx.set_node_attributes(G, {"b": {"demand": -4}})

+    pytest.raises(nx.NetworkXUnfeasible, nx.network_simplex, G)

+

+

+def test_negative_capacity_raise(simple_flow_graph):

+    G = simple_flow_graph

+    nx.set_edge_attributes(G, {("a", "b"): {"weight": 1}, ("b", "d"): {"capacity": -9}})

+    pytest.raises(nx.NetworkXUnfeasible, nx.network_simplex, G)

+

+

+def test_no_flow_satisfying_demands(simple_no_flow_graph):

+    G = simple_no_flow_graph

+    pytest.raises(nx.NetworkXUnfeasible, nx.network_simplex, G)

+

+

+def test_sum_demands_not_zero(simple_no_flow_graph):

+    G = simple_no_flow_graph

+    nx.set_node_attributes(G, {"t": {"demand": 4}})

+    pytest.raises(nx.NetworkXUnfeasible, nx.network_simplex, G)

+

+

+def test_google_or_tools_example():

+    """

+    https://developers.google.com/optimization/flow/mincostflow

+    """

+    G = nx.DiGraph()

+    start_nodes = [0, 0, 1, 1, 1, 2, 2, 3, 4]

+    end_nodes = [1, 2, 2, 3, 4, 3, 4, 4, 2]

+    capacities = [15, 8, 20, 4, 10, 15, 4, 20, 5]

+    unit_costs = [4, 4, 2, 2, 6, 1, 3, 2, 3]

+    supplies = [20, 0, 0, -5, -15]

+    answer = 150

+

+    for i in range(len(supplies)):

+        G.add_node(i, demand=(-1) * supplies[i])  # supplies are negative of demand

+

+    for i in range(len(start_nodes)):

+        G.add_edge(

+            start_nodes[i],

+            end_nodes[i],

+            weight=unit_costs[i],

+            capacity=capacities[i],

+        )

+

+    flowCost, flowDict = nx.network_simplex(G)

+    assert flowCost == answer

+    assert flowCost == get_flowcost_from_flowdict(G, flowDict)

+

+

+def test_google_or_tools_example2():

+    """

+    https://developers.google.com/optimization/flow/mincostflow

+    """

+    G = nx.DiGraph()

+    start_nodes = [0, 0, 1, 1, 1, 2, 2, 3, 4, 3]

+    end_nodes = [1, 2, 2, 3, 4, 3, 4, 4, 2, 5]

+    capacities = [15, 8, 20, 4, 10, 15, 4, 20, 5, 10]

+    unit_costs = [4, 4, 2, 2, 6, 1, 3, 2, 3, 4]

+    supplies = [23, 0, 0, -5, -15, -3]

+    answer = 183

+

+    for i in range(len(supplies)):

+        G.add_node(i, demand=(-1) * supplies[i])  # supplies are negative of demand

+

+    for i in range(len(start_nodes)):

+        G.add_edge(

+            start_nodes[i],

+            end_nodes[i],

+            weight=unit_costs[i],

+            capacity=capacities[i],

+        )

+

+    flowCost, flowDict = nx.network_simplex(G)

+    assert flowCost == answer

+    assert flowCost == get_flowcost_from_flowdict(G, flowDict)

+

+

+def test_large():

+    fname = os.path.join(os.path.dirname(__file__), "netgen-2.gpickle.bz2")

+    G = nx.read_gpickle(fname)

+    flowCost, flowDict = nx.network_simplex(G)

+    assert 6749969302 == flowCost

+    assert 6749969302 == nx.cost_of_flow(G, flowDict)

+

+

+def test_simple_digraph():

+    G = nx.DiGraph()

+    G.add_node("a", demand=-5)

+    G.add_node("d", demand=5)

+    G.add_edge("a", "b", weight=3, capacity=4)

+    G.add_edge("a", "c", weight=6, capacity=10)

+    G.add_edge("b", "d", weight=1, capacity=9)

+    G.add_edge("c", "d", weight=2, capacity=5)

+    flowCost, H = nx.network_simplex(G)

+    soln = {"a": {"b": 4, "c": 1}, "b": {"d": 4}, "c": {"d": 1}, "d": {}}

+    assert flowCost == 24

+    assert nx.min_cost_flow_cost(G) == 24

+    assert H == soln

+

+

+def test_negcycle_infcap():

+    G = nx.DiGraph()

+    G.add_node("s", demand=-5)

+    G.add_node("t", demand=5)

+    G.add_edge("s", "a", weight=1, capacity=3)

+    G.add_edge("a", "b", weight=3)

+    G.add_edge("c", "a", weight=-6)

+    G.add_edge("b", "d", weight=1)

+    G.add_edge("d", "c", weight=-2)

+    G.add_edge("d", "t", weight=1, capacity=3)

+    pytest.raises(nx.NetworkXUnfeasible, nx.network_simplex, G)

+

+

+def test_transshipment():

+    G = nx.DiGraph()

+    G.add_node("a", demand=1)

+    G.add_node("b", demand=-2)

+    G.add_node("c", demand=-2)

+    G.add_node("d", demand=3)

+    G.add_node("e", demand=-4)

+    G.add_node("f", demand=-4)

+    G.add_node("g", demand=3)

+    G.add_node("h", demand=2)

+    G.add_node("r", demand=3)

+    G.add_edge("a", "c", weight=3)

+    G.add_edge("r", "a", weight=2)

+    G.add_edge("b", "a", weight=9)

+    G.add_edge("r", "c", weight=0)

+    G.add_edge("b", "r", weight=-6)

+    G.add_edge("c", "d", weight=5)

+    G.add_edge("e", "r", weight=4)

+    G.add_edge("e", "f", weight=3)

+    G.add_edge("h", "b", weight=4)

+    G.add_edge("f", "d", weight=7)

+    G.add_edge("f", "h", weight=12)

+    G.add_edge("g", "d", weight=12)

+    G.add_edge("f", "g", weight=-1)

+    G.add_edge("h", "g", weight=-10)

+    flowCost, H = nx.network_simplex(G)

+    soln = {

+        "a": {"c": 0},

+        "b": {"a": 0, "r": 2},

+        "c": {"d": 3},

+        "d": {},

+        "e": {"r": 3, "f": 1},

+        "f": {"d": 0, "g": 3, "h": 2},

+        "g": {"d": 0},

+        "h": {"b": 0, "g": 0},

+        "r": {"a": 1, "c": 1},

+    }

+    assert flowCost == 41

+    assert H == soln

+

+

+def test_digraph1():

+    # From Bradley, S. P., Hax, A. C. and Magnanti, T. L. Applied

+    # Mathematical Programming. Addison-Wesley, 1977.

+    G = nx.DiGraph()

+    G.add_node(1, demand=-20)

+    G.add_node(4, demand=5)

+    G.add_node(5, demand=15)

+    G.add_edges_from(

+        [

+            (1, 2, {"capacity": 15, "weight": 4}),

+            (1, 3, {"capacity": 8, "weight": 4}),

+            (2, 3, {"weight": 2}),

+            (2, 4, {"capacity": 4, "weight": 2}),

+            (2, 5, {"capacity": 10, "weight": 6}),

+            (3, 4, {"capacity": 15, "weight": 1}),

+            (3, 5, {"capacity": 5, "weight": 3}),

+            (4, 5, {"weight": 2}),

+            (5, 3, {"capacity": 4, "weight": 1}),

+        ]

+    )

+    flowCost, H = nx.network_simplex(G)

+    soln = {

+        1: {2: 12, 3: 8},

+        2: {3: 8, 4: 4, 5: 0},

+        3: {4: 11, 5: 5},

+        4: {5: 10},

+        5: {3: 0},

+    }

+    assert flowCost == 150

+    assert nx.min_cost_flow_cost(G) == 150

+    assert H == soln

+

+

+def test_zero_capacity_edges():

+    """Address issue raised in ticket #617 by arv."""

+    G = nx.DiGraph()

+    G.add_edges_from(

+        [

+            (1, 2, {"capacity": 1, "weight": 1}),

+            (1, 5, {"capacity": 1, "weight": 1}),

+            (2, 3, {"capacity": 0, "weight": 1}),

+            (2, 5, {"capacity": 1, "weight": 1}),

+            (5, 3, {"capacity": 2, "weight": 1}),

+            (5, 4, {"capacity": 0, "weight": 1}),

+            (3, 4, {"capacity": 2, "weight": 1}),

+        ]

+    )

+    G.nodes[1]["demand"] = -1

+    G.nodes[2]["demand"] = -1

+    G.nodes[4]["demand"] = 2

+

+    flowCost, H = nx.network_simplex(G)

+    soln = {1: {2: 0, 5: 1}, 2: {3: 0, 5: 1}, 3: {4: 2}, 4: {}, 5: {3: 2, 4: 0}}

+    assert flowCost == 6

+    assert nx.min_cost_flow_cost(G) == 6

+    assert H == soln

+

+

+def test_digon():

+    """Check if digons are handled properly. Taken from ticket

+    #618 by arv."""

+    nodes = [(1, {}), (2, {"demand": -4}), (3, {"demand": 4})]

+    edges = [

+        (1, 2, {"capacity": 3, "weight": 600000}),

+        (2, 1, {"capacity": 2, "weight": 0}),

+        (2, 3, {"capacity": 5, "weight": 714285}),

+        (3, 2, {"capacity": 2, "weight": 0}),

+    ]

+    G = nx.DiGraph(edges)

+    G.add_nodes_from(nodes)

+    flowCost, H = nx.network_simplex(G)

+    soln = {1: {2: 0}, 2: {1: 0, 3: 4}, 3: {2: 0}}

+    assert flowCost == 2857140

+

+

+def test_deadend():

+    """Check if one-node cycles are handled properly. Taken from ticket

+    #2906 from @sshraven."""

+    G = nx.DiGraph()

+

+    G.add_nodes_from(range(5), demand=0)

+    G.nodes[4]["demand"] = -13

+    G.nodes[3]["demand"] = 13

+

+    G.add_edges_from([(0, 2), (0, 3), (2, 1)], capacity=20, weight=0.1)

+    pytest.raises(nx.NetworkXUnfeasible, nx.network_simplex, G)

+

+

+def test_infinite_capacity_neg_digon():

+    """An infinite capacity negative cost digon results in an unbounded

+    instance."""

+    nodes = [(1, {}), (2, {"demand": -4}), (3, {"demand": 4})]

+    edges = [

+        (1, 2, {"weight": -600}),

+        (2, 1, {"weight": 0}),

+        (2, 3, {"capacity": 5, "weight": 714285}),

+        (3, 2, {"capacity": 2, "weight": 0}),

+    ]

+    G = nx.DiGraph(edges)

+    G.add_nodes_from(nodes)

+    pytest.raises(nx.NetworkXUnbounded, nx.network_simplex, G)

+

+

+def test_multidigraph():

+    """Multidigraphs are acceptable."""

+    G = nx.MultiDiGraph()

+    G.add_weighted_edges_from([(1, 2, 1), (2, 3, 2)], weight="capacity")

+    flowCost, H = nx.network_simplex(G)

+    assert flowCost == 0

+    assert H == {1: {2: {0: 0}}, 2: {3: {0: 0}}, 3: {}}

+

+

+def test_negative_selfloops():

+    """Negative selfloops should cause an exception if uncapacitated and

+    always be saturated otherwise.

+    """

+    G = nx.DiGraph()

+    G.add_edge(1, 1, weight=-1)

+    pytest.raises(nx.NetworkXUnbounded, nx.network_simplex, G)

+

+    G[1][1]["capacity"] = 2

+    flowCost, H = nx.network_simplex(G)

+    assert flowCost == -2

+    assert H == {1: {1: 2}}

+

+    G = nx.MultiDiGraph()

+    G.add_edge(1, 1, "x", weight=-1)

+    G.add_edge(1, 1, "y", weight=1)

+    pytest.raises(nx.NetworkXUnbounded, nx.network_simplex, G)

+

+    G[1][1]["x"]["capacity"] = 2

+    flowCost, H = nx.network_simplex(G)

+    assert flowCost == -2

+    assert H == {1: {1: {"x": 2, "y": 0}}}

+

+

+def test_bone_shaped():

+    # From #1283

+    G = nx.DiGraph()

+    G.add_node(0, demand=-4)

+    G.add_node(1, demand=2)

+    G.add_node(2, demand=2)

+    G.add_node(3, demand=4)

+    G.add_node(4, demand=-2)

+    G.add_node(5, demand=-2)

+    G.add_edge(0, 1, capacity=4)

+    G.add_edge(0, 2, capacity=4)

+    G.add_edge(4, 3, capacity=4)

+    G.add_edge(5, 3, capacity=4)

+    G.add_edge(0, 3, capacity=0)

+    flowCost, H = nx.network_simplex(G)

+    assert flowCost == 0

+    assert H == {0: {1: 2, 2: 2, 3: 0}, 1: {}, 2: {}, 3: {}, 4: {3: 2}, 5: {3: 2}}

+

+

+def test_graphs_type_exceptions():

+    G = nx.Graph()

+    pytest.raises(nx.NetworkXNotImplemented, nx.network_simplex, G)

+    G = nx.MultiGraph()

+    pytest.raises(nx.NetworkXNotImplemented, nx.network_simplex, G)

+    G = nx.DiGraph()

+    pytest.raises(nx.NetworkXError, nx.network_simplex, G)