Skip to content

Commit 3e64667

Browse files
0/1 fixed variable linearization also for DAGs
Refactored and enhanced the logic for fixing edge variables to zero or one in AbstractPathModelDAG and related classes, improving reachability-based pruning and protection for safe paths. Updated product constraint encoding in kFlowDecomp, kLeastAbsErrors, and kMinPathError to respect fixed edge variables. Added new demo and test files for acyclic graphs and updated example usage to reflect solver and optimization changes.
1 parent 7ee8dd2 commit 3e64667

13 files changed

+1066
-138
lines changed

examples/cycles_demo.py

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -27,13 +27,13 @@ def test_min_flow_decomp(filename: str):
2727
weight_type=int,
2828
subset_constraints=graph.graph["constraints"], # try with and without
2929
optimization_options={
30-
"optimize_with_safe_sequences": False, # set to false to deactivate the safe sequences optimization
30+
"optimize_with_safe_sequences": True, # set to false to deactivate the safe sequences optimization
3131
"optimize_with_safe_sequences_allow_geq_constraints": True,
3232
# "optimize_with_safe_sequences_fix_via_bounds": True,
3333
"optimize_with_safe_sequences_fix_zero_edges": True,
3434
},
3535
solver_options={
36-
"external_solver": "gurobi", # we can try also "highs" at some point
36+
"external_solver": "highs", # we can try also "highs" at some point
3737
"time_limit": 300, # 300s = 5min, is it ok?
3838
"threads": 1
3939
},
@@ -187,12 +187,12 @@ def process_solution(model):
187187
print("avg_size_of_non_trivial_SCC:", solve_statistics['avg_size_of_non_trivial_SCC']) # size = number of edges
188188

189189
def main():
190-
test_min_flow_decomp(filename = "tests/cyclic_graphs/gt3.kmer15.(130000.132000).V23.E32.cyc100.graph")
190+
# test_min_flow_decomp(filename = "tests/cyclic_graphs/gt3.kmer15.(130000.132000).V23.E32.cyc100.graph")
191191
# test_min_flow_decomp(filename = "tests/cyclic_graphs/gt5.kmer27.(1300000.1400000).V809.E1091.mincyc1000.graph")
192192
test_min_flow_decomp(filename = "tests/cyclic_graphs/gt32.kmer63.(0.10000).V231.E336.mincyc1.e1.0.graph")
193193
# test_min_flow_decomp(filename = "tests/cyclic_graphs/gt4.kmer15.(0.10000).V1096.E1622.mincyc100.e1.0.graph")
194-
test_least_abs_errors(filename = "tests/cyclic_graphs/gt5.kmer27.(655000.660000).V18.E27.mincyc4.e0.75.graph")
195-
test_min_path_error(filename = "tests/cyclic_graphs/gt5.kmer27.(655000.660000).V18.E27.mincyc4.e0.75.graph")
194+
# test_least_abs_errors(filename = "tests/cyclic_graphs/gt5.kmer27.(655000.660000).V18.E27.mincyc4.e0.75.graph")
195+
# test_min_path_error(filename = "tests/cyclic_graphs/gt5.kmer27.(655000.660000).V18.E27.mincyc4.e0.75.graph")
196196

197197
if __name__ == "__main__":
198198
# Configure logging

