Improve safe sequence optimizations and reachability

Refactored safe sequence optimizations to fix edge variables directly via constraints and track zero/one assignments for use in product constraints. Added efficient node reachability queries to stDiGraph using SCC condensation, replacing networkx ancestor/descendant calls. Updated demo and test files to use new options and added a new cyclic graph test case.

Author: alexandrutomescu

examples/cycles_demo.py

@@ -27,14 +27,15 @@ def test_min_flow_decomp(filename: str):
         weight_type=int,
         subset_constraints=graph.graph["constraints"], # try with and without
         optimization_options={
-            "optimize_with_safe_sequences": True, # set to false to deactivate the safe sequences optimization
+            "optimize_with_safe_sequences": False, # set to false to deactivate the safe sequences optimization
             "optimize_with_safe_sequences_allow_geq_constraints": True,
-            "optimize_with_safe_sequences_fix_via_bounds": True,
+            # "optimize_with_safe_sequences_fix_via_bounds": True,
             "optimize_with_safe_sequences_fix_zero_edges": True,
         },
         solver_options={
             "external_solver": "gurobi", # we can try also "highs" at some point
             "time_limit": 300, # 300s = 5min, is it ok?
+            "threads": 1
         },
     )
     mfd_model.solve()
@@ -56,7 +57,7 @@ def test_least_abs_errors(filename):
             "optimize_with_safe_sequences_allow_geq_constraints": False,
         },
         solver_options={
-            "external_solver": "gurobi", # we can try also "highs" at some point
+            "external_solver": "highs", # we can try also "highs" at some point
             "time_limit": 300, # 300s = 5min, is it ok?
         },
         trusted_edges_for_safety_percentile=0, # we trust for safety edges whose weight in >= 0 percentile, that is, all edges
@@ -75,7 +76,7 @@ def test_least_abs_errors(filename):
             "optimize_with_safe_sequences_allow_geq_constraints": False,
         },
         solver_options={
-            "external_solver": "gurobi", # we can try also "highs" at some point
+            "external_solver": "highs", # we can try also "highs" at some point
             "time_limit": 300, # 300s = 5min, is it ok?
         },
         trusted_edges_for_safety_percentile=25, # we trust for safety edges whose weight in >= 25 percentile, remove this if not using the safety optimization
@@ -99,7 +100,7 @@ def test_min_path_error(filename):
             "optimize_with_safe_sequences_allow_geq_constraints": True,
         },
         solver_options={
-            "external_solver": "gurobi", # we can try also "highs" at some point
+            "external_solver": "highs", # we can try also "highs" at some point
             "time_limit": 300, # 300s = 5min, is it ok?
         },
     )
@@ -115,10 +116,9 @@ def test_min_path_error(filename):
         optimization_options={
             "optimize_with_safe_sequences": True, # set to false to deactivate the safe sequences optimization
             "optimize_with_safe_sequences_allow_geq_constraints": False,
-            "optimize_with_safe_sequences_fix_via_bounds": True,
         },
         solver_options={
-            "external_solver": "gurobi", # we can try also "highs" at some point
+            "external_solver": "highs", # we can try also "highs" at some point
             "time_limit": 300, # 300s = 5min, is it ok?
         },
         trusted_edges_for_safety_percentile=25, # we trust for safety edges whose weight in >= 25 percentile, remove this if not using the safety optimization
@@ -136,10 +136,9 @@ def test_min_path_error(filename):
         optimization_options={
             "optimize_with_safe_sequences": True, # set to false to deactivate the safe sequences optimization
             "optimize_with_safe_sequences_allow_geq_constraints": False,
-            "optimize_with_safe_sequences_fix_via_bounds": True,
         },
         solver_options={
-            "external_solver": "gurobi", # we can try also "highs" at some point
+            "external_solver": "highs", # we can try also "highs" at some point
             "time_limit": 300, # 300s = 5min, is it ok?
         },
     )
@@ -155,20 +154,20 @@ def process_solution(model):
     else:
         print("Model could not be solved.")
 
-    fp.utils.draw(
-            G=model.G,
-            filename= "solution.pdf",
-            flow_attr="flow",
-            paths=model.get_solution().get('walks', None),
-            weights=model.get_solution().get('weights', None),
-            draw_options={
-                "show_graph_edges": False,
-                "show_edge_weights": False,
-                "show_path_weights": False,
-                "show_path_weight_on_first_edge": True,
-                "pathwidth": 2,
-                # "style": "points",
-            })
+    # fp.utils.draw(
+    #         G=model.G,
+    #         filename= "solution.pdf",
+    #         flow_attr="flow",
+    #         paths=model.get_solution().get('walks', None),
+    #         weights=model.get_solution().get('weights', None),
+    #         draw_options={
+    #             "show_graph_edges": False,
+    #             "show_edge_weights": False,
+    #             "show_path_weights": False,
+    #             "show_path_weight_on_first_edge": True,
+    #             "pathwidth": 2,
+    #             # "style": "points",
+    #         })
 
     solve_statistics = model.solve_statistics
     print(solve_statistics)
