From b788ffc8248d6074d8c2735b0d1c200559844df5 Mon Sep 17 00:00:00 2001 From: Xuan Wang Date: Sun, 26 Apr 2026 22:59:16 -0400 Subject: [PATCH 1/2] fix: complete pinned_step body, restore driver-set search, with strict pcstg sufficiency check MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This PR completes Alex Gates' 2018 pinning-controllability design, fixes a long-standing bug in pinned_step, and restores two utility functions for minimum-driver-set search with a strict pcstg sufficiency criterion. The pinned_step bug fix ----------------------- The pcstg edge-loop scaffolding for limit-cycle attractors has been in place since Alex Gates' 2018 introduction of pinning controllability (98792ed): the iteration over attractor STG edges, the (attsource, attsink) pair, and the call-site form pinned_step(initial, pinned_binstate=attsink, pinned_var=...). What was missing was the body of pinned_step that would honor the pinned_binstate argument — pinned positions in the output were copied from initial instead of being written to pinned_binstate. For fixed-point attractors src and sink share the same pin pattern, so the missing-body path happened to produce correct results. For limit-cycle attractors where the pinned variable flips between cycle states, src and sink differ; the buggy body kept pin=src in the output, fragmenting the pcstg into one disjoint snapshot per cycle position. The cycle's states were split across snapshots, set(att) was never a subset of any single weakly-connected component, and fraction_pinned_configurations returned 0 for every such attractor. A related independent fix was made in 2019 (44a96d3, on python3-friendly) via a different mechanism — splitting pinned_step into two functions and changing the call site. That fix shipped on its branch but didn't make it to upstream when subsequent merges branched off pre-fix history. This commit completes Alex's original design: the body of pinned_step now writes pinned_binstate into the destination's pinned positions while still using the source pin pattern (from initial) to compute the unpinned step. The call site stays exactly as it was written in 2018. Changes ------- * pinned_step: body now reads pin from `initial` (as inputs to node update functions) and writes pin from `pinned_binstate`. Signature unchanged from Alex's original. Docstring rewritten to document all three arguments correctly (the original docstring had a phantom `n` arg from a copy-paste and never mentioned `pinned_binstate`). * pinning_controlled_state_transition_graph: rename loop variables attsource/attsink → src_pin/dst_pin. Same variables, just renamed — the old names suggested *attractor states* when they actually hold the *pin-bit projections* of those states. The rename is purely internal: these are local names inside one method's body, not part of any signature, attribute, or external API. An inline comment in the function explains the meaning and notes the rename. Restored search and Stage-2 design ---------------------------------- pinning_control_driver_nodes — minimum-driver-set search via exhaustive enumeration bracketed by ceil(log2(N_attractors)) lower bound and FVS upper bound. Discrete/Boolean analog of FVS-based open-loop control (Mochizuki & Fiedler 2013, JTB §7). The search uses a two-stage filter: * Stage 1 (necessary) — _signatures_distinguish_attractors. Cheap pre-check that the candidate's per-attractor signatures are distinct. Handles fixed-point, pin-constant cycle, and pin-flipping cycle attractors separately. * Stage 2 (strict sufficient) — attracting_components(pcstg) == [set(att)]. Each pcstg must have exactly one attracting SCC, equal to the target attractor's state set. Pure integer/set comparison — no floating-point exposure even at large pcstg sizes. The historically-tempting weaker criterion is WCC == 1 per pcstg. WCC=1 is necessary (no orphan basins) but not sufficient: it only verifies the pcstg is connected, not that the pcstg's attractor equals the intended target. For deterministic pcstg (fixed-point pinning) each WCC has exactly one attracting SCC, but that SCC may be a different attractor than the target if the pin pattern's projection collides with multiple original attractors and the unpinned dynamics converge to one of the others. Concrete example: Arabidopsis thaliana under pin = {AP3, UFO, AP1, LFY, WUS} — fpc = [1, 0, 1, 1, 1, 0, 1, 1, 1, 0] across the ten attractors, so three target attractors are not the pcstg's attractor under that pinning. Yet all ten pcstg are single-WCC. The size-6 sets Alex Gates' 2018 slide-3 reports for Thaliana pass the strict check; the WCC=1-only check would have admitted the size-5 set as a false positive. Drosophila and Yeast cell-cycle happen to be robust to the difference because their pin signatures distinguish all attractors uniquely. _signatures_distinguish_attractors — module-level helper used by the Stage-1 filter above. Three attractor types handled inline. A third orphaned function from a 2019 effort (full_control_driver_nodes) is deferred to a separate PR. The right name (configuration_perturbation_driver_nodes) needs lab consensus on the broader 4-quadrant control taxonomy (perturbation/pinning × attractor-level/configuration-level). Tests ----- 20 new tests in tests/test_pinning_control.py: - 18 tests on a 3-node BN fixture (A *= B; B *= A; C *= A and B) with all three attractor types simultaneously (length-2 cycle, two fixed points). Exercises: * bug surface (pin A or {A,B} — pin flips in cycle): pcstg single-WCC + mpcf=1.0; verifies the bug is fixed without tripping over fixed-point cases that incidentally worked pre-fix. * pin-constant control case (pin C — pin invariant in cycle): values unchanged from pre-fix; guards against regressing the fixed-point-equivalent case. * minimum driver-set search: returns {A,B} as the unique 2-element pinning control set (ceil(log2(3))=2 lower bound). * pinned_step unit tests: pin-unchanged (pinned_binstate matches initial[pinned_var]), pin-advanced (pinned_binstate differs), multi-pin, length-mismatch error. * _signatures_distinguish_attractors: four cases covering distinct fixed points, colliding fixed points, pin-constant cycle, and pin-flipping cycle. - 2 real-network regression tests on Arabidopsis thaliana (Chaos 2006), the network that motivated the strict Stage-2 criterion: * test_pinning_control_driver_nodes_thaliana_strict_criterion: asserts the search returns exactly Alex's three size-6 sets. * test_thaliana_size5_false_positive_rejected_by_strict_check: direct assertion on the strict criterion against the historical false-positive size-5 set; lighter-weight regression that survives even if the search algorithm changes shape. All 20 new tests pass; full CANA test suite (81 tests total) passes. Sanity check on real CC model: Fanconi anemia (length-2 cycle with CHKREC flipping), pinning CHKREC → pcf=1.0 (pcstg single WCC, correct full-control on the cycle). Lac Operon (fixed-point only) — unchanged. References ---------- - Fiedler, Mochizuki, Kurosawa, Saito (2013). Dynamics and Control at Feedback Vertex Sets I. J. Dyn. Diff. Equat. 25, 563-604. - Mochizuki, Fiedler, Kurosawa, Saito (2013). Dynamics and Control at Feedback Vertex Sets II. J. Theor. Biol. 335, 130-146. - Zañudo, Yang, Albert (2017). Structure-based control of complex networks with nonlinear dynamics. PNAS 114, 7234-7239. Co-Authored-By: Claude Opus 4.7 (1M context) --- cana/boolean_network.py | 225 ++++++++++++++++-- tests/test_pinning_control.py | 432 ++++++++++++++++++++++++++++++++++ 2 files changed, 636 insertions(+), 21 deletions(-) create mode 100644 tests/test_pinning_control.py diff --git a/cana/boolean_network.py b/cana/boolean_network.py index 0085e5c..4e61a21 100644 --- a/cana/boolean_network.py +++ b/cana/boolean_network.py @@ -13,6 +13,7 @@ # All rights reserved. # MIT license. from collections import defaultdict +from math import ceil, log2 try: import cStringIO.StringIO as StringIO # type: ignore @@ -51,6 +52,64 @@ from cana.utils import entropy, flip_binstate_bit_set, output_transitions +def _signatures_distinguish_attractors(candidate_nodes, bin_attractors): + """Necessary-condition pre-filter for pinning controllability. + + Returns ``True`` iff each attractor in ``bin_attractors`` has a + distinct *signature* on ``candidate_nodes``. The signature is the + tuple of values the candidate nodes take in the attractor: + + - **Fixed-point attractor** (length 1): signature is the tuple of + values at the single attractor state. + - **Limit-cycle attractor, pin constant within the cycle**: all + cycle states agree on the candidate nodes; the cycle is treated + as fixed-point-like with that single signature. + - **Limit-cycle attractor, pin flips within the cycle**: each + cycle state contributes its own signature; each must not + collide with another attractor's fixed-point-style signature. + + This is a *necessary* condition, not sufficient: two attractors + that pass here may still fail the downstream PCSTG-WCC check + (e.g., two flipping cycles whose per-state signatures happen to + overlap on the candidate nodes — rare in practice). Sufficiency + is verified in :func:`BooleanNetwork.pinning_control_driver_nodes`. + + Args: + candidate_nodes (list of int): node indices to check. + bin_attractors (list of list of str): attractors as lists of + binary state strings; ``bin_attractors[i][j][k]`` is the + value of node ``k`` at the ``j``-th state of attractor ``i``. + + Returns: + bool: ``True`` iff signatures distinguish all attractors; + ``False`` on any signature collision (including the trivial + ``len(candidate_nodes) == 0`` case where signatures are all + empty tuples). + """ + if len(candidate_nodes) == 0: + return False + fixed_signatures = set() + flipping_attractors = [] + for attr in bin_attractors: + sig = tuple(attr[0][node] for node in candidate_nodes) + is_pin_constant = all( + tuple(state[node] for node in candidate_nodes) == sig + for state in attr[1:] + ) + if len(attr) == 1 or is_pin_constant: + if sig in fixed_signatures: + return False + fixed_signatures.add(sig) + else: + flipping_attractors.append(attr) + for attr in flipping_attractors: + for state in attr: + sig = tuple(state[node] for node in candidate_nodes) + if sig in fixed_signatures: + return False + return True + + class BooleanNetwork: """ """ @@ -1050,6 +1109,99 @@ def attractor_driver_nodes(self, min_dvs=1, max_dvs=4, verbose=False): return attractor_controllers_found + def pinning_control_driver_nodes(self): + """Find minimum-size driver sets that achieve pinning control. + + A driver set ``D`` achieves *pinning control* if for every + attractor ``A`` of this Boolean network, pinning the nodes in + ``D`` to ``A``'s projection drives every initial configuration + to ``A``. Operationally, this requires the pinning-controlled + STG (``pcstg``) for each attractor to have exactly one + attracting strongly-connected component, and that SCC must + equal the target attractor's state set. + + The search starts at the information-theoretic lower bound + ``ceil(log2(N_attractors))`` (you need that many bits to + distinguish the attractors at all) and increments until at + least one valid driver set is found. All minimum-size valid + sets are returned. A two-stage filter is used: + + 1. *Necessary-condition pre-filter* + (:func:`_signatures_distinguish_attractors`): cheap + combinatorial check that the candidate's signatures + distinguish all attractors. Skips the expensive pcstg build + for obviously-invalid candidates. + 2. *Sufficiency check*: build pcstg per attractor; require + ``attracting_components(pcstg) == [set(att)]`` (a single + attracting SCC, equal to the target). The weaker + ``WCC == 1`` check used in earlier revisions only verifies + that the pcstg is connected — it does not verify that the + pcstg's unique attractor equals the intended target. For a + deterministic pcstg (fixed-point pinning) the pcstg always + has exactly one attracting SCC per WCC, but that SCC may + be some other state if the pin pattern happens to drive + the unpinned dynamics to a different attractor. The + Thaliana network exposed this in 2026; see test + ``test_thaliana_size5_false_positive_rejected_by_strict_check``. + + This is the discrete/Boolean analog of FVS-based open-loop + control proved for ODE systems in Mochizuki & Fiedler 2013 + ("Dynamics and Control at Feedback Vertex Sets II", JTB §7); + the Boolean version is supported by stable-motif theory + (Zañudo & Albert). FVS provides an upper bound on the driver- + set size; the search here may find smaller sets when the + discrete dynamics permit, since the FVS theorem requires a + guarantee for *all* nonlinearities while a specific Boolean + network may admit a smaller set. + + Returns: + list: minimum-size driver sets (as tuples of node indices) + that achieve pinning control. If no set up to size + ``Nnodes - 1`` works (degenerate networks), returns + ``[list(range(Nnodes))]`` as the trivial fallback. + + See also: + :func:`pinning_controlled_state_transition_graph`, + :func:`fraction_pinned_configurations`, + :func:`feedback_vertex_set_driver_nodes`, + :func:`_signatures_distinguish_attractors`. + """ + self._check_compute_variables(attractors=True) + if len(self._attractors) == 1: + return [] + lower_bound = ceil(log2(len(self._attractors))) + nodeids = list(range(self.Nnodes)) + bin_attractors = [ + [self.num2bin(state) for state in attr] for attr in self._attractors + ] + result = [] + for n_pin in range(lower_bound, self.Nnodes): + if result: + break + for pvs in itertools.combinations(nodeids, n_pin): + if not _signatures_distinguish_attractors( + list(pvs), bin_attractors + ): + continue + controlled = True + pcstg_dict = self.pinning_controlled_state_transition_graph( + list(pvs) + ) + for att, pcstg in pcstg_dict.items(): + # Strict check: pcstg must have exactly one + # attracting SCC, and that SCC must equal the + # target attractor's state set. Pure + # set/integer comparison — no floating point. + attracting = list(nx.attracting_components(pcstg)) + if len(attracting) != 1 or attracting[0] != set(att): + controlled = False + break + if controlled: + result.append(pvs) + if not result: + return [list(range(self.Nnodes))] + return result + def controlled_state_transition_graph(self, driver_nodes=[]): """Returns the Controlled State-Transition-Graph (CSTG). In practice, it copies the original STG, flips driver nodes (variables), and updates the CSTG. @@ -1124,6 +1276,16 @@ def pinning_controlled_state_transition_graph(self, driver_nodes=[]): pcstg_dict = {} for att in self._attractors: + # For each STG edge ``(s_src, s_dst)`` *inside* the attractor, + # ``src_pin`` and ``dst_pin`` are the projections of those + # states onto the pinned variables. For a fixed-point + # attractor the self-loop gives ``src_pin == dst_pin``; for a + # length-L cycle the L tuples have ``dst_pin`` rotated one + # cycle-step ahead of ``src_pin``. (These are the same loop + # variables previously named ``attsource`` and ``attsink``; + # renamed because the old names suggested *attractor states* + # when they actually hold the *pin-bit projections* of those + # states.) dn_attractor_transitions = [ tuple( "".join([self.num2bin(s)[dn] for dn in driver_nodes]) @@ -1136,12 +1298,12 @@ def pinning_controlled_state_transition_graph(self, driver_nodes=[]): self.bin2num( binstate_pinned_to_binstate( statenum_to_binstate(statenum, base=uncontrolled_system_size), - attsource, + src_pin, pinned_var=driver_nodes, ) ) for statenum in range(2**uncontrolled_system_size) - for attsource, attsink in dn_attractor_transitions + for src_pin, _dst_pin in dn_attractor_transitions ] pcstg = nx.DiGraph(name="STG: " + self.name) @@ -1155,19 +1317,27 @@ def pinning_controlled_state_transition_graph(self, driver_nodes=[]): pcstg.add_nodes_from((ps, {"label": ps}) for ps in pcstg_states) - for attsource, attsink in dn_attractor_transitions: + for src_pin, dst_pin in dn_attractor_transitions: for statenum in range(2**uncontrolled_system_size): initial = binstate_pinned_to_binstate( statenum_to_binstate(statenum, base=uncontrolled_system_size), - attsource, + src_pin, pinned_var=driver_nodes, ) + # ``pinned_step`` advances the unpinned variables + # using ``initial`` (which has ``src_pin`` at the + # pinned positions) and writes ``dst_pin`` at the + # pinned positions of the output. For fixed-point + # attractors ``src_pin == dst_pin`` so the pinned + # positions are unchanged; for cycles where the pin + # flips, this is what connects positions of the + # pcstg around the cycle. pcstg.add_edge( self.bin2num(initial), self.bin2num( self.pinned_step( initial, - pinned_binstate=attsink, + pinned_binstate=dst_pin, pinned_var=driver_nodes, ) ), @@ -1178,29 +1348,42 @@ def pinning_controlled_state_transition_graph(self, driver_nodes=[]): return pcstg_dict def pinned_step(self, initial, pinned_binstate, pinned_var): - """Steps the boolean network 1 step from the given initial input condition when the driver variables are pinned - to their controlled states. + """Advance the network one Boolean step under pinning control. + + Pinned variables are read as inputs to the node update + functions from ``initial`` (so the unpinned variables see the + *source* pin pattern when computing their next values), and + written as ``pinned_binstate`` in the output (the *destination* + pin pattern). For fixed-point pinning ``pinned_binstate`` is + the same pattern at every step; for limit-cycle pinning it + rotates one cycle position ahead of the source pin. Args: - initial (string) : the initial state. - n (int) : the number of steps. + initial (str) : the source binary state of length ``Nnodes``. + pinned_binstate (str) : destination values for the pinned + positions; must satisfy + ``len(pinned_binstate) == len(pinned_var)``. + pinned_var (list of int) : indices of the pinned variables. Returns: - (string) : The stepped binary state. + (str) : the next binary state, with pinned positions equal + to ``pinned_binstate`` and unpinned positions equal to + one Boolean step from ``initial`` (using the values in + ``initial`` — including the source pin — as inputs). + + See also: + :func:`pinning_controlled_state_transition_graph`. """ - # for every node: - # node input = breaks down initial by node input - # asks node to step with the input - # append output to list - # joins the results from each node output assert len(initial) == self.Nnodes + assert len(pinned_binstate) == len(pinned_var) + # Build a quick lookup so the comprehension is O(Nnodes) + # rather than O(Nnodes * |pinned_var|). + pin_map = dict(zip(pinned_var, pinned_binstate)) return "".join( - [ - str(node.step("".join(initial[j] for j in self.logic[i]["in"]))) - if not (i in pinned_var) - else initial[i] - for i, node in enumerate(self.nodes, start=0) - ] + pin_map[i] + if i in pin_map + else str(node.step("".join(initial[j] for j in self.logic[i]["in"]))) + for i, node in enumerate(self.nodes, start=0) ) def controlled_attractor_graph(self, driver_nodes=[]): diff --git a/tests/test_pinning_control.py b/tests/test_pinning_control.py new file mode 100644 index 0000000..d8158ca --- /dev/null +++ b/tests/test_pinning_control.py @@ -0,0 +1,432 @@ +"""Tests for pinning controllability: pcstg construction, mpcf, +and minimum driver-set search. + +Background +---------- +The pcstg edge-loop scaffolding for limit-cycle attractors has been +in place since Alex Gates' 2018 introduction of pinning +controllability (`98792ed`): the iteration over attractor STG edges, +the ``(attsource, attsink)`` pair, and the call-site form +``pinned_step(initial, pinned_binstate=attsink, pinned_var=...)``. +What was missing was the body of ``pinned_step`` that would honor +the ``pinned_binstate`` argument — pinned positions in the output +were copied from ``initial`` instead of being written to +``pinned_binstate``. For limit-cycle attractors where the pinned +variable flips, the pcstg fragmented into one disjoint snapshot per +cycle position; the cycle's states were split across snapshots, and +``fraction_pinned_configurations`` returned 0 for every such +attractor. + +A related independent fix was made in 2019 (`44a96d3`, on the +``python3-friendly`` branch) via a different mechanism — splitting +``pinned_step`` into two functions and changing the call site. That +fix shipped on its branch but didn't make it to upstream when +subsequent merges branched off pre-fix history. + +This test module covers the completion of Alex's original design: +the body of ``pinned_step`` now honors ``pinned_binstate``, plus two +restored utility functions for minimum driver-set search. + +Functions exercised: + +- ``pinned_step`` — advance the network one Boolean step under + pinning control (read pin from ``initial``, write pin from + ``pinned_binstate``, advance unpinned normally). +- ``pinning_controlled_state_transition_graph`` — pcstg should be a + single connected component for any controlled driver set. +- ``fraction_pinned_configurations`` / + ``mean_fraction_pinned_configurations`` — return correct values + for fixed-point and limit-cycle attractors. +- ``pinning_control_driver_nodes`` — minimum driver-set search. +- ``_signatures_distinguish_attractors`` — necessary-condition + pre-filter (handles fixed-point and limit-cycle attractor types). + +Test fixture +------------ +A 3-node Boolean network designed so that all three attractor types +coexist and a single pin choice exercises each scenario: + + A *= B + B *= A + C *= A and B + +Attractor structure (verified by ``BN.attractors(mode='stg')``): + + [[2, 4], [0], [7]] + = [{(010), (100)}, {(000)}, {(111)}] + = [length-2 cycle, fixed point, fixed point] + +Within the cycle, A and B each *flip* between cycle states (010 ↔ +100), while C stays at 0. So: + +- Pinning ``A`` (or ``B``) exposes the bug — pin flips in cycle. +- Pinning ``C`` is a control case — pin constant in cycle. +- Pinning ``{A, B}`` is the unique minimum pinning-control driver + set (3 attractors → ``ceil(log2(3)) = 2`` lower bound; only + ``{A, B}`` distinguishes all three attractor signatures). +""" + +import networkx as nx +import pytest + +from cana.boolean_network import BooleanNetwork + + +# ---------------------------------------------------------------- # +# Fixture # +# ---------------------------------------------------------------- # + + +def _build_two_fp_one_cycle(): + """Build the 3-node fixture and pre-compute attractors. + + Returns a ``BooleanNetwork`` with attractors already computed + via ``BN.attractors(mode='stg')`` so callers can use + ``BN._attractors`` directly. + """ + rules = "\n".join( + [ + "A *= B", + "B *= A", + "C *= A and B", + ] + ) + BN = BooleanNetwork.from_string_boolean(rules, keep_constants=True) + BN.attractors(mode="stg") + return BN + + +# Node indices for readability in tests. +A, B, C = 0, 1, 2 + +# Attractor positions in BN._attractors after _build_two_fp_one_cycle(). +ATT_CYCLE = 0 # [2, 4] = {(010), (100)} +ATT_FP_000 = 1 # [0] +ATT_FP_111 = 2 # [7] + + +# ---------------------------------------------------------------- # +# Sanity: the fixture has the expected attractor structure # +# ---------------------------------------------------------------- # + + +def test_fixture_has_expected_attractors(): + BN = _build_two_fp_one_cycle() + assert BN.Nnodes == 3 + # Order is implementation-defined but stable for this BN. + assert BN._attractors == [[2, 4], [0], [7]] + + +# ---------------------------------------------------------------- # +# pinned_step — unpinned advances; pin reads from initial, writes # +# from pinned_binstate # +# ---------------------------------------------------------------- # + + +def test_pinned_step_pin_unchanged_keeps_pin_at_source(): + """When ``pinned_binstate`` matches ``initial[pinned_var]`` (no pin + change), the pinned positions in the output equal those in + ``initial``; unpinned positions advance one Boolean step computed + from ``initial``.""" + BN = _build_two_fp_one_cycle() + # From '010' (A=0, B=1, C=0) with A pinned to its current value '0': + # B_next = A = 0; C_next = A and B = 0; A stays at '0'. + # Result: '000'. + assert BN.pinned_step("010", pinned_binstate="0", pinned_var=[A]) == "000" + # From '110' with A pinned to its current value '1': + # B_next = A = 1; C_next = A and B = 1; A stays at '1'. + # Result: '111'. + assert BN.pinned_step("110", pinned_binstate="1", pinned_var=[A]) == "111" + + +def test_pinned_step_pin_advanced_writes_destination_pin(): + """When ``pinned_binstate`` differs from ``initial[pinned_var]`` (pin + is being advanced — e.g., one cycle position to the next), the + pinned positions in the output equal ``pinned_binstate``; unpinned + positions still advance using values from ``initial`` (the source + pin pattern).""" + BN = _build_two_fp_one_cycle() + # From '010' (A=0, B=1, C=0) advancing A from '0' to '1': + # B_next = A_initial = 0; C_next = A_initial and B = 0; A → '1'. + # Result: '100'. + assert BN.pinned_step("010", pinned_binstate="1", pinned_var=[A]) == "100" + # From '100' (A=1, B=0, C=0) advancing A from '1' to '0': + # B_next = A_initial = 1; C_next = A_initial and B = 0; A → '0'. + # Result: '010'. + assert BN.pinned_step("100", pinned_binstate="0", pinned_var=[A]) == "010" + + +def test_pinned_step_multiple_pinned_vars(): + """Pinning multiple variables: pin pattern is written in the same + order as ``pinned_var``; unpinned positions advance from ``initial``.""" + BN = _build_two_fp_one_cycle() + # From '010' with {A, B} pinned to '10' (A→1, B→0): + # C_next = A_initial and B_initial = 0 and 1 = 0. + # A → '1', B → '0', C → '0'. Result: '100'. + assert BN.pinned_step("010", pinned_binstate="10", pinned_var=[A, B]) == "100" + + +def test_pinned_step_length_mismatch_raises(): + """``pinned_binstate`` length must match ``pinned_var`` length.""" + BN = _build_two_fp_one_cycle() + with pytest.raises(AssertionError): + BN.pinned_step("000", pinned_binstate="11", pinned_var=[A]) + + +# ---------------------------------------------------------------- # +# pcstg construction — single-WCC invariant under pinning control # +# ---------------------------------------------------------------- # + + +def test_pcstg_single_wcc_fixed_point_pin_A(): + """Fixed-point attractor: pcstg restricted to states with the + attractor's src_pin should still be one WCC.""" + BN = _build_two_fp_one_cycle() + pcstg_dict = BN.pinning_controlled_state_transition_graph([A]) + fp_pcstg = pcstg_dict[tuple(BN._attractors[ATT_FP_000])] + n_wccs = sum(1 for _ in nx.weakly_connected_components(fp_pcstg)) + # All 4 states with A=0 flow into (000) → single WCC. + assert n_wccs == 1 + assert len(fp_pcstg) == 4 + + +def test_pcstg_single_wcc_cycle_pin_flips_post_fix(): + """Length-2 cycle, pin flips around cycle. Post-fix the pcstg + should be a single WCC containing all 8 states (with edges that + cross between src_pin and dst_pin layers). + + Pre-fix this WCC count is 2 — the snapshot fragmentation bug. + Test asserts the post-fix invariant. + """ + BN = _build_two_fp_one_cycle() + pcstg_dict = BN.pinning_controlled_state_transition_graph([A]) + cyc_pcstg = pcstg_dict[tuple(BN._attractors[ATT_CYCLE])] + n_wccs = sum(1 for _ in nx.weakly_connected_components(cyc_pcstg)) + assert n_wccs == 1 + assert len(cyc_pcstg) == 8 + + +def test_pcstg_cycle_pin_constant_unaffected_by_fix(): + """Pin C: C is constant in the cycle, so pcstg construction + behaves identically pre-fix and post-fix. The cycle's pcstg + only contains the 4 states with C=0 (the cycle's pin value). + """ + BN = _build_two_fp_one_cycle() + pcstg_dict = BN.pinning_controlled_state_transition_graph([C]) + cyc_pcstg = pcstg_dict[tuple(BN._attractors[ATT_CYCLE])] + # 4 states with C=0; the cycle is one of three WCCs in this + # restricted pcstg (the others are FP basins for (000) and (110)). + assert len(cyc_pcstg) == 4 + n_wccs = sum(1 for _ in nx.weakly_connected_components(cyc_pcstg)) + assert n_wccs == 3 + # The cycle states must share a WCC. + cycle_states = set(BN._attractors[ATT_CYCLE]) + cycle_wcc = next( + wcc for wcc in nx.weakly_connected_components(cyc_pcstg) + if cycle_states <= wcc + ) + assert cycle_wcc == cycle_states + + +# ---------------------------------------------------------------- # +# fraction_pinned_configurations + mpcf # +# ---------------------------------------------------------------- # + + +def test_pcf_pin_A_post_fix(): + """Pin A. Pre-fix the cycle attractor returns 0 (bug). Post-fix + every attractor's basin reaches its target → all pcf = 1.0.""" + BN = _build_two_fp_one_cycle() + pcstg_dict = BN.pinning_controlled_state_transition_graph([A]) + pcf = BN.fraction_pinned_configurations(pcstg_dict) + # _attractors order: [cycle, fp000, fp111] + assert pcf == [1.0, 1.0, 1.0] + assert BN.mean_fraction_pinned_configurations(pcstg_dict) == 1.0 + + +def test_pcf_pin_C_unchanged(): + """Pin C: pin is constant in cycle, bug does not fire. Pre- and + post-fix should both give the same values.""" + BN = _build_two_fp_one_cycle() + pcstg_dict = BN.pinning_controlled_state_transition_graph([C]) + pcf = BN.fraction_pinned_configurations(pcstg_dict) + # Per-attractor: cycle has 2 states out of 4 in its WCC = 0.5; + # fixed points each have 1 of 4 = 0.25. + assert pcf == [0.5, 0.25, 0.25] + + +def test_pcf_pin_AB_post_fix(): + """Pin {A, B}: minimum pinning-control set. Both pinned vars + flip in cycle; post-fix gives full pinning control, mpcf = 1.0.""" + BN = _build_two_fp_one_cycle() + pcstg_dict = BN.pinning_controlled_state_transition_graph([A, B]) + pcf = BN.fraction_pinned_configurations(pcstg_dict) + assert pcf == [1.0, 1.0, 1.0] + + +def test_mpcf_consistency_with_fraction_pinned(): + """``mean_fraction_pinned_configurations`` is the mean of + ``fraction_pinned_configurations`` — keeps the contract simple.""" + BN = _build_two_fp_one_cycle() + pcstg_dict = BN.pinning_controlled_state_transition_graph([A]) + pcf = BN.fraction_pinned_configurations(pcstg_dict) + mpcf = BN.mean_fraction_pinned_configurations(pcstg_dict) + assert mpcf == pytest.approx(sum(pcf) / len(pcf)) + + +# ---------------------------------------------------------------- # +# pinning_control_driver_nodes — minimum driver-set search # +# ---------------------------------------------------------------- # + + +def test_pinning_control_driver_nodes_finds_unique_minimum(): + """For our 3-attractor fixture, the minimum pinning control set + is {A, B}. ``{A, C}`` and ``{B, C}`` fail the necessary-condition + pre-filter (signature collision between fp000 and the cycle). + """ + BN = _build_two_fp_one_cycle() + result = BN.pinning_control_driver_nodes() + # Expect exactly one minimum set: {A, B} (size 2 = ceil(log2(3))). + assert len(result) == 1 + assert set(result[0]) == {A, B} + + +def test_pinning_control_driver_nodes_size_above_log2_lower_bound(): + """The search starts at the info-theoretic lower bound + ``ceil(log2(N_attractors))`` and increments. For our fixture + (3 attractors), the lower bound is 2 and the answer matches it, + so size-1 candidates are not enumerated.""" + BN = _build_two_fp_one_cycle() + result = BN.pinning_control_driver_nodes() + for dvs in result: + assert len(dvs) == 2 + + +# ---------------------------------------------------------------- # +# Real-network regression: Arabidopsis thaliana (Chaos 2006) # +# ---------------------------------------------------------------- # +# +# Thaliana exposed a sufficiency-check bug: prior to this fix the +# search used ``WCC count == 1`` per pcstg, which is necessary but +# not sufficient — multiple original attractors can collapse into a +# single pcstg in which the target is reachable but not the unique +# attracting SCC. Returns a false-positive size-5 set +# ``{AP3, UFO, AP1, LFY, WUS}``. The strict criterion (unique +# attracting SCC equal to the target) recovers Alex Gates' 2018 +# slide result of three size-6 sets. +# +# Cost: ~15s. Sole network-scale regression test in this module. + + +def test_pinning_control_driver_nodes_thaliana_strict_criterion(): + """Search must apply the strict criterion (unique attracting SCC + equals target), not the weaker WCC==1. Returns Alex Gates' three + size-6 sets, not the historical size-5 false positive.""" + from cana.datasets.bio import THALIANA + + bn = THALIANA() + result = bn.pinning_control_driver_nodes() + names = [n.name for n in bn.nodes] + sets = {frozenset(names[i] for i in dvs) for dvs in result} + + # Alex Gates 2018 slide-3 result (three alternative size-6 sets, + # all sharing {UFO, WUS, AP3, AG} and extended by one of three pairs): + expected = { + frozenset({"UFO", "WUS", "AP3", "AG", "AP1", "LFY"}), + frozenset({"UFO", "WUS", "AP3", "AG", "TFL1", "EMF1"}), + frozenset({"UFO", "WUS", "AP3", "AG", "TFL1", "LFY"}), + } + assert sets == expected + assert all(len(dvs) == 6 for dvs in result) + + +def test_thaliana_size5_false_positive_rejected_by_strict_check(): + """Direct regression on the strict criterion: the historical + false-positive size-5 set ``{AP3, UFO, AP1, LFY, WUS}`` must fail + the strict per-attractor check. Targets at least one attractor + whose pcstg's unique attracting SCC differs from the target. + + Faster than the full search (~3s vs ~15s); exists so the + regression survives even if the search algorithm changes shape. + """ + from cana.datasets.bio import THALIANA + + bn = THALIANA() + names = [n.name for n in bn.nodes] + pin = sorted(names.index(n) for n in ["AP3", "UFO", "AP1", "LFY", "WUS"]) + pcstg_dict = bn.pinning_controlled_state_transition_graph(pin) + + # At least one attractor must fail the strict check. The pcstg's + # set of attracting SCCs must not equal {set(att)}. + failures = 0 + for att, pcstg in pcstg_dict.items(): + attracting = list(nx.attracting_components(pcstg)) + if not (len(attracting) == 1 and attracting[0] == set(att)): + failures += 1 + assert failures > 0, ( + "Strict criterion should reject the size-5 set on Thaliana; " + "if this assertion passes, the WCC==1-only check has been " + "reintroduced or the bug does not manifest on Thaliana anymore." + ) + + +# ---------------------------------------------------------------- # +# _signatures_distinguish_attractors — necessary-condition filter # +# ---------------------------------------------------------------- # + + +def test_signatures_distinguish_attractors_handles_fixed_point(): + """Two fixed points with distinct signatures on the candidate + nodes pass the necessary-condition test.""" + from cana.boolean_network import _signatures_distinguish_attractors + + # bin_attractors[i] is a list of binary state strings for attractor i. + fp1 = ["000"] + fp2 = ["111"] + assert _signatures_distinguish_attractors([0], [fp1, fp2]) + # Both signatures collapse to () when no nodes are pinned → + # cannot distinguish; return False. + assert not _signatures_distinguish_attractors([], [fp1, fp2]) + + +def test_signatures_distinguish_attractors_collision_fixed_point(): + """Two fixed points sharing a signature on the candidate nodes + cannot be distinguished by pinning.""" + from cana.boolean_network import _signatures_distinguish_attractors + + # Both attractors agree on bit 1 (B=0). + fp1 = ["000"] + fp2 = ["100"] + assert not _signatures_distinguish_attractors([1], [fp1, fp2]) + + +def test_signatures_distinguish_attractors_cycle_pin_constant(): + """A length-2 cycle whose pin is constant within the cycle is + treated as fixed-point-like: its signature must not collide with + any other attractor's signature on the candidate nodes. + + Setup: fixed point (000) has C=0; cycle ((010), (100)) has C + constant=0 throughout. On candidate ``[C]`` both attractors have + signature ``("0",)`` — collision. + """ + from cana.boolean_network import _signatures_distinguish_attractors + + fp = ["000"] + cycle_pin_constant = ["010", "100"] + assert not _signatures_distinguish_attractors([C], [fp, cycle_pin_constant]) + + +def test_signatures_distinguish_attractors_cycle_pin_flips(): + """A length-2 cycle whose pin flips within the cycle has its + per-state signatures checked against fixed-point signatures. + Per-state collision with a fixed-point fails the test. + + Setup: fixed point (010) has A=0; cycle ((010), (100)) has A=0 in + one state and A=1 in the other. On candidate ``[A]`` the cycle's + per-state signature ``("0",)`` collides with the fixed point's. + """ + from cana.boolean_network import _signatures_distinguish_attractors + + fp = ["010"] + cycle_pin_flips = ["010", "100"] + assert not _signatures_distinguish_attractors([A], [fp, cycle_pin_flips]) From 39295066ec71a859f990d5dce313b9b158202a7f Mon Sep 17 00:00:00 2001 From: Xuan Wang Date: Mon, 27 Apr 2026 19:32:43 -0400 Subject: [PATCH 2/2] Address code review feedback - Return [()] for single-attractor networks and tuple in fallback for consistent list-of-tuples return type - Exclude constant nodes from driver-set search (matching attractor_driver_nodes); add comment about modifying models if a constant node needs to be controllable - Bound search loop by len(nodeids) instead of Nnodes to avoid no-op iterations when constants are excluded - Clarify _signatures_distinguish_attractors docstring: collisions among flipping cycles are not checked (caught downstream) - Replace bare asserts in pinned_step with ValueError for clearer error messages - Look up attractors by content in tests instead of hard-coded indices to avoid dependence on implementation-defined ordering Co-Authored-By: Claude Opus 4.6 (1M context) --- cana/boolean_network.py | 54 +++++++++++++++++++++++++---------- tests/test_pinning_control.py | 33 ++++++++++++--------- 2 files changed, 59 insertions(+), 28 deletions(-) diff --git a/cana/boolean_network.py b/cana/boolean_network.py index 4e61a21..5630259 100644 --- a/cana/boolean_network.py +++ b/cana/boolean_network.py @@ -55,8 +55,8 @@ def _signatures_distinguish_attractors(candidate_nodes, bin_attractors): """Necessary-condition pre-filter for pinning controllability. - Returns ``True`` iff each attractor in ``bin_attractors`` has a - distinct *signature* on ``candidate_nodes``. The signature is the + Returns ``True`` iff no obvious signature collision is detected + among attractors on ``candidate_nodes``. The signature is the tuple of values the candidate nodes take in the attractor: - **Fixed-point attractor** (length 1): signature is the tuple of @@ -81,10 +81,13 @@ def _signatures_distinguish_attractors(candidate_nodes, bin_attractors): value of node ``k`` at the ``j``-th state of attractor ``i``. Returns: - bool: ``True`` iff signatures distinguish all attractors; - ``False`` on any signature collision (including the trivial - ``len(candidate_nodes) == 0`` case where signatures are all - empty tuples). + bool: ``True`` if no collision is found among fixed-point / + pin-constant signatures, and no flipping-cycle per-state + signature collides with a fixed-point signature. ``False`` + on any detected collision (including the trivial + ``len(candidate_nodes) == 0`` case). Note: collisions + *among* flipping cycles are not checked here — those are + caught by the downstream pcstg sufficiency check. """ if len(candidate_nodes) == 0: return False @@ -1155,10 +1158,12 @@ def pinning_control_driver_nodes(self): network may admit a smaller set. Returns: - list: minimum-size driver sets (as tuples of node indices) - that achieve pinning control. If no set up to size - ``Nnodes - 1`` works (degenerate networks), returns - ``[list(range(Nnodes))]`` as the trivial fallback. + list of tuple: minimum-size driver sets (as tuples of node + indices) that achieve pinning control. For a single- + attractor network, returns ``[()]`` (the empty driver set). + If no set up to size ``Nnodes - 1`` works (degenerate + networks), returns ``[tuple(range(Nnodes))]`` as the + trivial fallback. See also: :func:`pinning_controlled_state_transition_graph`, @@ -1168,14 +1173,25 @@ def pinning_control_driver_nodes(self): """ self._check_compute_variables(attractors=True) if len(self._attractors) == 1: - return [] + return [()] lower_bound = ceil(log2(len(self._attractors))) nodeids = list(range(self.Nnodes)) + # Exclude constant nodes: they cannot distinguish attractors + # and waste combinatorial search effort. If you need to treat + # a constant node as a controllable driver (e.g. toggling a + # stimulus), modify the model to make it non-constant before + # calling this function. + if self.keep_constants: + constant_nodeids = set(self.get_constants().keys()) + nodeids = [nodeid for nodeid in nodeids if nodeid not in constant_nodeids] bin_attractors = [ [self.num2bin(state) for state in attr] for attr in self._attractors ] result = [] - for n_pin in range(lower_bound, self.Nnodes): + max_pin = len(nodeids) + if lower_bound > max_pin: + return [tuple(range(self.Nnodes))] + for n_pin in range(lower_bound, max_pin + 1): if result: break for pvs in itertools.combinations(nodeids, n_pin): @@ -1199,7 +1215,7 @@ def pinning_control_driver_nodes(self): if controlled: result.append(pvs) if not result: - return [list(range(self.Nnodes))] + return [tuple(range(self.Nnodes))] return result def controlled_state_transition_graph(self, driver_nodes=[]): @@ -1374,8 +1390,16 @@ def pinned_step(self, initial, pinned_binstate, pinned_var): See also: :func:`pinning_controlled_state_transition_graph`. """ - assert len(initial) == self.Nnodes - assert len(pinned_binstate) == len(pinned_var) + if len(initial) != self.Nnodes: + raise ValueError( + "initial state length must equal Nnodes: " + "expected %d, got %d" % (self.Nnodes, len(initial)) + ) + if len(pinned_binstate) != len(pinned_var): + raise ValueError( + "pinned_binstate length must match pinned_var: " + "expected %d, got %d" % (len(pinned_var), len(pinned_binstate)) + ) # Build a quick lookup so the comprehension is O(Nnodes) # rather than O(Nnodes * |pinned_var|). pin_map = dict(zip(pinned_var, pinned_binstate)) diff --git a/tests/test_pinning_control.py b/tests/test_pinning_control.py index d8158ca..132e9b1 100644 --- a/tests/test_pinning_control.py +++ b/tests/test_pinning_control.py @@ -96,14 +96,18 @@ def _build_two_fp_one_cycle(): return BN +def _find_attractor(BN, states): + """Find the index of an attractor by its state set (order-independent).""" + target = set(states) + for i, attr in enumerate(BN._attractors): + if set(attr) == target: + return i + raise ValueError(f"Attractor {states} not found in {BN._attractors}") + + # Node indices for readability in tests. A, B, C = 0, 1, 2 -# Attractor positions in BN._attractors after _build_two_fp_one_cycle(). -ATT_CYCLE = 0 # [2, 4] = {(010), (100)} -ATT_FP_000 = 1 # [0] -ATT_FP_111 = 2 # [7] - # ---------------------------------------------------------------- # # Sanity: the fixture has the expected attractor structure # @@ -113,8 +117,11 @@ def _build_two_fp_one_cycle(): def test_fixture_has_expected_attractors(): BN = _build_two_fp_one_cycle() assert BN.Nnodes == 3 - # Order is implementation-defined but stable for this BN. - assert BN._attractors == [[2, 4], [0], [7]] + # Check attractor content, not ordering (order is implementation-defined). + att_sets = [set(a) for a in BN._attractors] + assert {2, 4} in att_sets # length-2 cycle + assert {0} in att_sets # fixed point 000 + assert {7} in att_sets # fixed point 111 # ---------------------------------------------------------------- # @@ -169,7 +176,7 @@ def test_pinned_step_multiple_pinned_vars(): def test_pinned_step_length_mismatch_raises(): """``pinned_binstate`` length must match ``pinned_var`` length.""" BN = _build_two_fp_one_cycle() - with pytest.raises(AssertionError): + with pytest.raises(ValueError): BN.pinned_step("000", pinned_binstate="11", pinned_var=[A]) @@ -183,7 +190,7 @@ def test_pcstg_single_wcc_fixed_point_pin_A(): attractor's src_pin should still be one WCC.""" BN = _build_two_fp_one_cycle() pcstg_dict = BN.pinning_controlled_state_transition_graph([A]) - fp_pcstg = pcstg_dict[tuple(BN._attractors[ATT_FP_000])] + fp_pcstg = pcstg_dict[tuple(BN._attractors[_find_attractor(BN, [0])])] n_wccs = sum(1 for _ in nx.weakly_connected_components(fp_pcstg)) # All 4 states with A=0 flow into (000) → single WCC. assert n_wccs == 1 @@ -200,7 +207,7 @@ def test_pcstg_single_wcc_cycle_pin_flips_post_fix(): """ BN = _build_two_fp_one_cycle() pcstg_dict = BN.pinning_controlled_state_transition_graph([A]) - cyc_pcstg = pcstg_dict[tuple(BN._attractors[ATT_CYCLE])] + cyc_pcstg = pcstg_dict[tuple(BN._attractors[_find_attractor(BN, [2, 4])])] n_wccs = sum(1 for _ in nx.weakly_connected_components(cyc_pcstg)) assert n_wccs == 1 assert len(cyc_pcstg) == 8 @@ -213,14 +220,14 @@ def test_pcstg_cycle_pin_constant_unaffected_by_fix(): """ BN = _build_two_fp_one_cycle() pcstg_dict = BN.pinning_controlled_state_transition_graph([C]) - cyc_pcstg = pcstg_dict[tuple(BN._attractors[ATT_CYCLE])] + cyc_pcstg = pcstg_dict[tuple(BN._attractors[_find_attractor(BN, [2, 4])])] # 4 states with C=0; the cycle is one of three WCCs in this # restricted pcstg (the others are FP basins for (000) and (110)). assert len(cyc_pcstg) == 4 n_wccs = sum(1 for _ in nx.weakly_connected_components(cyc_pcstg)) assert n_wccs == 3 # The cycle states must share a WCC. - cycle_states = set(BN._attractors[ATT_CYCLE]) + cycle_states = set(BN._attractors[_find_attractor(BN, [2, 4])]) cycle_wcc = next( wcc for wcc in nx.weakly_connected_components(cyc_pcstg) if cycle_states <= wcc @@ -239,7 +246,7 @@ def test_pcf_pin_A_post_fix(): BN = _build_two_fp_one_cycle() pcstg_dict = BN.pinning_controlled_state_transition_graph([A]) pcf = BN.fraction_pinned_configurations(pcstg_dict) - # _attractors order: [cycle, fp000, fp111] + # All attractors fully controlled regardless of ordering. assert pcf == [1.0, 1.0, 1.0] assert BN.mean_fraction_pinned_configurations(pcstg_dict) == 1.0