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.354 adapt_delay: int = 10 # periods before a code is re-optimized after a distribution shock55 cost_per_bit: float = 0.00256 value_per_task: float = 1.057 seed: int = 058 regime: str = "none"596061def zipf_dist(s: int, a: float) -> np.ndarray:62 w = 1.0 / np.arange(1, s + 1) ** a if a > 0 else np.ones(s)63 return w / w.sum()646566def entropy_bits(p: np.ndarray) -> float:67 p = p[p > 0]; return float(-(p * np.log2(p)).sum())686970def huffman_lengths(p: np.ndarray) -> np.ndarray:71 """Optimal prefix-code lengths for distribution p (Huffman)."""72 import heapq73 n = len(p); heap = [(float(p[i]), i, [i]) for i in range(n)]; heapq.heapify(heap); L = np.zeros(n)74 cnt = n75 while len(heap) > 1:76 a = heapq.heappop(heap); b = heapq.heappop(heap)77 for i in a[2] + b[2]:78 L[i] += 179 heapq.heappush(heap, (a[0] + b[0], cnt, a[2] + b[2])); cnt += 180 return L818283def build_network(P: Params, rng: np.random.Generator):84 pos = rng.random((P.N, 2))85 d = np.linalg.norm(pos[:, None, :] - pos[None, :, :], axis=2); np.fill_diagonal(d, np.inf)86 nbr = np.argsort(d, axis=1)[:, :P.degree]87 A = np.zeros((P.N, P.N), bool)88 for i in range(P.N):89 A[i, nbr[i]] = True; A[nbr[i], i] = True90 if P.network == "smallworld":91 for i in range(P.N):92 for j in np.where(A[i])[0]:93 if j > i and rng.random() < P.rewire:94 A[i, j] = A[j, i] = False; m = rng.integers(P.N)95 if m != i:96 A[i, m] = A[m, i] = True97 elif P.network == "random":98 A[:] = False; m = P.N * P.degree // 299 ii = rng.integers(P.N, size=m); jj = rng.integers(P.N, size=m); ok = ii != jj100 A[ii[ok], jj[ok]] = True; A[jj[ok], ii[ok]] = True101 return pos, A102103104def tree_structure(P: Params, pos: np.ndarray):105 """Assign agents to leaves of a b-ary tree by spatial sorting (modules are spatially contiguous)."""106 b = P.branching107 depth = max(1, math.ceil(math.log(P.N, b)))108 order = np.lexsort((pos[:, 1], np.floor(pos[:, 0] * b))) # column strips then y109 leaf = np.empty(P.N, int); leaf[order] = np.arange(P.N)110 # path of each agent: digits of its leaf index in base b, depth levels111 paths = np.array([[(l // b ** (depth - 1 - lvl)) % b for lvl in range(depth)] for l in leaf])112 return depth, paths113114115def lca_depth(paths: np.ndarray, i: np.ndarray, j: np.ndarray) -> np.ndarray:116 """Level of the lowest common ancestor (0 = root). Vectorized over task pairs."""117 eq = paths[i] == paths[j]118 first_diff = np.where(eq.all(axis=1), eq.shape[1], np.argmin(eq, axis=1))119 return first_diff120121122def simulate(P: Params) -> dict:123 rng = np.random.default_rng(P.seed)124 pos, A = build_network(P, rng)125 dialect = rng.random(P.N) # heterogeneous information: cultural distance126 p = zipf_dist(P.s, P.zipf); H = entropy_bits(p)127 code_len = huffman_lengths(p); L_code = float((p * code_len).sum()); L_raw = math.ceil(math.log2(P.s))128 q = p.copy()129 if P.shock_type == "distribution" and P.shock_period >= 0:130 # q = p permuted so that D_KL(q||p) is close to the target (search over cyclic shifts)131 best = None132 for sh in range(1, P.s):133 qq = np.roll(p, sh); kl = float((qq * np.log2(qq / p)).sum())134 if best is None or abs(kl - P.shock_kl_target) < abs(best[0] - P.shock_kl_target):135 best = (kl, qq)136 q = best[1]137 regime = P.regime138 if P.admin_capacity_ratio > 0:139 P.admin_bandwidth = max(1, int(round(P.admin_capacity_ratio * 2 * P.lam * P.N * (P.k - 1))))140 use_code = regime in ("norms", "records", "central", "federated")141 use_records = regime in ("records", "central", "federated")142 hier = regime in ("central", "federated")143 depth, paths = tree_structure(P, pos) if hier else (0, None)144 n_admin = 0145 if hier:146 n_admin = sum(P.branching ** lvl for lvl in range(depth)) # one administrator per internal node147 producers = np.ones(P.N, bool)148 if hier: # administrators are withdrawn from production (overhead): the first n_admin agents by index149 producers[: min(n_admin, P.N // 4)] = False150 # per-agent FIFO queues: cumulative message count served at `bandwidth` per period151 agent_backlog = np.zeros(P.N); admin_backlog = {} # admin id -> backlog152 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,153 tasks_pre=0, success_pre=0, tasks_post=0, success_post=0, bits_pre=0.0, bits_post=0.0, root_load=0.0)154 benefit = np.zeros(P.N)155 learn_cost = P.N * P.learn_bits if use_code else 0.0156 cur_code_len = code_len.copy(); shock_time = None157 for t in range(P.T):158 if P.shock_period == t and P.shock_type == "network":159 drop = rng.random(A.shape) < P.shock_link_removal; A &= ~drop; A &= A.T160 if P.shock_period == t and P.shock_type == "distribution":161 shock_time = t162 dist_now = q if (shock_time is not None and t >= shock_time) else p163 if use_code and shock_time is not None and t >= shock_time + P.adapt_delay:164 cur_code_len = huffman_lengths(q) # institution re-learns the code after adapt_delay165 # drain queues by one period of service166 agent_backlog = np.maximum(0, agent_backlog - P.bandwidth)167 for a_ in admin_backlog:168 admin_backlog[a_] = max(0.0, admin_backlog[a_] - P.admin_bandwidth)169 M = rng.poisson(P.lam * producers.sum())170 if M == 0:171 continue172 init = rng.choice(np.where(producers)[0], size=M)173 plans = rng.choice(P.s, size=M, p=dist_now)174 # participants: k-1 network neighbours of the initiator (random if none)175 parts = np.empty((M, P.k - 1), int)176 for m in range(M):177 nb = np.where(A[init[m]] & producers)[0]178 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)179 # bits per recipient180 if use_code:181 bits_r = cur_code_len[plans]182 else:183 bits_r = np.full(M, float(L_raw))184 # decoding errors185 if use_records and P.verify_records:186 err = np.zeros((M, P.k - 1), bool)187 elif use_code:188 err = rng.random((M, P.k - 1)) < P.eps_code189 else:190 dd = np.abs(dialect[init][:, None] - dialect[parts]); err = rng.random((M, P.k - 1)) < P.eps * dd191 task_fail_decode = err.any(axis=1)192 # messages and routing193 msgs_per_recip = 1 if use_records else 2 # ack round without records194 if not hier:195 n_msgs = msgs_per_recip * (P.k - 1)196 # load: each message is processed by its recipient; the ack by the initiator197 np.add.at(agent_backlog, parts.ravel(), 1.0)198 if msgs_per_recip == 2:199 np.add.at(agent_backlog, init, float(P.k - 1))200 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)201 root_load_t = 0.0202 else:203 # relay through the LCA of initiator and each participant204 lvl = np.stack([lca_depth(paths, init, parts[:, c]) for c in range(P.k - 1)], axis=1) # 0 = root205 hops = 2 * (depth - lvl) # messages up and down206 n_msgs = hops.sum(axis=1)207 if regime == "federated":208 pass # locality is captured by the LCA level: same-module tasks never reach the root209 else:210 lvl[:] = 0 # central: everything is relayed through the root211 hops = np.full_like(hops, 2 * depth); n_msgs = hops.sum(axis=1)212 # admin load: the LCA administrator processes 2 messages per recipient213 keys = []214 for c in range(P.k - 1):215 for m in range(M):216 key = (int(lvl[m, c]), tuple(paths[init[m]][: lvl[m, c]])); keys.append(key)217 admin_backlog[key] = admin_backlog.get(key, 0.0) + 2.0218 wait = np.zeros(M)219 idx = 0220 for c in range(P.k - 1):221 for m in range(M):222 wait[m] = max(wait[m], math.ceil(admin_backlog[keys[idx]] / P.admin_bandwidth)); idx += 1223 np.add.at(agent_backlog, parts.ravel(), 1.0)224 wait = np.maximum(wait, np.ceil(agent_backlog[parts] / P.bandwidth).max(axis=1))225 root_load_t = float(admin_backlog.get((0, ()), 0.0))226 task_fail_latency = wait > P.tau_max227 success = ~task_fail_decode & ~task_fail_latency228 # bits per task = codeword/raw bits to each recipient + 1-bit acknowledgements (if any) + record overhead229 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)230 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())231 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)232 stats["latency_sum"] += float(wait[success].sum()); stats["root_load"] += root_load_t233 if shock_time is None or t < shock_time:234 stats["tasks_pre"] += M; stats["success_pre"] += int(success.sum()); stats["bits_pre"] += float(bits.sum())235 else:236 stats["tasks_post"] += M; stats["success_post"] += int(success.sum()); stats["bits_post"] += float(bits.sum())237 np.add.at(benefit, init[success], P.value_per_task); np.add.at(benefit, parts[success].ravel(), P.value_per_task)238 tot_bits = stats["bits"] + learn_cost239 succ = max(stats["success"], 1)240 b = benefit[producers] if producers.any() else benefit241 g = float(np.abs(b[:, None] - b[None, :]).mean() / (2 * b.mean())) if b.mean() > 0 else float("nan")242 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,243 "success_rate": stats["success"] / max(stats["tasks"], 1), "msgs_per_task": stats["messages"] / max(stats["tasks"], 1),244 "bits_per_task": stats["bits"] / max(stats["tasks"], 1), "bits_per_success_incl_setup": tot_bits / succ, "bits_per_success": stats["bits"] / succ,245 "tasks_per_kbit": succ / (tot_bits / 1000.0), "latency_mean": stats["latency_sum"] / succ, "fail_latency_share": stats["fail_latency"] / max(stats["tasks"], 1),246 "fail_decode_share": stats["fail_decode"] / max(stats["tasks"], 1), "admin_overhead_share": (stats["admin_bits"] + learn_cost) / max(tot_bits, 1e-9),247 "welfare": float(benefit.sum() - P.cost_per_bit * tot_bits), "welfare_per_agent": float((benefit.sum() - P.cost_per_bit * tot_bits) / P.N),248 "gini_benefit": g, "root_load_mean": stats["root_load"] / P.T,249 "success_pre": stats["success_pre"] / max(stats["tasks_pre"], 1), "success_post": stats["success_post"] / max(stats["tasks_post"], 1),250 "bits_per_task_pre": stats["bits_pre"] / max(stats["tasks_pre"], 1), "bits_per_task_post": stats["bits_post"] / max(stats["tasks_post"], 1),251 "kl_q_p_bits": float((q * np.log2(q / p)).sum()) if P.shock_type == "distribution" else 0.0, "tasks": stats["tasks"]}252 return out253