L04 — Echo Chambers

"""L04 Echo Chambers: Deffuant bounded-confidence opinions on an adaptive network (Mesa 3).

Run from the repo root:  uv run python L04_echo_chambers/sim.py
"""
import os

os.environ.setdefault("OMP_NUM_THREADS", "4")

import time
import warnings
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import mesa
import networkx as nx
import numpy as np
from matplotlib.collections import LineCollection
from mesa.discrete_space import FixedAgent, Network
from scipy import stats

OUT = Path(__file__).parent / "output"


# ---------------------------------------------------------------- metrics (pure functions)
def major_clusters(x, gap=0.02, min_frac=0.02):
    """Opinion clusters: sorted opinions split wherever the gap > `gap`; count groups >= 2% of agents."""
    xs = np.sort(x)
    sizes = np.diff(np.r_[0, np.flatnonzero(np.diff(xs) > gap) + 1, len(xs)])
    return int((sizes >= min_frac * len(xs)).sum())


def bimodality(x):
    """Sarle's bimodality coefficient (skew^2+1)/(excess kurtosis + small-sample term); >5/9 hints at >1 mode."""
    n = len(x)
    if np.std(x) < 0.01:  # collapsed to one point: unimodal by definition (shape of 1e-4 noise is meaningless)
        return 0.0
    g, k = stats.skew(x, bias=False), stats.kurtosis(x, bias=False)
    return float((g**2 + 1) / (k + 3 * (n - 1) ** 2 / ((n - 2) * (n - 3))))


def opinion_assortativity(G, x):
    """Pearson correlation of opinions across the two ends of every tie (= nx numeric assortativity)."""
    e = np.array(G.edges(), dtype=int)
    a, b = np.r_[x[e[:, 0]], x[e[:, 1]]], np.r_[x[e[:, 1]], x[e[:, 0]]]
    return float(np.corrcoef(a, b)[0, 1]) if a.std() > 0.01 else np.nan  # undefined once opinions collapse


# ---------------------------------------------------------------- the agent-based model
class Person(FixedAgent):
    """One network node. State: opinion x in [0, 1]. Fixed to its cell; its ties are the cell's connections."""

    def __init__(self, model, opinion, cell):
        super().__init__(model)
        self.x = float(opinion)
        self.cell = cell

    def step(self):
        m = self.model
        friends = list(self.cell.connections.values())
        if not friends:  # every tie was rewired away: an isolate listens to no one
            return
        their_cell = self.random.choice(friends)
        other = their_cell.agents[0]
        d = other.x - self.x
        if abs(d) < m.epsilon:  # Deffuant: close enough -> both move toward each other
            self.x += m.mu * d
            other.x -= m.mu * d
        elif self.random.random() < m.rewire:  # too far -> "recommender" swaps the tie for a like-minded one
            new = m.recommend(self)
            if new is not None:
                m.grid.remove_connection(self.cell, their_cell)
                m.grid.add_connection(self.cell, new.cell)
                m.rewires += 1


