Refactor SCC expansion

biobalm/_sd_attractors/attractor_candidates.py

@@ -302,18 +302,22 @@ def compute_attractor_candidates(
 
                 # At this point, we know the candidate count increased and so we should
                 # try to bring it back down.
-                if sd.config["debug"]:
-                    print(f"[{node_id}] Optimizing partial retained set...")
-                optimized = asp_greedy_retained_set_optimization(
-                    sd,
-                    node_id,
-                    petri_net=pn_reduced,
-                    retained_set=retained_set,
-                    candidate_states=candidate_states,
-                    avoid_dnf=child_motifs_reduced,
-                )
-                retained_set = optimized[0]
-                candidate_states = optimized[1]
+                if (
+                    len(candidate_states)
+                    > sd.config["retained_set_optimization_threshold"]
+                ):
+                    if sd.config["debug"]:
+                        print(f"[{node_id}] Optimizing partial retained set...")
+                    optimized = asp_greedy_retained_set_optimization(
+                        sd,
+                        node_id,
+                        petri_net=pn_reduced,
+                        retained_set=retained_set,
+                        candidate_states=candidate_states,
+                        avoid_dnf=child_motifs_reduced,
+                    )
+                    retained_set = optimized[0]
+                    candidate_states = optimized[1]
 
     # Terminate if done.
     if len(candidate_states) == 0:

biobalm/interaction_graph_utils.py

@@ -14,7 +14,7 @@
 
 from typing import cast
 
-from biodivine_aeon import BooleanNetwork, RegulatoryGraph, SignType
+from biodivine_aeon import BooleanNetwork, RegulatoryGraph, SymbolicContext, SignType
 from networkx import DiGraph  # type: ignore
 
 
@@ -181,3 +181,34 @@ def source_SCCs(bn: BooleanNetwork) -> list[list[str]]:
             result.append(scc_names)
 
     return sorted(result)
+
+
+def source_nodes(
+    network: BooleanNetwork, ctx: SymbolicContext | None = None
+) -> list[str]:
+    """
+    Return the source nodes of a `BooleanNetwork`. That is, variables whose value
+    cannot change, but is not fixed to a `true`/`false` constant.
+
+    Note that this internally uses BDD translation to detect identity functions
+    semantically rather than syntactically. If you already have a `SymbolicContext`
+    for the given `network` available, you can supply it as the second argument.
+    """
+    if ctx is None:
+        ctx = SymbolicContext(network)
+
+    result: list[str] = []
+    for var in network.variable_names():
+        update_function = network.get_update_function(var)
+        if update_function is None:
+            # This is an input variable with unspecified update
+            # (this defaults to identity in most tools).
+            assert len(network.predecessors(var)) == 0
+            result.append(var)
+        else:
+            fn_bdd = ctx.mk_update_function(update_function)
+            var_bdd = ctx.mk_network_variable(var)
+            if fn_bdd == var_bdd:
+                result.append(network.get_variable_name(var))
+
+    return result

biobalm/succession_diagram.py

@@ -5,8 +5,9 @@
 if TYPE_CHECKING:
     from typing import Iterator
 
+import copy
 import networkx as nx  # type: ignore
-from biodivine_aeon import AsynchronousGraph, BooleanNetwork, VariableId, VertexSet
+from biodivine_aeon import AsynchronousGraph, BooleanNetwork, VertexSet
 
 # Attractor detection algorithms.
 from biobalm._sd_attractors.attractor_candidates import compute_attractor_candidates
@@ -173,7 +174,7 @@ def default_config() -> SuccessionDiagramConfiguration:
             "nfvs_size_threshold": 2_000,
             "pint_goal_size_limit": 8_192,
             "attractor_candidates_limit": 100_000,
-            "retained_set_optimization_threshold": 100,
+            "retained_set_optimization_threshold": 1_000,
         }
 
     @staticmethod
