From c1ed30bc641b6310a7bef67beb7f885db3665120 Mon Sep 17 00:00:00 2001
From: Peter Call <157658209+PlanetHopf@users.noreply.github.com>
Date: Thu, 9 Oct 2025 09:55:42 -0700
Subject: [PATCH] Create Multilayer_Network.ipynb

Just some code I've worked on
---
 Multilayer_Network.ipynb | 529 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 529 insertions(+)
 create mode 100644 Multilayer_Network.ipynb

diff --git a/Multilayer_Network.ipynb b/Multilayer_Network.ipynb
new file mode 100644
index 0000000..46d54f9
--- /dev/null
+++ b/Multilayer_Network.ipynb
@@ -0,0 +1,529 @@
+# =============================================================================
+# Lattice Boolean Networks — Bit-Packed States + Memoization
+# -----------------------------------------------------------------------------
+# - rows×cols lattice of identical Boolean networks (CellSpec)
+# - Neighbor directions: up/down/left/right
+# - Global state is an integer with N bits (N = rows*cols*P). Bit i stores node i.
+# - Two-phase coupling:
+#       Phase A: identity inputs of each cell are set from neighbor outputs (prev state)
+#       Phase B: truth-table nodes update synchronously using Phase A values
+#
+# Measurements (Borriello & Daniels style):
+#   - Attractors (fixed points and cycles)
+#   - Control kernels (CKs) via static pinning:
+#       CK = input nodes + minimal distinguishing (witness) nodes + small additional
+#
+# Performance:
+#   - step_int(state:int)->int is memoized with LRU cache
+#   - pinned transitions are separately memoized per pin-set
+#
+# Dependencies: numpy, matplotlib, pandas
+# =============================================================================
+
+import itertools
+import math
+import random
+from collections import Counter
+from functools import lru_cache
+from typing import Dict, List, Tuple, Optional, Set, Callable
+
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+
+
+# --------------------------- bit helpers --------------------------------
+
+def get_bit(x: int, i: int) -> int:
+    """Return bit i (0/1) from integer x where bit 0 is LSB."""
+    return (x >> i) & 1
+
+def set_bit(x: int, i: int, v: int) -> int:
+    """Set bit i in integer x to v (0/1)."""
+    mask = 1 << i
+    return (x | mask) if v else (x & ~mask)
+
+def ints_to_bits_list(x: int, idxs: List[int]) -> Tuple[int, ...]:
+    """Extract bits at positions in idxs into a tuple (kept in given order)."""
+    return tuple(get_bit(x, i) for i in idxs)
+
+
+# --------------------------- homogeneous cell spec ----------------------
+
+class CellSpec:
+    """
+    One Boolean micro-network to be cloned at every lattice site.
+
+    Reserved, distinct local indices:
+      - out_idx: output node (feeds neighbors)
+      - in_up, in_down, in_left, in_right: identity input nodes (set from neighbors, not computed)
+
+    All other nodes are random K-input Boolean functions with a bias in {0,1}.
+    """
+    def __init__(self, P: int, out_idx: int,
+                 in_up: int, in_down: int, in_left: int, in_right: int,
+                 K: int = 2, bias: float = 0.5, seed: Optional[int] = None):
+        assert P >= 5, "P must be >= 5 (out + 4 identity inputs)."
+        assert len({out_idx, in_up, in_down, in_left, in_right}) == 5, "Reserved indices must be distinct."
+        self.P = P
+        self.out_idx = out_idx
+        self.in_ports = {"up": in_up, "down": in_down, "left": in_left, "right": in_right}
+        self.K = K
+        self.bias = bias
+        self._rng = random.Random(seed)
+        self.specs = self._build_specs()
+
+    def _build_specs(self):
+        """Build per-node update specs (identity vs truth table)."""
+        specs = []
+        for i in range(self.P):
+            if i in self.in_ports.values():
+                specs.append({"type": "identity"})  # overwritten by neighbor outputs each step
+            else:
+                # choose K distinct inputs (from all but self i)
+                choices = [j for j in range(self.P) if j != i]
+                inputs = self._rng.sample(choices, k=min(self.K, len(choices)))
+                table = {}
+                for pat in itertools.product([0,1], repeat=len(inputs)):
+                    table[pat] = 1 if self._rng.random() < self.bias else 0
+                specs.append({"type": "truth", "inputs": inputs, "table": table})
+        return specs
+
+    def identity_indices(self) -> Set[int]:
+        return set(self.in_ports.values())
+
+
+# --------------------------- lattice network (bit-packed) ---------------
+
+class LatticeBooleanNetworkInt:
+    """
+    A rows×cols lattice; each cell uses the same CellSpec. Global state is an int of N bits.
+    Bit layout: node index i corresponds to bit i (0-based, LSB=0).
+    Flatten index for local (r,c,i_local) is: i_global = ((r*cols + c) * P) + i_local.
+    """
+    def __init__(self, rows: int, cols: int, cell: CellSpec, wrap: bool = True, seed: Optional[int] = None):
+        self.rows, self.cols, self.cell = rows, cols, cell
+        self.P = cell.P
+        self.N = rows * cols * self.P
+        self.wrap = wrap
+        self._rng = random.Random(seed)
+
+        # Memoized core transition:
+        self._memo_step = lru_cache(maxsize=200_000)(self._step_impl)
+        # Per-pin-set memoized transitions:
+        self._pinned_steps: Dict[Tuple[Tuple[int,int], ...], Callable[[int], int]] = {}
+
+    # -------------- indexing helpers --------------
+
+    def idx(self, r: int, c: int, i_local: int) -> int:
+        return ((r * self.cols) + c) * self.P + i_local
+
+    def neighbors(self, r: int, c: int) -> Dict[str, Tuple[int,int]]:
+        def mod(x, m): return (x % m) if self.wrap else min(max(x, 0), m-1)
+        return {
+            "up":    (mod(r-1, self.rows), c),
+            "down":  (mod(r+1, self.rows), c),
+            "left":  (r, mod(c-1, self.cols)),
+            "right": (r, mod(c+1, self.cols)),
+        }
+
+    def input_nodes_global(self) -> Set[int]:
+        s = set()
+        for r in range(self.rows):
+            for c in range(self.cols):
+                for i_local in self.cell.identity_indices():
+                    s.add(self.idx(r, c, i_local))
+        return s
+
+    # -------------- core transition (memoized) --------------
+
+    def step(self, x: int) -> int:
+        """Public transition with memoization."""
+        return self._memo_step(x)
+
+    def _step_impl(self, x: int) -> int:
+        """
+        Apply one full two-phase step to integer state x.
+        Phase A: set identity inputs for each cell from neighbor outputs (from x).
+        Phase B: compute truth nodes using the Phase A cell-local bits; identity inputs remain as set.
+        """
+        P = self.P
+        out_idx = self.cell.out_idx
+        in_ports = self.cell.in_ports
+        specs = self.cell.specs
+
+        # Phase A: for each cell, build a list of P bits with identity ports replaced by neighbor outputs.
+        phaseA_cells: Dict[Tuple[int,int], List[int]] = {}
+        for r in range(self.rows):
+            for c in range(self.cols):
+                # read current cell bits
+                base = ((r * self.cols) + c) * P
+                cell_bits = [get_bit(x, base + i) for i in range(P)]
+                # fetch neighbor outputs
+                nb = self.neighbors(r, c)
+                up_out    = get_bit(x, self.idx(*nb["up"],    out_idx))
+                down_out  = get_bit(x, self.idx(*nb["down"],  out_idx))
+                left_out  = get_bit(x, self.idx(*nb["left"],  out_idx))
+                right_out = get_bit(x, self.idx(*nb["right"], out_idx))
+                # set identity input bits
+                cell_bits[in_ports["up"]]    = up_out
+                cell_bits[in_ports["down"]]  = down_out
+                cell_bits[in_ports["left"]]  = left_out
+                cell_bits[in_ports["right"]] = right_out
+                phaseA_cells[(r,c)] = cell_bits
+
+        # Phase B: compute next bits
+        y = 0
+        for r in range(self.rows):
+            for c in range(self.cols):
+                st = phaseA_cells[(r,c)]
+                new_bits = st[:]  # identity inputs kept as set
+                for i_local, spec in enumerate(specs):
+                    if spec["type"] == "truth":
+                        inp = tuple(st[j] for j in spec["inputs"])
+                        new_bits[i_local] = spec["table"][inp]
+                # write bits into result integer
+                base = ((r * self.cols) + c) * P
+                for i_local, v in enumerate(new_bits):
+                    if v:
+                        y |= (1 << (base + i_local))
+        return y
+
+    def clear_caches(self):
+        self._memo_step.cache_clear()
+        self._pinned_steps.clear()
+
+    # -------------- pinned transitions (memoized per pin-set) --------------
+
+    def pinned_step(self, pins: Dict[int,int]) -> Callable[[int], int]:
+        """
+        Return a memoized transition function f(x) for this pin-set.
+        'pins' maps global indices -> fixed values {0,1}.
+        """
+        key = tuple(sorted(pins.items()))
+        if key in self._pinned_steps:
+            return self._pinned_steps[key]
+
+        @lru_cache(maxsize=100_000)
+        def f(x: int) -> int:
+            # enforce pins before
+            for i, v in pins.items():
+                x = set_bit(x, i, v)
+            # transition
+            y = self._memo_step(x)
+            # enforce pins after (static pinning)
+            for i, v in pins.items():
+                y = set_bit(y, i, v)
+            return y
+
+        self._pinned_steps[key] = f
+        return f
+
+    # -------------- attractors (on ints) --------------
+
+    def trajectory(self, x0: int, max_steps: int = 100000) -> Tuple[List[int], Tuple[int,int]]:
+        """
+        Run from x0 until a repeat: returns (visited_list, (mu, lam)).
+        mu = preperiod length, lam = cycle length.
+        """
+        seen = {}
+        path: List[int] = []
+        x = x0
+        for t in range(max_steps):
+            if x in seen:
+                mu = seen[x]
+                lam = t - mu
+                return path, (mu, lam)
+            seen[x] = t
+            path.append(x)
+            x = self.step(x)
+        return path, (len(path), 0)
+
+    def enumerate_attractors(self, exhaustive: bool = False, sample_size: int = 1024, seed: Optional[int] = None):
+        """
+        Returns dict rep -> {"type","cycle_len","states","basin_count"}.
+        Canonical rep for a cycle = minimal integer in the cycle.
+        """
+        rng = random.Random(seed)
+        reps: Dict[int, Dict] = {}
+        basin = Counter()
+
+        def canonical_cycle(cyc: List[int]) -> Tuple[int, str]:
+            if len(cyc) == 1:
+                return cyc[0], "fixed"
+            return min(cyc), "cycle"
+
+        if exhaustive:
+            total = 1 << self.N
+            for x in range(total):
+                visited, (mu, lam) = self.trajectory(x)
+                cyc = visited[mu:mu+lam] if lam > 0 else [visited[-1]]
+                rep, kind = canonical_cycle(cyc)
+                basin[rep] += 1
+                if rep not in reps:
+                    reps[rep] = {"type": kind, "cycle_len": len(cyc), "states": tuple(cyc)}
+        else:
+            for _ in range(sample_size):
+                x0 = rng.getrandbits(self.N)
+                visited, (mu, lam) = self.trajectory(x0)
+                cyc = visited[mu:mu+lam] if lam > 0 else [visited[-1]]
+                rep, kind = canonical_cycle(cyc)
+                basin[rep] += 1
+                if rep not in reps:
+                    reps[rep] = {"type": kind, "cycle_len": len(cyc), "states": tuple(cyc)}
+        for rep in reps:
+            reps[rep]["basin_count"] = basin[rep]
+        return reps
+
+    # -------------- CK helpers --------------
+
+    def attractor_signature(self, cyc_states: Tuple[int, ...]) -> Tuple[int, ...]:
+        """
+        Signature: per-bit 0/1 if constant across the cycle; -1 if varies.
+        """
+        sig = []
+        for i in range(self.N):
+            b0 = get_bit(cyc_states[0], i)
+            const = all(get_bit(s, i) == b0 for s in cyc_states)
+            sig.append(b0 if const else -1)
+        return tuple(sig)
+
+    def minimal_witness(self, sig_target: Tuple[int, ...], sig_others: List[Tuple[int, ...]], max_k: int = 6) -> Tuple[int, ...]:
+        """
+        Find a small witness set S of indices (const bits in target) such that
+        for every other signature o, o|_S != target|_S. Brute force up to max_k.
+        """
+        eligible = [i for i, v in enumerate(sig_target) if v in (0,1)]
+        for k in range(min(max_k, len(eligible)) + 1):
+            for S in itertools.combinations(eligible, k):
+                ok = True
+                for o in sig_others:
+                    if all(o[j] == sig_target[j] for j in S):
+                        ok = False
+                        break
+                if ok:
+                    return S
+        return ()
+
+    def converges_to_single_attractor(self, pins: Dict[int,int],
+                                      exhaustive: bool = False, sample_size: int = 1024,
+                                      seed: Optional[int] = None) -> bool:
+        """
+        Under static pins, do all (or sampled) initial states converge to a single attractor?
+        """
+        rng = random.Random(seed)
+        step_fn = self.pinned_step(pins)
+        reps = set()
+
+        def traj(x0: int, T: int = 100000) -> int:
+            seen = {}
+            path: List[int] = []
+            x = x0
+            for t in range(T):
+                if x in seen:
+                    mu = seen[x]
+                    lam = t - mu
+                    cyc = path[mu:mu+lam] if lam > 0 else [path[-1]]
+                    return min(cyc)  # canonical rep = minimal int in cycle
+                seen[x] = t
+                path.append(x)
+                x = step_fn(x)
+            return path[-1]
+
+        if exhaustive:
+            for x in range(1 << self.N):
+                reps.add(traj(x))
+                if len(reps) > 1:
+                    return False
+            return True
+        else:
+            for _ in range(sample_size):
+                x0 = rng.getrandbits(self.N)
+                reps.add(traj(x0))
+                if len(reps) > 1:
+                    return False
+            return True
+
+    def control_kernel_for(self, cyc_states: Tuple[int, ...], all_sigs: List[Tuple[int, ...]],
+                           exhaustive: bool = False, sample_size: int = 2048,
+                           seed: Optional[int] = None, max_witness_k: int = 6,
+                           max_additional_rounds: int = 3) -> Dict:
+        """
+        CK = input nodes + minimal witness + small additional nodes if needed.
+        Returns dict with sets and sizes.
+        """
+        rng = random.Random(seed)
+        sig_t = self.attractor_signature(cyc_states)
+        sig_others = [s for s in all_sigs if s is not sig_t]
+
+        # If no constant bits in the target cycle, static pinning cannot stabilize globally.
+        if all(v == -1 for v in sig_t):
+            return {"controllable": False, "reason": "cycle_all_bits_vary"}
+
+        # Round 0: pin all input nodes to their target values
+        input_nodes = self.input_nodes_global()
+        rep_state = cyc_states[0]
+        pins = {i: get_bit(rep_state, i) for i in input_nodes}
+
+        # Round 1: minimal witness (distinguishing set)
+        witness = self.minimal_witness(sig_t, sig_others, max_k=max_witness_k)
+        for j in witness:
+            pins[j] = sig_t[j]
+
+        if self.converges_to_single_attractor(pins, exhaustive, sample_size, seed):
+            return {
+                "controllable": True,
+                "input_nodes": tuple(sorted(input_nodes)),
+                "witness_nodes": tuple(sorted(witness)),
+                "additional_nodes": tuple(),
+                "CK_nodes": tuple(sorted(pins.keys()))
+            }
+
+        # Additional small rounds (greedy, 1-bit increments then pairs)
+        additional: List[int] = []
+        current = dict(pins)
+        for _ in range(max_additional_rounds):
+            eligible = [j for j, v in enumerate(sig_t) if v in (0,1) and j not in current]
+            found = False
+            # try single
+            for j in eligible:
+                trial = dict(current); trial[j] = sig_t[j]
+                if self.converges_to_single_attractor(trial, exhaustive, sample_size, seed):
+                    current = trial; additional.append(j); found = True
+                    break
+            if found:
+                # one more check; if good, stop
+                if self.converges_to_single_attractor(current, exhaustive, sample_size, seed):
+                    break
+            else:
+                # try pairs if still not found
+                for j1, j2 in itertools.combinations(eligible, 2):
+                    trial = dict(current); trial[j1] = sig_t[j1]; trial[j2] = sig_t[j2]
+                    if self.converges_to_single_attractor(trial, exhaustive, sample_size, seed):
+                        current = trial; additional.extend([j1, j2]); found = True
+                        break
+                if not found:
+                    break
+
+        success = self.converges_to_single_attractor(current, exhaustive, sample_size, seed)
+        return {
+            "controllable": success,
+            "input_nodes": tuple(sorted(input_nodes)),
+            "witness_nodes": tuple(sorted(witness)),
+            "additional_nodes": tuple(sorted(additional)),
+            "CK_nodes": tuple(sorted(current.keys())) if success else tuple(sorted(pins.keys()))
+        }
+
+
+# --------------------------- experiment runner -------------------------
+
+def run_experiment_bitpacked(rows=2, cols=2, P=5, K=2, bias=0.5, seed=1,
+                             exhaustive: bool = False, sample_size: int = 512,
+                             max_witness_k: int = 6, max_additional_rounds: int = 3,
+                             n_networks: int = 1):
+    """
+    Build multiple random lattices and measure:
+      r, r_fixed, r_cycle, ⟨|CK|⟩, ⟨|w(1)|⟩, ⟨additional⟩, union CK size, log2 r, etc.
+    Returns (df_network, df_attractor).
+    """
+    rng = random.Random(seed)
+    net_rows, attr_rows = [], []
+
+    for nid in range(n_networks):
+        # Reserve local indices: out=0, in_up=1, in_down=2, in_left=3, in_right=4
+        out_idx, in_up, in_down, in_left, in_right = 0, 1, 2, 3, 4
+        cell = CellSpec(P=P, out_idx=out_idx, in_up=in_up, in_down=in_down,
+                        in_left=in_left, in_right=in_right, K=K, bias=bias,
+                        seed=rng.randint(0, 1 << 30))
+        net = LatticeBooleanNetworkInt(rows, cols, cell, wrap=True, seed=rng.randint(0, 1 << 30))
+
+        # Attractors
+        attractors = net.enumerate_attractors(exhaustive=exhaustive, sample_size=sample_size,
+                                              seed=rng.randint(0, 1 << 30))
+        reps = list(attractors.keys())
+        r = len(reps)
+        r_fixed = sum(1 for rep in reps if attractors[rep]["type"] == "fixed")
+        r_cycle = r - r_fixed
+
+        # Signatures for distinguishing sets
+        sigs = [net.attractor_signature(attractors[rep]["states"]) for rep in reps]
+
+        # CK per attractor
+        ck_sizes, w1_sizes, add_sizes = [], [], []
+        union_nodes: Set[int] = set()
+
+        for rep, sig in zip(reps, sigs):
+            info = net.control_kernel_for(attractors[rep]["states"], sigs,
+                                          exhaustive=exhaustive, sample_size=max(sample_size, 1024),
+                                          seed=rng.randint(0, 1 << 30),
+                                          max_witness_k=max_witness_k,
+                                          max_additional_rounds=max_additional_rounds)
+            if info["controllable"]:
+                ck_sizes.append(len(info["CK_nodes"]))
+                w1_sizes.append(len(info["witness_nodes"]))
+                add_sizes.append(len(info["additional_nodes"]))
+                union_nodes.update(info["CK_nodes"])
+
+            attr_rows.append({
+                "network_id": nid,
+                "type": attractors[rep]["type"],
+                "cycle_len": attractors[rep]["cycle_len"],
+                "basin_count": attractors[rep]["basin_count"],
+                "CK_size": len(info["CK_nodes"]) if info["controllable"] else np.nan,
+                "w1_size": len(info["witness_nodes"]) if info["controllable"] else np.nan,
+                "add_size": len(info["additional_nodes"]) if info["controllable"] else np.nan,
+                "r": r
+            })
+
+        mean_ck = float(np.nanmean(ck_sizes)) if ck_sizes else float('nan')
+        mean_w1 = float(np.nanmean(w1_sizes)) if w1_sizes else float('nan')
+        mean_add = float(np.nanmean(add_sizes)) if add_sizes else float('nan')
+        net_rows.append({
+            "network_id": nid,
+            "rows": rows, "cols": cols, "P": P, "K": K, "bias": bias,
+            "N_total": net.N,
+            "r": r, "r_fixed": r_fixed, "r_cycle": r_cycle,
+            "mean_CK": mean_ck, "mean_w1": mean_w1, "mean_add": mean_add,
+            "union_CK": len(union_nodes),
+            "log2_r": math.log2(r) if r > 0 else float('nan'),
+        })
+
+    return pd.DataFrame(net_rows), pd.DataFrame(attr_rows)
+
+
+# --------------------------- example run -------------------------
+
+if __name__ == "__main__":
+    # Small defaults so it runs fast; scale sample_size or n_networks for stronger stats.
+    df_net, df_attr = run_experiment_bitpacked(
+        rows=3, cols=3, P=5, K=2, bias=0.5, seed=123,
+        exhaustive=False, sample_size=256,
+        max_witness_k=10, max_additional_rounds=5,
+        n_networks=10
+    )
+
+    print("Network summary:\n", df_net)
+    print("\nAttractor details (head):\n", df_attr.head())
+
+    # Plot: mean CK size vs log2(r)
+    plt.figure()
+    plt.scatter(df_net["log2_r"], df_net["mean_CK"])
+    for _, row in df_net.iterrows():
+        if not np.isnan(row["mean_CK"]):
+            plt.text(row["log2_r"], row["mean_CK"], f'n{int(row["network_id"])}')
+    plt.xlabel("log2(# attractors)")
+    plt.ylabel("Mean CK size")
+    plt.title("Mean Control Kernel Size vs log2(r)")
+    plt.show()
+
+    # Plot: witness vs additional contributions
+    plt.figure()
+    x = np.arange(len(df_net))
+    w = 0.35
+    plt.bar(x - w/2, df_net["mean_w1"], width=w, label="Mean |w(1)|")
+    plt.bar(x + w/2, df_net["mean_add"], width=w, label="Mean additional")
+    plt.xticks(x, [f'n{int(i)}' for i in df_net["network_id"]])
+    plt.ylabel("Size")
+    plt.title("CK contributions: witness vs additional")
+    plt.legend()
+    plt.show()