class EchoChamber(mesa.Model):
    """Deffuant et al. (2000) bounded confidence + opinion-driven rewiring (cf. Kozma & Barrat 2008).

    n agents on a G(n, m) random graph with mean degree k. Each step every agent (random order) talks
    to one random neighbour. |dx| < epsilon: both move mu*dx closer. Otherwise, with probability
    `rewire`, the agent drops that tie and follows someone the recommender finds with |dx| < epsilon.
    Edge count is conserved, so only *who* talks to whom changes.
    """

    def __init__(self, n=400, k=12, epsilon=0.2, rewire=0.0, mu=0.3, t_max=600,
                 collect_every=0, track_agents=False, rng=None):
        super().__init__(rng=rng)
        self.epsilon, self.rewire, self.mu, self.t_max = epsilon, rewire, mu, t_max
        self.collect_every = collect_every
        self.rewires, self._q = 0, (-1, 0.0)  # ties rewired so far; Louvain cache keyed on that count
        G = nx.gnm_random_graph(n, n * k // 2, seed=self.random)
        self.grid = Network(G, random=self.random, layout=None)  # layout=None: skip the spring layout
        self.people = list(Person.create_agents(self, n, self.rng.random(n), list(self.grid.all_cells)))
        self.datacollector = mesa.DataCollector(
            model_reporters={
                "clusters": self.n_clusters,
                "bimodality": lambda m: bimodality(m.opinions()),
                "modularity": self.modularity,
                "assortativity": lambda m: opinion_assortativity(m.grid.G, m.opinions()),
                "components": self.n_components,
                "t_end": "steps",
            },
            agent_reporters={"opinion": "x"} if track_agents else None,
        )
        if collect_every:
            self.datacollector.collect(self)

    def opinions(self):
        return np.array([p.x for p in self.people])

    def recommend(self, person):
        # ponytail: rejection sampling of 20 random strangers; misses when like-minded strangers are rare
        for _ in range(20):
            cand = self.random.choice(self.people)
            if (cand is not person and abs(cand.x - person.x) < self.epsilon
                    and cand.cell.coordinate not in person.cell.connections):
                return cand
        return None

    def n_clusters(self):
        return major_clusters(self.opinions())

    def modularity(self):
        """Louvain modularity Q of the tie structure alone (opinions not used). Recomputed only if ties changed."""
        if self._q[0] != self.rewires:
            G = self.grid.G
            self._q = (self.rewires, nx.community.modularity(G, nx.community.louvain_communities(G, seed=0)))
        return self._q[1]

    def n_components(self):
        return sum(len(c) >= 0.02 * len(self.people) for c in nx.connected_components(self.grid.G))

    def converged(self, tol=1e-3):
        x = self.opinions()
        e = np.array(self.grid.G.edges(), dtype=int)
        d = np.abs(x[e[:, 0]] - x[e[:, 1]])
        talking = d < self.epsilon
        return (d[talking] < tol).all() and (self.rewire == 0 or talking.all())

    def step(self):
        self.agents.shuffle_do("step")
        done = self.steps >= self.t_max or (self.steps % 10 == 0 and self.converged())
        if done or (self.collect_every and self.steps % self.collect_every == 0):
            self.datacollector.collect(self)
        self.running = not done


# ---------------------------------------------------------------- experiments
N, K, T_MAX = 400, 12, 600
EPS = [round(e, 3) for e in np.arange(0.05, 0.501, 0.025)]  # confidence bound epsilon
RW = [0.0, 0.05, 0.1, 0.2, 0.3, 0.5, 0.75, 1.0]              # rewiring probability r, dense near 0
SEEDS = list(range(5))  # same seeds in every cell = common random numbers (L03), sharper contrasts
METRICS = ["clusters", "modularity", "assortativity", "bimodality", "components", "t_end"]


def sweep():
    """Parallel eps x r x seed factorial with mesa.batch_run -> dict metric -> array [r, eps, seed]."""
    rows = mesa.batch_run(
        EchoChamber,
        parameters={"n": N, "k": K, "epsilon": EPS, "rewire": RW, "t_max": T_MAX},
        number_processes=4, rng=SEEDS, max_steps=T_MAX, data_collection_period=-1,
        display_progress=False,
    )
    Z = {m: np.full((len(RW), len(EPS), len(SEEDS)), np.nan) for m in METRICS}
    for row in rows:
        i, j = RW.index(row["rewire"]), EPS.index(row["epsilon"])
        for m in METRICS:
            Z[m][i, j, row["iteration"]] = row[m]
    return Z


def verify(n_seeds=10):
    """No rewiring: number of major opinion clusters ~ floor(1/(2 eps)) (Deffuant et al. 2000)."""
    eps = [0.4, 0.3, 0.2, 0.15, 0.1]
    rows = mesa.batch_run(
        EchoChamber, parameters={"n": N, "k": K, "epsilon": eps, "rewire": 0.0, "t_max": T_MAX},
        number_processes=4, rng=list(range(100, 100 + n_seeds)), max_steps=T_MAX,
        data_collection_period=-1, display_progress=False,
    )
    got = {e: np.array([r["clusters"] for r in rows if r["epsilon"] == e]) for e in eps}
    print("\nVERIFY  r = 0 (static network), %d seeds each, N = %d, <k> = %d" % (n_seeds, N, K))
    print("  eps    floor(1/2eps)   clusters mean +- sd   [min..max]")
    ok = True
    for e in eps:
        rule, c = int(1 / (2 * e) + 1e-9), got[e]
        hit = abs(c.mean() - rule) <= 1
        ok &= hit
        print("  %.2f   %d               %.2f +- %.2f          [%d..%d]  %s"
              % (e, rule, c.mean(), c.std(ddof=1), c.min(), c.max(), "ok" if hit else "MISS"))
    return eps, got, ok


def run_tracked(eps, r, seed, t_max=250, every=1, keep_edges=False):
    """One run with history every `every` steps: times t, opinions X[t, agent] (DataCollector agent reporter),
    model metrics ts, and (optionally) the edge list after every step."""
    m = EchoChamber(n=N, k=K, epsilon=eps, rewire=r, t_max=t_max, collect_every=every, track_agents=True, rng=seed)
    edges = [np.array(m.grid.G.edges())]
    while m.running:
        m.step()
        if keep_edges:
            edges.append(np.array(m.grid.G.edges()))
    X = m.datacollector.get_agent_vars_dataframe()["opinion"].unstack().to_numpy()
    ts = {k: np.array(v, dtype=float) for k, v in m.datacollector.model_vars.items()}
    return ts["t_end"].astype(int), X, ts, edges


# ---------------------------------------------------------------- figures
BG, FG, DIM = "#0b0e14", "#d8dee9", "#8a93a3"
CMAP = plt.get_cmap("Spectral_r")  # opinion 0 = blue ... 1 = red
plt.style.use("dark_background")
plt.rcParams.update({
    "figure.facecolor": BG, "axes.facecolor": BG, "savefig.facecolor": BG, "axes.edgecolor": "#4c5566",
    "axes.labelcolor": FG, "xtick.color": DIM, "ytick.color": DIM, "text.color": FG,
    "font.size": 10, "axes.titlesize": 11, "axes.grid": False,
})


def trajectories(ax, t, X, t_max, lw=0.5, alpha=0.35):
    """Deffuant 'spaghetti': one line per agent, coloured by its starting opinion; drawn flat after convergence."""
    if t[-1] < t_max:
        t, X = np.r_[t, t_max], np.vstack([X, X[-1]])
    segs = np.stack([np.broadcast_to(t, X.T.shape), X.T], axis=-1)
    ax.add_collection(LineCollection(segs, colors=CMAP(X[0]), lw=lw, alpha=alpha))
    ax.set(xlim=(0, t_max), ylim=(-0.02, 1.02))


def draw_network(ax, pos, e, x, eps, node_size=16):
    """Nodes coloured by opinion (with a soft glow); ties coloured if the two ends can still talk, grey if not."""
    col = CMAP((x[e[:, 0]] + x[e[:, 1]]) / 2)
    col[:, 3] = 0.28
    col[np.abs(x[e[:, 0]] - x[e[:, 1]]) >= eps] = (0.62, 0.64, 0.70, 0.10)
    ax.add_collection(LineCollection(pos[e], colors=col, lw=0.45, zorder=1))
    ax.scatter(*pos.T, s=node_size * 5, c=x, cmap=CMAP, vmin=0, vmax=1, alpha=0.12, edgecolors="none", zorder=2)
    ax.scatter(*pos.T, s=node_size, c=x, cmap=CMAP, vmin=0, vmax=1, edgecolors="none", zorder=3)
    ax.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1), aspect="equal")
    ax.set_axis_off()


