Example Code

import numpy as np
from itertools import product

class BooleanNetworkNP:
    """
    NumPy Boolean network with synchronous updates.
    State is a 1-D np.bool_ array s of length n.
    Update rules: list of n callables: f_i(s) -> bool (use &, |, ^, and ~ or 1- for NOT).
    """
    def __init__(self, n, update_funcs):
        self.n = int(n)
        if len(update_funcs) != n:
            raise ValueError(f"Expected {n} update funcs, got {len(update_funcs)}")
        self.f = list(update_funcs)

    @staticmethod
    def as_bool_array(x):
        a = np.asarray(x, dtype=bool)
        if a.ndim != 1:
            raise ValueError("State must be 1-D")
        return a

    @staticmethod
    def bits(s):
        return '(' + ','.join('1' if b else '0' for b in s) + ')'

    @staticmethod
    def bits_flat(s):
        return ''.join('1' if b else '0' for b in s)

    @staticmethod
    def lex_index1(s):
        return int(''.join('1' if v else '0' for v in s), 2) + 1

    @staticmethod
    def key_bytes(s):
        # compact, hashable key for dicts
        return np.packbits(s, bitorder='big').tobytes()

    def step(self, s):
        # ---- OPTIMIZED: bind attributes to locals; avoid repeated lookups in loop ----
        as_bool = self.as_bool_array
        funcs   = self.f
        n       = self.n

        s = as_bool(s)
        out = np.empty(n, dtype=bool)
        out_set = out.__setitem__       # localize setter
        # localize range for micro gain
        for i in range(n):
            out_set(i, bool(funcs[i](s)))
        return out

    def run_until_cycle_with_trace(self, s0, max_steps=100_000):
        """
        From initial state s0 (bool vector or list/tuple of 0/1), run until a state repeats.
        Returns:
          transient_len (int),
          cycle (tuple of states, each np.bool_ vector),
          trace (list of pairs: (input_state, output_state) for each step until the first repeat).
        """
        # ---- OPTIMIZED: bind helpers & step to locals ----
        as_bool    = self.as_bool_array
        step       = self.step
        key_bytes  = self.key_bytes
        arr_equal  = np.array_equal

        s = as_bool(s0)
        seen  = {}     # key_bytes(state) -> time index
        trace = []     # [(s_t, s_{t+1}), ...] until repeat
        t = 0
        while t < max_steps:
            k = key_bytes(s)
            if k in seen:
                # s is first repeated state; rebuild cycle starting from s
                cyc = [s.copy()]
                s2 = step(s)
                while not arr_equal(s2, s):
                    cyc.append(s2.copy())
                    s2 = step(s2)
                return seen[k], tuple(cyc), trace
            seen[k] = t
            nxt = step(s)
            trace.append((s.copy(), nxt.copy()))
            s = nxt
            t += 1
        raise RuntimeError("Cycle not found within max_steps")

# ---------- Canonicalization helpers for attractor identity ----------
def canonical_cycle_key(cycle_states):
    """
    Given a tuple/list of states (np.bool_ vectors) forming a cycle in order,
    return a canonical, hashable key (bytes) invariant to rotation (phase).
    """
    # Represent each state as its packed bytes
    byte_states = [np.packbits(s, bitorder='big').tobytes() for s in cycle_states]
    # Consider all rotations; choose lexicographically smallest rotation as canonical
    k = len(byte_states)
    rotations = [tuple(byte_states[i:]+byte_states[:i]) for i in range(k)]
    return min(rotations)

# ---------- Pretty printers for Option B ----------
def print_trajectory_until_cycle(bn, s0, attractor_registry):
    """
    For one initial state s0, print the full trajectory until first cycle,
    showing each step as: f(input) -> output, with both binary and y_i labels.
    Use attractor_registry (dict) to assign stable IDs to cycles across starts.
    """
    # ---- OPTIMIZED: one-time local binds ----
    as_bool = bn.as_bool_array
    bits    = bn.bits
    bits_f  = bn.bits_flat
    lex1    = bn.lex_index1

    mu, cycle, trace = bn.run_until_cycle_with_trace(s0)

    # Determine attractor ID (stable across starts)
    key = canonical_cycle_key(cycle)
    if key not in attractor_registry:
        attractor_registry[key] = len(attractor_registry) + 1
    att_id = attractor_registry[key]

    s0_arr  = as_bool(s0)
    print(f"\nStart {bits(s0_arr)}  [y_{lex1(s0_arr)}]")
    print("Trajectory (synchronous):")
    for t, (inp, outp) in enumerate(trace):
        print(
            f"  t={t:>3d}:  f{bits(inp)} [y_{lex1(inp)}]  ->  "
            f"{bits(outp)} [y_{lex1(outp)}]"
        )

    cyc_bits   = [bits_f(s) for s in cycle]
    cyc_labels = [f"y_{lex1(s)}" for s in cycle]
    print(f"Attractor #{att_id}: cycle length = {len(cycle)}")
    print("  cycle (order): " + " -> ".join(f"{lab}({b})" for lab, b in zip(cyc_labels, cyc_bits)))

    if mu > 0:
        print(f"  transient length: {mu}")
    else:
        print("  transient length: 0 (started inside the cycle)")

def sweep_all_initial_states(bn, order='lex'):
    """
    Run trajectories from all 2^n initial states.
    order='lex' means lexicographic order (000..0, 000..1, ..., 111..1).
    """
    # ---- OPTIMIZED: bind utilities once for the whole sweep ----
    as_bool = bn.as_bool_array
    product_ = product  # local alias

    attractor_registry = {}  # canonical_key -> attractor_id
    states = product_((0,1), repeat=bn.n) if order == 'lex' else product_((0,1), repeat=bn.n)
    for tup in states:
        # np.fromiter is marginally faster than np.asarray for tiny tuples; either is fine
        s0 = np.fromiter(tup, dtype=bool, count=bn.n)
        print_trajectory_until_cycle(bn, s0, attractor_registry)

# ============================
# Minimal example (edit freely)
# ============================

# Define rules using indices s[0], s[1], ...
# AND: &   OR: |   XOR: ^   NOT: ~s[i] or (1 - s[i])
rules = [
    lambda s: s[1] ^ s[3],        # x1(t+1) = x2 XOR x4
    lambda s: s[0] | s[2],        # x2(t+1) = x1 OR  x3
    lambda s: s[0] & s[1],        # x3(t+1) = x1 AND x2
    lambda s: (1 - s[3]) ^ s[2],  # x4(t+1) = NOT x4 XOR x3
]
bn = BooleanNetworkNP(n=4, update_funcs=rules)

# Run Option B for EVERY initial state (all 2^n starts)
# (For larger n, this prints a lot; start with n<=6 or so)
sweep_all_initial_states(bn)