1"""SIMULATED — agent-based coordination model for H3 (pre-registration §1, framework Theorems 1–5).23Agents live on a spatial interaction network. Each period, tasks arrive; a task needs k participants to hold the same4plan x in S (|S| = s), x drawn from a Zipf-like distribution p. The initiator must communicate x to the k-1 others.5Communication is counted in MESSAGES and BITS; each agent (and each administrator) can process at most `bandwidth`6messages per period (limited attention) — excess messages queue (FIFO), producing latency; a task whose last message7is served more than `tau_max` periods after arrival FAILS (coordination failure). Without a shared code, decoding a8message across agents of different "dialects" fails with probability eps * |dialect_i - dialect_j| (heterogeneous9information / trust); with a shared code the error is eps_code; with records it is ~0 (verification against the record).1011Five institutional regimes are IMPOSED by the experimenter (nothing emerges):12 none : direct messages, log2(s) bits per recipient, dialect errors, one acknowledgement round (2 msgs/recipient)13 norms : shared prefix code optimal for p (bits = codeword length), lower error; learning cost per agent; after a14 distribution shock the stale code costs KL(q||p) extra bits until re-learned (adapt_delay periods)15 records : norms + a written record per task (admin overhead bits) that recipients verify (error -> 0, 1 msg each)16 central : hierarchy (b-ary tree of administrators); messages relayed via the lowest common ancestor; admin nodes17 have their own bandwidth; per-task admin overhead; administrators are agents withdrawn from production18 federated : hierarchy with local modules: tasks whose participants share a module are settled by the module admin,19 cross-module tasks go to the root; same overhead as central20All quantities below are measured, not assumed. Shannon quantities (H(p), code lengths, KL) are computed only because21p and the code are explicitly defined in the simulation.22"""23from __future__ import annotations24from dataclasses import dataclass, asdict25import math26import numpy as np272829@dataclass30class Params:31 N: int = 20032 k: int = 3 # participants per task33 lam: float = 0.5 # tasks per agent per period34 T: int = 40 # periods35 s: int = 64 # plan set size36 zipf: float = 1.0 # 0 -> uniform plans (H = log2 s); larger -> concentrated37 network: str = "lattice" # lattice | smallworld | random38 rewire: float = 0.139 degree: int = 640 bandwidth: int = 8 # messages an agent can process per period41 admin_bandwidth: int = 60 # messages an administrator can process per period (used when admin_capacity_ratio <= 0)42 admin_capacity_ratio: float = 0.0 # if > 0: admin_bandwidth = ratio * expected relayed messages per period (2*lam*N*(k-1))43 verify_records: bool = True # ablation: records without verification (overhead only)44 tau_max: int = 3 # periods allowed before a task fails45 eps: float = 0.25 # per-unit-dialect-distance decoding error (no code)46 eps_code: float = 0.02 # decoding error with a shared code47 learn_bits: float = 200.0 # bits each agent spends to learn a code (fixed cost, once)48 admin_bits: float = 4.0 # per-task administrative overhead bits (records/central/federated)49 branching: int = 8 # tree branching factor for central/federated50 shock_period: int = -1 # period at which a shock occurs (-1: none)51 shock_type: str = "none" # none | distribution | network52 shock_kl_target: float = 1.0 # for distribution shocks: how far q is from p (approx bits)53 shock_link_removal: float = 0.554 split_period: int = 20 # pre/post accounting window boundary used in ALL runs (shocked or not), for difference-in-differences55 adapt_delay: int = 10 # periods before a code is re-optimized after a distribution shock56 cost_per_bit: float = 0.00257 value_per_task: float = 1.058 seed: int = 059 regime: str = "none"606162def zipf_dist(s: int, a: float) -> np.ndarray:63 w = 1.0 / np.arange(1, s + 1) ** a if a > 0 else np.ones(s)64 return w / w.sum()656667def entropy_bits(p: np.ndarray) -> float:68 p = p[p > 0]; return float(-(p * np.log2(p)).sum())697071def shannon_lengths(p: np.ndarray) -> np.ndarray:72 """Shannon code lengths ceil(-log2 p(x)) — a prefix code with these lengths exists by Kraft's inequality; this is the73 code used in Theorem 3 of research/framework.md (so that the measured stale-code penalty is comparable to D_KL)."""74 return np.ceil(-np.log2(np.maximum(p, 1e-300)))757677def huffman_lengths(p: np.ndarray) -> np.ndarray:78 """Optimal prefix-code lengths for distribution p (Huffman); kept for reference, not used by default."""79 import heapq80 n = len(p); heap = [(float(p[i]), i, [i]) for i in range(n)]; heapq.heapify(heap); L = np.zeros(n)81 cnt = n82 while len(heap) > 1:83 a = heapq.heappop(heap); b = heapq.heappop(heap)84 for i in a[2] + b[2]:85 L[i] += 186 heapq.heappush(heap, (a[0] + b[0], cnt, a[2] + b[2])); cnt += 187 return L888990def build_network(P: Params, rng: np.random.Generator):91 pos = rng.random((P.N, 2))92 d = np.linalg.norm(pos[:, None, :] - pos[None, :, :], axis=2); np.fill_diagonal(d, np.inf)93 nbr = np.argsort(d, axis=1)[:, :P.degree]94 A = np.zeros((P.N, P.N), bool)95 for i in range(P.N):96 A[i, nbr[i]] = True; A[nbr[i], i] = True97 if P.network == "smallworld":98 for i in range(P.N):99 for j in np.where(A[i])[0]:100 if j > i and rng.random() < P.rewire:101 A[i, j] = A[j, i] = False; m = rng.integers(P.N)102 if m != i:103 A[i, m] = A[m, i] = True104 elif P.network == "random":105 A[:] = False; m = P.N * P.degree // 2106 ii = rng.integers(P.N, size=m); jj = rng.integers(P.N, size=m); ok = ii != jj107 A[ii[ok], jj[ok]] = True; A[jj[ok], ii[ok]] = True108 return pos, A109110111def tree_structure(P: Params, pos: np.ndarray):112 """Assign agents to leaves of a b-ary tree by spatial sorting (modules are spatially contiguous)."""113 b = P.branching114 depth = max(1, math.ceil(math.log(P.N, b)))115 order = np.lexsort((pos[:, 1], np.floor(pos[:, 0] * b))) # column strips then y116 leaf = np.empty(P.N, int); leaf[order] = np.arange(P.N)117 # path of each agent: digits of its leaf index in base b, depth levels118 paths = np.array([[(l // b ** (depth - 1 - lvl)) % b for lvl in range(depth)] for l in leaf])119 return depth, paths120121122def lca_depth(paths: np.ndarray, i: np.ndarray, j: np.ndarray) -> np.ndarray:123 """Level of the lowest common ancestor (0 = root). Vectorized over task pairs."""124 eq = paths[i] == paths[j]125 first_diff = np.where(eq.all(axis=1), eq.shape[1], np.argmin(eq, axis=1))126 return first_diff127128129def simulate(P: Params) -> dict:130 rng = np.random.default_rng(P.seed)131 pos, A = build_network(P, rng)132 dialect = 0.5 * pos[:, 0] + 0.5 * rng.random(P.N) # heterogeneous information: cultural distance, partly spatial (geography)133 p = zipf_dist(P.s, P.zipf); H = entropy_bits(p)134 code_len = shannon_lengths(p); L_code = float((p * code_len).sum()); L_raw = math.ceil(math.log2(P.s))135 q = p.copy()136 if P.shock_type == "distribution" and P.shock_period >= 0:137 # q = p permuted so that D_KL(q||p) is close to the target (search over cyclic shifts)138 best = None139 for sh in range(1, P.s):140 qq = np.roll(p, sh); kl = float((qq * np.log2(qq / p)).sum())141 if best is None or abs(kl - P.shock_kl_target) < abs(best[0] - P.shock_kl_target):142 best = (kl, qq)143 q = best[1]144 regime = P.regime145 if P.admin_capacity_ratio > 0:146 P.admin_bandwidth = max(1, int(round(P.admin_capacity_ratio * 2 * P.lam * P.N * (P.k - 1))))147 use_code = regime in ("norms", "records", "central", "federated")148 use_records = regime in ("records", "central", "federated")149 hier = regime in ("central", "federated")150 depth, paths = tree_structure(P, pos) if hier else (0, None)151 n_admin = 0152 if hier:153 n_admin = sum(P.branching ** lvl for lvl in range(depth)) # one administrator per internal node154 producers = np.ones(P.N, bool)155 if hier: # administrators are withdrawn from production (overhead): the first n_admin agents by index156 producers[: min(n_admin, P.N // 4)] = False157 # per-agent FIFO queues: cumulative message count served at `bandwidth` per period158 agent_backlog = np.zeros(P.N); admin_backlog = {} # admin id -> backlog159 stats = dict(tasks=0, success=0, fail_latency=0, fail_decode=0, messages=0, bits=0.0, admin_bits=0.0, latency_sum=0.0,160 tasks_pre=0, success_pre=0, tasks_post=0, success_post=0, bits_pre=0.0, bits_post=0.0, root_load=0.0)161 benefit = np.zeros(P.N)162 learn_cost = P.N * P.learn_bits if use_code else 0.0163 cur_code_len = code_len.copy(); shock_time = None164 for t in range(P.T):165 if P.shock_period == t and P.shock_type == "network":166 drop = np.triu(rng.random(A.shape) < P.shock_link_removal, 1); drop = drop | drop.T; A &= ~drop # symmetric removal of a share of links167 if P.shock_period == t and P.shock_type == "distribution":168 shock_time = t169 dist_now = q if (shock_time is not None and t >= shock_time) else p170 if use_code and shock_time is not None and t >= shock_time + P.adapt_delay:171 cur_code_len = shannon_lengths(q) # institution re-learns the code after adapt_delay172 # drain queues by one period of service173 agent_backlog = np.maximum(0, agent_backlog - P.bandwidth)174 for a_ in admin_backlog:175 admin_backlog[a_] = max(0.0, admin_backlog[a_] - P.admin_bandwidth)176 M = rng.poisson(P.lam * producers.sum())177 if M == 0:178 continue179 init = rng.choice(np.where(producers)[0], size=M)180 plans = rng.choice(P.s, size=M, p=dist_now)181 # participants: k-1 network neighbours of the initiator (random if none)182 parts = np.empty((M, P.k - 1), int)183 for m in range(M):184 nb = np.where(A[init[m]] & producers)[0]185 parts[m] = rng.choice(nb, size=P.k - 1, replace=len(nb) < P.k - 1) if len(nb) > 0 else rng.choice(np.where(producers)[0], size=P.k - 1)186 # bits per recipient187 if use_code:188 bits_r = cur_code_len[plans]189 else:190 bits_r = np.full(M, float(L_raw))191 # decoding errors192 if use_records and P.verify_records:193 err = np.zeros((M, P.k - 1), bool)194 elif use_code:195 err = rng.random((M, P.k - 1)) < P.eps_code196 else:197 dd = np.abs(dialect[init][:, None] - dialect[parts]); err = rng.random((M, P.k - 1)) < P.eps * dd198 task_fail_decode = err.any(axis=1)199 # messages and routing200 msgs_per_recip = 1 if use_records else 2 # ack round without records201 if not hier:202 n_msgs = msgs_per_recip * (P.k - 1)203 # load: each message is processed by its recipient; the ack by the initiator204 np.add.at(agent_backlog, parts.ravel(), 1.0)205 if msgs_per_recip == 2:206 np.add.at(agent_backlog, init, float(P.k - 1))207 wait = np.maximum(np.ceil(agent_backlog[parts] / P.bandwidth).max(axis=1), np.ceil(agent_backlog[init] / P.bandwidth) if msgs_per_recip == 2 else 0)208 root_load_t = 0.0209 else:210 # relay through the LCA of initiator and each participant211 lvl = np.stack([lca_depth(paths, init, parts[:, c]) for c in range(P.k - 1)], axis=1) # 0 = root212 hops = 2 * (depth - lvl) # messages up and down213 n_msgs = hops.sum(axis=1)214 if regime == "federated":215 pass # locality is captured by the LCA level: same-module tasks never reach the root216 else:217 lvl[:] = 0 # central: everything is relayed through the root218 hops = np.full_like(hops, 2 * depth); n_msgs = hops.sum(axis=1)219 # admin load: the LCA administrator processes 2 messages per recipient220 keys = []221 for c in range(P.k - 1):222 for m in range(M):223 key = (int(lvl[m, c]), tuple(paths[init[m]][: lvl[m, c]])); keys.append(key)224 admin_backlog[key] = admin_backlog.get(key, 0.0) + 2.0225 wait = np.zeros(M)226 idx = 0227 for c in range(P.k - 1):228 for m in range(M):229 wait[m] = max(wait[m], math.ceil(admin_backlog[keys[idx]] / P.admin_bandwidth)); idx += 1230 np.add.at(agent_backlog, parts.ravel(), 1.0)231 wait = np.maximum(wait, np.ceil(agent_backlog[parts] / P.bandwidth).max(axis=1))232 root_load_t = float(admin_backlog.get((0, ()), 0.0))233 task_fail_latency = wait > P.tau_max234 success = ~task_fail_decode & ~task_fail_latency235 # bits per task = codeword/raw bits to each recipient + 1-bit acknowledgements (if any) + record overhead236 bits = bits_r * (P.k - 1) + (1.0 * (P.k - 1) if msgs_per_recip == 2 else 0.0) + (P.admin_bits if use_records else 0.0)237 stats["tasks"] += M; stats["success"] += int(success.sum()); stats["fail_decode"] += int(task_fail_decode.sum()); stats["fail_latency"] += int((task_fail_latency & ~task_fail_decode).sum())238 stats["messages"] += int(np.sum(n_msgs) if np.ndim(n_msgs) else M * n_msgs); stats["bits"] += float(bits.sum()); stats["admin_bits"] += float(M * P.admin_bits if use_records else 0.0)239 stats["latency_sum"] += float(wait[success].sum()); stats["root_load"] += root_load_t240 if t < P.split_period:241 stats["tasks_pre"] += M; stats["success_pre"] += int(success.sum()); stats["bits_pre"] += float(bits.sum())242 else:243 stats["tasks_post"] += M; stats["success_post"] += int(success.sum()); stats["bits_post"] += float(bits.sum())244 np.add.at(benefit, init[success], P.value_per_task); np.add.at(benefit, parts[success].ravel(), P.value_per_task)245 tot_bits = stats["bits"] + learn_cost246 succ = max(stats["success"], 1)247 b = benefit[producers] if producers.any() else benefit248 g = float(np.abs(b[:, None] - b[None, :]).mean() / (2 * b.mean())) if b.mean() > 0 else float("nan")249 out = {**asdict(P), "H_p_bits": H, "L_code_bits": L_code, "L_raw_bits": float(L_raw), "n_admin": int((~producers).sum()), "depth": depth,250 "success_rate": stats["success"] / max(stats["tasks"], 1), "msgs_per_task": stats["messages"] / max(stats["tasks"], 1),251 "bits_per_task": stats["bits"] / max(stats["tasks"], 1), "bits_per_success_incl_setup": tot_bits / succ, "bits_per_success": stats["bits"] / succ,252 "tasks_per_kbit": succ / (tot_bits / 1000.0), "latency_mean": stats["latency_sum"] / succ, "fail_latency_share": stats["fail_latency"] / max(stats["tasks"], 1),253 "fail_decode_share": stats["fail_decode"] / max(stats["tasks"], 1), "admin_overhead_share": (stats["admin_bits"] + learn_cost) / max(tot_bits, 1e-9),254 "welfare": float(benefit.sum() - P.cost_per_bit * tot_bits), "welfare_per_agent": float((benefit.sum() - P.cost_per_bit * tot_bits) / P.N),255 "gini_benefit": g, "root_load_mean": stats["root_load"] / P.T,256 "success_pre": stats["success_pre"] / max(stats["tasks_pre"], 1), "success_post": stats["success_post"] / max(stats["tasks_post"], 1),257 "bits_per_task_pre": stats["bits_pre"] / max(stats["tasks_pre"], 1), "bits_per_task_post": stats["bits_post"] / max(stats["tasks_post"], 1),258 "kl_q_p_bits": float((q * np.log2(q / p)).sum()) if P.shock_type == "distribution" else 0.0, "tasks": stats["tasks"]}259 return out260