def layouts(edges, frames):
    """Force-directed layout continued frame-to-frame from the previous positions, so the picture morphs smoothly."""
    pos, P = None, {}
    for f in frames:
        G = nx.empty_graph(N)
        G.add_edges_from(edges[f])
        # gravity pulls disconnected chambers back toward the centre, so they stay big on screen
        pos = nx.spring_layout(G, pos=pos, method="energy", gravity=5.0, iterations=100 if pos is None else 10, seed=1)
        P[f] = np.array([pos[i] for i in range(N)])
    return P


def plot_verify(eps_v, got, Z):
    fig, ax = plt.subplots(figsize=(7.5, 4.8))
    e = np.linspace(0.045, 0.5, 800)
    ax.plot(e, np.floor(1 / (2 * e)), color="#ffb454", lw=2, label=r"$\lfloor 1/(2\varepsilon) \rfloor$  (Deffuant et al. 2000)")
    ax.plot(e, 1 / (2 * e), color="#ffb454", lw=1, ls=":", alpha=0.7, label=r"$1/(2\varepsilon)$")
    c0 = Z["clusters"][0]
    ax.plot(EPS, c0.mean(1), "o-", color="#59c2ff", ms=4, lw=1, alpha=0.8,
            label="sweep, r = 0 (mean of %d seeds)" % c0.shape[1])
    for ev in eps_v:
        c = got[ev]
        ax.errorbar(ev, c.mean(), yerr=c.std(ddof=1), fmt="s", color="#aad94c", ms=7, capsize=4, zorder=5)
    ax.plot([], [], "s", color="#aad94c", label="verification, r = 0 (mean +- sd, %d seeds)" % len(got[eps_v[0]]))
    ax.set(xlabel=r"confidence bound $\varepsilon$ (opinion units)", ylabel="major opinion clusters (>= 2% of agents)",
           title="Verification: static network reproduces Deffuant's cluster law", ylim=(0, 11), xlim=(0.04, 0.52))
    ax.legend(frameon=False, fontsize=9)
    fig.tight_layout()
    fig.savefig(OUT / "verify_clusters.png", dpi=150)
    plt.close(fig)


