Made FES configurable

src/macrogen/etc/default.yaml

@@ -17,6 +17,17 @@ bibliography: https://raw.githubusercontent.com/faustedition/faust-gen-html/mast
 half_interval_correction: 182.5    # if we only have a start or end date, the other limit is max. this many days away
 render_node_limit: 1000
 
+## Options for solving the FES
+fes_method: baharev   # ip, baharev: exact methods; eades: inexact, list of two: select by fes_threshold
+fes_threshold: 64     # if two fes_methods, number of edges above which to select the second one
+solvers:              # the first installed solver is chosen. see log message for values.
+  - GUROBI
+  - auto              # auto = let CVXPY choose
+solver_options:
+  all:
+    max_iters: 10000
+
+
 ## Other data
 namespaces:
   f: http://www.faustedition.net/ns

src/macrogen/fes.py

@@ -157,11 +157,37 @@ def __init__(self, graph: nx.DiGraph):
         self.weights = np.array(weights)
         self.edges = edges
         self.m = len(self.edges)
-
+        self.solver_args = {}
         self.solution_vector = None
         self.solution = None
         self.objective = None
         self.iterations = None
+        self._load_solver_args()
+
+    def _load_solver_args(self):
+        # get solver settings from config
+        solvers: List[str] = config.solvers
+        installed = cp.installed_solvers()
+        index = 0
+        while index < len(solvers):
+            if solvers[index] not in installed:
+                del solvers[index]
+            else:
+                index += 1
+        if solvers and solvers[0]:
+            solver = solvers[0]
+        else:
+            solver = None
+        options = {}
+        if 'all' in config.solver_options:
+            options.update(config.solver_options['all'])
+        if solver and solver in config.solver_options:
+            options.update(config.solver_options[solver])
+        if solver:
+            options['solver'] = solver
+        self.solver_args = options
+        logger.info('configured solver: %s, options: %s (installed solvers: %s)',
+                    solver, options, ', '.join(installed))
 
     def edge_vector(self, edges: Iterable[Tuple[V, V]]) -> np.ndarray:
         """
@@ -211,21 +237,29 @@ def solve(self):
             cycle_vectors = [self.edge_vector(nx.utils.pairwise(cycle)) for cycle in simple_cycles]
             constraints = [cp.sum(a * y) >= 1 for a in cycle_vectors]
             problem = cp.Problem(objective, constraints)
-            resolution = problem.solve(solver=cp.GUROBI, verbose=False, max_iters=10000)
-            logger.debug("Solved optimization problem with %d constraints: %s -> %s (%s)", len(constraints), resolution,
-                         problem.solution.status, problem.solver_stats)
+            resolution = problem.solve(**self.solver_args)
+            if problem.status != 'optimal':
+                logger.warning('Optimization solution is %s. Try solver != %s?', problem.status,
+                               problem.solver_stats.solver_name)
+            logger.debug("Solved optimization problem with %d constraints: %s -> %s (%g + %g seconds, %d iterations, solver %s)",
+                         len(constraints), resolution, problem.solution.status,
+                         problem.solver_stats.solve_time, problem.solver_stats.setup_time or 0,
+                         problem.solver_stats.num_iters, problem.solver_stats.solver_name)
             current_solution = np.abs(y.value) >= 0.5
             S = self.edges_for_vector(current_solution)
             logger.debug('Iteration %d, resolution: %s, %d feedback edges', iteration, resolution, len(S))
             # S, the feedback edge set calculated using the constraint subset, can be an incomplete solution
             # (i.e. cycles remain after removing S from the graph). So lets compare this with the upper bound
             # from the heuristic
-
             lower_bound = max(lower_bound, objective.value)
             if lower_bound == upper_bound:
                 logger.info('upper == lower bound == %g, optimal solution found', lower_bound)
                 break  # y.value is the optimal solution
 
+            if resolution > upper_bound:
+                logger.error('Solution %g > upper bound %g!', resolution, upper_bound)
+                break
+
             Gi = self.graph.copy()
             Gi.remove_edges_from(S)
             if nx.is_directed_acyclic_graph(Gi):

src/macrogen/graph.py

@@ -3,7 +3,7 @@
 """
 
 import csv
-from collections import defaultdict, Counter
+from collections import defaultdict, Counter, Sequence
 from dataclasses import dataclass
 from datetime import date, timedelta
 from pathlib import Path
@@ -141,31 +141,35 @@ def expand_edges(graph: nx.MultiDiGraph, edges: Iterable[Tuple[Any, Any]]) -> Ge
             yield u, v, key, atlas[key]
 
 
-def feedback_arcs(graph: nx.MultiDiGraph, method='baharev', auto_threshold=64):
+def feedback_arcs(graph: nx.MultiDiGraph, method=None):
     """
     Calculates the feedback arc set using the given method and returns a
     list of edges in the form (u, v, key, data)
 
     Args:
         graph: NetworkX DiGraph
-        method: 'eades' (approximation, fast) or 'ip' (exact, exponential), or 'auto'
+        method: 'eades', 'baharev', or 'ip'; if None, look at config
     """
-    if method == 'auto':
-        method = 'eades' if len(graph.edges) > auto_threshold else 'ip'
+    if method is None:
+        method = config.fes_method
+    if isinstance(method, Sequence) and not isinstance(method, str):
+        try:
+            threshold = config.fes_threshold
+        except AttributeError:
+            threshold = 64
+        method = method[0] if len(graph.edges > threshold) else method[1]
+
+    logger.debug('Calculating MFAS for a %d-node graph using %s, may take a while', graph.number_of_nodes(), method)
     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))
     elif method == 'baharev':
         solver = FES_Baharev(graph)
         fes = solver.solve()
         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))