Enhance safe sequence handling and new drawing style

Added a 'points' style option to the graph drawing utility for improved visualization. Improved handling of edge multiplicities in safe sequence constraints, now supporting multiple traversals for SCC edges and validating multiplicity for non-SCC edges. Introduced a new test case in safe_seq_cycles.py to demonstrate safe sequence and incompatible sequence computation, and updated cycles_demo.py to enable additional output and drawing.

Author: alexandrutomescu

examples/cycles_demo.py

@@ -4,7 +4,7 @@
 def test_min_flow_decomp(filename: str):
     graph = fp.graphutils.read_graphs(filename)[0]
     print("graph id", graph.graph["id"])
-    # print("subset_constraints", graph.graph["constraints"])
+    print("subset_constraints", graph.graph["constraints"])
     # fp.utils.draw(
     #         G=graph,
     #         filename=filename + ".pdf",
@@ -16,6 +16,7 @@ def test_min_flow_decomp(filename: str):
     #         "show_path_weights": False,
     #         "show_path_weight_on_first_edge": True,
     #         "pathwidth": 2,
+    #         "style": "points",
     #     })
 
     print(graph.graph["n"], graph.graph["m"], graph.graph["w"])
@@ -154,6 +155,21 @@ 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",
+            })
+
     solve_statistics = model.solve_statistics
     print(solve_statistics)
     print("node_number:", solve_statistics['node_number'])

examples/safe_seq_cycles.py

@@ -114,6 +114,72 @@ def test2():
     lae_model.solve()
     process_solution(graph, lae_model)
 
+def test6():
+    graph = nx.DiGraph()
+    graph.add_edge("s", "a", flow=10)   
+    graph.add_edge("a", "a", flow=10)
+    graph.add_edge("a", "b", flow=10)
+    graph.add_edge("a", "f", flow=10)
+    graph.add_edge("a", "h", flow=10)
+    graph.add_edge("b", "c", flow=10)
+    graph.add_edge("c", "d", flow=10)
+    graph.add_edge("d", "b", flow=10)
+    # graph.add_edge("d", "h", flow=10)
+    graph.add_edge("c", "t", flow=10)
+    graph.add_edge("s", "e", flow=10)
+    graph.add_edge("e", "f", flow=10)
+    graph.add_edge("f", "g", flow=10)
+    graph.add_edge("g", "h", flow=10)
+    graph.add_edge("h", "f", flow=10)
+    graph.add_edge("g", "t", flow=10)
+
+    stDiGraph = fp.stDiGraph(graph)
+    X = set(stDiGraph.edges())
+    safe_seqs = safety.maximal_safe_sequences_via_dominators(stDiGraph, X)
+    for seq in safe_seqs:
+        print("Safe sequence:", seq)
+
+    # For every edge of stDiGraph, compute the length of the longest safe seq using it
+    longest_safe_seq_length = dict()
+    for edge in stDiGraph.edges():
+        longest_safe_seq_length[edge] = max(len(seq) for seq in safe_seqs if edge in seq)
+    print("Longest safe seq lengths:", longest_safe_seq_length)
+
+    incompatible_sequences = stDiGraph.get_longest_incompatible_sequences(safe_seqs)
+    for seq in incompatible_sequences:
+        print("Incompatible sequence:", seq)
+
+    fp.graphutils.draw(
+        G=stDiGraph,
+        filename="safe_sequences_example.pdf",
+        subpath_constraints=incompatible_sequences,
+        flow_attr="flow",
+        draw_options={
+            "show_graph_edges": True,
+            "show_edge_weights": False,
+            "show_path_weights": False,
+            "show_path_weight_on_first_edge": True,
+            "pathwidth": 2,
+            "style": "points",
+        }
+    )
+
+    kmpe_model = fp.kMinPathErrorCycles(
+        G=graph,
+        flow_attr="flow",
+        weight_type=float,
+        optimization_options={
+            "optimize_with_safe_sequences": True, # set to false to deactivate the safe sequences optimization
+        },
+        solver_options={
+            "external_solver": "gurobi", # we can try also "highs" at some point
+            "time_limit": 300, # 300s = 5min, is it ok?
+        },
+    )
+    kmpe_model.solve()
+    assert(kmpe_model.is_solved())
+    assert(kmpe_model.is_valid_solution())
+
 def process_solution(graph, model: fp.kLeastAbsErrors):
     if model.is_solved():
         print(model.get_solution())
@@ -135,11 +201,12 @@ def process_solution(graph, model: fp.kLeastAbsErrors):
         print("Model could not be solved.")
 
 def main():
-    test1()
-    test3()
-    test2()
-    test4()
-    test5()
+    # test1()
+    # test3()
+    # test2()
+    # test4()
+    # test5()
+    test6()
 
 if __name__ == "__main__":
     # Configure logging

flowpaths/abstractwalkmodeldigraph.py

@@ -6,6 +6,7 @@
 import time
 import copy
 from abc import ABC, abstractmethod
+from collections import Counter
 
 class AbstractWalkModelDiGraph(ABC):
     # storing some defaults
