Improve SCC edge handling

Refactored SCC edge checks to use a unified is_scc_edge method, updated edge upper bounds for non-SCC edges, and improved logging for initialization and solving steps. Enhanced solve statistics reporting in MinFlowDecompCycles and clarified time tracking for ILP and model phases. Also adjusted demo script to print new statistics and commented out some test calls for clarity.

Author: alexandrutomescu

examples/cycles_demo.py

@@ -29,7 +29,7 @@ def test_min_flow_decomp(filename: str):
             "optimize_with_safe_sequences": True, # 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_zero_edges": True,
+            "optimize_with_safe_sequences_fix_zero_edges": False,
         },
         solver_options={
             "external_solver": "highs", # we can try also "highs" at some point
@@ -146,15 +146,16 @@ def process_solution(model):
     print("graph_width:", solve_statistics['graph_width']) # the the minimum number of s-t walks needed to cover all edges
     print("model_status:", solve_statistics['model_status'])
     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
+    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
     print("number_of_nontrivial_SCCs:", solve_statistics['number_of_nontrivial_SCCs']) # trivial = at least one edge
     print("size_of_largest_SCC:", solve_statistics['size_of_largest_SCC']) # size = number of edges
     print("avg_size_of_non_trivial_SCC:", solve_statistics['avg_size_of_non_trivial_SCC']) # size = number of edges
 
 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/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_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")
+    # 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")
 
 if __name__ == "__main__":
     # Configure logging

flowpaths/abstractwalkmodeldigraph.py

@@ -122,6 +122,11 @@ def __init__(
                 if edge not in self.edge_upper_bounds:
                     utils.logger.critical(f"{__name__}: Missing max_edge_repetition in max_edge_repetition_dict for edge {edge}")
                     raise ValueError(f"Missing max_edge_repetition for edge {edge}")
+        # We set to 1 in edge_upper_bounds if the edge is not inside an SCC of self.G,
+        # because these edges cannot be traversed more than 1 times by any walk
+        for edge in self.G.edges():
+            if not self.G.is_scc_edge(edge[0], edge[1]):
+                self.edge_upper_bounds[edge] = 1
 
         self.subset_constraints = copy.deepcopy(subset_constraints)
         if self.subset_constraints is not None:
@@ -404,7 +409,7 @@ def _apply_safety_optimizations(self):
                 # print("Fixing variables for safe list #", i)
                 # iterate over the edges in the safe list to fix variables to 1
                 for u, v in self.walks_to_fix[i]:
-                    if self.G._is_scc_edge(u, v):
+                    if self.G.is_scc_edge(u, v):
                         if self.optimize_with_safe_sequences_allow_geq_constraints:
                             # Raise LB via bounds only when enabled; else add constraint
                             if self.optimize_with_safe_sequences_fix_via_bounds:
@@ -605,7 +610,7 @@ def solve(self) -> bool:
         # self.write_model(f"model-{self.id}.lp")
         start_time = time.perf_counter()
         self.solver.optimize()
-        self.solve_statistics[f"solve_time"] = time.perf_counter() - start_time
+        self.solve_statistics[f"solve_time_ilp"] = time.perf_counter() - start_time
         self.solve_statistics[f"model_status"] = self.solver.get_model_status()
         self.solve_statistics[f"number_of_nontrivial_SCCs"] = self.G.get_number_of_nontrivial_SCCs()
         self.solve_statistics[f"avg_size_of_non_trivial_SCC"] = self.G.get_avg_size_of_non_trivial_SCC()

flowpaths/kflowdecompcycles.py

@@ -95,7 +95,8 @@ def __init__(
             - If the graph contains edges with negative flow values.
             - ValueError: If `flow_attr_origin` is not `node` or `edge`.
         """
-    
+        utils.logger.info(f"{__name__}: START initializing with graph id = {utils.fpid(G)}, k = {k}")
+
         # Handling node-weighted graphs
         self.flow_attr_origin = flow_attr_origin
         if self.flow_attr_origin == "node":
@@ -164,32 +165,46 @@ def __init__(
         self.optimization_options["trusted_edges_for_safety"] = self.G.get_non_zero_flow_edges(flow_attr=self.flow_attr, edges_to_ignore=self.edges_to_ignore)
 
         # Call the constructor of the parent class AbstractPathModelDAG
+        # Build per-edge repetition upper bounds: use the edge flow when available,
+        # otherwise fall back to self.w_max (e.g., for source/sink helper edges).
+        self.edge_upper_bounds_dict = {
+            (u, v): (data[self.flow_attr] if self.flow_attr in data else self.w_max)
+            for u, v, data in self.G.edges(data=True)
+        }
         super().__init__(
             G=self.G,
             k=self.k,
             # max_edge_repetition=self.w_max,
-            max_edge_repetition_dict=self.G.compute_edge_max_reachable_value(flow_attr=self.flow_attr),
+            max_edge_repetition_dict=self.edge_upper_bounds_dict,
             subset_constraints=self.subset_constraints,
             subset_constraints_coverage=self.subset_constraints_coverage,
             optimization_options=self.optimization_options,
             solver_options=solver_options,
             solve_statistics=self.solve_statistics
         )
 
+        utils.logger.debug(f"{__name__}: START create_solver_and_walks()")
         # This method is called from the super class AbstractWalkModelDiGraph
         self.create_solver_and_walks()
+        utils.logger.debug(f"{__name__}: END create_solver_and_walks()")
 
+        utils.logger.debug(f"{__name__}: START encoding flow decomposition")
         # This method is called from the current class 
         self._encode_flow_decomposition()
+        utils.logger.debug(f"{__name__}: END encoding flow decomposition")
 
+        utils.logger.debug(f"{__name__}: START encoding given weights")
         # This method is called from the current class
         self._encode_given_weights()
+        utils.logger.debug(f"{__name__}: END encoding given weights")
 
-        utils.logger.info(f"{__name__}: initialized with graph id = {utils.fpid(G)}, k = {self.k}")
+        utils.logger.info(f"{__name__}: END initialized with graph id = {utils.fpid(G)}, k = {self.k}")
 
     def _encode_flow_decomposition(self):
 
+        # print("edge_upper_bounds", self.edge_upper_bounds)
         # pi vars 
+        # edge_ubs = [self.edge_upper_bounds[(u, v)] for (u, v, i) in self.edge_indexes]
         self.pi_vars = self.solver.add_variables(
             self.edge_indexes,
             name_prefix="pi",

flowpaths/minflowdecompcycles.py

@@ -165,6 +165,7 @@ def __init__(
         self.solver_options = solver_options
         self.time_limit = self.solver_options.get("time_limit", sw.SolverWrapper.time_limit)
         self.solve_time_start = None
+        self.solve_time_ilp_total = 0
 
         self.solve_statistics = {}
         self._solution = None
@@ -200,6 +201,7 @@ def solve(self) -> bool:
             This overloads the `solve()` method from `AbstractWalkModelDiGraph` class.
         """
         self.solve_time_start = time.perf_counter()
+        utils.logger.info(f"{__name__}: starting to solve the MinFlowDecompCycles model for graph id = {utils.fpid(self.G)}")
 
         if self.optimization_options.get("optimize_with_given_weights", MinFlowDecompCycles.optimize_with_given_weights):            
             self._solve_with_given_weights()
@@ -213,10 +215,9 @@ def solve(self) -> bool:
                 if len(self._given_weights_model.get_solution(remove_empty_walks=True)["walks"]) == i:
                     fd_model = self._given_weights_model
 
-            fd_solver_options = copy.deepcopy(self.solver_options)
-            fd_solver_options["time_limit"] = self.time_limit - self.solve_time_elapsed
-
             if fd_model is None:
+                fd_solver_options = copy.deepcopy(self.solver_options)
+                fd_solver_options["time_limit"] = self.time_limit - self.solve_time_elapsed
                 fd_model = kflowdecompcycles.kFlowDecompCycles(
                     G=self.G,
                     flow_attr=self.flow_attr,
@@ -231,12 +232,17 @@ def solve(self) -> bool:
                     solver_options=fd_solver_options,
                 )
                 fd_model.solve()
-                self.solve_statistics = fd_model.solve_statistics
 
-                # If the previous run exceeded the time limit, 
-                # we still stop the search, even if we might have managed to solve it
-                if self.solve_time_elapsed > self.time_limit:
-                    return False
+            self.solve_statistics = fd_model.solve_statistics
+            self.solve_time_ilp_total += self.solve_statistics.get("solve_time_ilp", 0)
+            self.solve_statistics["solve_time"] = self.solve_time_elapsed
+            self.solve_statistics["solve_time_ilp"] = self.solve_time_ilp_total
+            self.solve_statistics["min_gen_set_solve_time"] = self._mingenset_model.solve_statistics.get("total_solve_time", 0) if self._mingenset_model is not None else 0
+
+            # If the previous run exceeded the time limit, 
+            # we still stop the search, even if we might have managed to solve it
+            if self.solve_time_elapsed > self.time_limit:
+                return False
 
             if fd_model.is_solved():
                 self._solution = fd_model.get_solution(remove_empty_walks=True)
@@ -245,17 +251,11 @@ def solve(self) -> bool:
                     self._solution["_walks_internal"] = self._solution["walks"]
                     self._solution["walks"] = self.G_internal.get_condensed_paths(self._solution["walks"])
                 self.set_solved()
-                self.solve_statistics = fd_model.solve_statistics
-                # we overwrite solve_time with the total time MFD has taken
-                self.solve_statistics["solve_time"] = time.perf_counter() - self.solve_time_start
-                self.solve_statistics["min_gen_set_solve_time"] = self._mingenset_model.solve_statistics.get("total_solve_time", 0) if self._mingenset_model is not None else 0
                 self.fd_model = fd_model
                 return True
             elif fd_model.solver.get_model_status() == sw.SolverWrapper.infeasible_status:
-                self.solve_statistics = fd_model.solve_statistics
                 utils.logger.info(f"{__name__}: model is infeasible for k = {i}")
             else:
-                self.solve_statistics = fd_model.solve_statistics
                 # If the model is not solved and the status is not infeasible,
                 # it means that the solver stopped because of an unexpected termination,
                 # thus we cannot conclude that the model is infeasible.

flowpaths/stdigraph.py

@@ -181,7 +181,7 @@ def _edge_to_condensation_node(self, u, v) -> str:
         - `ValueError` if the edge (u,v) is not inside an SCC.
         """
 
-        if not self._is_scc_edge(u, v):
+        if not self.is_scc_edge(u, v):
             utils.logger.error(f"{__name__}: Edge ({u},{v}) is not an edge inside an SCC.")
             raise ValueError(f"Edge ({u},{v}) is not an edge inside an SCC.")
 
@@ -237,7 +237,7 @@ def get_width(self, edges_to_ignore: list = None) -> int:
         for u, v in (edges_to_ignore or []):
             # If (u,v) is an edge between different SCCs
             # Then the corresponding edge to ignore is between the two SCCs
-            if not self._is_scc_edge(u, v):
+            if not self.is_scc_edge(u, v):
                 edge_multiplicity[self._edge_to_condensation_edge(u, v)] -= 1
             else:
                 # (u,v) is an edge within the same SCC
@@ -286,9 +286,9 @@ def get_width(self, edges_to_ignore: list = None) -> int:
 
         return width
 
-    def _is_scc_edge(self, u, v) -> bool:
+    def is_scc_edge(self, u, v) -> bool:
         """
-        Returns True if (u,v) is an edge inside an SCCs, False otherwise.
+        Returns True if (u,v) is an edge inside an SCC of self, False otherwise.
         """
 
         # Check if (u,v) is an edge of the graph