def plot_phase(Z):
    de, ri = EPS[1] - EPS[0], np.arange(len(RW))  # r grid is non-uniform: plot rows by index, label with r
    ext = [EPS[0] - de / 2, EPS[-1] + de / 2, -0.5, len(RW) - 0.5]
    C, Q = Z["clusters"].mean(2), Z["modularity"].mean(2)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)  # "mean of empty slice"
        A = np.nanmean(Z["assortativity"], axis=2)  # NaN (grey) where every seed reached consensus
    panels = [
        (C, "(a) opinion clusters  ->  regime", "magma", matplotlib.colors.LogNorm(1, 10), "major clusters"),
        (Q, "(b) network modularity Q (Louvain)", "viridis", matplotlib.colors.Normalize(0.2, 0.9), "Q"),
        (A, "(c) opinion assortativity of ties", "cividis", matplotlib.colors.Normalize(0, 1), "Pearson r across ties"),
        (Z["t_end"].mean(2), "(d) steps to converge (cap %d)" % T_MAX, "inferno",
         matplotlib.colors.Normalize(0, T_MAX), "steps"),
    ]
    fig, axes = plt.subplots(2, 2, figsize=(14, 9.5))
    for ax, (D, title, cm, norm, cb) in zip(axes.flat, panels):
        cmo = plt.get_cmap(cm).with_extremes(bad="#2a2f3a")
        im = ax.imshow(D, origin="lower", extent=ext, aspect="auto", cmap=cmo, norm=norm, interpolation="nearest")
        fig.colorbar(im, ax=ax, label=cb, pad=0.02)
        ax.set(title=title, xlabel=r"confidence bound $\varepsilon$", ylabel="rewiring probability r",
               yticks=ri, yticklabels=[f"{r:g}" for r in RW])
        ax.contour(EPS, ri, C, levels=[1.5, 2.5], colors="white", linewidths=1.2)
        ax.contour(EPS, ri, Z["components"].mean(2), levels=[1.5], colors="#5ccfe6", linewidths=1.6, linestyles="--")
    a = axes[0, 0]
    for x_, txt in [(0.41, "CONSENSUS"), (0.215, "POLARISATION"), (0.095, "FRAGMENTATION")]:
        a.text(x_, len(RW) / 2, txt, ha="center", va="center", rotation=90, fontsize=11, weight="bold", color="white", alpha=0.9)
    fig.suptitle("Phase diagram: %d agents, <k> = %d, %d seeds per cell, %d runs via mesa.batch_run (4 processes)"
                 % (N, K, Z["clusters"].shape[2], Z["clusters"].size), fontsize=12)
    fig.text(0.5, 0.01, "white lines: mean clusters = 1.5 and 2.5 (regime borders).   cyan dashed: mean network "
             "components = 1.5 (upper left of it the network has shattered into sealed echo chambers).\n(c) grey = "
             "opinions collapsed to one value, assortativity undefined; yellow specks in the consensus zone are tiny "
             "extremist cliques (< 2% of agents) the recommender sealed off.", ha="center", va="bottom", fontsize=9,
             color=DIM)
    fig.tight_layout(rect=(0, 0.04, 1, 1))
    fig.savefig(OUT / "phase_diagram.png", dpi=130)
    plt.close(fig)


