This repository has been archived by the owner on Sep 2, 2023. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 5
/
grb_lop.py
165 lines (149 loc) · 6.57 KB
/
grb_lop.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
# Copyright (C) 2014, 2015 University of Vienna
# All rights reserved.
# BSD license.
# Author: Ali Baharev <ali.baharev@gmail.com>
from __future__ import print_function
from itertools import combinations, product
from time import time
from benchmarks import gen_benchmark_digraphs
from grb_simplifier import iteratively_remove_runs_and_bypasses
from grb_pcm import get_orig_edges_map, feasible_solution
from mfes import noncopy_split_to_nontrivial_sccs
from py3compat import irange
from utils import double_check, solve_ilp, info_short
from networkx.algorithms.dag import topological_sort
################################################################################
#
# The code is admittedly ugly: This is a port of the CPLEX OPL script to Python.
# The linear ordering problem (LOP) formulation is ineffective for sparse
# graphs. The goal of this script is to demonstrate this.
#
################################################################################
def main():
for g_input in gen_benchmark_digraphs():
# Giving up on Problem 10 after 2 hours:
# 0 0 11.00000 0 402 12.00000 11.00000 8.33% - 7981s
if g_input.graph['name'] == 'Problem 10 (opt=12)':
continue
solve_problem(g_input)
# FIXME solve_triangle_inequalities assumes that we have a single SCC when it
# sets the feasible solution. Problems can arise if a feasible solution is
# passed in but the graph has more than one SCC.
def solve_problem(g_orig, stats=None, feasible_sol=None):
'Returns: [torn edges], cost.'
elims, cost = [ ], 0
for sc in noncopy_split_to_nontrivial_sccs(g_orig.copy()): # <- Copy passed!
partial_elims, partial_cost = \
solve_triangle_inequalities(sc, stats, feasible_sol)
elims.extend(partial_elims)
cost += partial_cost
double_check(g_orig, cost, elims, is_labeled=True)
print('Input graph')
info_short(g_orig)
return elims, cost
def solve_triangle_inequalities(g, stats, feasible_sol):
iteratively_remove_runs_and_bypasses(g)
# The above simplification removes edges, but the d['orig_edges'] still
# refers to these edges. As a result, the double_check calls in this module
# and in mfes will fail in the cost calculation, when trying to access
# those not present edges to figure out their edge weights. So we have to
# replace d['orig_edges'] appropriately and undo this at the end of this
# function. Additionally, as runs are removed, new edges are introduced
# that are not in the original graph. It is much simpler to just relabel
# the graph and then undo it on the elimination order.
origedges_map = get_orig_edges_map(g)
#
model, y = build_lp(g, feasible_sol)
start = time()
success = solve_ilp(model, stats)
end = time()
print('Overall solution time: {0:0.1f} s'.format(end-start))
assert success, 'Solver failures are not handled at the moment...'
cost = int(round(model.getObjective().getValue()))
elim_order = recover_order(model, y, len(g))
elims = get_torn_edges(g, elim_order)
# Done. Now undo the d['orig_edges'] mess.
final_elims = [ ]
for edge in elims:
final_elims.extend(origedges_map[edge])
double_check(g, cost, elims, is_labeled=True)
return final_elims, cost
def build_lp(g, feasible_sol):
from gurobipy import LinExpr, GRB, Model #, setParam
# We introduce an integer index per node ID. Sometimes we will use this
# index, and sometimes the node ID; this makes the code a bit messy.
i_nodeid = {i: n for i, n in enumerate(g)}
n_nodes = len(i_nodeid)
model = Model()
# The lower half of the n^2 binary variables, the rest: y_{j,i} = 1-y_{i,j}
y = { (i, j) : model.addVar(vtype=GRB.BINARY)
for i, j in combinations(irange(n_nodes), 2) }
model.update()
# The triangle inequalities
for i, j, k in combinations(irange(n_nodes), 3):
lhs = LinExpr([( 1,y[(i,j)]), ( 1,y[(j,k)]), (-1,y[(i,k)])])
model.addConstr(lhs, GRB.LESS_EQUAL, 1.0)
lhs = LinExpr([(-1,y[(i,j)]), (-1,y[(j,k)]), ( 1,y[(i,k)])])
model.addConstr(lhs, GRB.LESS_EQUAL, 0.0)
# The objective
shift = 0
c = [[0]*n_nodes for _ in irange(n_nodes)]
indices = ((i, j) for i, j in product(irange(n_nodes), repeat=2)
if g.has_edge(i_nodeid[j], i_nodeid[i]))
for j, k in indices:
w = g[i_nodeid[k]][i_nodeid[j]]['weight']
if k < j:
c[k][j] += w
elif k > j:
c[j][k] -= w
shift += w
obj = [(c[i][j], y[(i,j)]) for i, j in combinations(irange(n_nodes), 2)]
model.setObjective(LinExpr(obj), GRB.MINIMIZE)
model.setAttr('ObjCon', shift)
set_start_point(g, model, y, i_nodeid, feasible_sol)
model.update()
return model, y
def set_start_point(g_orig, model, y, i_nodeid, feasible_sol):
# FIXME It assumes that we only have a single SCC
elims, _cost = feasible_solution(g_orig) if feasible_sol is None else feasible_sol
g = g_orig.copy()
g.remove_edges_from(elims)
order = {n: i for i, n in enumerate(topological_sort(g))}
for i, j in combinations(irange(len(i_nodeid)), 2):
pos_a = order[i_nodeid[i]]
pos_b = order[i_nodeid[j]]
y[(i,j)].start = 0 if pos_a < pos_b else 1
def recover_order(model, y, n_nodes):
# Port of the OPL script: Admittedly ugly and inefficient
# order[j in 1..n] = sum(i in 1..j-1) x[i,j] + sum(k in j+1..n) (1-x[j,k])+1;
order = [-1]*n_nodes
for j in irange(n_nodes): # counting the number of nodes preceeding j
s = 0
for i in irange(0, j):
s += int(round(y[(i,j)].X))
for i in irange(j+1, n_nodes):
s += int(round((1-y[(j,i)].X)))
order[j] = s
# elim_order[k in 1..n] = sum(i in 1..n) i*(order[i]==k);
print(order)
elim_order = [-1]*n_nodes
for k in irange(n_nodes): # the permutation realizing `order`
elim_order[k] = sum(i if order[i]==k else 0 for i in irange(n_nodes))
print(elim_order)
return elim_order
def get_torn_edges(g, elim_order):
# Torn edges are the back edges (pointing backwards)
i_nodeid = {i: n for i, n in enumerate(g)}
keys = {i_nodeid[n]: i for i, n in enumerate(elim_order)}
n_edges = 0
edges = [ ]
for pos, i in enumerate(elim_order):
n = i_nodeid[i]
for nbr in g[n]:
if keys[nbr] > pos: # apparently the elim_order should be reversed?
n_edges +=1
edges.append((n, nbr))
print(n_edges)
return edges
if __name__ == '__main__':
main()