examples/mfd_demo.py

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
import flowpaths as fp
2+
import networkx as nx
3+
4+
def test_min_flow_decomp(filename: str):
5+
graph = fp.graphutils.read_graphs(filename)[0]
6+
print("graph id", graph.graph["id"])
7+
print("subset_constraints", graph.graph["constraints"])
8+
# fp.utils.draw(
9+
# G=graph,
10+
# filename=filename + ".pdf",
11+
# flow_attr="flow",
12+
# subpath_constraints=graph.graph["constraints"],
13+
# draw_options={
14+
# "show_graph_edges": True,
15+
# "show_edge_weights": True,
16+
# "show_path_weights": False,
17+
# "show_path_weight_on_first_edge": True,
18+
# "pathwidth": 2,
19+
# "style": "points",
20+
# })
21+
22+
print(graph.graph["n"], graph.graph["m"], graph.graph["w"])
23+
24+
mfd_model = fp.kFlowDecomp(
25+
G=graph,
26+
flow_attr="flow",
27+
k = 15,
28+
weight_type=int,
29+
optimization_options={
30+
"optimize_with_safe_sequences": True, # set to false to deactivate the safe sequences optimization
31+
"optimize_with_safe_paths": False, # set to false to deactivate the safe paths optimization
32+
"optimize_with_flow_safe_paths": False,
33+
"optimize_with_safe_zero_edges": True,
34+
"optimize_with_greedy": False,
35+
},
36+
solver_options={
37+
"external_solver": "highs", # we can try also "highs" at some point
38+
"time_limit": 300, # 300s = 5min, is it ok?
39+
"threads": 1
40+
},
41+
)
42+
mfd_model.solve()
43+
process_solution(mfd_model)
44+
45+
def process_solution(model):
46+
if model.is_solved():
47+
solution = model.get_solution()
48+
print("solution paths:", solution['paths'])
49+
print("solution weights:", solution['weights'])
50+
print("model.is_valid_solution()", model.is_valid_solution()) # Keep this to verify the solution
51+
else:
52+
print("Model could not be solved.")
53+
54+
# fp.utils.draw(
55+
# G=model.G,
56+
# filename= "solution.pdf",
57+
# flow_attr="flow",
58+
# paths=model.get_solution().get('walks', None),
59+
# weights=model.get_solution().get('weights', None),
60+
# draw_options={
61+
# "show_graph_edges": False,
62+
# "show_edge_weights": False,
63+
# "show_path_weights": False,
64+
# "show_path_weight_on_first_edge": True,
65+
# "pathwidth": 2,
66+
# # "style": "points",
67+
# })
68+
69+
solve_statistics = model.solve_statistics
70+
print(solve_statistics)
71+
# print("node_number:", solve_statistics['node_number'])
72+
# print("edge_number:", solve_statistics['edge_number'])
73+
# print("safe_sequences_time:", solve_statistics.get('safe_sequences_time', 0)) # the time to compute safe sequences. use get(), as this is not set if not using safe sequences
74+
# print("edge_variables_total:", solve_statistics['edge_variables_total']) # number of edges * number of solution walks in the last iteration
75+
# print("edge_variables=1:", solve_statistics['edge_variables=1'])
76+
# print("edge_variables>=1:", solve_statistics['edge_variables>=1'])
77+
# print("edge_variables=0:", solve_statistics['edge_variables=0'])
78+
# print("graph_width:", solve_statistics['graph_width']) # the the minimum number of s-t walks needed to cover all edges
79+
# print("model_status:", solve_statistics['model_status'])
80+
# print("solve_time:", solve_statistics['solve_time']) # time taken by the ILP for a given k, or by MFD to iterate through k and do small internal things
81+
# print("solve_time_ilp:", solve_statistics['solve_time_ilp']) # time taken by the ILP for a given k, or by MFD to iterate through k and do small internal things
82+
83+
def main():
84+
# test_min_flow_decomp(filename = "tests/acyclic_graphs/gt15.kmer21.(288000.294000).V390.E560.acyc.graph")
85+
test_min_flow_decomp(filename = "tests/acyclic_graphs/gt15.kmer21.(612000.618000).V89.E128.acyc.graph")
86+
87+
if __name__ == "__main__":
88+
# Configure logging
89+
fp.utils.configure_logging(
90+
level=fp.utils.logging.INFO,
91+
log_to_console=True,
92+
)
93+
main()

flowpaths/abstractpathmodeldag.py

