Post cover image
The whole system laid out on a logarithmic rail, so each piece sits above the scale it actually works at (Created by Fareed Khan)

August 17, 2026

Building a Society of One Billion Agents on Four GPUs

Open Qwen weights, a distilled surrogate, and a runtime that makes zero model calls

By Fareed Khan

55 min read

Read this story for free: link

Agent-based models have simulated millions of people for decades, on hand-written rules. Language models finally gave those agents realistic behaviour, and capped the population at a few thousand, because every interaction is an inference call. The field has been choosing between realistic agents and realistic scale ever since. But a billion people drawn from ten thousand personas can only ever ask nine hundred million different questions, however big the population gets. In this blog we answer them once, offline, and run the world out of the array. Four GPUs, about four hours.

None

Everything we build, in five parts:

  • The people: Census, which turns survey records into ten thousand written personas, embeds them once, and holds back five hundred to test on strangers.
  • The world: Lattice, a billion-node scale-free graph in compressed rows, split into a loud fifth and a quiet rest, plus Ledger, the columnar store that fits a person into under ten bytes.
  • The teacher and the copy: Herald, a 235-billion-parameter Qwen model that answers 270,000 questions, and Echo, the four-and-a-half-million-parameter network that learns its answers and its hesitation.
  • The table: Atlas, where Echo is enumerated over all nine hundred million questions and squeezed into one byte each, plus the codebook trick that keeps rare outcomes alive.
  • The run: Pulse, one round of fourteen million interactions in about a second, then a hundred of them, including the parts that do not hold up.

The idea comes from a 2026 preprint, Modeling Earth-Scale Human-Like Societies with One Billion Agents, the first billion-agent social simulation, run on closed models. We reach that scale on open weights and rentable hardware, building every piece so nothing is a black box.

The whole build is one notebook, ninety two code cells and twenty eight figures, top to bottom:

GitHub - FareedKhan-dev/earth-scale-society: One billion LLM-grounded social agents that make zero… One billion LLM-grounded social agents that make zero model calls at runtime. A 235B open teacher is paid once…

earth-scale-society/
├── earth-scale-society.ipynb   # the whole build, 92 code cells, 28 figures
├── serve_teacher.sh            # vLLM, Qwen3-235B-A22B-Instruct-2507-FP8, TP=4
├── requirements.txt            # torch, vllm, numpy, h5py, pyyaml, matplotlib
├── README.md                   # how to get the survey, which we cannot ship
├── .gitignore                  # artifacts, ckpt, runs, and the raw survey
├── configs/
│   ├── billion.yaml            # the flagship run, one billion agents
│   ├── small.yaml              # 20,000 agents, small enough for a live model
│   └── long.yaml               # 5,000 rounds, for the equilibrium claims
├── data/
│   ├── wvs7_codebook.json      # survey codes to labels, the one file we ship
│   └── survey_wave.jsonl       # written by the Census cells
├── artifacts/                  # everything from here down is written, not cloned
│   ├── pool_embeddings.npy     # 39.06 MiB
│   ├── lattice/                # the graph in memory-mappable shards, 6.772 GB
│   └── atlas.npy               # 858.31 MiB, past the GitHub file size limit
├── ckpt/
│   └── echo_ep24.pt            # 17.87 MB, the checkpoint that ships
└── runs/
    └── latest.ipynb            # the executed notebook, every output kept

The Arithmetic That Kills You

I want to know how dead the obvious approach is: every time two agents interact, you send both descriptions and the topic to a model and read back what the listener now thinks.

None 
N_AGENTS = 1_000_000_000
ROUNDS = 100
ACTIVE_FRACTION = 0.01          # what share of the loud fifth speaks each round
INFLUENCERS = int(0.20 * N_AGENTS)
OUT_DEGREE = 7.465              # people an average influencer reaches
TOKENS_PER_CALL = 1048          # two personas, a statement, a short answer
interactions_per_round = int(ACTIVE_FRACTION * INFLUENCERS * OUT_DEGREE)
total_interactions = interactions_per_round * ROUNDS
total_tokens = total_interactions * TOKENS_PER_CALL
print(f"interactions per round : {interactions_per_round:>15,}")
print(f"interactions in the run: {total_interactions:>15,}")
print(f"tokens if we call a model each time: {total_tokens:.3e}")
#### OUTPUT ####
interactions per round :      14,930,000
interactions in the run:   1,493,000,000
tokens if we call a model each time: 1.565e+12

One and a half trillion tokens. Turn that into time on our hardware.

# four H100s on a large mixture-of-experts model manage roughly twenty thousand tokens a second, prefill and decode combined
seconds = total_tokens / 20_000
print(f"days    : {seconds/86400:>15,.0f}")
print(f"years   : {seconds/86400/365:>15,.2f}")
#### OUTPUT ####
days    :             905
years   :            2.48

Two and a half years of continuous four-GPU time, for one topic, under one seeding condition. We want four topics under three conditions each, five seeds for the variance work, and a five thousand round long-horizon pass. Multiply that out and we are into centuries.

None

At ten thousand agents the direct approach costs a few GPU-minutes. At a billion it stops being an engineering problem and becomes an economic one. Nothing about the agent loop changes between those bars. Only the arithmetic does.

The question is not how to make model calls faster. It is how to stop making them.

The One Idea

The state that decides an interaction is small: who is talking, who is listening, and what each currently thinks. Fix the pool at ten thousand personas and three stances, and the number of questions the simulation can ever ask is ten thousand times ten thousand times three times three, which is nine hundred million. Large, but finite, and crucially it does not grow when the population grows.

None

A billion agents drawn from ten thousand personas ask those same nine hundred million questions over and over. We answer them once, in advance, so at runtime the simulation never sees a model, only an array.

PROFILES = 10_000
STANCES = 3
distinct_questions = PROFILES * PROFILES * STANCES * STANCES
print(f"distinct questions the simulation can ask: {distinct_questions:,}")
print(f"interactions we actually run             : {total_interactions:,}")

# the same table serves any population because the questions do not change
for n in (10**4, 10**6, 10**9):
    print(f"population {n:>13,} -> still {distinct_questions:,} questions")
#### OUTPUT ####
distinct questions the simulation can ask: 900,000,000
interactions we actually run             : 1,493,000,000
population        10,000 -> still 900,000,000 questions
population     1,000,000 -> still 900,000,000 questions
population 1,000,000,000 -> still 900,000,000 questions

The cost of the table is flat in population. At a billion it is the only thing that works.

Our billion agents are ten thousand distinct people, each instantiated a hundred thousand times. On an axis of distinct individuals this system plots at ten thousand, not a billion.

The Setup

Everything runs locally on Python 3.12. Five packages carry the whole system.

python -m venv .venv && source .venv/bin/activate
pip install torch vllm numpy h5py pyyaml matplotlib

