Merge pull request #123 from jcrozum/symbolic-fallback

Add a fallback to the fully symbolic attractor detection if the NFVS search fails

biobalm/_sd_algorithms/expand_source_blocks.py

@@ -201,10 +201,16 @@ def expand_source_blocks(
                     # MAAs in the problematic nodes while using the nice properties of the expansion to
                     # still disprove MAAs in the remaining nodes. If we used `seeds`, the expansion could
                     # just get stuck on this node and the "partial" results wouldn't be usable.
-                    block_sd_candidates = block_sd.node_attractor_candidates(
-                        block_sd.root(), compute=True
-                    )
-                    if len(block_sd_candidates) == 0:
+                    is_clean = False
+                    try:
+                        block_sd_candidates = block_sd.node_attractor_candidates(
+                            block_sd.root(), compute=True
+                        )
+                        is_clean = len(block_sd_candidates) == 0
+                    except RuntimeError:
+                        is_clean = False
+
+                    if is_clean:
                         if sd.config["debug"]:
                             print(
                                 f" > [{node}] Found clean block with no MAAs ({len(block_nodes)}): {block_nodes}"
@@ -217,7 +223,7 @@ def expand_source_blocks(
                     else:
                         if sd.config["debug"]:
                             print(
-                                f"[{node}] > Found {len(block_sd_candidates)} MAA cnadidates in a block. Delaying expansion."
+                                f"[{node}] > Found MAA candidates in a block (or failed candidate search). Delaying expansion."
                             )
                 if not clean_block_found:
                     # If all blocks have MAAs, we expand all successors.

biobalm/_sd_attractors/attractor_symbolic.py

@@ -6,9 +6,104 @@
     from biobalm.succession_diagram import SuccessionDiagram
     from biodivine_aeon import VariableId
 
-from biodivine_aeon import AsynchronousGraph, ColoredVertexSet, VertexSet
+from biodivine_aeon import (
+    AsynchronousGraph,
+    ColoredVertexSet,
+    VertexSet,
+    Attractors,
+    Reachability,
+)
 from biobalm.symbolic_utils import state_list_to_bdd
 from biobalm.types import BooleanSpace
+import biodivine_aeon
+
+
+def symbolic_attractor_fallback(
+    sd: SuccessionDiagram,
+    node_id: int,
+) -> tuple[list[BooleanSpace], list[VertexSet]]:
+    """
+    In case the attractor candidates cannot be computed, this fallback method
+    can be used to attempt fully symbolic attractor computation for the given node.
+
+    If the normal method works, then it is usually faster, but in the rare cases
+    where it fails, this can sometimes help to solve the few outlier nodes that
+    would otherwise block the computation.
+    """
+
+    old_log_level = biodivine_aeon.LOG_LEVEL
+
+    if sd.config["debug"]:
+        print(f"[{node_id}] > Start symbolic fallback.")
+        biodivine_aeon.LOG_LEVEL = biodivine_aeon.LOG_ESSENTIAL
+
+    node_data = sd.node_data(node_id)
+    node_space = node_data["space"]
+
+    candidates = sd.symbolic.mk_subspace(node_space)
+    if node_data["expanded"]:
+        # This node has successors that we should exclude from the attractor search.
+        for s in sd.node_successors(node_id, compute=False):
+            s_space = sd.node_data(s)["space"]
+            candidates = candidates.minus(sd.symbolic.mk_subspace(s_space))
+
+    fixed_variables = list(node_space.keys())
+    free_variables = [
+        v for v in sd.network.variable_names() if v not in fixed_variables
+    ]
+
+    if sd.config["debug"]:
+        print(f"[{node_id}] > Initial attractor candidates: {candidates}")
+
+    candidates = Attractors.transition_guided_reduction(
+        sd.symbolic, candidates, free_variables
+    )
+
+    if not sd.node_is_minimal(node_id) and not candidates.is_empty():
+        # This like it could be a motif-avoidant attractor. We better investigate this further,
+        # because this could get complicated...
+        avoid = sd.symbolic.mk_empty_colored_vertices()
+        if node_data["expanded"]:
+            for s in sd.node_successors(node_id):
+                s_space = sd.node_data(s)["space"]
+                avoid = avoid.union(sd.symbolic.mk_subspace(s_space))
+
+        if sd.config["debug"]:
+            print(
+                f"[{node_id}] > Reduction was ineffective. Start computing avoid set: {avoid}"
+            )
+
+        avoid = Reachability.reach_bwd(sd.symbolic, avoid)
+
+        if sd.config["debug"]:
+            print(f"[{node_id}] > Avoid set: {avoid}")
+
+        candidates = candidates.minus(avoid)
+    else:
+        if sd.config["debug"]:
+            print(
+                f"[{node_id}] > Node is minimal and the reduction did not finish with an empty set. This is fine."
+            )
+
+    if sd.config["debug"]:
+        print(f"[{node_id}] > Attractor candidates after reduction: {candidates}")
+
+    attractors = Attractors.xie_beerel(sd.symbolic, candidates)
+
+    biodivine_aeon.LOG_LEVEL = old_log_level
+
+    result_seeds: list[BooleanSpace] = []
+    result_sets: list[VertexSet] = []
+    for attr in attractors:
+        attr_vertices = attr.vertices()
+        attr_seed = next(attr_vertices.items()).to_dict()
+        attr_seed_named = {
+            sd.network.get_variable_name(k): v for (k, v) in attr_seed.items()
+        }
+        result_seeds.append(cast(BooleanSpace, attr_seed_named))
+        result_sets.append(attr_vertices)
+
+    return (result_seeds, result_sets)
 
 
 def compute_attractors_symbolic(

biobalm/succession_diagram.py

@@ -11,7 +11,10 @@
 
 # Attractor detection algorithms.
 from biobalm._sd_attractors.attractor_candidates import compute_attractor_candidates
-from biobalm._sd_attractors.attractor_symbolic import compute_attractors_symbolic
+from biobalm._sd_attractors.attractor_symbolic import (
+    compute_attractors_symbolic,
+    symbolic_attractor_fallback,
+)
 
 # SD expansion algorithms/heuristics.
 from biobalm._sd_algorithms.expand_attractor_seeds import expand_attractor_seeds
@@ -780,6 +783,7 @@ def node_attractor_seeds(
         self,
         node_id: int,
         compute: bool = False,
+        symbolic_fallback: bool = False,
     ) -> list[BooleanSpace]:
         """
         Return the list of attractor seed states for the given `node_id`.
@@ -797,6 +801,11 @@ def node_attractor_seeds(
             The ID of the node.
         compute: bool
             Whether to compute the attractor seeds if they are not already known.
+        symbolic_fallback: bool
+            If active, the method will attempt to compute the attractor seeds fully
+            symbolically if the default NFVS-based method fails. However, note that
+            the program can become unresponsive if the symbolic encoding of the
+            result grows to be too large. [Default: False]
 
         Returns
         -------
@@ -810,32 +819,48 @@ def node_attractor_seeds(
         if seeds is None and not compute:
             raise KeyError(f"Attractor seeds not computed for node {node_id}.")
 
-        if seeds is None:
-            candidates = self.node_attractor_candidates(node_id, compute=True)
-            # Typically, this should be done when computing the candidates, but just in case
-            # something illegal happended... if we can show that the current candidate set
-            # is optimal, we just keep it and don't compute the attractors symbolically.
-            node_is_pseudo_minimal = (not node["expanded"]) or self.node_is_minimal(
-                node_id
-            )
-            if len(candidates) == 0 or (
-                node_is_pseudo_minimal and len(candidates) == 1
-            ):
-                node["attractor_seeds"] = candidates
-                seeds = candidates
-            else:
-                result = compute_attractors_symbolic(
-                    self, node_id, candidate_states=candidates, seeds_only=True
+        try:
+            if seeds is None:
+                candidates = self.node_attractor_candidates(node_id, compute=True)
+                # Typically, this should be done when computing the candidates, but just in case
+                # something illegal happended... if we can show that the current candidate set
+                # is optimal, we just keep it and don't compute the attractors symbolically.
+                node_is_pseudo_minimal = (not node["expanded"]) or self.node_is_minimal(
+                    node_id
                 )
-                node["attractor_seeds"] = result[0]
-                # At this point, attractor_sets could be `None`, but that
-                # is valid, as long as we actually compute them later when
-                # they are needed.
-                node["attractor_sets"] = result[1]
-                seeds = result[0]
-
-            # Release memory once attractor seeds are known. We might need these
-            # for attractor set computation later, but only if the seeds are not empty.
+                if len(candidates) == 0 or (
+                    node_is_pseudo_minimal and len(candidates) == 1
+                ):
+                    node["attractor_seeds"] = candidates
+                    seeds = candidates
+                else:
+                    result = compute_attractors_symbolic(
+                        self, node_id, candidate_states=candidates, seeds_only=True
+                    )
+                    node["attractor_seeds"] = result[0]
+                    # At this point, attractor_sets could be `None`, but that
+                    # is valid, as long as we actually compute them later when
+                    # they are needed.
+                    node["attractor_sets"] = result[1]
+                    seeds = result[0]
+
+                # Release memory once attractor seeds are known. We might need these
+                # for attractor set computation later, but only if the seeds are not empty.
+                if len(seeds) == 0:
+                    node["percolated_network"] = None
+                    node["percolated_nfvs"] = None
+                    node["percolated_petri_net"] = None
+        except RuntimeError as e:
+            if not symbolic_fallback:
+                raise e
+
+            # The NFVS method failed, likely because the candidate set was too large.
+            # We can still try to fix this and compute the attractors symbolically.
+            (seeds, sets) = symbolic_attractor_fallback(self, node_id)
+
+            node["attractor_seeds"] = seeds
+            node["attractor_sets"] = sets
+
             if len(seeds) == 0:
                 node["percolated_network"] = None
                 node["percolated_nfvs"] = None

tests/succession_diagram_test.py

@@ -6,6 +6,7 @@
 import biobalm.succession_diagram
 from biobalm.succession_diagram import SuccessionDiagram
 from biobalm.types import BooleanSpace
+from biobalm._sd_attractors.attractor_symbolic import symbolic_attractor_fallback
 
 
 class SuccessionDiagramTest(unittest.TestCase):
@@ -171,13 +172,13 @@ def test_expansion_comparisons(network_file: str):
         # succession_diagram is too large for this test.
         return
 
+    assert sd_bfs.is_isomorphic(sd_dfs)
+
     # Normal minimal trap space expansion.
     sd_min = SuccessionDiagram(bn)
     assert sd_min.expand_minimal_spaces(size_limit=NODE_LIMIT)
 
-    assert sd_bfs.is_isomorphic(sd_dfs)
     assert sd_min.is_subgraph(sd_bfs)
-    assert sd_min.is_subgraph(sd_dfs)
 
     # Expand the first node fully, and then expand the rest
     # until minimal trap spaces are found.
@@ -188,7 +189,6 @@ def test_expansion_comparisons(network_file: str):
             assert sd_min_larger.expand_minimal_spaces(node_id=node_id)
 
     assert sd_min_larger.is_subgraph(sd_bfs)
-    assert sd_min_larger.is_subgraph(sd_dfs)
     assert len(sd_min_larger) >= len(sd_min)
 
     # Normal block expansion with size limit.
@@ -200,7 +200,6 @@ def test_expansion_comparisons(network_file: str):
     )
 
     assert sd_block.is_subgraph(sd_bfs)
-    assert sd_block.is_subgraph(sd_dfs)
 
     block_trap = [
         sd_block.node_data(i)["space"] for i in sd_block.minimal_trap_spaces()
@@ -211,6 +210,15 @@ def test_expansion_comparisons(network_file: str):
     for t in min_trap:
         assert t in block_trap
 
+    # Block expansion with symbolic fallback.
+    sd_config = SuccessionDiagram.default_config()
+    sd_fallback = SuccessionDiagram(bn, sd_config)
+    assert sd_fallback.expand_bfs()
+    for n in sd_fallback.node_ids():
+        symbolic = symbolic_attractor_fallback(sd_fallback, n)
+        normal = sd_fallback.node_attractor_seeds(n, compute=True)
+        assert len(normal) == len(symbolic[0])
+
     sd_attr = SuccessionDiagram(bn)
     assert sd_attr.expand_attractor_seeds(size_limit=NODE_LIMIT)