Lines changed: 104 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -195,6 +195,9 @@ def __init__(
195195
self.encode_path_length = encode_path_length
196196
self.edge_position_vars = {}
197197

198+
self.edges_set_to_zero = {}
199+
self.edges_set_to_one = {}
200+
198201
self.solver_options = solver_options
199202
if self.solver_options is None:
200203
self.solver_options = {}
@@ -288,6 +291,8 @@ def create_solver_and_paths(self):
288291

289292
self._encode_paths()
290293

294+
self._apply_safety_optimizations_fix_zero_edges()
295+
291296
def _encode_paths(self):
292297

293298
# Encodes the paths in the graph by creating variables for edges and subpaths.
@@ -447,49 +452,108 @@ def _encode_paths(self):
447452
name=f"path_length_constr_i={i}"
448453
)
449454

450-
########################################
451-
# #
452-
# Fixing variables based on safe lists #
453-
# #
454-
########################################
455+
def _apply_safety_optimizations(self):
455456

456457
if self.safe_lists is not None:
457-
paths_to_fix = self._get_paths_to_fix_from_safe_lists()
458-
459-
if not self.optimize_with_safety_as_subpath_constraints:
460-
# iterating over safe lists
461-
for i in range(min(len(paths_to_fix), self.k)):
462-
# print("Fixing variables for safe list #", i)
463-
# iterate over the edges in the safe list to fix variables to 1
464-
for u, v in paths_to_fix[i]:
465-
self.solver.add_constraint(
466-
self.edge_vars[(u, v, i)] == 1,
467-
name=f"safe_list_u={u}_v={v}_i={i}",
468-
)
458+
self.paths_to_fix = self._get_paths_to_fix_from_safe_lists()
459+
460+
if not self.optimize_with_safety_as_subpath_constraints:
461+
# iterating over safe lists
462+
for i in range(min(len(self.paths_to_fix), self.k)):
463+
# print("Fixing variables for safe list #", i)
464+
# iterate over the edges in the safe list to fix variables to 1
465+
for u, v in self.paths_to_fix[i]:
466+
self.solver.add_constraint(
467+
self.edge_vars[(u, v, i)] == 1,
468+
name=f"safe_list_u={u}_v={v}_i={i}",
469+
)
470+
self.edges_set_to_one[(u, v, i)] = True
471+
472+
self._apply_safety_optimizations_fix_zero_edges()
473+
474+
def _apply_safety_optimizations_fix_zero_edges(self):
475+
"""
476+
Prune layer-edge variables to zero using safe-walk reachability while
477+
preserving edges that can be part of the walk or its connectors.
478+
479+
For each walk i in `walks_to_fix` we build a protection set of edges that
480+
must not be fixed to 0 for layer i:
481+
1) Protect all edges that appear in the walk itself.
482+
2) Whole-walk reachability: let first_node be the first node of the walk
483+
and last_node the last node. Protect any edge (u,v) such that
484+
- u is reachable (forward) from last_node, OR
485+
- v can reach (backward) first_node.
486+
3) Gap-bridging between consecutive edges: for every pair of consecutive
487+
edges whose endpoints do not match (a gap), let
488+
- current_last = end node of the first edge, and
489+
- current_start = start node of the next edge.
490+
Protect any edge (u,v) such that
491+
- u is reachable (forward) from current_last, AND
492+
- v can reach (backward) current_start.
493+
494+
All remaining edges (u,v) not in the protection set are fixed to 0 in
495+
layer i.
496+
497+
Notes:
498+
- Requires `self.paths_to_fix` already computed and `self.edge_vars` created.
499+
"""
500+
if not hasattr(self, "paths_to_fix") or self.paths_to_fix is None:
501+
return
502+
503+
fixed_zero_count = 0
504+
# Ensure we don't go beyond k layers
505+
for i in range(min(len(self.paths_to_fix), self.k)):
506+
path = self.paths_to_fix[i]
507+
if not path or len(path) == 0:
508+
continue
509+
510+
# Build the set of edges that should NOT be fixed to 0 for this layer i
511+
# Start by protecting all edges in the path itself
512+
protected_edges = set((u, v) for (u, v) in path if self.G.has_edge(u, v))
513+
514+
# Also protect edges that are reachable from the last node of the path
515+
# or that can reach the first node of the path
516+
first_node = path[0][0]
517+
last_node = path[-1][1]
518+
for (u, v) in self.G.edges:
519+
if (u in self.G.reachable_nodes_from[last_node]) or (v in self.G.nodes_reaching(first_node)):
520+
protected_edges.add((u, v))
521+
522+
# Collect pairs of non-contiguous consecutive edges (gaps)
523+
gap_pairs = []
524+
for idx in range(len(path) - 1):
525+
end_prev = path[idx][1]
526+
start_next = path[idx + 1][0]
527+
# We consider all consecutive edges as gap pairs, because there could be a cycle
528+
# formed between them (this is not the case in DAGs)
529+
if end_prev != start_next:
530+
gap_pairs.append((end_prev, start_next))
531+
532+
# For each gap, add edges that can lie on some path bridging the gap
533+
for (current_last, current_start) in gap_pairs:
534+
for (u, v) in self.G.edges:
535+
if (u in self.G.nodes_reachable(current_last)) and (v in self.G.nodes_reaching(current_start)):
536+
# if (u in reachable_from_last) and (v in can_reach_start):
537+
protected_edges.add((u, v))
538+
539+
# Now fix every other edge to 0 for this layer i
540+
for (u, v) in self.G.edges:
541+
if (u, v) in protected_edges:
542+
continue
543+
# Queue zero-fix for batch bounds update
544+
# self.solver.queue_fix_variable(self.edge_vars[(u, v, i)], int(0))
545+
self.solver.add_constraint(
546+
self.edge_vars[(u, v, i)] == 0,
547+
name=f"i={i}_u={u}_v={v}_fix0",
548+
)
549+
self.edges_set_to_zero[(u, v, i)] = True
550+
fixed_zero_count += 1
551+
552+
if fixed_zero_count:
553+
# Accumulate into solve statistics
554+
self.solve_statistics["edge_variables=0"] = self.solve_statistics.get("edge_variables=0", 0) + fixed_zero_count
555+
utils.logger.debug(f"{__name__}: Fixed {fixed_zero_count} edge variables to 0 via reachability pruning.")
469556

