diff --git a/src/macrogen/etc/default.yaml b/src/macrogen/etc/default.yaml index c6fb13d..8a3884d 100644 --- a/src/macrogen/etc/default.yaml +++ b/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 diff --git a/src/macrogen/fes.py b/src/macrogen/fes.py index d72bdc0..b5dd38b 100644 --- a/src/macrogen/fes.py +++ b/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): diff --git a/src/macrogen/graph.py b/src/macrogen/graph.py index 5f6ebba..805b849 100644 --- a/src/macrogen/graph.py +++ b/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,20 +141,26 @@ 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': @@ -162,10 +168,8 @@ def feedback_arcs(graph: nx.MultiDiGraph, method='baharev', auto_threshold=64): 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))