@@ -407,32 +408,42 @@ def _apply_safety_optimizations(self):
 
         # Otherwise, we fix variables using the walks to fix
         if self.optimize_with_safe_sequences:
-            # iterating over safe lists
+            # Iterate over walks to fix (up to k layers) and enforce per-edge multiplicity
             for i in range(min(len(self.walks_to_fix), self.k)):
-                # 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]:
+                walk = self.walks_to_fix[i]
+                if not walk:
+                    continue
+
+                # Count multiplicities of each edge in this safe walk
+                edge_multiplicities = Counter(walk)  # keys are (u,v), values are number of occurrences
+
+                for (u, v), m in edge_multiplicities.items():
                     if self.G.is_scc_edge(u, v):
+                        # Inside an SCC we allow multiple traversals; enforce lower bound = m
                         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:
-                                self.solver.queue_set_var_lower_bound(self.edge_vars[(u, v, i)], 1)
+                                # Tighten variable lower bound directly
+                                self.solver.queue_set_var_lower_bound(self.edge_vars[(u, v, i)], m)
                             else:
+                                # Add inequality constraint x >= m
                                 self.solver.add_constraint(
-                                    self.edge_vars[(u, v, i)] >= 1,
-                                    name=f"safe_list_u={u}_v={v}_i={i}",
+                                    self.edge_vars[(u, v, i)] >= m,
+                                    name=f"safe_list_u={u}_v={v}_i={i}_geq{m}",
                                 )
                             self.solve_statistics["edge_variables>=1"] += 1
+                        # If allow_geq_constraints is False, fall through without enforcing (original logic only acted when True)
                     else:
-                        # Instead of adding an equality constraint x==1 we tighten bounds directly.
+                        # Edge not in SCC: x == 1 either via direct fix (bounds) or equality constraint
+                        if m != 1:
+                            utils.logger.critical(f"{__name__}: Unexpected multiplicity {m} != 1 for non-SCC edge ({u},{v})")
+                            raise ValueError(f"Unexpected multiplicity {m} != 1 for non-SCC edge ({u},{v})")
                         if self.optimize_with_safe_sequences_fix_via_bounds:
-                            # Queue a fix to value 1
-                            self.solver.queue_fix_variable(self.edge_vars[(u, v, i)], int(1))
+                            # Since UB=1, fixing is safe and sets both LB & UB to m
+                            self.solver.queue_fix_variable(self.edge_vars[(u, v, i)], 1)
                         else:
-                            # Keep old behaviour via equality constraint
                             self.solver.add_constraint(
                                 self.edge_vars[(u, v, i)] == 1,
-                                name=f"safe_list_u={u}_v={v}_i={i}",
+                                name=f"safe_list_u={u}_v={v}_i={i}_eq{m}",
                             )
                         self.solve_statistics["edge_variables=1"] += 1
 

flowpaths/utils/graphutils.py

@@ -368,7 +368,8 @@ def draw(
             "show_node_weights": False,
             "show_path_weights": False,
             "show_path_weight_on_first_edge": True,
-            "pathwidth": 3.0
+            "pathwidth": 3.0,
+            "style": "default",
         },
         ):
         """
@@ -437,6 +438,10 @@ def draw(
             
                 The width of the path to be drawn. Default is `3.0`.
 
+            - `style`: str
+
+                The style of the drawing. Available options: `default`, `points`.
+
         """
 
         if len(paths) != len(weights) and len(weights) > 0:
@@ -447,8 +452,16 @@ def draw(
 
             dot = gv.Digraph(format="pdf")
             dot.graph_attr["rankdir"] = "LR"  # Display the graph in landscape mode
-            dot.node_attr["shape"] = "rectangle"  # Rectangle nodes
-            dot.node_attr["style"] = "rounded"  # Rounded rectangle nodes
+            
+            style = draw_options.get("style", "default")
+            if style == "default":
+                dot.node_attr["shape"] = "rectangle"  # Rectangle nodes
+                dot.node_attr["style"] = "rounded"  # Rounded rectangle nodes
+            elif style == "points":
+                dot.node_attr["shape"] = "point"  # Point nodes
+                dot.node_attr["style"] = "filled"  # Filled point nodes
+                # dot.node_attr['label'] = '' 
+                dot.node_attr['width'] = '0.1' 
 
             colors = [
                 "red",
@@ -483,17 +496,17 @@ def draw(
                     elif node in additional_ends:
                         color = "red"
                         penwidth = "2.0"
-                    
+
                     if draw_options.get("show_node_weights", False) and flow_attr is not None and flow_attr in G.nodes[node]:
+                        label = f"{G.nodes[node][flow_attr]}\\n{node}" if style != "points" else ""
                         dot.node(
                             name=str(node),
-                            #label=f"{{{node} | {G.nodes[node][flow_attr]}}}",
-                            label=f"{G.nodes[node][flow_attr]}\\n{node}",
-                            #xlabel=str(G.nodes[node][flow_attr]),
+                            label=label,
                             shape="record",
                             color=color, 
                             penwidth=penwidth)
                     else:
+                        label = str(node) if style != "points" else ""
                         dot.node(
                             name=str(node), 
                             label=str(node),