Add a fallback to the fully symbolic attractor detection if the NFVS …

…method 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