@@ -189,7 +188,8 @@ def process_solution(model):
 
 def main():
     test_min_flow_decomp(filename = "tests/cyclic_graphs/gt3.kmer15.(130000.132000).V23.E32.cyc100.graph")
-    test_min_flow_decomp(filename = "tests/cyclic_graphs/gt5.kmer27.(1300000.1400000).V809.E1091.mincyc1000.graph")
+    # test_min_flow_decomp(filename = "tests/cyclic_graphs/gt5.kmer27.(1300000.1400000).V809.E1091.mincyc1000.graph")
+    test_min_flow_decomp(filename = "tests/cyclic_graphs/gt32.kmer63.(0.10000).V231.E336.mincyc1.e1.0.graph")
     # test_min_flow_decomp(filename = "tests/cyclic_graphs/gt4.kmer15.(0.10000).V1096.E1622.mincyc100.e1.0.graph")
     test_least_abs_errors(filename = "tests/cyclic_graphs/gt5.kmer27.(655000.660000).V18.E27.mincyc4.e0.75.graph")
     test_min_path_error(filename = "tests/cyclic_graphs/gt5.kmer27.(655000.660000).V18.E27.mincyc4.e0.75.graph")

flowpaths/abstractwalkmodeldigraph.py

@@ -13,12 +13,11 @@ class AbstractWalkModelDiGraph(ABC):
     optimize_with_safe_sequences = True
     optimize_with_safe_sequences_allow_geq_constraints = True
     optimize_with_safe_sequences_fix_via_bounds = False
-    optimize_with_safe_sequences_fix_zero_edges = False
+    optimize_with_safe_sequences_fix_zero_edges = True
     # TODO: optimize_with_subset_constraints_as_safe_sequences = True
     optimize_with_safety_as_subset_constraints = False
     optimize_with_max_safe_antichain_as_subset_constraints = False
     allow_empty_walks = False