470-
if self.optimize_with_safe_zero_edges:
471-
# get the endpoints of the longest safe path in the sequence
472-
first_node, last_node = (
473-
safetypathcovers.get_endpoints_of_longest_safe_path_in(paths_to_fix[i])
474-
)
475-
# get the reachable nodes from the last node
476-
reachable_nodes = self.G.reachable_nodes_from[last_node]
477-
# get the backwards reachable nodes from the first node
478-
reachable_nodes_reverse = self.G.reachable_nodes_rev_from[first_node]
479-
# get the edges in the path
480-
path_edges = set((u, v) for (u, v) in paths_to_fix[i])
481-
482-
for u, v in self.G.base_graph.edges():
483-
if (
484-
(u, v) not in path_edges
485-
and u not in reachable_nodes
486-
and v not in reachable_nodes_reverse
487-
):
488-
# print(f"Adding zero constraint for edge ({u}, {v}) in path {i}")
489-
self.solver.add_constraint(
490-
self.edge_vars[(u, v, i)] == 0,
491-
name=f"safe_list_zero_edge_u={u}_v={v}_i={i}",
492-
)
493557

494558

495559
def _get_paths_to_fix_from_safe_lists(self) -> list:

flowpaths/abstractwalkmodeldigraph.py

Lines changed: 1 addition & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -462,8 +462,7 @@ def _apply_safety_optimizations_fix_zero_edges(self):
462462
and last_node the last node. Protect any edge (u,v) such that
463463
- u is reachable (forward) from last_node, OR
464464
- v can reach (backward) first_node.
465-
3) Gap-bridging between consecutive edges: for every pair of consecutive
466-
edges whose endpoints do not match (a gap), let
465+
3) Gap-bridging between consecutive edges: for every pair of consecutive, let
467466
- current_last = end node of the first edge, and
468467
- current_start = start node of the next edge.
469468
Protect any edge (u,v) such that
@@ -475,11 +474,6 @@ def _apply_safety_optimizations_fix_zero_edges(self):
475474
476475
Notes:
477476
- Requires `self.walks_to_fix` already computed and `self.edge_vars` created.
478-
- Reachability is computed with networkx `descendants` (forward) and
479-
`ancestors` (backward), including the seed node itself in each set.
480-
- If a walk has no gaps, only the whole-walk reachability protection (2)
481-
applies. If the graph structure makes many edges reachable, little or no
482-
pruning may occur for that walk.
483477
"""
484478
if not hasattr(self, "walks_to_fix") or self.walks_to_fix is None:
485479
return
@@ -513,10 +507,6 @@ def _apply_safety_optimizations_fix_zero_edges(self):
513507
if True or end_prev != start_next:
514508
gap_pairs.append((end_prev, start_next))
515509

516-
# If there are no gaps, do not prune anything for this walk/layer
517-
if not gap_pairs:
518-
continue
519-
520510
# For each gap, add edges that can lie on some path bridging the gap
521511
for (current_last, current_start) in gap_pairs:
522512
for (u, v) in self.G.edges:

flowpaths/kflowdecomp.py

Lines changed: 22 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -117,6 +117,8 @@ def __init__(
117117
- ValueError: If `flow_attr_origin` is not "node" or "edge".
118118
"""
119119