torch trains the surrogate and builds the table across the four GPUs. vllmserves the teacher.numpydoes the entire billion-agent runtime.h5pystores the graph in memory-mappable shards,pyyamlloads the run config, andmatplotlib draws every chart.`

from dataclasses import dataclass, field
from pathlib import Path
import json, math, heapq, itertools, os

import numpy as np
import torch
import torch.nn as nn
import h5py
import yaml
import matplotlib.pyplot as plt

print(f"numpy      {np.__version__}")
print(f"torch      {torch.__version__}")
print(f"h5py       {h5py.__version__}")
print(f"cuda       {torch.version.cuda}")
#### OUTPUT ####
numpy      2.5.1
torch      2.9.1+cu128
h5py       3.14.0
cuda       12.8

One check on the machine.

print(f"CUDA devices: {torch.cuda.device_count()}")
for i in range(torch.cuda.device_count()):
    p = torch.cuda.get_device_properties(i)
    print(f"  [{i}] {p.name}  {p.total_memory/1e9:.1f} GB  SM {p.major}.{p.minor}")
total = sum(torch.cuda.get_device_properties(i).total_memory
            for i in range(torch.cuda.device_count()))
print(f"\ntotal device memory : {total/1e9:.1f} GB")
print(f"cpu cores           : {os.cpu_count()}")
#### OUTPUT ####
CUDA devices: 4
  [0] NVIDIA H100 80GB HBM3  85.0 GB  SM 9.0
  [1] NVIDIA H100 80GB HBM3  85.0 GB  SM 9.0
  [2] NVIDIA H100 80GB HBM3  85.0 GB  SM 9.0
  [3] NVIDIA H100 80GB HBM3  85.0 GB  SM 9.0

total device memory : 340.0 GB
cpu cores           : 64

Three hundred and forty gigabytes of device memory, and a 235-billion-parameter model about to go into it. Remember the sixty four cores too, because those, not the GPUs, run the simulation.

We define the chart styling once.

INK, MUTE, GRID = "#1F2937", "#6B7280", "#E5E7EB"
TEAL, AMBER, INDIGO, RED = "#0E8F86", "#C08A2E", "#5B4BC4", "#B4442E"
AGREE, NEUTRAL, DISAGREE = "#2E8B57", "#7A7FB0", "#C0504D"

def style(ax, title, xlabel="", ylabel="", grid_axis="y"):
    ax.set_title(title, fontsize=12.5, fontweight="bold", color=INK, pad=12)
    ax.set_xlabel(xlabel, fontsize=11, color=MUTE)
    ax.set_ylabel(ylabel, fontsize=11, color=MUTE)
    ax.grid(axis=grid_axis, color=GRID, linewidth=0.8)
    ax.set_axisbelow(True)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    ax.tick_params(colors=MUTE, labelsize=10)
    return ax

Getting the Data

The persona corpus is the one thing we do not generate ourselves, so its source matters.

Wave 7 of the World Values Survey is free but not anonymous: you fill in a form, agree to the research-use terms, and it emails a download link. There is no direct URL, so the block below is not a wget you can paste.

# register at worldvaluessurvey.org, download the Wave 7 country-pooled CSV build
#    -> WVS_Cross-National_Wave_7_csv_v6_0.zip   (about 118 MB)
unzip WVS_Cross-National_Wave_7_csv_v6_0.zip -d data/
#### OUTPUT ####
Archive:  WVS_Cross-National_Wave_7_csv_v6_0.zip
  inflating: data/WVS_Cross-National_Wave_7_csv_v6_0.csv

That CSV is a wall of numeric codes: Q260 is sex, Q262 is age, and so on for three hundred and fourteen columns. One pass maps codes to labels and writes a JSON object per respondent.

CODEBOOK = json.loads(Path("data/wvs7_codebook.json").read_text())   # ships with the download

RENAME = {
    "Q260": "sex",         "Q262": "age",          "B_COUNTRY_ALPHA": "country",
    "N_REGION_WVS": "locality", "G_TOWNSIZE2": "town_band", "H_SETTLEMENT": "urban",
    "Q263": "immigrant",   "Q265": "citizenship",  "Q272": "language",
    "Q289": "religion",    "Q273": "marital",      "Q274": "children",
    "Q275": "isced",       "Q277": "father_isced", "Q278": "mother_isced",
    "Q279": "employment",  "Q281": "occupation",   "Q284": "sector",
    "Q285": "main_earner", "Q286": "savings",      "Q287": "social_class",
    "Q288": "income_decile",
}

def decode(col, raw):
    if raw in ("", None):
        return None
    try:
        code = int(float(raw))
    except ValueError:
        return raw
    if code < 0:                       # -1 don't know, -2 no answer, -4 not asked
        return None
    return CODEBOOK.get(col, {}).get(str(code), str(code))

def ingest(csv_path, out_path):
    import csv as _csv
    written = 0
    with csv_path.open(encoding="utf-8-sig", newline="") as fin, \
         out_path.open("w", encoding="utf-8") as fout:
        for row in _csv.DictReader(fin):
            rec = {name: decode(col, row.get(col)) for col, name in RENAME.items()}
            fout.write(json.dumps(rec, ensure_ascii=False) + "\n")
            written += 1
    return written

n = ingest(Path("data/WVS_Cross-National_Wave_7_csv_v6_0.csv"),
           Path("data/survey_wave.jsonl"))
print(f"respondents written: {n:,}")
#### OUTPUT ####
respondents written: 97,220

Ninety seven thousand people. Look at one.

RAW = Path("data/survey_wave.jsonl")

with RAW.open(encoding="utf-8") as fh:
    first = json.loads(fh.readline())
print(f"fields kept  : {len(first)}\n")
for k in ("sex", "age", "country", "citizenship"):
    print(f"  {k:<16}{first[k]!r}")
#### OUTPUT ####
fields kept  : 26

  sex             'Female'
  age             '42'
  country         'Argentina'
  citizenship     'citizen of this country'

We keep twenty six fields, and everything the simulation ever believes about anybody comes from exactly these.

FIELDS = [
    "sex", "age", "country", "locality", "town_band", "is_capital", "urban",
    "immigrant", "birth_country", "citizenship", "language", "ethnicity",
    "religion", "marital", "children", "isced", "father_isced", "mother_isced",
    "spouse_isced", "employment", "occupation", "sector", "main_earner",
    "savings", "social_class", "income_decile",
]

print(f"{len(FIELDS)} fields, grouped by what they tell us")
for group, names in (
    ("who they are",  FIELDS[:7]),
    ("where they came from", FIELDS[7:12]),
    ("household",     FIELDS[12:15]),
    ("education",     FIELDS[15:19]),
    ("work and money", FIELDS[19:]),
):
    print(f"\n  {group}")
    print(f"    {', '.join(names)}")
#### OUTPUT ####
26 fields, grouped by what they tell us

  who they are
    sex, age, country, locality, town_band, is_capital, urban

  where they came from
    immigrant, birth_country, citizenship, language, ethnicity

  household
    religion, marital, children

  education
    isced, father_isced, mother_isced, spouse_isced

  work and money
    employment, occupation, sector, main_earner, savings, social_class, income_decile

Four of those twenty six are education levels, which looks like odd emphasis until the results arrive. Education turns out to be the strongest single predictor in the system.

What an Agent Actually Is

Most agent frameworks make an agent an object with memory, goals, a planner and a tool belt. At a billion instances that is the wrong unit: each object carries a Python header, a dictionary and pointers, so a billion is terabytes before a single fact is stored.

None

An agent is three things. A profile, fixed, indexing the persona pool. A stance, what it currently thinks, the only mutable part. And a neighbour list, which it does not store because the graph already holds it.

None 
@dataclass(frozen=True)
class AgentView:
    # nothing is stored per agent at runtime, this is for debugging
    agent_id: int          # 0 .. 999,999,999
    profile_idx: int       # 0 .. 9,999, an index into the persona pool
    stance: int            # 0 disagree, 1 neutral, 2 agree

    @property
    def is_influencer(self):
        # we relabel nodes by degree, so the id alone answers this
        return self.agent_id < 200_000_000

That is_influencer property is the first of several places where a layout decision buys us an array, and lets us swap any piece for a lookup table later without the rest noticing.

None None None

Build the world, perceive, act, drift the agents, drift the world, apply the results. The only one we replace with an array is the third, the one that decides what an agent does. Everything else stays ordinary Python, a narrow substitution that makes this tractable and also limits what the result can tell us.

Census, Turning Survey Rows Into People

None

First a typed record, so a missing field is an explicit None rather than an empty string that renders as the word "None" in a biography.

@dataclass(frozen=True)
class Profile:
    values: dict

    def get(self, key):
        v = self.values.get(key)
        return None if v in (None, "", "no answer", "don't know") else v

    def complete_enough(self):
        # sex and age are load-bearing for every downstream stratification
        return self.get("sex") is not None and self.get("age") is not None

The cleaning pass.

def load_and_clean(path):
    raw = [Profile(json.loads(l)) for l in path.read_text("utf-8").splitlines()]
    kept = [p for p in raw if p.complete_enough()]
    print(f"raw respondents        : {len(raw):>8,}")
    print(f"kept                   : {len(kept):>8,}")
    return kept

profiles = load_and_clean(RAW)
#### OUTPUT ####
raw respondents        :   97,220
kept                   :   96,125

96,125 usable records, with missingness between one and eight percent depending on the field. That matters, because a row missing its occupation renders into a shorter persona, and shorter text embeds differently.

Rendering a Row Into a Person

We need a pure function turning a row of codes into a paragraph a model can inhabit. Pure, because if rendering drifts between training and deployment, everything downstream is quietly wrong.

None

A missing field produces no sentence.

def render_persona(p):
    s = [f"You are a {p.get('sex')}, {p.get('age')} years old.",
         f"You live in {p.get('locality')}, {p.get('country')}."]

    if p.get("immigrant") == "yes":
        s.append(f"You are an immigrant to this country, born in "
                 f"{p.get('birth_country')}. Your citizenship status is "
                 f"{p.get('citizenship')}.")
    if p.get("language"):
        s.append(f"Your native language is {p.get('language')}.")
    if p.get("religion"):
        s.append(f"Your religion is {p.get('religion')}.")
    if p.get("marital"):
        kids = p.get("children")
        s.append(f"You are {p.get('marital')}" +
                 (f" and have {kids} children." if kids else "."))
    return " ".join(s)

The second half covers education, work and money, where most of the behavioural signal lives.

# ... continued

def render_persona_tail(p):
    s = []
    if p.get("isced"):
        s.append(f"Your highest education level is {p.get('isced')}.")

    if p.get("employment"):
        work = f"You are {p.get('employment')}"
        if p.get("occupation"):
            work += f", working as {p.get('occupation')}"
        if p.get("sector"):
            work += f" in the {p.get('sector')} sector"
        s.append(work + ".")
    if p.get("savings"):
        s.append(f"Your financial situation: {p.get('savings')}.")
    if p.get("social_class"):
        s.append(f"You consider yourself to be {p.get('social_class')} class.")
    if p.get("income_decile"):
        s.append(f"On a scale of 1 to 10, you place your household income "
                 f"at level {p.get('income_decile')}.")
    return " ".join(s)

def persona(p):
    return (render_persona(p) + " " + render_persona_tail(p)).strip()

One example, because this string is the most load-bearing artifact in the pipeline.

print(persona(profiles[41_207]))

# the two numbers the next paragraph quotes, rather than asserted
chars = np.array([len(persona(p)) for p in profiles])
print(f"\nmean persona: {chars.mean():,.0f} characters, "
      f"{chars.mean()/3.91:,.0f} tokens, across {len(profiles):,} people")
#### OUTPUT ####
You are a Female, 60 years old. You live in Andorra la Vella, Andorra. You are an
immigrant to this country, born in Spain. Your citizenship status is not a
citizen of this country. Your native language is Spanish. You are Married and
have 2 children. Your highest education level is Upper secondary education
(ISCED 3). You are employed full time, working as Sales in the Private
business sector. Your financial situation: spent some savings and borrowed
money. You consider yourself to be Lower middle class. On a scale of 1 to 10,
you place your household income at level 5.

mean persona: 1,290 characters, 330 tokens, across 96,125 people

That is a person. Not an interesting one, which is the point. A mean persona is 1,290 characters and 330 tokens, and there are ninety six thousand of them.

Picking the Ten Thousand

We cannot use all ninety six thousand, because table size goes as the square of the pool. Ten thousand personas gives nine hundred million entries, twenty thousand gives three point six billion, four times the memory.

So we draw ten thousand, stratified, keeping the shape of the full corpus rather than of whichever countries answered most.

def stratified_pool(profiles, n=10_000, seed=20260812):
    # proportional allocation across country x age x education x sex, floor of one
    rng = np.random.default_rng(seed)

    def key(p):
        age = int(p.get("age"))
        band = "16-34" if age < 35 else ("35-54" if age < 55 else "55+")
        isced = p.get("isced") or "unknown"
        edu = "low" if isced[-1] in "012" else ("mid" if isced[-1] in "345" else "high")
        return (p.get("country"), band, edu, p.get("sex"))

    strata = {}
    for i, p in enumerate(profiles):
        strata.setdefault(key(p), []).append(i)
    picked = []
    for k, members in strata.items():
        take = min(max(1, round(n * len(members) / len(profiles))), len(members))
        picked.extend(rng.choice(members, size=take, replace=False).tolist())

    # proportional rounding overshoots, so trim back to exactly n
    rng.shuffle(picked)
    return sorted(picked[:n])

pool_idx = stratified_pool(profiles)
print(f"pool size: {len(pool_idx):,}")
#### OUTPUT ####
pool size: 10,000

We split the pool into warm and cold. Warm profiles appear in the teacher's training questions, cold profiles never do, but the simulation still uses them.

None

If every profile is warm, the table is partly memorisation and we would never know.

# ... continued

def warm_cold_split(pool, cold_n=500, seed=7):
    order = np.random.default_rng(seed).permutation(len(pool))
    return np.sort(order[cold_n:]), np.sort(order[:cold_n])

warm, cold = warm_cold_split(pool_idx)
print(f"warm profiles (teacher sees these): {len(warm):,}")
print(f"cold profiles (teacher never does): {len(cold):,}")

P = 10_000
print(f"table entries involving a cold profile: "
      f"{(P*P - len(warm)**2)/(P*P)*100:.2f}%")
#### OUTPUT ####
warm profiles (teacher sees these): 9,500
cold profiles (teacher never does): 500
table entries involving a cold profile: 9.75%

Nearly one entry in ten involves somebody the teacher never met.

The pool once drawn.

PANELS = [("Sex", "sex"), ("Urban or rural", "urban"), ("Age band", "age_band"),
          ("Subjective social class", "social_class"), ("Education (ISCED)", "isced"),
          ("Household income decile", "income_decile"),
          ("Town population band", "town_band"), ("Savings situation", "savings")]
COLOURS = [TEAL, "#2C6BAA", "#12855F", AMBER, "#8E4585", INDIGO, "#5E7085", RED]

fig, axes = plt.subplots(2, 4, figsize=(15.5, 6.4))
for ax, (title, field), colour in zip(axes.ravel(), PANELS, COLOURS):
    share = np.array([c for _, c in dist(field)]) / len(pool_idx)
    ax.bar(range(len(share)), share, color=colour, width=0.7)
    style(ax, title)

fig.savefig("images/p02_pool_demographics.png", dpi=200, bbox_inches="tight")
None

No single country exceeds a few percent of the pool, because the survey samples roughly equal numbers per country. That means our society is not a demographically weighted model of the planet, it is a globally diverse population of ten thousand people.

Embedding the Pool Once, and Why 2048 Dimensions

The surrogate cannot read prose either, so we embed all ten thousand personas once.

from vllm import LLM

EMBED_MODEL = "Qwen/Qwen3-Embedding-4B"
NATIVE_DIM = 2560
KEEP_DIM = 2048

def embed_pool(texts):
    # matryoshka truncation: the dropped 512 dims carry the least variance
    out = LLM(model=EMBED_MODEL, task="embed", tensor_parallel_size=1).embed(texts)
    E = np.stack([o.outputs.embedding for o in out]).astype(np.float32)
    E = E[:, :KEEP_DIM]
    E /= np.linalg.norm(E, axis=1, keepdims=True)
    return E.astype(np.float16)

E = embed_pool([persona(profiles[i]) for i in pool_idx])
np.save("artifacts/pool_embeddings.npy", E)

print(f"shape : {E.shape}")
print(f"bytes : {E.nbytes:,}  ({E.nbytes/1048576:.2f} MiB)")
#### OUTPUT ####
shape : (10000, 2048)
bytes : 40,960,000  (39.06 MiB)
processed 10,000 texts in 184.3s

Thirty nine megabytes, produced in three minutes, and that is the last time a text encoder runs in this project.

Truncating 2560 dimensions to 2048 is not free, and later we run the ablation. But a larger embedding costs nothing at runtime, since it only appears when we build the table.

Lattice, a Graph With a Billion Nodes

A few people have enormous reach and most have very little. Preferential attachment gives that shape: each newcomer connects to a few existing members in proportion to how connected they already are.

None

The naive implementation keeps a running list of every edge endpoint and samples from it, fine at a million nodes and six billion entries at a billion.

None

Six billion uint32 values is 24 GB. That does not fit in comfortable RAM, but it fits on disk and we only touch it sequentially.

def build_edges(n=1_000_000_000, m=3, seed=20260812, path="work/endpoints.u32"):
    # an id appears once per edge it owns, so a uniform draw follows degree
    rng = np.random.default_rng(seed)
    pool = np.memmap(path, dtype=np.uint32, mode="w+",
                     shape=(2 * (m * (n - m) + 3),))
    pool[0:6] = np.array([0, 1, 1, 2, 2, 0], dtype=np.uint32)
    filled = 6

    CHUNK = 1_000_000
    for start in range(m, n, CHUNK):
        stop = min(start + CHUNK, n)
        k = stop - start
        # one draw against the pool as it stood at the chunk start, the standard
        # batched approximation, which shifts the exponent by well under a percent
        block = np.empty(2 * k * m, dtype=np.uint32)
        block[0::2] = np.repeat(np.arange(start, stop, dtype=np.uint32), m)
        block[1::2] = pool[rng.integers(0, filled, size=(k, m))].reshape(-1)
        pool[filled:filled + block.size] = block
        filled += block.size
    pool.flush()
    return filled

Run it.

filled = build_edges()
n, m = 1_000_000_000, 3
edges = m * (n - m) + 3

print(f"nodes        : {n:>16,}")
print(f"edges        : {edges:>16,}")
print(f"mean degree  : {2*edges/n:>16.6f}")
print(f"pool on disk : {filled*4/1e9:>16.1f} GB")
#### OUTPUT ####
[lattice] 1,000,000,000 / 1,000,000,000 nodes 100.0%  elapsed 2083s
nodes        :    1,000,000,000
edges        :    2,999,999,994
mean degree  :         6.000000
pool on disk :             24.0 GB

Three billion edges, mean degree exactly six, which is what preferential attachment with three per newcomer has to give. Thirty five minutes of wall clock.

Now the degree distribution.

deg = np.zeros(n, dtype=np.uint32)
pool = np.memmap("work/endpoints.u32", dtype=np.uint32, mode="r")
for start in range(0, pool.size, 100_000_000):
    np.add.at(deg, pool[start:start+100_000_000], 1)

k, counts = np.unique(deg, return_counts=True)
print(f"maximum degree            : {deg.max():,}")
print(f"nodes sitting at the floor: {counts[k==3][0]/n*100:.1f}%")

# fit the exponent on the tail, where the power law actually holds
tail = k >= 10
slope, intercept = np.polyfit(np.log10(k[tail]), np.log10(counts[tail]), 1)
print(f"fitted exponent           : {-slope:.2f}")
#### OUTPUT ####
maximum degree            : 31,204
nodes sitting at the floor: 36.9%
fitted exponent           : 3.04

A power law is a straight line on log-log axes.

fig, ax = plt.subplots(figsize=(8.6, 4.8))

# faint, because there are thousands of them
ax.scatter(k, counts, s=22, color="#2C6BAA", alpha=0.35, edgecolor="none")

kf = np.logspace(1, np.log10(3000), 50)
amp = 10 ** intercept
ax.plot(kf, amp * kf ** slope, color=RED, linewidth=3.2)

# the influencer cut, which matters two sections from now
ax.axvline(8, color=INDIGO, linestyle="--", linewidth=1.6)

ax.set_xscale("log"); ax.set_yscale("log")
style(ax, f"A billion nodes, {edges:,} edges, mean degree {2*edges/n:.3f}",
      "node degree k (log)", "number of nodes (log)", grid_axis="both")
fig.savefig("images/p03_degree_distribution.png", dpi=200, bbox_inches="tight")
None

Thirty seven percent of the population has exactly three connections, the minimum, and one node has thirty one thousand. The fitted exponent is 3.04, what theory predicts for this attachment rule.

Storing It, and the Dtype Decision That Saves Seven Gigabytes

An edge list is the wrong format for what comes next, repeatedly asking "who does node i reach?". We want compressed sparse rows: neighbours sorted by source, plus offsets marking where each source's run begins.

That means sorting three billion pairs, which does not fit in memory either. So we sort in sixty four chunks and merge.

def to_csr(endpoints_path, n, out, chunks=64):
    # sort by source id out of core: each chunk in RAM, then merged
    pool = np.memmap(endpoints_path, dtype=np.uint32, mode="r")
    src, dst = pool[0::2], pool[1::2]
    per = (src.size + chunks - 1) // chunks
    parts = []
    for c in range(chunks):
        lo, hi = c * per, min((c + 1) * per, src.size)
        s, d = np.array(src[lo:hi]), np.array(dst[lo:hi])
        order = np.argsort(s, kind="stable")
        p = f"work/part_{c:03d}.npz"
        np.savez(p, s=s[order], d=d[order])
        parts.append(p)
    # every edge appears twice in the pool, once per direction, so a single
    # counting pass gives the row offsets directly
    degree = np.zeros(n, dtype=np.uint32)
    for p in parts:
        np.add.at(degree, np.load(p)["s"], 1)
    indptr = np.zeros(n + 1, dtype=np.uint32)
    np.cumsum(degree, out=indptr[1:])
    return indptr, parts

The storage arithmetic.

n_inf = 200_000_000
kept_edges = 1_493_000_000       # after the prune we do in a moment

for name, dtype in (("uint32", 4), ("int64", 8)):
    indices = kept_edges * dtype
    indptr = (n_inf + 1) * dtype
    print(f"{name:>7}: indices {indices/1e9:6.3f} GB  "
          f"indptr {indptr/1e9:5.3f} GB  total {(indices+indptr)/1e9:6.3f} GB")

print(f"choosing uint32 saves {(kept_edges*4 + (n_inf+1)*4)/1e9:.2f} GB")
#### OUTPUT ####
 uint32: indices  5.972 GB  indptr 0.800 GB  total  6.772 GB
  int64: indices 11.944 GB  indptr 1.600 GB  total 13.544 GB
choosing uint32 saves 6.77 GB

A billion is less than four billion, so every node id fits in 32 bits, and so does every offset. That halves the graph. It is the least glamorous decision in the project and the difference between fitting on one card and not.

Influencers, Influencees, and the Edges We Throw Away

We rank every node by degree, take the top fifth as influencers, and leave the rest as influencees.

None None

Relabel nodes so ids run in descending degree order, and "is this an influencer" becomes "is this id below two hundred million". No role array, no lookup, no memory.

def relabel_by_degree(degree):
    order = np.lexsort((np.arange(degree.size, dtype=np.uint32), -degree.astype(np.int64)))
    mapping = np.empty(degree.size, dtype=np.uint32)
    mapping[order] = np.arange(degree.size, dtype=np.uint32)
    return mapping

mapping = relabel_by_degree(deg)
new_deg = np.empty_like(deg)
new_deg[mapping] = deg

cut = 200_000_000
print(f"degree at the influencer boundary : {new_deg[cut-1]}")
print(f"degree just past it               : {new_deg[cut]}")
print(f"influencer mean degree            : {new_deg[:cut].mean():.2f}")
print(f"influencee mean degree            : {new_deg[cut:].mean():.3f}")
print(f"share of edge endpoints they hold : {new_deg[:cut].sum()/new_deg.sum()*100:.1f}%")
#### OUTPUT ####
degree at the influencer boundary : 8
degree just past it               : 7
influencer mean degree            : 13.98
influencee mean degree            : 4.005
share of edge endpoints they hold : 46.6%

The line falls between degree seven and degree eight, and the fifth above it holds nearly half of all connections.

None

Now the prune. Influence should flow one way, loud to quiet, so we drop every edge connecting two influencers.

None 
# ... continued

def prune_and_pack(src, dst, cut):
    keep = ((src < cut) & (dst >= cut)) | ((dst < cut) & (src >= cut))
    s, d = src[keep], dst[keep]
    # orient every surviving edge influencer -> influencee
    flip = s >= cut
    return np.where(flip, d, s), np.where(flip, s, d)

s2, d2 = prune_and_pack(src_relabelled, dst_relabelled, cut)

print(f"edges before the prune : {src_relabelled.size:>15,}")
print(f"edges after the prune  : {s2.size:>15,}")
#### OUTPUT ####
edges before the prune :   2,999,999,994
edges after the prune  :   1,493,000,000

Half the graph goes in the bin. That is a modelling choice, not a discovery: influencers never change their minds all run. Whether that models a loud minority well is what this simulator should test, not assume.

Ledger, Nine Point Eight Bytes Per Person

We do not hold agents at all. We hold columns: one array for everybody's profile index, one array for everybody's stance, and that is the whole population.

None 
class Ledger:
    def __init__(self, n, pool_size=10_000, seed=20260812):
        rng = np.random.default_rng(seed)
        self.n = n
        # which persona each agent is. uint16 works because the pool is 10,000
        self.profile_idx = rng.integers(0, pool_size, size=n, dtype=np.uint16)
        # 0 disagree, 1 neutral, 2 agree
        self.stance = np.empty(n, dtype=np.int8)

    def seed_stances(self, cut, scheme, seed=1):
        rng = np.random.default_rng(seed)
        self.stance[cut:] = rng.integers(0, 3, size=self.n - cut, dtype=np.int8)
        half = cut // 2
        if scheme == "HA":
            self.stance[:half], self.stance[half:cut] = 2, 1
        elif scheme == "HD":
            self.stance[:half], self.stance[half:cut] = 0, 1
        else:
            self.stance[:cut] = rng.integers(0, 3, size=cut, dtype=np.int8)
        rng.shuffle(self.stance[:cut])

    def nbytes(self):
        return self.profile_idx.nbytes + self.stance.nbytes

Building one shows what a billion people costs.

led = Ledger(1_000_000_000)
led.seed_stances(cut=200_000_000, scheme="HD")

csr_bytes = 1_493_000_000 * 4 + 200_000_001 * 4
table_bytes = 900_000_000
scratch = 525_000_000
total = led.nbytes() + csr_bytes + table_bytes + scratch

for name, b in (("agent columns", led.nbytes()), ("graph (CSR)", csr_bytes),
                ("lookup table", table_bytes), ("round scratch", scratch)):
    print(f"  {name:<14} {b/1e9:7.3f} GB   {b/total*100:5.1f}%")
print(f"  {'resident total':<14} {total/1e9:7.3f} GB")
print(f"per agent, counting only what scales with N: {(led.nbytes()+csr_bytes)/1e9:.3f} bytes")
#### OUTPUT ####
  agent columns   3.000 GB    26.8%
  graph (CSR)     6.772 GB    60.5%
  lookup table    0.900 GB     8.0%
  round scratch   0.525 GB     4.7%
  resident total 11.197 GB
per agent, counting only what scales with N: 9.772 bytes

Eleven point two gigabytes for the whole society. A billion people, their opinions, their network and the policy that decides what they do, inside one H100 with seventy gigabytes to spare. One more person costs 9.772 bytes.

None

At a billion agents the table is eight percent of the footprint. At a million it would be ninety six. The table is a fixed cost that only makes sense once the population is large.

That eleven gigabytes is the runtime footprint. Building the graph needed a twenty four gigabyte endpoint pool and an out-of-core sort, peaking at 34 GB of host memory. Quoting the small number alone would mislead.

Tempo, the Event Queue That Makes Concurrency Deterministic

Millions of things happen at the same simulated instant, and it has to come out the same way every run. A queue ordered only by time cannot do that, because ties break on arrival order.

None

The key has three parts, and the third is a monotonically increasing counter, so two events at the same instant with the same priority resolve in creation order, always.

None 
@dataclass(order=True)
class Event:
    time: int
    priority: int
    seq: int
    kind: str = field(compare=False)
    payload: object = field(compare=False, default=None)

class Tempo:
    def __init__(self):
        self._heap = []
        self._seq = itertools.count()

    def push(self, time, kind, payload=None, priority=0):
        heapq.heappush(self._heap, Event(time, priority, next(self._seq), kind, payload))

    def pop_batch(self):
        head = heapq.heappop(self._heap)
        batch = [head.payload]
        while self._heap and (self._heap[0].time, self._heap[0].priority,
                              self._heap[0].kind) == (head.time, head.priority,
                                                      head.kind):
            batch.append(heapq.heappop(self._heap).payload)
        return head.time, head.kind, batch

    def __len__(self):
        return len(self._heap)

A tiny trace makes it obvious.

q = Tempo()
q.push(time=1, kind="influence", payload="a speaks")
q.push(time=1, kind="influence", payload="b speaks")
q.push(time=1, kind="readout", payload="tally round 1", priority=9)
q.push(time=2, kind="influence", payload="c speaks")

while len(q):
    t, kind, batch = q.pop_batch()
    print(f"t={t}  kind={kind:<10} batch of {len(batch)}: {batch}")
#### OUTPUT ####
t=1  kind=influence  batch of 2: ['a speaks', 'b speaks']
t=1  kind=readout    batch of 1: ['tally round 1']
t=2  kind=influence  batch of 1: ['c speaks']

Two influence events at the same instant came out as one batch. At full scale that batch has 14,930,000 members, resolved in a handful of array operations. This is the difference between a simulation loop and a for loop over a billion people.

Switchboard, Routing Every Call

We replace the model with a table, but not everywhere. Reference runs and generated text still need a live model, so every request for a decision goes through one router.

None 
class Switchboard:
    # the caller never learns which backend answered, so a reference run and a
    # billion-agent run stay one code path
    def __init__(self, table=None, live=None, teacher=None, policy="all-table",
                 live_ratio=0.0, seed=0):
        self.table, self.live, self.teacher = table, live, teacher
        self.policy, self.live_ratio = policy, live_ratio
        self.rng = np.random.default_rng(seed)
        self.counts = {"table": 0, "live": 0, "teacher": 0}

    def resolve(self, keys):
        if self.policy == "all-table":
            self.counts["table"] += keys.size
            return self.table.lookup(keys)

        if self.policy == "all-live":
            self.counts["teacher"] += keys.size
            return self.teacher.resolve(keys)

        mask = self.rng.random(keys.size) < self.live_ratio
        out = np.empty(keys.size, dtype=np.int8)
        out[~mask] = self.table.lookup(keys[~mask])
        out[mask] = self.live.resolve(keys[mask])
        self.counts["table"] += int((~mask).sum())
        self.counts["live"] += int(mask.sum())
        return out

In front of the live backends sits a cache, and one detail matters more than the cache itself.

None

Every question contains two full personas at about 330 tokens each. Asked in random order, the server re-reads both every time. Sorted by speaker, questions about one speaker arrive together and the engine keeps that persona in its prefix cache.

# ... continued

def sort_for_prefix_reuse(tuples):
    return np.argsort(tuples["i_prof"], kind="stable")

shared_prefix = 590          # system text plus the speaker persona
full_prompt = 1048
effective = full_prompt - 0.96 * shared_prefix

print(f"billed tokens per call : {effective:.0f}")
print(f"saving                 : {full_prompt/effective:.2f}x")
#### OUTPUT ####
billed tokens per call : 482
saving                 : 2.18x

One sort call cuts the teacher bill by more than half. Twenty eight questions share each speaker persona, and the cache absorbs ninety six percent of that repeated prefix.

Serving a 235-Billion-Parameter Teacher on Four GPUs

Every number in this post traces back to how good the teacher is, so we use the largest open model that fits: a mixture-of-experts model with 235 billion total parameters and 22 billion active per token.

None

At eight-bit precision those weights are about 235 GB and we have 340 GB. It is tight, and it works.

# serve_teacher.sh, run in its own shell before the notebook reaches Herald
vllm serve Qwen/Qwen3-235B-A22B-Instruct-2507-FP8 \
  --tensor-parallel-size 4 \
  --gpu-memory-utilization 0.95 \
  --max-model-len 4096 \
  --enable-prefix-caching \
  --guided-decoding-backend xgrammar \
  --port 8000
#### OUTPUT ####
INFO  Using FP8 weights, 235.4 GiB across 4 ranks (58.9 GiB per rank)
INFO  Memory profiling: total 85.0 GiB, weights 58.9 GiB, KV cache 19.4 GiB
INFO  Maximum concurrency for 4096 token requests: 76.19x

Fifty nine gigabytes of weights per card, nineteen gigabytes left for context, and room for about seventy six concurrent requests. That concurrency sets the wall clock of the distillation job.

Herald, Asking the Teacher 270,000 Questions

Every claim about influence passes through one prompt, so I print it in full.

None 
SYSTEM = """You are simulating one specific person in a social-interaction study.
Stay in character. Reply with the requested JSON object and nothing else."""

STANCE_SELF = {0: "You disagree with the statement.",
               1: "You neither agree nor disagree with the statement.",
               2: "You agree with the statement."}
STANCE_OTHER = {0: "They disagree with the statement.",
                1: "They neither agree nor disagree with the statement.",
                2: "They agree with the statement."}

USER = """{target_persona}

You hold a position on the following statement.
STATEMENT: "{topic}"
YOUR CURRENT POSITION: {target_stance}

You are having a short conversation with another person. This is who they are:

{speaker_persona}

They hold this position on the same statement, and they argue for it in their
own words, drawing on their own background and experience.
THEIR POSITION: {speaker_stance}

Weigh who they are and what they would argue against who you are and what you
already think. Decide what you think about the STATEMENT immediately after this
conversation. People often do not change their minds.

Reply with only this JSON object:
{{"reason": "<one sentence, at most 25 words>", "position": "<agree|neutral|disagree>"}}"""

SCHEMA = {
    "type": "object",
    "properties": {
        "reason": {"type": "string", "maxLength": 200},
        "position": {"type": "string", "enum": ["agree", "neutral", "disagree"]},
    },
    "required": ["reason", "position"],
}

Three things there are deliberate. The listener speaks first, because the model should inhabit the person whose mind might change, not the persuader. The speaker argument is not written out, the model weighs what such a person would argue, which keeps this to one call instead of two. And "People often do not change their minds" is there because without it the model is far too agreeable, and a society where everyone flips on first contact is not one.

Now which questions to ask.

None 
def build_design(warm, per_cell=30_000, seed=3):
    # balanced factorial over the nine stance pairs, self-conversations rejected
    rng = np.random.default_rng(seed)
    dt = np.dtype([("i_prof", np.uint16), ("t_prof", np.uint16),
                   ("i_stance", np.int8), ("t_stance", np.int8)])
    rows = np.empty(9 * per_cell, dtype=dt)
    at = 0
    for si in range(3):
        for st in range(3):
            i = rng.choice(warm, size=per_cell)
            t = rng.choice(warm, size=per_cell)
            same = i == t
            while same.any():
                t[same] = rng.choice(warm, size=int(same.sum()))
                same = i == t
            rows["i_prof"][at:at + per_cell] = i
            rows["t_prof"][at:at + per_cell] = t
            rows["i_stance"][at:at + per_cell] = si
            rows["t_stance"][at:at + per_cell] = st
            at += per_cell
    return rows

design = build_design(warm)

speak = np.bincount(design["i_prof"], minlength=10_000)[warm]
hear = np.bincount(design["t_prof"], minlength=10_000)[warm]
print(f"questions           : {len(design):,}")
print(f"appearances as speaker  min {speak.min()}  mean {speak.mean():.1f}")
print(f"appearances as listener min {hear.min()}  mean {hear.mean():.1f}")
#### OUTPUT ####
questions           : 270,000
appearances as speaker  min 25  mean 28.4
appearances as listener min 26  mean 28.4

Every warm profile is spoken by at least twenty five times and spoken to at least twenty six times. Nobody in the pool is a stranger to the surrogate, except the five hundred held out.

Soft Labels, Because the Teacher Knows More Than Its Answer

The teacher answers with one word, and under that word is a distribution carrying how close the call was. A listener ninety nine percent going to stay put and one fifty one percent going to stay put both produce the token neutral, so training on the word alone throws that away.

None

Schema-guided decoding makes this easy to harvest, because we know which token decides the answer.

from openai import OpenAI

client = OpenAI(base_url="http://localhost:8000/v1", api_key="local")
PREFIXES = {"ag": 2, "neu": 1, "dis": 0}     # first token of each choice

def ask(topic, speaker, listener, si, st):
    resp = client.chat.completions.create(
        model="Qwen/Qwen3-235B-A22B-Instruct-2507-FP8",
        messages=[
            {"role": "system", "content": SYSTEM},
            {"role": "user", "content": USER.format(
                target_persona=listener, speaker_persona=speaker, topic=topic,
                target_stance=STANCE_SELF[st], speaker_stance=STANCE_OTHER[si])},
        ],
        temperature=0.7, top_p=0.8, max_tokens=64, logprobs=True, top_logprobs=20,
        extra_body={"guided_json": SCHEMA, "top_k": 20},
    )
    text = resp.choices[0].message.content
    hard = {"disagree": 0, "neutral": 1, "agree": 2}[json.loads(text)["position"]]
    # find the token right after the opening quote of the position value
    toks = resp.choices[0].logprobs.content
    idx = next(i for i, t in enumerate(toks)
               if t.token.strip().strip('"').lower()[:3] in PREFIXES)
    soft = np.zeros(3, dtype=np.float32)
    for alt in toks[idx].top_logprobs:
        key = alt.token.strip().strip('"').lower()[:3]
        if key in PREFIXES:
            soft[PREFIXES[key]] += math.exp(alt.logprob)
    soft /= soft.sum()
    return hard, soft

One question, and what comes back.

topic = "AI automation will lead to mass unemployment"
hard, soft = ask(topic,
                 speaker=persona(profiles[pool_idx[812]]),
                 listener=persona(profiles[pool_idx[41]]),
                 si=2, st=1)

print(f"hard label : {['disagree','neutral','agree'][hard]}")
print(f"soft label : disagree {soft[0]:.3f}  neutral {soft[1]:.3f}  agree {soft[2]:.3f}")
#### OUTPUT ####
hard label : neutral
soft label : disagree 0.041  neutral 0.532  agree 0.427

The soft label says this listener was almost persuaded, a coin flip from agreeing. The word alone would teach the surrogate that the pairing is stable. The distribution teaches it that the pairing is on a knife edge.

Before trusting that distribution, check it. If the logprob triple is well calibrated, sampling the same question repeatedly should reproduce it.

def audit(n_questions=10_000, repeats=8):
    devs = []
    for q in rng.choice(len(design), n_questions, replace=False):
        _, soft = ask(topic, *render_pair(design[q]))
        draws = [ask(topic, *render_pair(design[q]))[0] for _ in range(repeats)]
        devs.append(np.abs(np.bincount(draws, minlength=3) / repeats - soft).mean())
    return float(np.mean(devs))

mad = audit()
print(f"mean absolute deviation, logprob vs 8-sample frequency: {mad:.4f}")
print("threshold for trusting the logprob triple: 0.0300")
#### OUTPUT ####
mean absolute deviation, logprob vs 8-sample frequency: 0.0190
threshold for trusting the logprob triple: 0.0300

0.019 against a threshold of 0.03, so one call is a fair stand-in for eight. That check saved seven eighths of the labelling budget.

Now the whole job.

# ... continued

def run_job(design, topic, workers=64):
    order = sort_for_prefix_reuse(design)
    hard = np.empty(len(design), dtype=np.int8)
    soft = np.empty((len(design), 3), dtype=np.float16)
    with ThreadPoolExecutor(max_workers=workers) as pool:
        for k, (h, s) in enumerate(pool.map(lambda q: ask_with_retry(topic, q),
                                            design[order])):
            hard[order[k]], soft[order[k]] = h, s
    return hard, soft
#### OUTPUT ####
[herald] 270,000 / 270,000 100.0% |  9,420 tok/s  |  elapsed 2h 59m
[herald] schema violations: 0
[herald] empty completions retried: 312 (0.12%)
[herald] label mix: disagree 31.4%  neutral 38.9%  agree 29.7%
[herald] churn rate (final stance differs from initial): 26.3%

Three hours, zero schema violations, and a churn rate of 26.3 percent. That last number is the most important statistic about this teacher, because the surrogate has to match it.

When a matrix is small enough to read, print the numbers rather than making the reader guess at a colour.

fig, axes = plt.subplots(1, 3, figsize=(14.2, 4.4))
cmap = plt.matplotlib.colors.LinearSegmentedColormap.from_list("t", ["#FFFFFF", AMBER])

for ax, start in zip(axes, range(3)):
    m = T3[:, start, :]                      # rows: speaker stance, cols: outcome
    ax.imshow(m, cmap=cmap, vmin=0, vmax=m.max())
    for i in range(3):
        for j in range(3):
            ax.text(j, i, f"{m[i, j]:,}", ha="center", va="center", fontsize=11,
                    color=INK if m[i, j] < m.max() * 0.6 else "white")
    ax.set_xticks(range(3), STANCES, fontsize=9.5)
    ax.set_yticks(range(3), STANCES, fontsize=9.5)
    ax.set_title(f"target starts {STANCES[start]}", fontsize=11,
                 fontweight="bold", color=INK, pad=10)

fig.suptitle(f"All {len(design):,} teacher decisions, and the neutral column that swallows everything",
             fontsize=13, fontweight="bold", color=INK, y=1.04)
fig.tight_layout()
fig.savefig("images/p07_teacher_transitions.png", dpi=200, bbox_inches="tight")
None

Two things in there matter. Mass stays in the column the target started in unless the influencer pushes, and the far corners are nearly empty: a listener who disagrees almost never jumps straight to agreeing, they go through neutral. That pattern is not something we built in. It is what the teacher does, and it drives every result later.

Echo, Four and a Half Million Parameters Standing In for 235 Billion

We need something small enough to evaluate nine hundred million times that behaves the same way. The input is fixed and tiny: two profile embeddings we already have, plus two stances as one-hot vectors of length three.

None None 
def make_features(E, rows):
    # embeddings are looked up, never recomputed, so this is cheap
    ei = torch.from_numpy(E[rows["i_prof"]].astype(np.float32))
    et = torch.from_numpy(E[rows["t_prof"]].astype(np.float32))
    hot = lambda k: torch.nn.functional.one_hot(
        torch.from_numpy(rows[k].astype(np.int64)), 3).float()
    return torch.cat([ei, et, hot("i_stance"), hot("t_stance")], dim=1)

print(f"feature shape: {tuple(make_features(E, design[:4]).shape)}")
#### OUTPUT ####
feature shape: (4, 4102)

The network on top of that is three linear layers.

None None 
class Echo(nn.Module):
    # small on purpose: 900 million evaluations, and it must fit in cache
    def __init__(self, d_in=4102, h1=1024, h2=256, n_classes=3, p_drop=0.1):
        super().__init__()
        self.fc1, self.ln1 = nn.Linear(d_in, h1), nn.LayerNorm(h1)
        self.fc2, self.ln2 = nn.Linear(h1, h2), nn.LayerNorm(h2)
        self.head = nn.Linear(h2, n_classes)
        self.drop, self.act = nn.Dropout(p_drop), nn.GELU()

    def forward(self, x):
        x = self.drop(self.act(self.ln1(self.fc1(x))))
        x = self.drop(self.act(self.ln2(self.fc2(x))))
        return self.head(x)

echo = Echo()
n_params = sum(p.numel() for p in echo.parameters())
print(f"parameters : {n_params:,}  ({n_params*4/1e6:.2f} MB in fp32)")
print(f"ratio      : {235_000_000_000/n_params:,.0f}x smaller")
#### OUTPUT ####
parameters : 4,467,203  (17.87 MB in fp32)
ratio      : 52,606x smaller

Four and a half million parameters standing in for two hundred and thirty five billion. Fifty two thousand times smaller, seventeen megabytes on disk.

Echo answers one question over a fixed population of ten thousand people. We are not compressing the model. We are compressing one narrow behaviour of the model, and narrow behaviours are small.

Training It on the Hesitation

The loss is where the soft labels earn their keep.

None 
def distillation_loss(logits, hard, soft, alpha=0.3):
    # at alpha = 1 this is ordinary classification, and the copy turns
    # overconfident in exactly the cases the teacher was not
    ce = torch.nn.functional.cross_entropy(logits, hard)
    kl = torch.nn.functional.kl_div(torch.nn.functional.log_softmax(logits, dim=-1),
                                    soft, reduction="batchmean")
    return alpha * ce + (1 - alpha) * kl
# ... continued

def train(E, rows, hard, soft, epochs=30, batch=4096, lr=3e-4, seed=0):
    torch.manual_seed(seed)
    model = Echo().cuda()
    opt = torch.optim.AdamW(model.parameters(), lr=lr, weight_decay=0.01)
    sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=epochs)
    idx = np.random.default_rng(seed).permutation(len(rows))
    tr, va, te = idx[:240_000], idx[240_000:255_000], idx[255_000:]
    for ep in range(1, epochs + 1):
        model.train()
        for start in range(0, len(tr), batch):
            b = tr[start:start + batch]
            logits = model(make_features(E, rows[b]).cuda())
            loss = distillation_loss(logits,
                torch.from_numpy(hard[b].astype(np.int64)).cuda(),
                torch.from_numpy(soft[b].astype(np.float32)).cuda())
            opt.zero_grad(); loss.backward(); opt.step()
        sched.step()
        if ep % 3 == 0:
            torch.save(model.state_dict(), f"ckpt/echo_ep{ep:02d}.pt")
    return model, (tr, va, te)

The one that matters is the churn gap.

None 
# the fidelity metrics, used from here to the end

def churn_rate(pred, initial):
    # a property of the whole distribution of outcomes, which is why
    # per-sample accuracy cannot see it
    return float((pred != initial).mean())

def evaluate(model, E, rows, hard, soft, split):
    model.eval()
    with torch.no_grad():
        q = torch.softmax(model(make_features(E, rows[split]).cuda()),
                          dim=-1).cpu().numpy()
    pred = q.argmax(1)
    return dict(f1=macro_f1(pred, hard[split]),
                gap=abs(churn_rate(pred, rows[split]["t_stance"])
                        - churn_rate(hard[split], rows[split]["t_stance"])) * 100,
                tvd=per_cell_tvd(q, hard[split], rows[split]),
                ece=expected_calibration_error(q, hard[split], bins=15))
#### OUTPUT ####
[echo] epoch  3  train 0.681  val 0.694  F1 0.7412  gap 9.81  tvd 0.0614  ece 0.0721
[echo] epoch 18  train 0.518  val 0.569  F1 0.7724  gap 8.31  tvd 0.0472  ece 0.0413
[echo] epoch 24  train 0.503  val 0.574  F1 0.7698  gap 3.18  tvd 0.0298  ece 0.0219
[echo] epoch 30  train 0.495  val 0.586  F1 0.7665  gap 5.11  tvd 0.0367  ece 0.0344
[echo] total training time: 58s on one H100

F1 peaks at epoch 18 with 0.7724. The churn gap there is 8.31, the third worst number in the entire run. Epoch 24 gives up 0.0026 of F1 and cuts the gap to 3.18.

None None

A checkpoint can be better at guessing individual answers and worse at reproducing how much the population moves, and the second is what ships.

Choosing the Checkpoint That Ships

Shortlist by accuracy so we do not ship something broken, then choose inside the shortlist on what we care about.

None None 
def composite(m):
    # lower is better, and every term is a property of the output
    # distribution rather than of individual predictions
    return m["tvd"] + abs(m["gap"]) / 100.0 + 0.5 * m["ece"]

def choose(metrics, top_n=5):
    shortlist = sorted(metrics, key=lambda e: -metrics[e]["f1"])[:top_n]
    best = min(shortlist, key=lambda e: composite(metrics[e]))
    for e in shortlist:
        m = metrics[e]
        print(f"  epoch {e:>2}  F1 {m['f1']:.4f}  gap {m['gap']:5.2f}  "
              f"J {composite(m):.4f}" + ("  <- ships" if e == best else ""))
    return best

best = choose(metrics)
f1_best = max(metrics, key=lambda e: metrics[e]["f1"])
print(f"F1 given up : {metrics[f1_best]['f1'] - metrics[best]['f1']:.4f}")
print(f"gap improved: {metrics[f1_best]['gap']/metrics[best]['gap']:.1f}x")
#### OUTPUT ####
  epoch 18  F1 0.7724  gap  8.31  J 0.1510
  epoch 21  F1 0.7710  gap  6.47  J 0.1286
  epoch 24  F1 0.7698  gap  3.18  J 0.0726  <- ships
  epoch 27  F1 0.7681  gap  4.06  J 0.0874
  epoch 15  F1 0.7679  gap  5.83  J 0.1213
F1 given up : 0.0026
gap improved: 2.6x

We give up 0.0026 of accuracy and get a copy of the teacher 2.6 times closer on the statistic the simulation actually depends on. Select on the aggregate you care about, not the metric printed by default.

Five Surrogates, One Question

Does the architecture matter, or would anything do? We train five on the same labels, scored the same way.

# the same loop, pointed at five architectures instead of one

ARCHS = {
    "logistic regression": LogisticOnAttributes(),   # raw survey codes, no embedding
    "gradient boosting":   GradientBoostOnAttributes(),
    "Echo MLP":            Echo(),
    "Transformer":         TinyEncoder(d=2048, layers=2, heads=8),
    "Qwen3-0.6B tuned":    GenerativeSurrogate("Qwen/Qwen3-0.6B"),
}

for name, m in ARCHS.items():
    r = train_and_score(m, E, design, hard, soft)
    print(f"{name:<22} F1 {r['f1']:.4f}   gap {r['gap']:5.2f}   "
          f"dead channels {r['dead']:>2}/27")
#### OUTPUT ####
logistic regression    F1 0.7461   gap 18.62   dead channels  9/27
gradient boosting      F1 0.7688   gap 11.04   dead channels  5/27
Echo MLP               F1 0.7724   gap  8.31   dead channels  0/27
Transformer            F1 0.7702   gap 10.47   dead channels  4/27
Qwen3-0.6B tuned       F1 0.7715   gap  9.55   dead channels  2/27
METRICS = [("f1", "macro F1 (higher is better)", "{:.4f}", (0.735, 0.782)),
           ("gap", "churn gap in points (lower is better)", "{:.2f}", None),
           ("dead", "annihilated channels of 27 (lower is better)", "{:.0f}", None)]

fig, axes = plt.subplots(1, 3, figsize=(16.4, 5.4))
for ax, (key, title, fmt, fixed) in zip(axes, METRICS):
    vals = [results[n][key] for n in ARCHS]
    lo, hi = fixed or (0, max(vals) * 1.30)
    ax.set_ylim(lo, hi)                       # set the limits BEFORE placing labels
    ax.bar(range(len(vals)), vals, color=[COLOUR[n] for n in ARCHS], width=0.62)
    for x, v in enumerate(vals):
        ax.text(x, v + (hi - lo) * 0.028, fmt.format(v), ha="center", color=INK)
    style(ax, title)

fig.savefig("images/p11_bakeoff.png", dpi=200, bbox_inches="tight")
None

The F1 column spans 0.026. The gap column spans 10.3 points. Picking on accuracy means choosing between models that are essentially tied while ignoring a fourfold difference in what decides your results.

Dead channels counts how many of the twenty seven possible stance transitions the surrogate assigns essentially zero probability to when the teacher does not. Logistic regression kills nine.

Atlas, Freezing a Model Into an Array

We enumerate every possible question, answered in advance.

None

Nine hundred million forward passes done naively is worse than it sounds, because it builds a 4102-dimensional vector nine hundred million times.

print(f"materialising every feature vector: {900_000_000 * 4102 * 4 / 1e12:.2f} TB")
#### OUTPUT ####
materialising every feature vector: 14.77 TB

Fourteen terabytes of intermediate, so we do not do that. The first layer is linear and our input is a concatenation, which makes it separable.

None

Split the first weight matrix across the speaker embedding, the listener embedding, and the stances, and layer one for any question is the sum of three precomputed pieces.

def cast(model, E, codebook, out_path, tile_i=256, tile_t=2048):
    # fc1 is separable across the concatenation, so each half is computed
    # once for all 10,000 profiles instead of once per pair
    W1, b1 = model.fc1.weight.data, model.fc1.bias.data
    Ea = torch.from_numpy(E.astype(np.float32)).cuda()
    U = Ea @ W1[:, :2048].T                     # [10000, 1024] speaker halves
    V = Ea @ W1[:, 2048:4096].T                 # listener halves
    onehot9 = torch.cat([torch.eye(3).repeat_interleave(3, 0),
                         torch.eye(3).repeat(3, 1)], 1).cuda()
    G = onehot9 @ W1[:, 4096:].T + b1           # stance halves
    T = np.lib.format.open_memmap(out_path, mode="w+", dtype=np.uint8,
                                  shape=(10_000, 10_000, 3, 3))

    for i0 in range(0, 10_000, tile_i):
        for t0 in range(0, 10_000, tile_t):
            for k in range(9):
                h = model.act(model.ln1(
                    U[i0:i0+tile_i, None, :] + V[None, t0:t0+tile_t, :] + G[k]))
                h = model.act(model.ln2(model.fc2(h)))
                T[i0:i0+tile_i, t0:t0+tile_t, k // 3, k % 3] = codebook.encode(
                    torch.softmax(model.head(h), dim=-1)).cpu().numpy()
    T.flush()
    return T

The saving is the difference between three minutes and an afternoon.

naive_flops = 900_000_000 * 4102 * 1024 * 2
smart_flops = (2 * 10_000 * 2048 * 1024 * 2) + (900_000_000 * 1024 * 2)

print(f"saving               : {naive_flops/smart_flops:.0f}x")
#### OUTPUT ####
saving               : 3923x

Layer one gets almost four thousand times cheaper, because it stops being nine hundred million matrix multiplies and becomes twenty thousand plus a very large number of additions.

Argmax or Sample, and the Eighteen Channels We Nearly Deleted

This is the decision I got wrong first.

The obvious way to store an answer in one byte is the winning class. Three classes fits in two bits, so the whole table packs into 225 megabytes.

None

If the teacher leaves a listener disagreeing 97 percent of the time and neutral 3 percent, argmax records disagree, full stop.

teacher = teacher_transition_probs()          # [3,3,3], from the labelled set

argmax_alive = codebook_alive = 0
for si in range(3):
    for st in range(3):
        for sf in range(3):
            # anything the teacher gives real mass to is a channel the
            # simulation can travel, however rarely
            if teacher[si, st, sf] < 1e-4:
                continue
            # an argmax table stores one outcome per row, so a row can only
            # ever emit its own winner
            argmax_alive += int(teacher[si, st].argmax() == sf)
            codebook_alive += 1

frozen = [STANCES[st] for st in range(3)
          if all(teacher[si, st].argmax() == st for si in range(3))]

print(f"transitions the teacher actually uses : {codebook_alive}/27")
print(f"transitions an argmax table can produce: {argmax_alive}/27")
print(f"channels destroyed                     : {codebook_alive-argmax_alive}")
print(f"stances that become absorbing          : {frozen}")
#### OUTPUT ####
transitions the teacher actually uses : 27/27
transitions an argmax table can produce: 9/27
channels destroyed                     : 18
stances that become absorbing          : ['agree', 'neutral', 'disagree']

Eighteen of twenty seven, and it is worse than the count makes it sound. A stay rate of 73.7 percent means the most likely outcome of every single row is that the listener does not move. An argmax table stores exactly that, so all three stances become absorbing at once and the population freezes on the first round. The 26.3 percent of interactions that actually change somebody live entirely in the mass that argmax throws away.

The fix costs nothing. Instead of which class won, we store which distribution it was, using a codebook.

None None 
class Codebook:
    # 256 points on the 2-simplex, one byte each, and the byte names a
    # distribution rather than a decision
    def __init__(self, centroids):
        self.C = centroids.astype(np.float16)
        self.Ccum = np.cumsum(centroids, axis=1)[:, :2].astype(np.float16)

    @classmethod
    def fit(cls, samples, n=256, n_lattice=32, seed=0):
        # 32 centroids are frozen on a lattice because plain k-means chases
        # mass, and the corners hold little mass and a lot of consequence
        lattice = simplex_lattice(n_lattice)
        fitted = kmeans(samples, k=n - n_lattice, seed=seed)
        return cls(np.vstack([lattice, fitted]))

    def encode(self, q):
        C = torch.from_numpy(self.C.astype(np.float32)).to(q.device)
        # L1 on the simplex
        return (q.unsqueeze(-2) - C).abs().sum(-1).argmin(-1).to(torch.uint8)

    def sample(self, codes, u):
        cum = self.Ccum[codes].astype(np.float32)
        return ((u >= cum[:, 0]).astype(np.int8) +
                (u >= cum[:, 1]).astype(np.int8))
posteriors = sample_posteriors(echo, E, n=10_000_000, seed=5)
cb = Codebook.fit(posteriors)

err = np.abs(posteriors - cb.C[cb.encode(posteriors)].astype(np.float32)).sum(1) / 2
print(f"codebook bytes   : {cb.C.nbytes + cb.Ccum.nbytes:,}")
print(f"quantisation TVD : mean {err.mean():.4f}   p99 {np.quantile(err,0.99):.4f}")
print(f"dead channels    : {count_dead_channels(cb, echo)}/27")
#### OUTPUT ####
codebook bytes   : 2,560
quantisation TVD : mean 0.0041   p99 0.0163
dead channels    : 0/27
None

Mean error of four thousandths, a codebook that fits in two and a half kilobytes, and every transition channel intact. Same one byte per entry as the argmax table. Only the meaning of the byte changed.

Now run it over all nine hundred million.

T = cast(echo, E, cb, "artifacts/atlas.npy")

print(f"entries: {T.size:,}")
print(f"bytes  : {T.nbytes:,}  ({T.nbytes/1048576:.2f} MiB)")
#### OUTPUT ####
[atlas] tile 40/40 (i0=9984)  |  elapsed 192s
entries: 900,000,000
bytes  : 900,000,000  (858.31 MiB)
None

Nine hundred million answers in 858 mebibytes, built in three minutes and twelve seconds. The teacher is now a file.

The Lookup

At runtime this is four index arithmetic operations, one memory read, one codebook read, one comparison.

None 
class Atlas:
    def __init__(self, path, codebook):
        self.T = np.load(path, mmap_mode="r").reshape(-1)   # flat for speed
        self.cb = codebook

    @staticmethod
    def key(pi, pt, si, st):
        # int64 because the product overflows 32 bits
        return ((pi.astype(np.int64) * 10_000 + pt) * 3 + si) * 3 + st

    def lookup(self, keys, u):
        codes = self.T[keys]                # THE only policy call in the system
        return self.cb.sample(codes, u)

atlas = Atlas("artifacts/atlas.npy", cb)

k = atlas.key(np.array([812]), np.array([41]), np.array([2]), np.array([1]))
code = atlas.T[k]
print(f"distribution: {cb.C[code].astype(np.float32)[0].round(3)}")
print(f"teacher said: [0.041 0.532 0.427]")
#### OUTPUT ####
distribution: [0.043 0.529 0.428]
teacher said: [0.041 0.532 0.427]

That is the whole system. The question we asked the 235-billion-parameter teacher now costs one array index, and agrees to three decimal places.

One Round of a Billion People

A round is: pick who speaks, find who they reach, look up what happens, draw an outcome, write it back, count.

None 
def one_round(led, indptr, indices, atlas, tally, round_id, run_seed, cut=200_000_000, active=0.01):
    rng = np.random.Generator(np.random.Philox(key=run_seed, counter=round_id))
    n_active = int(active * cut)
    speakers = np.unique(rng.integers(0, cut, size=int(1.01 * n_active), dtype=np.uint32))[:n_active]
    counts = (indptr[speakers + 1] - indptr[speakers]).astype(np.int64)
    src = np.repeat(speakers, counts)
    targets = ragged_gather(indices, indptr[speakers], counts)
    keys = atlas.key(led.profile_idx[src], led.profile_idx[targets], led.stance[src], led.stance[targets])
    u = rng.random(targets.size, dtype=np.float32)
    new = atlas.lookup(keys, u)
    # sorted so duplicate writes always resolve the same way
    order = np.argsort(targets, kind="stable")
    uniq, first = np.unique(targets[order], return_index=True)
    last = np.append(first[1:], targets.size) - 1
    winners = order[last]
    old = led.stance[uniq]
    led.stance[uniq] = new[winners]
    tally.update(old, new[winners], led.stance[src[winners]])

    return {"speakers": speakers.size, "interactions": targets.size, "reached": uniq.size, "changed": int((old != new[winners]).sum())}

Run exactly one.

stats = one_round(led, indptr, indices, atlas, tally, round_id=0, run_seed=20260812)
for k, v in stats.items():
    print(f"{k:<14}: {v:>12,}")
print(f"{'change rate':<14}: {stats['changed']/stats['reached']*100:>11.2f}%")
#### OUTPUT ####
[pulse] round   0 | sample 0.02s gather 0.18s key 0.12s lookup 0.25s
                    draw 0.10s sort 0.14s scatter 0.20s tally 0.03s | 1.04s
speakers      :    2,000,000
interactions  :   14,930,000
reached       :   14,317,884
changed       :    3,761,504
change rate   :       26.27%

Fourteen point nine million interactions in 1.04 seconds, on CPU, no GPU. And the change rate is 26.27 percent, against the teacher's measured 26.3, on the first round the surrogate runs.

The counter-based generator means round seventeen draws the same numbers whether you ran the first sixteen or jumped straight to it, and whether the machine had four cores or sixty four.

The sort before the scatter handles a listener hearing from two people in one round. Without it, numpy fancy-index assignment resolves duplicates in memory order and the run stops being deterministic. With it, the highest-index writer wins.

Running All One Hundred Rounds

def simulate(cfg):
    led = Ledger(cfg.n_agents)
    led.seed_stances(cut=cfg.cut, scheme=cfg.seeding)
    atlas = Atlas(cfg.table_path, Codebook.load(cfg.codebook_path))
    indptr, indices = load_csr(cfg.graph_path)
    chron = Chronicle(rounds=cfg.rounds)

    q = Tempo()
    for r in range(cfg.rounds):
        q.push(time=r, kind="influence", payload=r)
        q.push(time=r, kind="readout", payload=r, priority=9)

    while len(q):
        t, kind, _ = q.pop_batch()
        if kind == "influence":
            one_round(led, indptr, indices, atlas, chron.tally, round_id=t, run_seed=cfg.seed, cut=cfg.cut)
        else:
            chron.record(t, led.stance[cfg.cut:])
    return chron

chron = simulate(load_config("configs/billion.yaml"))
#### OUTPUT ####
[sim] population 1,000,000,000 | influencers 200,000,000 | seeding HD
[sim] topic: "AI automation will lead to mass unemployment"
[sim] round   0 | agree 33.333% neutral 33.333% disagree 33.333% | 1.04s
[sim] round 100 | agree 23.473% neutral 43.214% disagree 33.313% | 1.04s
[sim] 100 rounds in 104.2s  (1.042 s/round)
[sim] 1,493,000,000 interactions resolved, 0 model calls
[sim] peak resident memory: 11.4 GB

One hundred rounds over eight hundred million listeners in one minute forty four seconds, and not a single model call.

frame = chron.frame()                 # [round, disagree%, neutral%, agree%]
rounds = frame[:, 0]
INFLUENCEES = 800_000_000

fig, ax = plt.subplots(figsize=(9.4, 4.8))
for col, colour, label in ((3, AGREE, "agree"), (2, NEUTRAL, "neutral"), (1, DISAGREE, "disagree")):
    share = frame[:, col]
    millions = share * INFLUENCEES / 100 / 1e6
    ax.plot(rounds, millions, color=colour, linewidth=2.8, label=label)
    ax.annotate(f"{share[-1]:.2f}%", (rounds[-1], millions[-1]), color=colour)

ax.legend(frameon=False, ncol=3)
style(ax, "One hundred rounds over eight hundred million people, in 1 minute 44 seconds", "round", "influencees (millions)")
fig.savefig("images/p13_trajectory_1b.png", dpi=200, bbox_inches="tight")
None

Read the numbers rather than the curve. Agree fell 9.860 points, neutral rose 9.881, and disagree moved 0.020.

We seeded half the influencers to disagree. Ten percent of the population left agree over a hundred rounds. Almost none arrived at disagree. They all stopped at neutral.

That is not a bug. It is the transition structure the teacher gave us, at population scale.

What the Seeding Does

What if we seed the loud fifth the other way, or not at all?

None 
for scheme in ("HA", "HD", "UN"):
    c = simulate(load_config("configs/billion.yaml", seeding=scheme))
    a, n, d = c.final_shares()
    print(f"{scheme}: agree {a:6.2f}%  neutral {n:6.2f}%  disagree {d:6.2f}%   (started 33.33 / 33.33 / 33.33)")
#### OUTPUT ####
HA: agree  42.65%  neutral  36.20%  disagree  21.15%   (started 33.33 / 33.33 / 33.33)
HD: agree  23.47%  neutral  43.21%  disagree  33.31%   (started 33.33 / 33.33 / 33.33)
UN: agree  36.55%  neutral  37.51%  disagree  25.94%   (started 33.33 / 33.33 / 33.33)
None

Seeding the loud fifth to agree moves 9.3 points into agree. Seeding them to disagree moves nothing into disagree at all, it empties agree into neutral instead. Even the uniform case drifts toward agree.

The population has a prior, and on this topic it leans toward agreeing. Pushing with the prior converts people. Pushing against it only makes them undecided.

Four Topics, Four Destinations

If that is a real effect and not an artifact, it should change with the topic. We rebuild the table for three more statements and rerun.

TOPICS = [
    "AI automation will lead to mass unemployment",
    "The Earth is flat",
    "Humans will establish a city on Mars within fifty years",
    "Short-form video is shortening human attention spans",
]

for topic in TOPICS:
    c = simulate(load_config("configs/billion.yaml", topic=topic, seeding="HA"))
    a, n, d = c.final_shares()
    print(f"{topic[:44]:<46} {a:6.2f} {n:6.2f} {d:6.2f}")
#### OUTPUT ####
AI automation will lead to mass unemployment    42.65  36.20  21.15
The Earth is flat                              18.42  33.71  47.87
Humans will establish a city on Mars within f   30.94  44.83  24.23
Short-form video is shortening human attentio   46.11  33.02  20.87
None

Identical machinery, identical seeding, opposite outcomes. The flat-Earth population moves to disagree even though every influencer was seeded to agree or stay neutral, because the prior against it is overwhelming. Mars, where nobody knows, drifts neutral. Short-form video, which most people already suspect, goes furthest toward agree.

Influence is not a force applied to a blank population. It is a force applied against a prior, and the prior wins more often than not.

Why Nobody Flips Straight to the Other Side

Every result so far has neutral in the middle of it.

None 
tr = chron.tally.matrix()        # [from, to], summed over the whole run
changes = tr.sum() - np.trace(tr)
cross = tr[0, 2] + tr[2, 0]
print(f"total stance changes      : {changes:>14,.0f}")
print(f"straight across the middle: {cross:>14,.0f}   ({cross/changes*100:.2f}%)")
print(f"through neutral           : {changes-cross:>14,.0f}   ({(changes-cross)/changes*100:.2f}%)")
#### OUTPUT ####
total stance changes      :    412,984,117
straight across the middle:      1,280,251   (0.31%)
through neutral           :    411,703,866   (99.69%)

Three tenths of one percent go straight from one pole to the other. Everything else stops in the middle first.

This is why the argmax decision mattered. Neutral is less a stance than a corridor, and a format that closes off the low-probability entrances does not slightly degrade the dynamics. It removes the only route people take.

Who Persuades and Who Bends

The table maps every pair of people to an outcome, so marginalising it against the run tells us who changes minds and who is easy to change.

# marginalise the table against what the run actually did

def persuasion_and_resistance(atlas, cb, led, observed_pairs):
    persuade, suscept = np.zeros(10_000), np.zeros(10_000)
    for pi, pt, si, st, w in observed_pairs:
        q = cb.C[atlas.T[atlas.key(pi, pt, si, st)]].astype(np.float32)
        p_change = 1.0 - q[st]
        persuade[pi] += w * p_change
        suscept[pt] += w * p_change
    return persuade / weights_by_speaker, suscept / weights_by_listener
persuade, suscept = persuasion_and_resistance(atlas, cb, led, chron.observed_pairs())

print(f"{'education':<11}{'persuasion':>12}{'susceptibility':>16}")
for level in range(0, 9, 2):
    m = [i for i in range(10_000) if profiles[pool_idx[i]].get("isced", "").endswith(str(level))]
    print(f"ISCED {level:<6}{persuade[m].mean()*100:>10.1f}%{suscept[m].mean()*100:>15.1f}%")
print(f"{'spread':<11}{PERSUASION[-1]/PERSUASION[0]:>10.2f}x{SUSCEPTIBILITY[0]/SUSCEPTIBILITY[-1]:>15.2f}x")
#### OUTPUT ####
education   persuasion   susceptibility
ISCED 0          8.2%           17.1%
ISCED 4         11.9%           10.1%
ISCED 8         13.7%            2.4%
spread          1.67x            7.13x
None

Both go the way you would guess, but not by the same amount, and the difference is the interesting part. Lowest education to highest multiplies persuasive power by 1.67 and divides susceptibility by 7.13.

Education buys about four times more resistance than it buys reach. A society modelled this way fails to converge not because educated people talk more persuasively, but because they stop listening.

Income adds a second axis.

fig, ax = plt.subplots(figsize=(9.6, 5.4))
im = ax.imshow(grid, cmap="RdYlBu_r", aspect="auto")

for i in range(9):
    for j in range(10):
        ax.text(j, i, f"{grid[i, j]:.1f}", ha="center", va="center", color=INK)

ax.set_yticks(range(9))
ax.set_yticklabels([f"ISCED {i}" for i in range(9)])
ax.set_xlabel("household income decile", color=MUTE)
fig.savefig("images/p21_edu_income.png", dpi=200, bbox_inches="tight")

# read the corners and both gradients off the grid rather than off the picture
print(f"ISCED 0, decile  1 : {grid[0, 0]:>5.1f}%")
print(f"ISCED 8, decile 10 : {grid[-1, -1]:>5.1f}%")
print(f"education spread {grid.mean(1)[-1] - grid.mean(1)[0]:.1f} points, "
      f"income spread {grid.mean(0)[-1] - grid.mean(0)[0]:.1f} points")
#### OUTPUT ####
ISCED 0, decile  1 :   7.4%
ISCED 8, decile 10 :  16.2%
education spread 5.1 points, income spread 1.8 points
None

The gradient runs corner to corner, from 7.4 percent at the top left to 16.2 at the bottom right. But read across a row and then down a column: the education gradient is about three times steeper than the income gradient.

Does the Language Change the Society

Everything so far ran in English, and the whole system sits downstream of one model reading prose. We rebuild the teacher labels in Spanish and Chinese, holding personas, stance pairs and prompt structure identical, and measure churn.

for topic in ("The Earth is flat", "A city on Mars", "Short-form video"):
    rates = {}
    for lang in ("en", "es", "zh"):
        hard, _ = run_job(design, topic, language=lang)
        rates[lang] = churn_rate(hard, design["t_stance"]) * 100
    spread = max(rates.values()) - min(rates.values())
    print(f"{topic:<24} {rates['en']:>8.1f} {rates['es']:>8.1f} {rates['zh']:>8.1f} {spread:>8.1f}")
#### OUTPUT ####
topic                     English  Spanish  Chinese   spread
The Earth is flat            45.9     47.2     44.1     3.1
A city on Mars               22.8     21.9     23.4     1.5
Short-form video             29.4     25.1     30.6     5.5
None

Up to 5.5 points of difference, and the sign flips between topics. Spanish beats Chinese on the flat-Earth statement and loses on short-form video.

A quantitative claim from this system is a claim about a model reading one language about one topic. Not a claim about people. Anyone using this for research should run both and report both.

How Much of the Teacher Survived the Compression

Take a population small enough to run with the live model in the loop, rerun it with the table, and compare.

for ratio in (0.0, 0.25, 0.5, 0.75, 1.0):
    c = simulate(load_config("configs/small.yaml", n_agents=20_000, policy="weighted", live_ratio=ratio))
    print(f"live ratio {ratio:>4}  ->  terminal MAD vs full-model run: {c.mad_against_reference():.4f}")
#### OUTPUT ####
live ratio  0.0  ->  terminal MAD vs full-model run: 0.0139
live ratio 0.25  ->  terminal MAD vs full-model run: 0.0106
live ratio  0.5  ->  terminal MAD vs full-model run: 0.0074
live ratio 0.75  ->  terminal MAD vs full-model run: 0.0038
live ratio  1.0  ->  terminal MAD vs full-model run: 0.0000
None

The curves keep their shape at every substitution level and are strictly ordered, so the difference is systematic bias, not noise. How big:

ref = simulate(load_config("configs/small.yaml", policy="all-live"))
for name, cfg in (("codebook table", "sampled"), ("argmax table", "argmax")):
    run = simulate(load_config("configs/small.yaml", table_mode=cfg))
    for stance in ("disagree", "neutral"):
        moved, moved_ref = run.delta(stance), ref.delta(stance)
        print(f"{name:<16}{moved:>+8.0f} agents  ({moved/moved_ref*100:.1f}%)")
#### OUTPUT ####
                       disagree decline   neutral rise
full model                  -548 agents     +542 agents
codebook table              -519 (94.7%)    +516 (95.3%)
argmax table                -446 (81.4%)    +451 (83.2%)

The codebook table reproduces about 95 percent of the teacher dynamics. The argmax table reproduces about 82 percent. Thirteen points of fidelity, free, from deciding what a byte means.

None

The errors are small everywhere and not uniformly signed. The surrogate leans slightly toward neutral, the familiar pull any distilled classifier feels toward its majority class. About two points per cell is the price of the compression, which is the shipped checkpoint's mean per-cell total variation of 0.0298.

Every agree curve also dips to a minimum around round 37 and then recovers. It appears even in the pure-model run, so it belongs to the dynamics rather than the compression: agents retreat under pressure, then return once enough neighbours have moved.

Over a long horizon the gap between the two storage choices stops being a percentage and becomes a difference in kind.

None 
for mode in ("codebook", "argmax"):
    run = simulate(load_config("configs/long.yaml", table_mode=mode, rounds=5000))
    ever = run.ever_changed[run.initial_stance == 0].mean()
    print(f"{mode:<12} P(ever changes) = {ever:.4f}")
#### OUTPUT ####
over 5,000 rounds, P(an agent that started disagreeing ever changes)
  full model, extrapolated   0.9951
  codebook table             0.9938
  argmax table               0.0000

Zero. Not small, zero. The argmax table produces a society in which one third of the population is frozen from the first round to the last.

Running It Five Times, and What Low Variance Does Not Prove

A simulation you cannot reproduce is a rumour. Five seeds, everything else held.

None 
runs = [simulate(load_config("configs/billion.yaml", seed=seed)) for seed in range(5)]
for seed, c in enumerate(runs):
    print(f"seed {seed}: agree {c.final_shares()[0]:.4f}%")
print(f"peak CV across rounds: {peak_cv(runs)['agree']:.4f}%")
#### OUTPUT ####
seed 0: agree 23.4731%
seed 4: agree 23.4733%
peak CV across rounds: 0.0041%

Most write-ups would stop there. I run two controls first.

CONTROLS = [("pure table, 1e9 agents, 5 seeds", dict(n_agents=10**9, policy="all-table")),
            ("pure table, 1e6 agents, 5 seeds", dict(n_agents=10**6, policy="all-table")),
            ("mixed 50% live, 1e6, 5 runs",     dict(n_agents=10**6, policy="weighted", live_ratio=0.5))]
for label, kw in CONTROLS:
    runs = [simulate(load_config("configs/billion.yaml", seed=s, **kw)) for s in range(5)]
    cv = peak_cv(runs)
    print(f"{label:<34}{cv['agree']:>7.4f}%{cv['neutral']:>12.4f}%{cv['disagree']:>12.4f}%")
#### OUTPUT ####
condition                          agree      neutral    disagree
pure table, 1e9 agents, 5 seeds   0.0041%     0.0057%     0.0036%
pure table, 1e6 agents, 5 seeds   0.0129%     0.0183%     0.0121%
mixed 50% live, 1e6, 5 runs       0.1512%     0.3218%     0.1744%

Hold scale at a million and change only whether a live model is in the loop, otherwise the third row is uninterpretable. Live inference adds roughly 15 to 20 times the variance of a frozen table at matched scale.

None None

Seed-to-seed variation peaks at 0.0057 percent. The churn gap against the teacher is 3.18 percent, a factor of about 550.

So the reproducibility number is not evidence that the answer is right. It is evidence that we froze the policy. The worry is that every run agrees and every run is off the teacher by the same three points.

Was the Billion Necessary

If the answer at a million matches the answer at a billion, the billion is a demonstration of machinery, not a scientific requirement.

for n in (10**5, 10**6, 10**7, 10**8, 10**9):
    runs = [simulate(load_config("configs/billion.yaml", n_agents=n, seed=s)) for s in range(5)]
    a = np.array([r.final_shares()[0] for r in runs]) / 100
    print(f"N = {n:>13,}  agree {a.mean():.4f}  +/- {a.std():.5f}")
#### OUTPUT ####
N =       100,000  agree 0.2361  +/- 0.00210
N =     1,000,000  agree 0.2352  +/- 0.00071
N = 1,000,000,000  agree 0.2347  +/- 0.00004
None

Flat to three decimal places above a million agents. For this topic and this network model, a million agents would have given the same scientific answer. The billion buys a fiftyfold tighter run-to-run spread and the ability to resolve rare subpopulations.

Coverage, or Who Never Gets Spoken To

None None 
exposures = (ROUNDS * ACTIVE_FRACTION * INFLUENCERS * OUT_DEGREE) / N_IEE
print(f"expected exposures per listener over 100 rounds: {exposures:.3f}")
print(f"P(never spoken to at all)  : {math.exp(-exposures)*100:.1f}%")
print(f"rounds for 99% coverage    : {math.ceil(math.log(100)/(exposures/ROUNDS))}")
#### OUTPUT ####
expected exposures per listener over 100 rounds: 1.866
P(never spoken to at all)  : 15.5%
rounds for 99% coverage    : 247

Fifteen and a half percent of the population is never contacted at all, and reaching almost everyone takes 247 rounds. A hundred-round run is an early transient, not an equilibrium, so every equilibrium claim here comes from the long runs.

Five Thousand Rounds

So we run those.

#### OUTPUT ####
[sim] long horizon: 4 topics x 3 seedings = 12 configurations, 5,000 rounds
[sim] all 12 configurations in 17m 26s (GPU path: 4m 08s)

99% of total change complete by round:  402 (fastest)  661 (slowest)
first round where every trajectory is flat to 0.01pp: 900

Every configuration settles, none collapses to a single stance, and the flat point is around round 900, nine times longer than the runs we have been reading.

Arena, Do These Agents Behave Like People

None of this means anything unless the agents behave like people, so we put them in economic games with known human results.

None 
TRUSTOR = """{persona}

You are in a one-shot, anonymous economic game with a stranger.

1. You begin with $10 and send an integer amount N between $0 and $10.
2. Whatever you send is tripled to $3N, and they choose how much to send back.

Reply with only this JSON object:
{{"reason": "<one sentence>", "amount": <integer 0 to 10>}}"""

def run_trust(model, personas, n=50_000):
    return np.array([ask_amount(model, TRUSTOR.format(persona=p)) for p in personas[:n]])
#### OUTPUT ####
backbone           sent   returned($9)  returned($21)
Qwen3-235B-A22B     4.10        3.80          9.40
Qwen3-14B           2.90        3.40          8.60
Qwen3-4B            2.30        3.00          7.10

human meta-analysis: trustors send 50% of endowment ($5.00)
                     trustees return 37% of what they received
our 235B agents:     send 41.0%, return 42.2% of a $9 receipt
None

Those two human numbers come from a meta-analysis of 162 trust-game replications.

sends, returns_at_9 = run_trust("Qwen3-235B-A22B", pool_idx, n=50_000)[:2]
by_class = group_mean(sends, key=lambda p: p.get("social_class"))
by_edu = group_mean(returns_at_9, key=lambda p: edu_band(p.get("isced")))
for (cls, s), (band, r) in zip(by_class.items(), by_edu.items()):
    print(f"{cls:<18}{s:.2f}      {band:<18}{r:>9.2f}")
#### OUTPUT ####
social class      sent      education band     returned($9)
Lower             3.40      ISCED 0-2                  3.81
Lower middle      3.90      ISCED 6-8                  4.12
Upper middle      4.40      rural                      3.94
Upper             4.90      urban                      4.06

Trust rises with social class and reciprocity with education, the same direction human studies find. The levels are wrong and the ordering is right, which is the strongest claim this evidence supports.

The Anchor That Moved More Than the Demographics

One control first. Our prompt contains a worked example of the JSON format, and that example contains a number.

None 
for anchor in (None, 2, 5, 8):
    example = "" if anchor is None else f'\n\nExample: {{"reason": "...", "amount": {anchor}}}'
    sends = run_trust(model, personas, prompt_suffix=example)
    print(f"anchor {str(anchor):<5} -> mean sent ${sends.mean():.2f}")
#### OUTPUT ####
anchor None  -> mean sent $3.60
anchor 2     -> mean sent $3.10
anchor 5     -> mean sent $4.10
anchor 8     -> mean sent $4.80
None

The worked example moves the mean by $1.70. The entire social class gradient is $1.50.

Every absolute number above is conditional on one arbitrary choice in a prompt. The orderings survive the anchor. The levels do not.

The Scaling Law, Read Honestly

One more result, where I disagree with the obvious reading of my own chart. We vary the number of agents and measure the young-old gap in sending.

for n in TRUST_SCALE_N:
    young = run_trust(model, sample_by_age(personas, "16-34", n))
    old = run_trust(model, sample_by_age(personas, "55+", n))
    half = 1.96 * math.sqrt(young.var()/len(young) + old.var()/len(old))
    print(f"{n:>6,}          ${young.mean()-old.mean():.2f}             ${half:.2f}")
#### OUTPUT ####
     N   young-old gap   95% CI half-width
    50          $0.71             $0.79
   250          $0.72             $0.33
 5,000          $0.62             $0.07
None

It is tempting to call this a demographic effect that sharpens with scale. It is not. The gap is flat and if anything slightly declining. What narrows is the confidence band, elevenfold, from 0.79 to 0.07.

What grows with scale is our ability to detect the effect, not the effect. At fifty agents the band contains zero. At five thousand it does not.

The Ultimatum Game

Chosen because the rational answer and the human answer are so far apart.

None 
single = run_ultimatum(model, personas, rounds=1, pairs=300)
print(f"single-round reference: mean offer {single.mean_offer:.1f} of 100  (human experiments: 40 to 50)")
print("game theory says      : offer 1, accept anything")
iterated = run_ultimatum(model, personas, rounds=10, pairs=300)
for r in (1, 3, 5, 7, 10):
    o, rej = iterated.offer[r-1], iterated.rejects[r-1]
    print(f"{r:>5}   {o:>10.1f}   {rej:>17}   {rej/300*100:>9.1f}%")
#### OUTPUT ####
single-round reference: mean offer 40.6 of 100  (human experiments: 40 to 50)
game theory says      : offer 1, accept anything

round   mean offer   rejections (of 300)   reject rate
    1        39.8              0                0.0%
    7        28.0             58               19.3%
   10        27.4              7                2.3%

Two quantities on very different scales, so this one needs a second y axis with each axis label coloured to match its line, otherwise nobody can tell which scale is which.

rounds = np.arange(1, 11)
fig, ax = plt.subplots(figsize=(9.4, 4.6))
ax.plot(rounds, iterated.offer, color="#2C6BAA", marker="o", linewidth=2.8)
ax.axhspan(40, 50, color="#12855F", alpha=0.10)          # the human fair band
ax.axhline(1, color=RED, linestyle=":", linewidth=1.6)
ax.set_ylabel("mean offer (of 100)", color="#2C6BAA", fontsize=11)
ax2 = ax.twinx()
ax2.plot(rounds, np.array(iterated.rejects) / 300 * 100, color=AMBER, marker="s", linewidth=2.6)
ax2.set_ylabel("rejection rate (%)", color=AMBER, fontsize=11)
fig.savefig("images/p26_ultimatum.png", dpi=200, bbox_inches="tight")
None

The single-round mean of 40.6 out of 100 sits inside the range the original ultimatum experiments and a fifteen-society replication both found. Over ten rounds proposers get greedier, responders punish them, and offers settle where responders tolerate them.

That is a negotiation. It emerges from personas and a prompt, with no bargaining logic anywhere in our code.

Does the Shape of the Network Change the Society

We take a thousand agents, wire them six ways, and run a richer interaction where agents carry a continuous stance and a confidence.

None

Five are generated: preferential-attachment and random graphs, each sparse and dense, plus a spanning tree. The sixth is empirical, a thousand-node subgraph of a public friendship network from a live-streaming platform, taken by a degree-prioritised search from the highest-degree node. That rule walks toward the dense core, so the sample is denser than its parent and any density conclusion inherits that choice.

for name in ("BA-sparse", "ER-sparse", "Tree", "BA-dense", "ER-dense", "Empirical"):
    g = build_substrate(name, n=1000)
    r = run_topology(g, days=20, model="Qwen/Qwen3-4B-Instruct-2507")
    print(f"{name:<13}{r.spread:>7.3f}{r.support:>10.0f}{r.skeptic:>10.0f}{g.mean_degree:>9.1f}{g.clustering:>13.2f}{r.messages_per_agent:>17.1f}")
#### OUTPUT ####
substrate     spread   support%  skeptic%  mean k   clustering   messages/agent
BA-sparse      0.487        34        38      6.0        0.03              2.4
ER-sparse      0.491        36        37      6.0        0.01              2.4
Tree           0.469        35        34      2.0        0.00              0.8
BA-dense       0.389        15        51     53.2        0.13             21.3
ER-dense       0.384        12        63     54.2        0.05             21.7
Empirical      0.331         6        55     54.3        0.31             21.7

initialisation baseline: 40% support, 40% skeptic, 20% neutral
None

Two clean groups. Sparse graphs keep their spread near 0.49, close to where they started. Dense graphs compress toward a shared, more skeptical position, the empirical graph furthest.

The obvious conclusion is that density drives convergence. But look at the last column before believing it. A dense-graph agent receives 21.7 messages over the run against 2.4 on a sparse graph, a ninefold difference in exposure sitting on top of the structural one.

So we run the control.

# cap fan-out so dense-graph agents get the sparse message budget
for name in ("BA-dense", "ER-dense", "Empirical"):
    c = run_topology(name, max_fanout=6)
    print(f"{name:<12} capped spread {c.spread:.3f}  (uncapped {UNCAPPED[name]:.3f})")
#### OUTPUT ####
BA-dense     capped spread 0.462  (uncapped 0.389)
ER-dense     capped spread 0.471  (uncapped 0.384)
Empirical    capped spread 0.458  (uncapped 0.331)

Cap the message volume and most of the effect disappears.

Dense networks expose people to more messages, and exposure is doing most of the work. The empirical graph does not fully return to baseline, so structure still matters, but the first reading of that chart would have been mostly wrong.

Free Discussion With Memory

One last experiment with the influencer constraints removed: 150 agents, anybody may talk to anybody, everyone remembers their last eight conversations and writes a summary each night.

None 
forum = run_forum(model, personas[:150], days=7, window=8, topic="Humans will establish a city on Mars within fifty years")
print(f"cosine diversity      {forum.diversity[0]:.3f} -> {forum.diversity[-1]:.3f}   ({(forum.diversity[-1]/forum.diversity[0]-1)*100:+.1f}%)")
for i, ratio in enumerate(forum.dispersion_ratio, 1):
    print(f"within-cluster spread cluster {i}: {ratio:.1f}x wider")
#### OUTPUT ####
[forum] day 1  conversations 412  mean cosine diversity 0.196
[forum] day 7  conversations 401  mean cosine diversity 0.128

cosine diversity      0.203 -> 0.128   (-36.9%)
within-cluster spread cluster 1: 3.0x wider   cluster 2: 2.2x wider
None

Two things happen at once. The population converges, diversity falling 37 percent, while individual reasoning diverges, the spread inside each opinion group tripling.

I would not have predicted it, and it only appears when agents have memory, which the billion-agent runtime does not.

Where the Time and the Tokens Went

One round is a sequence of steps laid end to end.

STEP_COLOURS = ["#5E7085", "#2C6BAA", "#2C6BAA", INDIGO, INDIGO, AMBER, "#12855F", "#5E7085"]
fig, axes = plt.subplots(1, 2, figsize=(14.4, 4.6))
names, vals = list(STAGES), list(STAGES.values())
axes[0].barh(range(len(names)), vals, height=0.64,
             color=[TEAL, TEAL, AMBER, AMBER, AMBER, INDIGO])
axes[0].set_yticks(range(len(names))); axes[0].set_yticklabels(names)
axes[0].invert_yaxis(); axes[0].set_xlim(0, 210)
style(axes[0], f"The billion-agent run is {vals[-1]/sum(vals)*100:.2f}% of the whole build",
      "minutes", grid_axis="x")
left = 0.0
for (step, secs), colour in zip(ROUND_STEPS.items(), STEP_COLOURS):
    axes[1].barh([0], [secs], left=[left], color=colour, height=0.5,
                 edgecolor="white", label=step)
    left += secs
axes[1].set_yticks([]); axes[1].set_xlim(0, 1.15)
axes[1].legend(frameon=False, fontsize=9, ncol=4, loc="lower center")
style(axes[1], f"One round of 14,930,000 interactions, in {left:.2f} seconds",
      "seconds", grid_axis="x")
fig.savefig("images/p29_wallclock.png", dpi=200, bbox_inches="tight")
None 
total = sum(STAGES.values())
print(f"{'stage':<30}{'wall clock':>12}{'share':>10}")
for name, minutes in STAGES.items():
    hms = (f"{int(minutes//60)}h {int(minutes%60):02d}m {int(minutes%1*60):02d}s"
           if minutes >= 60 else f"{int(minutes):>2}m {int(minutes%1*60):02d}s")
    print(f"{name:<30}{hms:>12}{minutes/total*100:>9.1f}%")
#### OUTPUT ####
stage                          wall clock     share
Census (clean, render, embed)      9m 44s      3.9%
Herald (teacher labelling)     2h 59m 00s     71.6%
Pulse (100 rounds, 1e9 agents)     1m 44s      0.7%
total                          4h 09m 38s

The billion-agent simulation is seven tenths of one percent of the build.

And the tokens.

None None 
fixed = 270_000 * 1048
per_agent = (1_493_000_000 / 1_000_000_000) * 1048
print(f"crossover                 : {fixed/per_agent:,.0f} agents")
for n in (1e3, 1e6, 1e9):
    print(f"  N={n:>10.0e}  live {n*per_agent:.3e}  ratio {n*per_agent/fixed:>10,.2f}x")
#### OUTPUT ####
crossover                 : 180,844 agents
  N=     1e+03  live 1.565e+06  ratio       0.01x
  N=     1e+06  live 1.565e+09  ratio       5.53x
  N=     1e+09  live 1.565e+12  ratio    5529.63x
None

Below one hundred and eighty thousand agents this pipeline is a net loss, and you should just call the model.

It is better above a threshold, and below it considerably worse.

What This Proves and What It Does Not

None

A billion-agent social simulation runs in eleven gigabytes and one hundred and four seconds on commodity hardware, with open weights, for a fixed cost that does not grow with population.

What is not proven is anything about people, for six reasons.

  1. Ten thousand people, not a billion.
  2. The runtime is not a language model. It is a stochastic automaton whose transition rule a language model estimated. Every social claim inherits that model's biases, uncorrectable mid-run.
  3. There is no memory, and no message. No history, no content, no fatigue. Repeated exposure is memoryless, which the persuasion literature calls wrong.
  4. The constants are choices. A fifth are influencers, one percent speak per round, influencers never change, a hundred rounds. None is derived, and the last one mattered a lot.
  5. The prompt is doing measurable work. A worked example moved trust further than the whole social class gradient.
  6. Fidelity is to one teacher, not to reality. Every fidelity number here measures agreement with one Qwen model. Humans enter only in the Arena comparison, where the levels were wrong even where the orderings were right.

None of that makes the build wrong. It makes the claims precise.

Running It End to End

None 
# the teacher runs in its own shell, the Herald cells talk to it over HTTP
bash serve_teacher.sh

# population 1000000000 and the topic list both come from here
cat configs/billion.yaml

# then the notebook, top to bottom, no timeout because Lattice takes 55 minutes
# and Herald takes three hours
jupyter nbconvert --to notebook --execute earth-scale-society.ipynb \
  --output runs/latest.ipynb \
  --ExecutePreprocessor.timeout=-1
#### OUTPUT ####
[census]  96,125 records -> 10,000 personas (9,500 warm / 500 cold), 39.06 MiB
[lattice] 1e9 nodes, 2,999,999,994 edges -> 1,493,000,000 kept, 6.772 GB
[herald]  270,000 questions in 2h 59m, churn 26.3%, 0 schema violations
[echo]    30 epochs in 58s, shipped epoch 24 (F1 0.7698, gap 3.18)
[atlas]   900,000,000 entries -> 858.31 MiB, 0/27 dead channels
[sim]     100 rounds, 1,493,000,000 interactions, 0 model calls, 104.2s
[readout] agree 23.473%  neutral 43.214%  disagree 33.313%

Where This Leaves Us

The state that decides a social interaction is small and finite, so it can be enumerated in advance, and once it has been the population becomes free.

The teacher runs once. The surrogate is small because it knows one narrow thing. The runtime is a numpy program reading an array, which is why a billion people fit inside one GPU.

Train on the teacher's hesitation and not just its answer. Choose the checkpoint on the aggregate you care about. Store a distribution rather than a decision, because the byte costs the same and rare outcomes drive long-run behaviour. Hold out a slice of your population to measure what happens to strangers. Run the control that separates your variable from the one beside it.

This system reproduces population-level dynamics faithfully enough to study, at any scale, on a fixed budget. Its agents are ten thousand people wearing a hundred thousand name tags each.

If you build one thing from this post, do not build a bigger simulator. Build the table. Then check whether the scale you paid for changed the answer.

The full code:

GitHub - FareedKhan-dev/earth-scale-society: One billion LLM-grounded social agents that make zero… One billion LLM-grounded social agents that make zero model calls at runtime. A 235B open teacher is paid once…

Thanks for reading this far.