-    allow_empty_walks = False
 
     def __init__(
         self,
@@ -149,6 +148,9 @@ def __init__(
             self.solver_options = {}
         self.threads = self.solver_options.get("threads", sw.SolverWrapper.threads)
 
+        self.edges_set_to_zero = {}
+        self.edges_set_to_one = {}
+
         # optimizations
         if optimization_options is None:
             optimization_options = {}
@@ -445,6 +447,7 @@ def _apply_safety_optimizations(self):
                                 self.edge_vars[(u, v, i)] == 1,
                                 name=f"safe_list_u={u}_v={v}_i={i}_eq{m}",
                             )
+                        self.edges_set_to_one[(u, v, i)] = True
                         self.solve_statistics["edge_variables=1"] += 1
 
     def _apply_safety_optimizations_fix_zero_edges(self):
@@ -496,18 +499,18 @@ def _apply_safety_optimizations_fix_zero_edges(self):
             # or that can reach the first node of the walk
             first_node = walk[0][0]
             last_node = walk[-1][1]
-            reachable_from_last_walk = nx.descendants(self.G, last_node) | {last_node}
-            can_reach_first_walk = nx.ancestors(self.G, first_node) | {first_node}
             for (u, v) in self.G.edges:
-                if (u in reachable_from_last_walk) or (v in can_reach_first_walk):
+                if (u in self.G.nodes_reachable(last_node)) or (v in self.G.nodes_reaching(first_node)):
                     protected_edges.add((u, v))
 
             # Collect pairs of non-contiguous consecutive edges (gaps)
             gap_pairs = []
             for idx in range(len(walk) - 1):
                 end_prev = walk[idx][1]
                 start_next = walk[idx + 1][0]
-                if end_prev != start_next:
+                # We consider all consecutive edges as gap pairs, because there could be a cycle
+                # formed between them (this is not the case in DAGs)
+                if True or end_prev != start_next:
                     gap_pairs.append((end_prev, start_next))
 
             # If there are no gaps, do not prune anything for this walk/layer
@@ -516,19 +519,22 @@ def _apply_safety_optimizations_fix_zero_edges(self):
 
             # For each gap, add edges that can lie on some path bridging the gap
             for (current_last, current_start) in gap_pairs:
-                reachable_from_last = nx.descendants(self.G, current_last) | {current_last}
-                can_reach_start = nx.ancestors(self.G, current_start) | {current_start}
-
                 for (u, v) in self.G.edges:
-                    if (u in reachable_from_last) and (v in can_reach_start):
+                    if (u in self.G.nodes_reachable(current_last)) and (v in self.G.nodes_reaching(current_start)):
+                    # if (u in reachable_from_last) and (v in can_reach_start):
                         protected_edges.add((u, v))
 
             # Now fix every other edge to 0 for this layer i
             for (u, v) in self.G.edges:
                 if (u, v) in protected_edges:
                     continue
                 # Queue zero-fix for batch bounds update
-                self.solver.queue_fix_variable(self.edge_vars[(u, v, i)], int(0))
+                # self.solver.queue_fix_variable(self.edge_vars[(u, v, i)], int(0))
+                self.solver.add_constraint(
+                    self.edge_vars[(u, v, i)] == 0,
+                    name=f"i={i}_u={u}_v={v}_fix0",
+                )
+                self.edges_set_to_zero[(u, v, i)] = True
                 fixed_zero_count += 1
 
         if fixed_zero_count:

flowpaths/kflowdecompcycles.py

@@ -228,6 +228,20 @@ def _encode_flow_decomposition(self):
             # We encode that edge_vars[(u,v,i)] * self.path_weights_vars[(i)] = self.pi_vars[(u,v,i)],
             # assuming self.w_max is a bound for self.path_weights_vars[(i)]
             for i in range(self.k):
+                if (u, v, i) in self.edges_set_to_zero:
+                    self.solver.add_constraint(
+                            self.pi_vars[(u, v, i)] == 0,
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
+                if (u, v, i) in self.edges_set_to_one:
+                    self.solver.add_constraint(
+                            self.pi_vars[(u, v, i)] == self.path_weights_vars[(i)],
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
                 self.solver.add_integer_continuous_product_constraint(
                     integer_var=self.edge_vars[(u, v, i)],
                     continuous_var=self.path_weights_vars[(i)],

flowpaths/kleastabserrorscycles.py

@@ -288,6 +288,20 @@ def _encode_leastabserrors_decomposition(self):
             # We encode that edge_vars[(u,v,i)] * self.path_weights_vars[(i)] = self.pi_vars[(u,v,i)],
             # assuming self.w_max is a bound for self.path_weights_vars[(i)]
             for i in range(self.k):
+                if (u, v, i) in self.edges_set_to_zero:
+                    self.solver.add_constraint(
+                            self.pi_vars[(u, v, i)] == 0,
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
+                if (u, v, i) in self.edges_set_to_one:
+                    self.solver.add_constraint(
+                            self.pi_vars[(u, v, i)] == self.path_weights_vars[(i)],
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
                 self.solver.add_integer_continuous_product_constraint(
                     integer_var=self.edge_vars[(u, v, i)],
                     continuous_var=self.path_weights_vars[(i)],

flowpaths/kminpatherrorcycles.py

@@ -326,6 +326,20 @@ def _encode_minpatherror_decomposition(self):
             # We encode that edge_vars[(u,v,i)] * self.path_weights_vars[(i)] = self.pi_vars[(u,v,i)],
             # assuming self.w_max is a bound for self.path_weights_vars[(i)]
             for i in range(self.k):
+                if (u, v, i) in self.edges_set_to_zero:
+                    self.solver.add_constraint(
+                            self.pi_vars[(u, v, i)] == 0,
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
+                if (u, v, i) in self.edges_set_to_one:
+                    self.solver.add_constraint(
+                            self.pi_vars[(u, v, i)] == self.path_weights_vars[(i)],
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
                 self.solver.add_integer_continuous_product_constraint(
                     integer_var=self.edge_vars[(u, v, i)],
                     continuous_var=self.path_weights_vars[(i)],
@@ -338,6 +352,20 @@ def _encode_minpatherror_decomposition(self):
             # We encode that edge_vars[(u,v,i)] * self.path_slacks_vars[(i)] = self.gamma_vars[(u,v,i)],
             # assuming self.w_max is a bound for self.path_slacks_vars[(i)]
             for i in range(self.k):
+                if (u, v, i) in self.edges_set_to_zero:
+                    self.solver.add_constraint(
+                            self.gamma_vars[(u, v, i)] == 0,
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
+                if (u, v, i) in self.edges_set_to_one:
+                    self.solver.add_constraint(
+                            self.gamma_vars[(u, v, i)] == self.path_slacks_vars[(i)],
+                            name=f"i={i}_u={u}_v={v}_10b",
+                        )
+                    continue
+
                 self.solver.add_integer_continuous_product_constraint(
                     integer_var=self.edge_vars[(u, v, i)],
                     continuous_var=self.path_slacks_vars[i],

flowpaths/stdigraph.py

@@ -4,7 +4,7 @@
 from flowpaths.stdag import stDAG
 import flowpaths.utils as utils
 from flowpaths.abstractsourcesinkgraph import AbstractSourceSinkGraph
-from typing import Dict, Tuple, Optional
+from typing import Dict, Tuple, Optional, Set
 
 
 class stDiGraph(AbstractSourceSinkGraph):
@@ -59,6 +59,27 @@ def _post_build(self):
             raise ValueError("The graph passed to stDiGraph must have at least one sink, or at least one node in `additional_ends`.")
         self.condensation_width = None
         self._build_condensation_expanded()
+        # Build indices and caches used by reachability queries
+        C: nx.DiGraph = self._condensation
+        mapping = C.graph["mapping"]  # original node -> condensation node (int)
+        # Per-SCC indices for fast unions
+        # - edges by SCC (kept for other utilities in this class)
+        self._edges_by_tail_scc: Dict[int, Set[Tuple[str, str]]] = {c: set() for c in C.nodes()}
+        self._edges_by_head_scc: Dict[int, Set[Tuple[str, str]]] = {c: set() for c in C.nodes()}
+        for a, b in self.edges():
+            ca = mapping[a]
+            cb = mapping[b]
+            self._edges_by_tail_scc[ca].add((a, b))
+            self._edges_by_head_scc[cb].add((a, b))
+
+        # - nodes by SCC (for fast node reachability queries)
+        self._nodes_by_scc: Dict[int, Set[str]] = {c: set() for c in C.nodes()}
+        for n in self.nodes():
+            self._nodes_by_scc[mapping[n]].add(n)
+
+        # Per-node caches for reachability queries (now returning nodes)
+        self._nodes_reachable_from_node_cache: Dict[str, Set[str]] = {}
+        self._nodes_reaching_node_cache: Dict[str, Set[str]] = {}
 
     def _expanded(self, v: int) -> str:
 
@@ -466,4 +487,80 @@ def compute_edge_max_reachable_value(self, flow_attr: str) -> Dict[Tuple[str, st
             result[(u, v)] = max(edge_weight[(u, v)], max_desc[cv], max_anc[cu])
 
         return result
+
+    def nodes_reachable(self, node: str) -> Set[str]:
+        """Return the set of nodes reachable from ``node`` (including itself).
+
+        The result is cached per query node. Reachability is computed on the SCC
+        condensation DAG: for the SCC containing ``node``, take all SCCs reachable
+        in the condensation (including itself) and return the union of original
+        nodes lying in any of those SCCs.
+
+        Parameters
+        ----------
+        node: str
+            The node ``v`` in this graph from which to evaluate forward reachability.
+
+        Returns
+        -------
+        Set[str]
+            All nodes ``a`` such that there exists a path from ``node`` to ``a``.
+        """
+        if node not in self.nodes():
+            utils.logger.error(f"{__name__}: Node {node} is not in the graph.")
+            raise ValueError(f"Node {node} is not in the graph.")
+        if node in self._nodes_reachable_from_node_cache:
+            return self._nodes_reachable_from_node_cache[node]
+
+        C: nx.DiGraph = self._condensation
+        mapping = C.graph["mapping"]
+        cv = mapping[node]
+
+        # All SCCs reachable from cv (descendants) plus itself
+        reachable_sccs = set(nx.descendants(C, cv)) | {cv}
+
+        result: Set[str] = set()
+        for c in reachable_sccs:
+            result |= self._nodes_by_scc.get(c, set())
+
+        self._nodes_reachable_from_node_cache[node] = result
+        return result
+
+    def nodes_reaching(self, node: str) -> Set[str]:
+        """Return the set of nodes that can reach ``node`` (including itself).
+
+        The result is cached per query node. Reachability is computed on the SCC
+        condensation DAG: for the SCC containing ``node``, take all SCCs that can
+        reach it (ancestors, including itself) and return the union of original
+        nodes lying in any of those SCCs.
+
+        Parameters
+        ----------
+        node: str
+            The node ``u`` in this graph to evaluate backward reachability to ``u``.
+
+        Returns
+        -------
+        Set[str]
+            All nodes ``a`` such that there exists a path from ``a`` to ``node``.
+        """
+        if node not in self.nodes():
+            utils.logger.error(f"{__name__}: Node {node} is not in the graph.")
+            raise ValueError(f"Node {node} is not in the graph.")
+        if node in self._nodes_reaching_node_cache:
+            return self._nodes_reaching_node_cache[node]
+
+        C: nx.DiGraph = self._condensation
+        mapping = C.graph["mapping"]
+        cu = mapping[node]
+
+        # All SCCs that can reach cu (ancestors) plus itself
+        ancestor_sccs = set(nx.ancestors(C, cu)) | {cu}
+
+        result: Set[str] = set()
+        for c in ancestor_sccs:
+            result |= self._nodes_by_scc.get(c, set())
+
+        self._nodes_reaching_node_cache[node] = result
+        return result