@@ -1124,47 +1125,16 @@ def component_subdiagrams(
         if node_id is None:
             node_id = self.root()
 
-        reference_bn = percolate_network(
-            self.network,
-            self.node_data(node_id)["space"],
-            remove_constants=True,
-        )
+        reference_bn = self.node_percolated_network(node_id, compute=True)
 
         source_scc_list = source_SCCs(reference_bn)
 
         for component_variables in source_scc_list:
-            new_bn = BooleanNetwork(component_variables)
-
-            # Build a mapping between the old and new network variables.
-            id_map: dict[VariableId, VariableId] = {}
-            for var in component_variables:
-                old_id = reference_bn.find_variable(var)
-                assert old_id is not None
-                new_id = new_bn.find_variable(var)
-                assert new_id is not None
-                id_map[old_id] = new_id
-
-            # Copy regulations that are in the source component.
-            for reg in reference_bn.regulations():
-                if reg["source"] in id_map and reg["target"] in id_map:
-                    new_bn.add_regulation(
-                        {
-                            "source": reference_bn.get_variable_name(reg["source"]),
-                            "target": reference_bn.get_variable_name(reg["target"]),
-                            "essential": reg["essential"],
-                            "sign": reg["sign"],
-                        }
-                    )
-
-            # Copy update functions from the source component after translating them
-            # to the new IDs.
-            for var_id in id_map.keys():
-                old_function = reference_bn.get_update_function(var_id)
-                assert old_function is not None
-                new_function = old_function.rename_all(new_bn, variables=id_map)
-                new_bn.set_update_function(id_map[var_id], new_function)
-
-            yield SuccessionDiagram(new_bn)
+            remaining_variables = [
+                v for v in reference_bn.variable_names() if v not in component_variables
+            ]
+            new_bn = reference_bn.drop(remaining_variables)
+            yield SuccessionDiagram(new_bn, copy.copy(self.config))
 
     def build(self):
         """
@@ -1436,11 +1406,7 @@ def _ensure_node(self, parent_id: int | None, stable_motif: BooleanSpace) -> int
         assert child_id is not None
 
         if parent_id is not None:
-            # TODO: It seems that there are some networks where the same child
-            # can be reached through multiple stable motifs. Not sure how to
-            # approach these... but this is probably good enough for now.
-            self.dag.add_edge(parent_id, child_id, motif=stable_motif)  # type: ignore
-            self._update_node_depth(child_id, parent_id)
+            self._ensure_edge(parent_id, child_id, stable_motif)
 
         # Compute the percolated petri net here, because we know the parent node ID
         # and using its already percolated petri net helps a lot.
@@ -1450,3 +1416,10 @@ def _ensure_node(self, parent_id: int | None, stable_motif: BooleanSpace) -> int
         self.node_percolated_petri_net(child_id, compute=True, parent_id=parent_id)
 
         return child_id
+
+    def _ensure_edge(self, parent_id: int, child_id: int, stable_motif: BooleanSpace):
+        # TODO: It seems that there are some networks where the same child
+        # can be reached through multiple stable motifs. Not sure how to
+        # approach these... but this is probably good enough for now.
+        self.dag.add_edge(parent_id, child_id, motif=stable_motif)  # type: ignore
+        self._update_node_depth(child_id, parent_id)

pyproject.toml

@@ -6,7 +6,7 @@ readme = "README.md"
 license = { file = "LICENSE" }
 requires-python = ">=3.11"
 dependencies = [
-    'biodivine_aeon ==1.0.0a6',
+    'biodivine_aeon ==1.0.0a8',
     'clingo >=5.6.2',
     'networkx >=2.8.8',    
 ]

requirements.txt

@@ -1,4 +1,4 @@
-biodivine_aeon==1.0.0a6
+biodivine_aeon==1.0.0a8
 clingo==5.6.2
 networkx==2.8.8
 pypint[pint]==1.6.2

tests/interaction_graph_utils_test.py

@@ -1,7 +1,8 @@
-from biodivine_aeon import BooleanNetwork
+from biodivine_aeon import BooleanNetwork, AsynchronousGraph
 from networkx import DiGraph  # type:ignore
 
-from biobalm.interaction_graph_utils import feedback_vertex_set
+from biobalm.interaction_graph_utils import feedback_vertex_set, source_nodes
+from biobalm.space_utils import percolate_network, percolate_space
 
 # There should be a negative cycle between b_1 and b_2,
 # a positive cycle between d_1 and d_2, and a negative cycle
@@ -41,6 +42,32 @@
 )
 
 
+def test_source_nodes():
+    bn = BooleanNetwork.from_bnet(
+        """targets,factors
+    constant1_1, (constant1_1 | !constant1_1)
+    constant1_0, (constant1_0 & !constant1_0)
+    constant2_1, true
+    constant2_0, false
+    source, source
+    oscillator, !oscillator
+    source_after_perc, source_after_perc & constant1_1
+    after_perc_0, after_perc_0 & constant1_0"""
+    ).infer_valid_graph()
+    graph = AsynchronousGraph(bn)
+
+    sources = source_nodes(bn)
+
+    assert sources == ["source"]
+
+    perc_space = percolate_space(graph, {})
+    perc_bn = percolate_network(bn, perc_space, symbolic_network=graph)
+
+    sources = source_nodes(perc_bn)
+
+    assert sources == ["source", "source_after_perc"]
+
+
 def test_fvs():
     fvs = feedback_vertex_set(CYCLES_BN)
     nfvs = feedback_vertex_set(CYCLES_BN, parity="negative")

tests/source_SCC_test.py

@@ -1,38 +1,8 @@
-from biodivine_aeon import AsynchronousGraph, BooleanNetwork
+from biodivine_aeon import BooleanNetwork
 
 from biobalm import SuccessionDiagram
-from biobalm._sd_algorithms.expand_source_SCCs import (
-    expand_source_SCCs,
-    find_source_nodes,
-    find_subnetwork_sd,
-)
-from biobalm.space_utils import percolate_network, percolate_space
-
-
-def test_find_source_nodes():
-    bn = BooleanNetwork.from_bnet(
-        """targets,factors
-    constant1_1, (constant1_1 | !constant1_1)
-    constant1_0, (constant1_0 & !constant1_0)
-    constant2_1, true
-    constant2_0, false
-    source, source
-    oscillator, !oscillator
-    source_after_perc, source_after_perc & constant1_1
-    after_perc_0, after_perc_0 & constant1_0"""
-    ).infer_valid_graph()
-    graph = AsynchronousGraph(bn)
-
-    source_nodes = find_source_nodes(bn)
-
-    assert source_nodes == ["source"]
-
-    perc_space = percolate_space(graph, {})
-    perc_bn = percolate_network(bn, perc_space, symbolic_network=graph)
-
-    source_nodes = find_source_nodes(perc_bn)
-
-    assert source_nodes == ["source", "source_after_perc"]
+from biobalm._sd_algorithms.expand_source_SCCs import expand_source_SCCs
+from biobalm.space_utils import percolate_network
 
 
 def test_perc_and_remove_constants_from_bn():
@@ -62,30 +32,6 @@ def test_perc_and_remove_constants_from_bn():
     assert clean_bnet == clean_bnet2
 
 
-def test_find_scc_sd():
-    bn = BooleanNetwork.from_bnet(
-        """targets,factors
-A, B
-B, A | C"""
-    )
-
-    # This "simulates" what would happen in the SCC expansion algorithm.
-    bn = percolate_network(bn, {"C": 0}, remove_constants=True)
-
-    scc_sd, _ = find_subnetwork_sd(
-        SuccessionDiagram(bn),
-        expander=SuccessionDiagram.expand_bfs,
-        check_maa=True,
-    )
-
-    for node in scc_sd.node_ids():
-        print(scc_sd.node_data(node)["space"])
-
-    assert scc_sd.dag.nodes[0]["space"] == {}
-    assert scc_sd.dag.nodes[1]["space"] == {"A": 0, "B": 0}
-    assert scc_sd.dag.nodes[2]["space"] == {"A": 1, "B": 1}
-
-
 def expansion(bn: BooleanNetwork):
     sd = SuccessionDiagram(bn)
     fully_expanded = expand_source_SCCs(sd, check_maa=False)
@@ -145,7 +91,7 @@ def test_expansion():
 
 def attractor_search(bn: BooleanNetwork):
     sd = SuccessionDiagram(bn)
-    fully_expanded = expand_source_SCCs(sd)
+    fully_expanded = sd.expand_scc()
     assert fully_expanded
 
     attractor_count = 0
@@ -300,7 +246,7 @@ def test_attractor_search():
     )
     n, size, depth, att, maa, min = attractor_search(bn)
     assert n == 6
-    assert size == 4
+    assert size == 3
     assert depth == 2
     assert att == 4
     assert maa == 3
@@ -318,7 +264,7 @@ def test_attractor_search():
     )
     n, size, depth, att, maa, min = attractor_search(bn)
     assert n == 6
-    assert size == 7
+    assert size == 6
     assert depth == 3
     assert att == 5
     assert maa == 2
@@ -343,7 +289,7 @@ def test_attractor_search():
     )
     n, size, depth, att, maa, min = attractor_search(bn)
     assert n == 13
-    assert size == 41
+    assert size == 28
     assert depth == 6
     assert att == 28
     assert maa == 14
@@ -363,6 +309,37 @@ def test_attractor_search():
     assert maa == 0
     assert min == 2
 
+    # interesting example uncovered during testing (random-nk2/n20_29.bnet)
+
+    bn = BooleanNetwork.from_bnet(
+        """targets, factors
+    n0, (n12 & !n5) | (n12 & n5)
+    n1, (!n1 & !n15) | (n1 & !n15) | (n1 & n15)
+    n2, (n3 & !n15) | (n3 & n15)
+    n3, (!n17 & n3) | (n17 & !n3)
+    n4, (n0 & !n18) | (n0 & n18)
+    n5, (!n13 & !n16) | (!n13 & n16) | (n13 & n16)
+    n6, true
+    n7, (!n10 & !n1) | (n10 & !n1) | (n10 & n1)
+    n8, true
+    n9, (!n0 & !n19) | (!n0 & n19) | (n0 & !n19)
+    n10, (!n2 & !n17) | (!n2 & n17)
+    n11, (n16 & !n1)
+    n12, (!n3 & n12) | (n3 & n12)
+    n13, (!n8 & n3) | (n8 & !n3)
+    n14, (!n3 & !n15) | (n3 & !n15)
+    n15, (!n14 & !n15) | (!n14 & n15)
+    n16, (!n15 & !n11) | (n15 & !n11)
+    n17, (n2 & !n11) | (n2 & n11)
+    n18, (!n7 & !n18) | (!n7 & n18)
+    n19, (!n5 & !n12) | (!n5 & n12)"""
+    )
+    n, size, depth, att, maa, min = attractor_search(bn)
+    assert size == 15
+    assert min == 6
+    assert att == 6
+    assert maa == 0
+
 
 def test_isomorph():
     path = "models/bbm-bnet-inputs-true/005.bnet"
@@ -372,7 +349,7 @@ def test_isomorph():
     sd_bfs.expand_bfs()
 
     sd_scc = SuccessionDiagram(bn)
-    expand_source_SCCs(sd_scc)
+    sd_scc.expand_scc()
 
     assert [sd_bfs.node_data(id)["space"] for id in sd_bfs.node_ids()] == [
         sd_scc.node_data(id)["space"] for id in sd_scc.node_ids()