Own implementation of the Eades heuristic

setup.py

@@ -27,6 +27,7 @@
                 'colorlog',
                 'tqdm',
                 'dataclasses',
+                'cvxpy'
         ],
         entry_points={
             'console_scripts': ['macrogen=macrogen.main:main']

src/macrogen/etc/logging.yaml

@@ -18,8 +18,10 @@ root:
 loggers:
   macrogen.visualize:
     level: INFO
-#  macrogen.graph:
-#    level: DEBUG
+  macrogen.fes:
+    level: DEBUG
+  macrogen.graph:
+    level: DEBUG
 #  macrogen.uris:
 #    level: INFO
 #  macrogen.datings:

src/macrogen/fes.py

@@ -0,0 +1,92 @@
+import itertools
+from typing import Any, Tuple, List, Generator
+from .config import config
+import networkx as nx
+
+logger = config.getLogger(__name__)
+
+def _exhaust_sinks(g: nx.DiGraph, sink: bool = True):
+    """
+    Produces all sinks until there are no more.
+
+    Warning: This modifies the graph g
+    """
+    sink_method = g.out_degree if sink else g.in_degree
+    while True:
+        sinks = [u for (u, d) in sink_method() if d == 0]
+        if sinks:
+            yield from sinks
+            g.remove_nodes_from(sinks)
+        else:
+            return
+
+
+def _exhaust_sources(g: nx.DiGraph):
+    """
+    Produces all sources until there are no more.
+
+    Warning: This modifies the given graph
+    """
+    return _exhaust_sinks(g, False)
+
+
+def eades(graph: nx.DiGraph, double_check=True) -> List[Tuple[Any, Any]]:
+    """
+    Fast heuristic for the minimum feedback arc set.
+
+    Eades’ heuristic creates an ordering of all nodes of the given graph,
+    such that each edge can be classified into *forward* or *backward* edges.
+    The heuristic tries to minimize the sum of the weights (`weight` attribute)
+    of the backward edges. It always produces an acyclic graph, however it can
+    produce more conflicting edges than the minimal solution.
+
+    Args:
+        graph: a directed graph, may be a multigraph.
+        double_check: check whether we’ve _really_ produced an acyclic graph
+
+    Returns:
+        a list of edges, removal of which guarantees a
+
+    References:
+        **Eades, P., Lin, X. and Smyth, W. F.** (1993). A fast and effective
+        heuristic for the feedback arc set problem. *Information Processing
+        Letters*, **47**\ (6): 319–23
+        doi:\ `10.1016/0020-0190(93)90079-O. <https://doi.org/10.1016/0020-0190(93)90079-O.>`__
+        http://www.sciencedirect.com/science/article/pii/002001909390079O
+        (accessed 27 July 2018).
+    """
+    g = graph.copy()
+    logger.info('Internal eades calculation for a graph with %d nodes and %d edges', g.number_of_nodes(), g.number_of_edges())
+    g.remove_edges_from(list(g.selfloop_edges()))
+    start = []
+    end = []
+    while g:
+        for v in _exhaust_sinks(g):
+            end.insert(0, v)
+        for v in _exhaust_sources(g):
+            start.append(v)
+        if g:
+            u = max(g.nodes, key=lambda v: g.out_degree(v, weight='weight') - g.in_degree(v, weight='weight'))
+            start.append(u)
+            g.remove_node(u)
+    ordering = start + end
+    logger.debug('Internal ordering: %s', ordering)
+    pos = dict(zip(ordering, itertools.count()))
+    feedback_edges = list(graph.selfloop_edges())
+    for u, v in graph.edges():
+        if pos[u] > pos[v]:
+            feedback_edges.append((u, v))
+    logger.info('Found %d feedback edges', len(feedback_edges))
+
+    if double_check:
+        check = graph.copy()
+        check.remove_edges_from(feedback_edges)
+        if not nx.is_directed_acyclic_graph(check):
+            logger.error('double-check: graph is not a dag!')
+            cycles = nx.simple_cycles()
+            counter_example = next(cycles)
+            logger.error('Counterexample cycle: %s', counter_example)
+        else:
+            logger.info('double-check: Graph is acyclic.')
+
+    return feedback_edges

src/macrogen/graph.py

@@ -2,13 +2,12 @@
 Functions to build the graphs and perform their analyses.
 """
 
-
 import csv
 from collections import defaultdict, Counter
 from dataclasses import dataclass
 from datetime import date, timedelta
 from pathlib import Path
-from typing import List, Callable, Any, Dict, Tuple, Union
+from typing import List, Callable, Any, Dict, Tuple, Union, Iterable, Generator
 
 import networkx as nx
 
@@ -17,6 +16,7 @@
 from .igraph_wrapper import to_igraph, nx_edges
 from .uris import Reference, Inscription, Witness, AmbiguousRef
 from .config import config
+from .fes import eades
 
 logger = config.getLogger(__name__)
 
@@ -123,6 +123,24 @@ def remove_edges(source: nx.MultiDiGraph, predicate: Callable[[Any, Any, Dict[st
     # return nx.restricted_view(source, source.nodes, [(u,v,k) for u,v,k,attr in source.edges if predicate(u,v,attr)])
 
 
+def expand_edges(graph: nx.MultiDiGraph, edges: Iterable[Tuple[Any, Any]]) -> Generator[
+    Tuple[Any, Any, int, dict], None, None]:
+    """
+    Expands a 'simple' edge list (of node pairs) to the corresponding full edge list, including keys and data.
+    Args:
+        graph: the graph with the edges
+        edges: edge list, a list of (u, v) node tuples
+
+    Returns:
+        all edges from the multigraph that are between any node pair from edges as tuple (u, v, key, attrs)
+
+    """
+    for u, v in edges:
+        atlas = graph[u][v]
+        for key in atlas:
+            yield u, v, key, atlas[key]
+
+
 def feedback_arcs(graph: nx.MultiDiGraph, method='auto', auto_threshold=64):
     """
     Calculates the feedback arc set using the given method and returns a
@@ -134,11 +152,17 @@ def feedback_arcs(graph: nx.MultiDiGraph, method='auto', auto_threshold=64):
     """
     if method == 'auto':
         method = 'eades' if len(graph.edges) > auto_threshold else 'ip'
-    logger.debug('Calculating MFAS for a %d-node graph using %s, may take a while', graph.number_of_nodes(), method)
-    igraph = to_igraph(graph)
-    iedges = igraph.es[igraph.feedback_arc_set(method=method, weights='weight')]
-    logger.debug('%d edges to remove', len(iedges))
-    return list(nx_edges(iedges, keys=True, data=True))
+    if method == 'eades':
+        logger.debug('Calculating MFAS for a %d-node graph using internal Eades, may take a while',
+                     graph.number_of_nodes())
+        fes = eades(graph)
+        return list(expand_edges(graph, fes))
+    else:
+        logger.debug('Calculating MFAS for a %d-node graph using %s, may take a while', graph.number_of_nodes(), method)
+        igraph = to_igraph(graph)
+        iedges = igraph.es[igraph.feedback_arc_set(method=method, weights='weight')]
+        logger.debug('%d edges to remove', len(iedges))
+        return list(nx_edges(iedges, keys=True, data=True))
 
 
 def mark_edges_to_delete(graph: nx.MultiDiGraph, edges: List[Tuple[Any, Any, int, Any]]):
@@ -282,6 +306,7 @@ def year_stats(self):
         years = [node.avg_year for node in self.base.nodes if hasattr(node, 'avg_year') and node.avg_year is not None]
         return Counter(years)
 
+
 def resolve_ambiguities(graph: nx.MultiDiGraph):
     """
     Replaces ambiguous refs with the referenced nodes, possibly duplicating edges.
@@ -293,12 +318,12 @@ def resolve_ambiguities(graph: nx.MultiDiGraph):
     for ambiguity in ambiguities:
         for u, _, k, attr in list(graph.in_edges(ambiguity, keys=True, data=True)):
             for witness in ambiguity.witnesses:
-                attr['from_ambiguity']=ambiguity
+                attr['from_ambiguity'] = ambiguity
                 graph.add_edge(u, witness, k, **attr)
             graph.remove_edge(u, ambiguity, k)
         for _, v, k, attr in list(graph.out_edges(ambiguity, keys=True, data=True)):
             for witness in ambiguity.witnesses:
-                attr['from_ambiguity']=ambiguity
+                attr['from_ambiguity'] = ambiguity
                 graph.add_edge(witness, v, k, **attr)
             graph.remove_edge(ambiguity, v, k)
         graph.remove_node(ambiguity)
@@ -332,8 +357,6 @@ def datings_from_inscriptions(base: nx.MultiDiGraph):
                 base.add_edge(witness, d, copy=(d, i, k), **attr)
 
 
-
-
 def adopt_orphans(graph: nx.MultiDiGraph):
     """
     Introduces auxilliary edges to witnesses that are referenced by an inscription or ambiguous ref, but are not
@@ -362,6 +385,7 @@ def add_inscription_links(base: nx.MultiDiGraph):
         if isinstance(node, Inscription):
             base.add_edge(node, node.witness, kind='inscription', source=BiblSource('faust://model/inscription'))
 
+
 def add_missing_wits(working: nx.MultiDiGraph):
     """
     Add known witnesses that are not in the graph yet.
@@ -414,8 +438,8 @@ def macrogenesis_graphs() -> MacrogenesisInfo:
 
     if not nx.is_directed_acyclic_graph(result_graph):
         logger.error('After removing %d conflicting edges, the graph is still not a DAG!', len(all_conflicting_edges))
-        cycles = list(nx.simple_cycles(result_graph))
-        logger.error('It contains %d simple cycles', len(cycles))
+        cycles = nx.simple_cycles(result_graph)
+        logger.error('Counterexample cycle: %s.', next(cycles))
     else:
         logger.info('Double-checking removed edges ...')
         for u, v, k, attr in sorted(all_conflicting_edges, key=lambda edge: edge[3].get('weight', 1), reverse=True):

tests/test_fes.py

@@ -0,0 +1,28 @@
+import pytest
+import networkx as nx
+
+from macrogen.fes import _exhaust_sources, _exhaust_sinks, eades
+
+@pytest.fixture
+def graph1():
+    """
+           ←
+    1 → 2 → 3 → 4 → 5
+
+    """
+    G = nx.DiGraph()
+    G.add_path([1,2,3,4,5])
+    G.add_edge(3,2)
+    return G
+
+
+def test_all_sinks(graph1):
+    assert list(_exhaust_sinks(graph1)) == [5, 4]
+
+
+def test_all_sources(graph1):
+    assert list(_exhaust_sources(graph1)) == [1]
+
+
+def test_eades(graph1):
+    assert list(eades(graph1)) == [(3,2)]