def plot_spaghetti(runs, eps_list, r_list, t_max):
    fig, axes = plt.subplots(len(r_list), len(eps_list), figsize=(15, 7.2), sharex=True, sharey=True)
    for i, r in enumerate(r_list):
        for j, eps in enumerate(eps_list):
            t, X, ts = runs[eps, r]
            ax = axes[i, j]
            trajectories(ax, t, X, t_max)
            c, parts = int(ts["clusters"][-1]), int(ts["components"][-1])
            ax.set_title(r"$\varepsilon$ = %.2f, r = %.1f:  %d cluster%s, Q = %.2f, %d network part%s"
                         % (eps, r, c, "s"[: c > 1], ts["modularity"][-1], parts, "s"[: parts > 1]), fontsize=10)
            if i == len(r_list) - 1:
                ax.set_xlabel("time (steps; 1 step = every agent talks once)")
        axes[i, 0].set_ylabel(("static network" if r == 0 else "recommender rewiring") + "\nopinion x")
    fig.suptitle("Deffuant trajectories: epsilon sets how many camps form; rewiring decides whether they still talk",
                 fontsize=12)
    fig.tight_layout()
    fig.savefig(OUT / "spaghetti.png", dpi=140)
    plt.close(fig)


def animate(t, X, ts, edges, eps, r, path, hold=10):
    """One GIF frame per DataCollector snapshot: t[i] = step, X[i] = opinions, edges[t[i]] = ties at that step."""
    from matplotlib.animation import FuncAnimation, PillowWriter

    T = t[-1]
    P = layouts(edges, t)  # ponytail: O(N^2) force layout per frame dominates; fewer frames or igraph if it drags
    fig = plt.figure(figsize=(12, 6.2), dpi=80)
    gs = fig.add_gridspec(2, 2, width_ratios=[1.05, 1], height_ratios=[1.5, 1], left=0.01, right=0.97,
                          top=0.87, bottom=0.09, wspace=0.12, hspace=0.3)
    axn, axs, axm = fig.add_subplot(gs[:, 0]), fig.add_subplot(gs[0, 1]), fig.add_subplot(gs[1, 1])
    axs.set(xlim=(0, T), ylim=(-0.02, 1.02), ylabel="opinion x")
    axm.set(xlim=(0, T), ylim=(-0.05, 1.05), xlabel="time (steps)")
    lc = LineCollection([], colors=CMAP(X[0]), lw=0.5, alpha=0.35)
    axs.add_collection(lc)
    lq, = axm.plot([], [], color="#aad94c", lw=2, label="modularity Q")
    la, = axm.plot([], [], color="#59c2ff", lw=2, label="opinion assortativity")
    axm.axhline(ts["modularity"][0], color="#aad94c", lw=0.8, ls=":", alpha=0.6)
    axm.legend(frameon=False, loc="center right", fontsize=9)
    head = fig.text(0.02, 0.95, "", fontsize=13, weight="bold")
    sub = fig.text(0.02, 0.91, "", fontsize=10, color=DIM)

    def update(i):
        f = t[i]
        axn.clear()
        draw_network(axn, P[f], edges[f], X[i], eps, node_size=12)
        lc.set_segments(np.stack([np.broadcast_to(t[: i + 1], (N, i + 1)), X[: i + 1].T], axis=-1))
        lq.set_data(t[: i + 1], ts["modularity"][: i + 1])
        la.set_data(t[: i + 1], ts["assortativity"][: i + 1])
        disc = np.mean(np.abs(X[i][edges[f][:, 0]] - X[i][edges[f][:, 1]]) >= eps)
        head.set_text("Echo chambers forming   t = %d" % f)
        sub.set_text(r"Deffuant $\varepsilon$ = %.2f, recommender r = %.1f, %d agents  |  ties across the divide (grey): %4.1f%%"
                     % (eps, r, N, 100 * disc))
        return ()

    anim = FuncAnimation(fig, update, frames=list(range(len(t))) + [len(t) - 1] * hold, blit=False)
    anim.save(path, writer=PillowWriter(fps=10))
    plt.close(fig)
    return P


