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 # (set after producers are known: total administrative budget = ratio x expected relayed messages per period)146 use_code = regime in ("norms", "records", "central", "federated")147 use_records = regime in ("records", "central", "federated")148 hier = regime in ("central", "federated")149 depth, paths = tree_structure(P, pos) if hier else (0, None)150 n_admin = 0151 if hier:152 n_admin = sum(P.branching ** lvl for lvl in range(depth)) # one administrator per internal node153 producers = np.ones(P.N, bool)154 if hier: # administrators are withdrawn from production (overhead): the first n_admin agents by index155 producers[: min(n_admin, P.N // 4)] = False156 # administrative budget (A-005): total = ratio x expected relayed load of PRODUCING agents; central: the root holds the157 # whole budget (only it relays); federated: the budget is split equally among all administrators158 # Allocation rule (A-005): the SAME total budget for central and federated (ratio x expected relayed load of a159 # centralized system, 2*depth*lam*N_prod*(k-1) messages per period); each administrator receives the share of the budget equal to its share of160 # the load measured in a 5-period warm-up with unconstrained capacity ("capacity follows load"). Thus the ratio is the161 # system-wide capacity/load of the CENTRAL design, and federation's advantage can only come from relaying less.162 admin_share = {}; per_admin_default = 1.0; total_budget = 0.0163 if hier:164 # expected relayed messages per period under the CENTRAL design: 2*depth hops per recipient165 total_budget = P.admin_capacity_ratio * 2 * depth * P.lam * producers.sum() * (P.k - 1) if P.admin_capacity_ratio > 0 else P.admin_bandwidth * n_admin166 wrng = np.random.default_rng(P.seed + 10_000); loads = {}167 for _ in range(5):168 Mw = wrng.poisson(P.lam * producers.sum()); iw = wrng.choice(np.where(producers)[0], size=max(Mw, 1))169 for m in range(len(iw)):170 nb = np.where(A[iw[m]] & producers)[0]171 pw_ = wrng.choice(nb, size=P.k - 1, replace=len(nb) < P.k - 1) if len(nb) > 0 else wrng.choice(np.where(producers)[0], size=P.k - 1)172 for part in pw_:173 l0 = int(lca_depth(paths, np.array([iw[m]]), np.array([part]))[0]) if regime == "federated" else 0174 for who in (iw[m], part):175 for level in range(l0, depth):176 key = (level, tuple(paths[who][:level])); loads[key] = loads.get(key, 0.0) + 1.0177 tot = sum(loads.values()); admin_share = {k_: v / tot for k_, v in loads.items()}178 def admin_bw(key):179 return max(1.0, total_budget * admin_share.get(key, 0.0)) if hier else 1.0180 # per-agent FIFO queues: cumulative message count served at `bandwidth` per period181 agent_backlog = np.zeros(P.N); admin_backlog = {} # admin id -> backlog182 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,183 tasks_pre=0, success_pre=0, tasks_post=0, success_post=0, bits_pre=0.0, bits_post=0.0, root_load=0.0)184 benefit = np.zeros(P.N)185 learn_cost = P.N * P.learn_bits if use_code else 0.0186 cur_code_len = code_len.copy(); shock_time = None187 for t in range(P.T):188 if P.shock_period == t and P.shock_type == "network":189 drop = np.triu(rng.random(A.shape) < P.shock_link_removal, 1); drop = drop | drop.T; A &= ~drop # symmetric removal of a share of links190 if P.shock_period == t and P.shock_type == "distribution":191 shock_time = t192 dist_now = q if (shock_time is not None and t >= shock_time) else p193 if use_code and shock_time is not None and t >= shock_time + P.adapt_delay:194 cur_code_len = shannon_lengths(q) # institution re-learns the code after adapt_delay195 # drain queues by one period of service196 agent_backlog = np.maximum(0, agent_backlog - P.bandwidth)197 for a_ in admin_backlog:198 admin_backlog[a_] = max(0.0, admin_backlog[a_] - admin_bw(a_))199 M = rng.poisson(P.lam * producers.sum())200 if M == 0:201 continue202 init = rng.choice(np.where(producers)[0], size=M)203 plans = rng.choice(P.s, size=M, p=dist_now)204 # participants: k-1 network neighbours of the initiator (random if none)205 parts = np.empty((M, P.k - 1), int)206 for m in range(M):207 nb = np.where(A[init[m]] & producers)[0]208 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)209 # bits per recipient210 if use_code:211 bits_r = cur_code_len[plans]212 else:213 bits_r = np.full(M, float(L_raw))214 # decoding errors215 if use_records and P.verify_records:216 err = np.zeros((M, P.k - 1), bool)217 elif use_code:218 err = rng.random((M, P.k - 1)) < P.eps_code219 else:220 dd = np.abs(dialect[init][:, None] - dialect[parts]); err = rng.random((M, P.k - 1)) < P.eps * dd221 task_fail_decode = err.any(axis=1)222 # messages and routing223 msgs_per_recip = 1 if use_records else 2 # ack round without records224 if not hier:225 n_msgs = msgs_per_recip * (P.k - 1)226 # load: each message is processed by its recipient; the ack by the initiator227 np.add.at(agent_backlog, parts.ravel(), 1.0)228 if msgs_per_recip == 2:229 np.add.at(agent_backlog, init, float(P.k - 1))230 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)231 root_load_t = 0.0232 else:233 # relay through the LCA of initiator and each participant234 lvl = np.stack([lca_depth(paths, init, parts[:, c]) for c in range(P.k - 1)], axis=1) # 0 = root235 hops = 2 * (depth - lvl) # messages up and down236 n_msgs = hops.sum(axis=1)237 if regime == "federated":238 pass # locality is captured by the LCA level: same-module tasks never reach the root239 else:240 lvl[:] = 0 # central: everything is relayed through the root241 hops = np.full_like(hops, 2 * depth); n_msgs = hops.sum(axis=1)242 # admin load (A-005): every administrator on the path leaf -> LCA -> leaf processes one message per hop243 wait = np.zeros(M); hop_bits = np.zeros(M)244 for c in range(P.k - 1):245 for m in range(M):246 l0 = int(lvl[m, c]); keys_m = []247 for who in (init[m], parts[m, c]):248 for level in range(l0, depth): # administrators at levels l0..depth-1 on this leaf's path249 keys_m.append((level, tuple(paths[who][:level])))250 for key in keys_m:251 admin_backlog[key] = admin_backlog.get(key, 0.0) + 1.0252 for key in keys_m:253 wait[m] = max(wait[m], math.ceil(admin_backlog[key] / admin_bw(key)))254 hop_bits[m] += hops[m, c] # bits are transmitted on every hop255 np.add.at(agent_backlog, parts.ravel(), 1.0)256 wait = np.maximum(wait, np.ceil(agent_backlog[parts] / P.bandwidth).max(axis=1))257 root_load_t = float(admin_backlog.get((0, ()), 0.0))258 task_fail_latency = wait > P.tau_max259 success = ~task_fail_decode & ~task_fail_latency260 # bits per task = codeword/raw bits to each recipient + 1-bit acknowledgements (if any) + record overhead261 if hier:262 bits = bits_r * hop_bits + (P.admin_bits if use_records else 0.0) # codeword repeated on each relay hop263 else:264 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)265 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())266 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)267 stats["latency_sum"] += float(wait[success].sum()); stats["root_load"] += root_load_t268 if t < P.split_period:269 stats["tasks_pre"] += M; stats["success_pre"] += int(success.sum()); stats["bits_pre"] += float(bits.sum())270 else:271 stats["tasks_post"] += M; stats["success_post"] += int(success.sum()); stats["bits_post"] += float(bits.sum())272 np.add.at(benefit, init[success], P.value_per_task); np.add.at(benefit, parts[success].ravel(), P.value_per_task)273 tot_bits = stats["bits"] + learn_cost274 succ = max(stats["success"], 1)275 b = benefit[producers] if producers.any() else benefit276 g = float(np.abs(b[:, None] - b[None, :]).mean() / (2 * b.mean())) if b.mean() > 0 else float("nan")277 out = {**asdict(P), "admin_budget_total": (total_budget if hier else 0.0), "n_producers": int(producers.sum()), "H_p_bits": H, "L_code_bits": L_code, "L_raw_bits": float(L_raw), "n_admin": int((~producers).sum()), "depth": depth,278 "success_rate": stats["success"] / max(stats["tasks"], 1), "msgs_per_task": stats["messages"] / max(stats["tasks"], 1),279 "bits_per_task": stats["bits"] / max(stats["tasks"], 1), "bits_per_success_incl_setup": tot_bits / succ, "bits_per_success": stats["bits"] / succ,280 "tasks_per_kbit": succ / (tot_bits / 1000.0), "latency_mean": stats["latency_sum"] / succ, "fail_latency_share": stats["fail_latency"] / max(stats["tasks"], 1),281 "fail_decode_share": stats["fail_decode"] / max(stats["tasks"], 1), "admin_overhead_share": (stats["admin_bits"] + learn_cost) / max(tot_bits, 1e-9),282 "welfare": float(benefit.sum() - P.cost_per_bit * tot_bits), "welfare_per_agent": float((benefit.sum() - P.cost_per_bit * tot_bits) / P.N),283 "gini_benefit": g, "root_load_mean": stats["root_load"] / P.T,284 "success_pre": stats["success_pre"] / max(stats["tasks_pre"], 1), "success_post": stats["success_post"] / max(stats["tasks_post"], 1),285 "bits_per_task_pre": stats["bits_pre"] / max(stats["tasks_pre"], 1), "bits_per_task_post": stats["bits_post"] / max(stats["tasks_post"], 1),286 "kl_q_p_bits": float((q * np.log2(q / p)).sum()) if P.shock_type == "distribution" else 0.0, "tasks": stats["tasks"]}287 return out288