120+
utils.logger.info(f"{__name__}: START initializing with graph id = {utils.fpid(G)}, k = {k}")
121+
120122
# Handling node-weighted graphs
121123
self.flow_attr_origin = flow_attr_origin
122124
if self.flow_attr_origin == "node":
@@ -253,7 +255,7 @@ def __init__(
253255
else:
254256
self._encode_flow_decomposition_with_given_weights()
255257

256-
utils.logger.info(f"{__name__}: initialized with graph id = {utils.fpid(G)}, k = {self.k}")
258+
utils.logger.info(f"{__name__}: END initialized with graph id = {utils.fpid(G)}, k = {self.k}")
257259

258260
def _encode_flow_decomposition(self):
259261

@@ -290,14 +292,25 @@ def _encode_flow_decomposition(self):
290292
# We encode that edge_vars[(u,v,i)] * self.path_weights_vars[(i)] = self.pi_vars[(u,v,i)],
291293
# assuming self.w_max is a bound for self.path_weights_vars[(i)]
292294
for i in range(self.k):
293-
self.solver.add_binary_continuous_product_constraint(
294-
binary_var=self.edge_vars[(u, v, i)],
295-
continuous_var=self.path_weights_vars[(i)],
296-
product_var=self.pi_vars[(u, v, i)],
297-
lb=0,
298-
ub=self.w_max,
299-
name=f"10_u={u}_v={v}_i={i}",
300-
)
295+
if (u, v, i) in self.edges_set_to_zero:
296+
self.solver.add_constraint(
297+
self.pi_vars[(u, v, i)] == 0,
298+
name=f"i={i}_u={u}_v={v}_10b",
299+
)
300+
elif (u, v, i) in self.edges_set_to_one:
301+
self.solver.add_constraint(
302+
self.pi_vars[(u, v, i)] == self.path_weights_vars[(i)],
303+
name=f"i={i}_u={u}_v={v}_10b",
304+
)
305+
else:
306+
self.solver.add_binary_continuous_product_constraint(
307+
binary_var=self.edge_vars[(u, v, i)],
308+
continuous_var=self.path_weights_vars[(i)],
309+
product_var=self.pi_vars[(u, v, i)],
310+
lb=0,
311+
ub=self.w_max,
312+
name=f"10_u={u}_v={v}_i={i}",
313+
)
301314

302315
self.solver.add_constraint(
303316
self.solver.quicksum(self.pi_vars[(u, v, i)] for i in range(self.k)) == f_u_v,

0 commit comments

Comments
 (0)