def plot_hero(t, X, edges, P, eps, r, i_mid=7):
    i_show = int(np.argmax(np.abs(X - X[-1]).max(axis=1) < 0.005))  # crop the flat tail of the spaghetti
    fig = plt.figure(figsize=(18, 9))
    gs = fig.add_gridspec(2, 3, height_ratios=[3, 1], left=0.04, right=0.98, top=0.86, bottom=0.08,
                          hspace=0.18, wspace=0.02)
    for j, (i, label) in enumerate([(0, "one random social network"), (i_mid, "the recommender sorts ties"),
                                    (len(t) - 1, "sealed echo chambers")]):
        ax = fig.add_subplot(gs[0, j])
        draw_network(ax, P[t[i]], edges[t[i]], X[i], eps, node_size=20)
        ax.set_title("t = %d: %s" % (t[i], label), fontsize=13, color=DIM)
    a2 = fig.add_subplot(gs[1, :])
    trajectories(a2, t, X, t[i_show] + 10, lw=0.6, alpha=0.4)
    for i in (0, i_mid):
        a2.axvline(t[i], color=DIM, lw=0.8, ls=":")
    a2.set(xlabel="time (steps)", ylabel="opinion x")
    fig.text(0.5, 0.95, "ECHO CHAMBERS", ha="center", fontsize=24, weight="bold")
    fig.text(0.5, 0.905, r"Deffuant bounded confidence ($\varepsilon$ = %.2f) + a like-minded recommender (r = %.1f) "
             "split one network of %d Mesa agents into camps that no longer talk" % (eps, r, N),
             ha="center", fontsize=12, color=DIM)
    fig.savefig(OUT / "hero.png", dpi=110)
    plt.close(fig)


# ---------------------------------------------------------------- main
def main():
    OUT.mkdir(exist_ok=True)
    t0 = time.perf_counter()

    eps_v, got, ok = verify()
    print("  [%.0f s]" % (time.perf_counter() - t0))

    Z = sweep()
    print("\nSWEEP  %d runs (eps x r x seeds = %d x %d x %d)  [%.0f s]"
          % (Z["clusters"].size, len(EPS), len(RW), len(SEEDS), time.perf_counter() - t0))
    C = Z["clusters"]
    for i, r in enumerate(RW):
        cons = [e for j, e in enumerate(EPS) if (C[i, j:] == 1).all()]
        q = Z["modularity"][i, EPS.index(0.2)].mean()
        split = (Z["components"][i] == C[i]).mean()
        print("  r = %-4g: consensus for all eps >= %.3f | Q(eps=0.2) = %.2f | network parts == opinion camps in %3.0f%% of runs"
              % (r, cons[0], q, 100 * split))
    plot_verify(eps_v, got, Z)
    plot_phase(Z)

    runs, eps_list, r_list = {}, [0.4, 0.2, 0.1], [0.0, 0.5]
    for r in r_list:
        for e in eps_list:
            t, X, ts, _ = run_tracked(e, r, seed=7, every=2)  # every 2nd step: Louvain is the costly reporter
            runs[e, r] = (t, X, ts)
    plot_spaghetti(runs, eps_list, r_list, t_max=250)
    print("\nspaghetti done  [%.0f s]" % (time.perf_counter() - t0))

    eps, r = 0.15, 0.3  # seed 6 = three similar-sized camps; illustrative, the sweep carries the statistics
    t, X, ts, edges = run_tracked(eps, r, seed=6, every=3, keep_edges=True)  # snapshot every 3 steps = GIF frames
    drift = abs(X[-1].mean() - X[0].mean())
    print("SHOWCASE eps=%.2f r=%.1f: %d steps, %d clusters, %d network parts, Q = %.2f, assortativity = %.2f"
          % (eps, r, t[-1], ts["clusters"][-1], ts["components"][-1], ts["modularity"][-1],
             ts["assortativity"][-1]))
    print("  mean opinion drift |<x>_T - <x>_0| = %.1e (pairwise update conserves the mean exactly)" % drift)
    P = animate(t, X, ts, edges, eps, r, OUT / "echo_chambers.gif")
    plot_hero(t, X, edges, P, eps, r)
    print("  gif + hero done  [%.0f s], gif = %.1f MB"
          % (time.perf_counter() - t0, (OUT / "echo_chambers.gif").stat().st_size / 1e6))

    assert drift < 1e-10, "Deffuant update must conserve the mean opinion"
    assert ok, "cluster count strays > 1 from floor(1/(2 eps))"
    print("\nall checks passed, total %.0f s" % (time.perf_counter() - t0))


if __name__ == "__main__":
    main()