Graph Learning

Notebooks Project Live marimo ↗

One graph, the whole course: a running project on an X-like retweet network

This is a read-only view of the notebook source: text and code, no outputs. To run it, open the live marimo server when the workstation tunnel is up, or run uv run marimo edit notebooks/project/retweet_thread.py.

import warnings

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", message=".*Level Zero.*")

import marimo as mo
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
import scipy.sparse as sp
import torch
import torch.nn as nn
import torch.nn.functional as F
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import StratifiedKFold

import json
import os

from gl.nb import reveal
from gl.retweet import FEATURE_NAMES, WorldConfig, make_world, stealthy

torch.set_num_threads(4)
plt.rcParams.update({"figure.dpi": 100, "axes.spines.top": False, "axes.spines.right": False})

One graph, the whole course: a running project on an X-like retweet network

CS224W companion · the thread that ties the 21 sessions to one research question

Every lecture notebook ends with "connecting to your research (SNA on X/Twitter)". This notebook makes that concrete: one network, carried through the course, with one question that gets sharper each time:

Can we find coordinated accounts from the structure of who amplifies whom, when their individual behaviour looks normal?

The data is a synthetic X world (gl.retweet), because real X data has cost about $42,000 a month since 2023 (Lecture 17 frontier). Its mechanisms are explicit and its ground truth is known, which is exactly what you want for checking a method before trusting it on real data: users in four communities with heavy-tailed activity, a follow graph from a degree-corrected SBM, time-decayed retweet cascades, and three coordinated rings of 10 accounts that amplify each other and a few target accounts. In the default "stealthy" setting the rings retweet with human-like delays and volume and add camouflage retweets of random content, so per-account statistics barely separate them.

Swap in real data at the end (gl.retweet.load_retweets_csv).

world_kind = mo.ui.dropdown(
    {"stealthy rings (human-like timing)": "stealthy", "naive rings (retweet within minutes)": "naive"},
    value="stealthy rings (human-like timing)",
    label="bots",
)
world_seed = mo.ui.number(0, 99, value=0, label="seed")
mo.hstack([world_kind, world_seed], justify="start")

world = make_world(int(world_seed.value), stealthy() if world_kind.value == "stealthy" else WorldConfig())
users = world.users
is_bot = users.is_bot.to_numpy()
community = users.community.to_numpy()
G = world.retweet_graph()  # u -> v: u retweeted v
n = G.number_of_nodes()
A = nx.to_scipy_sparse_array(G, nodelist=range(n), weight=None, format="csr").astype(float)
U = ((A + A.T) > 0).astype(float)  # undirected: any retweet either way
M = A.multiply(A.T)  # mutual: both retweeted each other
X_feat = world.features()
mo.md(
    f"**World:** {n} accounts ({int(is_bot.sum())} in rings), {len(world.follows):,} follows, {len(world.posts):,} posts, "
    f"{len(world.retweets):,} retweets → a directed retweet graph with {G.number_of_edges():,} edges "
    f"({int(M.sum() // 2):,} mutual pairs)."
)

L01 · Which graph? Representation choices change what you can see

The same log of retweets gives several graphs: directed or not, weighted or not, users only or users + hashtags. A common shortcut is to fold the user–hashtag bipartite graph into a user–user graph ("used the same hashtag").

Predict first. The folded graph will have a high clustering coefficient. Is that evidence of social closure? (This is the L01 proof workshop.)

def _fold_clustering(uh, rng=None):
    tags = uh.groupby("hashtag").user.apply(list)
    if rng is not None:  # degree-preserving shuffle of the bipartite graph: permute the user column
        uh = uh.assign(user=rng.permutation(uh.user.to_numpy()))
        tags = uh.groupby("hashtag").user.apply(lambda s: list(set(s)))
    g = nx.Graph()
    for us in tags:
        us = list(us)[:150]  # cap huge hashtags so the fold stays small
        g.add_edges_from((a, b) for i, a in enumerate(us) for b in us[i + 1 :])
    return nx.average_clustering(g), g.number_of_edges()

_uh = world.user_hashtag()
_obs = _fold_clustering(_uh)
_null = [_fold_clustering(_uh, np.random.default_rng(s)) for s in range(3)]
_und = nx.average_clustering(world.retweet_graph().to_undirected())
fold_table = pd.DataFrame(
    [
        {"graph": "retweet graph (undirected)", "avg clustering": _und, "edges": world.retweet_graph().number_of_edges()},
        {"graph": "user–hashtag fold", "avg clustering": _obs[0], "edges": _obs[1]},
        {"graph": "fold of a randomised bipartite graph (mean of 3)", "avg clustering": np.mean([c for c, _ in _null]),
         "edges": int(np.mean([e for _, e in _null]))},
    ]
).round(3)
fold_table
_c = fold_table["avg clustering"].tolist() · 6 lines
_c = fold_table["avg clustering"].tolist()
mo.md(rf"""
The fold's clustering ({_c[1]:.2f}) is far above the retweet graph's ({_c[0]:.2f}), and a fold of a **randomised**
bipartite graph with the same degrees is also high ({_c[2]:.2f}). Most of the fold's clustering is the artefact of folding (each hashtag
becomes a clique), not social closure. Compare against the folded null model before reading triangles in projected networks.
""")

L02–L04 · From embeddings to GNNs: who belongs to which community?

Task 1 (easy, homophilous): predict each account's community from 10% labelled accounts. Three models: per-account features only (logistic regression), a spectral embedding of the retweet graph (L02, Laplacian eigenmaps), and a 2-layer GCN written by hand (L03).

Predict first. Rank the three. How much does the graph add over the features?

def gcn_norm(Ad): · 53 lines
def gcn_norm(Ad):
    Ad = Ad + sp.eye(Ad.shape[0])
    d = np.asarray(Ad.sum(1)).ravel()
    Dm = sp.diags(1 / np.sqrt(d))
    Ah = (Dm @ Ad @ Dm).tocoo()
    return torch.sparse_coo_tensor(np.vstack([Ah.row, Ah.col]), Ah.data, Ah.shape, dtype=torch.float32).coalesce()

A_hat = gcn_norm(U)

class GCN(nn.Module):
    def __init__(self, d_in, d_out, h=32):
        super().__init__()
        self.l1, self.l2 = nn.Linear(d_in, h), nn.Linear(h, d_out)

    def forward(self, x):
        x = F.relu(torch.sparse.mm(A_hat, self.l1(x)))
        return torch.sparse.mm(A_hat, self.l2(F.dropout(x, 0.5, self.training)))

def train_gcn(x, y, train, epochs=200, seed=0, weight=None):
    torch.manual_seed(seed)
    m = GCN(x.shape[1], int(y.max()) + 1)
    opt = torch.optim.Adam(m.parameters(), lr=0.01, weight_decay=5e-4)
    xt, yt = torch.tensor(x), torch.tensor(y)
    w = None if weight is None else torch.tensor(weight, dtype=torch.float32)
    for _ in range(epochs):
        m.train()
        opt.zero_grad()
        F.cross_entropy(m(xt)[train], yt[train], weight=w).backward()
        opt.step()
    m.eval()
    with torch.no_grad():
        return F.softmax(m(xt), 1).numpy()

def spectral(Ad, k=16):
    L = sp.csgraph.laplacian(Ad, normed=True)
    vals, vecs = sp.linalg.eigsh(L.astype(float), k=k + 1, which="SM", tol=1e-4, maxiter=5000)
    return vecs[:, np.argsort(vals)[1:]].astype(np.float32)

_rng = np.random.default_rng(0)
_train = np.zeros(n, bool)
_train[_rng.choice(n, n // 10, replace=False)] = True
_test = ~_train
_Xs = (X_feat - X_feat.mean(0)) / X_feat.std(0)
_emb = spectral(U)
_emb = (_emb - _emb.mean(0)) / _emb.std(0)  # same scale as the features, or the L2 penalty ignores it
_rows = []
for _name, _X in [("features only (LR)", _Xs), ("spectral embedding (LR)", _emb), ("features + spectral (LR)", np.hstack([_Xs, _emb]))]:
    _p = LogisticRegression(max_iter=3000).fit(_X[_train], community[_train]).predict(_X[_test])
    _rows.append({"model": _name, "test accuracy": (_p == community[_test]).mean()})
_p = train_gcn(_Xs, community, torch.tensor(_train)).argmax(1)
_rows.append({"model": "GCN on features (2 layers)", "test accuracy": (_p[_test] == community[_test]).mean()})
comm_table = pd.DataFrame(_rows).round(3)
comm_table
_a = dict(zip(comm_table.model, comm_table["test accuracy"], strict=True)) · 7 lines
_a = dict(zip(comm_table.model, comm_table["test accuracy"], strict=True))
mo.md(rf"""
Features alone reach {_a["features only (LR)"]:.2f}: individual behaviour says little about camp. The graph says a lot:
the spectral embedding gets {_a["spectral embedding (LR)"]:.2f} and the GCN {_a["GCN on features (2 layers)"]:.2f}, because
retweets mostly stay inside a community (**homophily**, L02). This is the easy task, the one GNN benchmarks are made of.
Coordinated rings are the hard one, because they are built to look like members of the community they target.
""")

L05 · Honest splits: predicting future retweet edges

Task 2: link prediction. Score candidate pairs with Adamic–Adar (L20) on the undirected retweet graph and measure ROC AUC (L05 proof workshop). Two protocols: random edge split (hide 10% of edges at random) vs. temporal split (graph of days 0–20, predict new edges in days 20–30).

Predict first. Which protocol gives the higher AUC, and why should you trust the lower one?

def _aa_auc(train_g, pos, rng):
    neg = []
    while len(neg) < len(pos):
        a, b = rng.integers(n, size=2)
        if a != b and not train_g.has_edge(a, b):
            neg.append((int(a), int(b)))
    pairs = list(pos) + neg
    scores = [s for _, _, s in nx.adamic_adar_index(train_g, pairs)]
    return roc_auc_score([1] * len(pos) + [0] * len(neg), scores)

_rng = np.random.default_rng(1)
_full = world.retweet_graph().to_undirected()
_edges = list(_full.edges())
_rng.shuffle(_edges)
_k = len(_edges) // 10
_train_rand = nx.Graph(_edges[_k:])
_train_rand.add_nodes_from(range(n))
_auc_rand = _aa_auc(_train_rand, _edges[:_k], _rng)
_cut = 20 * 24
_past = world.retweet_graph(until=_cut).to_undirected()
_new = [e for e in _full.edges() if not _past.has_edge(*e)]
_new = [_new[i] for i in _rng.choice(len(_new), min(len(_new), _k), replace=False)]
_auc_time = _aa_auc(_past, _new, _rng)
split_table = pd.DataFrame(
    [{"protocol": "random 10% of edges", "Adamic–Adar AUC": _auc_rand},
     {"protocol": "temporal: days 0–20 → new edges in days 20–30", "Adamic–Adar AUC": _auc_time}]
).round(3)
split_table
_r, _t = split_table["Adamic–Adar AUC"] · 6 lines
_r, _t = split_table["Adamic–Adar AUC"]
mo.md(rf"""
Random split {_r:.3f} vs temporal {_t:.3f}. A random split hides edges whose endpoints already interacted *after* the
hidden edge appeared, so the "past" contains information from the future. The temporal protocol asks the question you
actually care about: who will start amplifying whom. Report it, and expect it to be {"lower" if _t < _r else "similar"}.
""")

L06–L07 · The research question: rings that look normal

Task 3: detect ring members (30 of ~1,230 accounts). Per-account features were designed to overlap. What separates rings is structure: they retweet each other (reciprocity) and form dense groups (triangles in the mutual graph).

Lecture 6 says message passing cannot count triangles (1-WL is blind to them), and Lecture 7 says closed walks can: \((M^3)_{vv}=2\,t(v)\) in the mutual graph \(M\) (L07 proof workshop). We compare, all with logistic regression and 5-fold CV:

  1. per-account features;
  2. + neighbourhood means of the features (what a 1-layer mean-aggregation GNN sees);
  3. + reciprocity and mutual triangles \(\tfrac12(M^3)_{vv}\) (structural counts);
  4. + both.

Predict first. Does neighbourhood averaging help find bots? Do the two structural counts?

_Xs = (X_feat - X_feat.mean(0)) / X_feat.std(0)
_deg = np.asarray(U.sum(1)).ravel()
_mean_nb = (U @ _Xs) / np.maximum(_deg, 1)[:, None]
_out = np.asarray(A.sum(1)).ravel()
recip = np.asarray(M.sum(1)).ravel() / np.maximum(_out, 1)
tri = 0.5 * (M @ M @ M).diagonal()
_struct = np.column_stack([recip, np.log1p(tri)])
feature_sets = {
    "1 · per-account features": _Xs,
    "2 · + neighbour means (≈ 1-layer mean GNN)": np.hstack([_Xs, _mean_nb]),
    "3 · + reciprocity + mutual triangles": np.hstack([_Xs, _struct]),
    "4 · all of the above": np.hstack([_Xs, _mean_nb, _struct]),
}

def cv_auc(X, y, seeds=(0, 1, 2)):
    """Mean ± std of per-fold test AUCs over 3 seeds × 5 folds (per-fold, not pooled: pooling
    out-of-fold scores from different models biases AUC)."""
    out = []
    for s in seeds:
        for tr, te in StratifiedKFold(5, shuffle=True, random_state=s).split(X, y):
            clf = LogisticRegression(max_iter=4000, class_weight="balanced").fit(X[tr], y[tr])
            out.append(roc_auc_score(y[te], clf.predict_proba(X[te])[:, 1]))
    return np.mean(out), np.std(out)

bot_table = pd.DataFrame(
    [{"features": k, "bot AUC (mean ± std, 3 seeds × 5 folds)": "{:.3f} ± {:.3f}".format(*cv_auc(v, is_bot))}
     for k, v in feature_sets.items()]
)
_fig, _ax = plt.subplots(1, 2, figsize=(8, 2.6))
for _a, (_v, _lab) in zip(_ax, [(recip, "reciprocity"), (np.log1p(tri), "log(1 + mutual triangles)")], strict=True):
    _a.hist(_v[~is_bot], bins=30, alpha=0.6, density=True, label="others")
    _a.hist(_v[is_bot], bins=30, alpha=0.6, density=True, label="ring members")
    _a.set_title(_lab, fontsize=9)
_ax[0].legend(fontsize=7)
plt.tight_layout()
mo.vstack([bot_table, _fig])

_g = nx.from_scipy_sparse_array(M)
_g.remove_nodes_from([v for v in list(_g) if _g.degree(v) == 0])
_pos = nx.spring_layout(_g, seed=1, k=0.15, iterations=60)
_fig, _ax = plt.subplots(figsize=(7.5, 4.2))
nx.draw_networkx_edges(_g, _pos, ax=_ax, alpha=0.25, width=0.6)
_nodes = list(_g)
nx.draw_networkx_nodes(_g, _pos, nodelist=_nodes, ax=_ax, node_size=[22 if is_bot[v] else 8 for v in _nodes],
                       node_color=["#c0392b" if is_bot[v] else "#7f8c8d" for v in _nodes], linewidths=0)
_ax.set_title(f"mutual-retweet graph: {_g.number_of_nodes()} accounts with a mutual tie; ring members in red", fontsize=9)
_ax.axis("off")
_fig
_v = [float(s.split(" ± ")[0]) for s in bot_table.iloc[:, 1]] · 16 lines
_v = [float(s.split(" ± ")[0]) for s in bot_table.iloc[:, 1]]
_better = "neighbour means" if _v[1] > _v[2] else "structural counts"
_share = (np.asarray(M.sum(1)).ravel()[is_bot] > 0).mean()
mo.md(rf"""
**This run:** features {_v[0]:.3f} → + neighbour means {_v[1]:.3f} → + structural counts {_v[2]:.3f} → all {_v[3]:.3f}.
The bigger single gain comes from the **{_better}**, and the two together are best.

*Why neighbour means help* (what a 1-layer mean GNN sees): a ring member's neighbours are other ring members and a handful of
amplified targets, so the *average neighbour* is unusual even when the account itself is not. The ring is visible through
**whom it touches**. *Why the counts help*: rings are reciprocal and dense, and those are **counts** of dyads and cycles, which
mean aggregation cannot produce and 1-WL cannot see. Their limit here: only {_share:.0%} of ring members have any mutual tie in
this world (ring posts are few), and counts cannot see the rest.

The practical content of L06 + L07: the two signals are complementary, so a detector should get both, either as features or with
a subgraph-aware model. Before concluding "the GNN found no bots", add the counts it cannot compute.
""")

L11 · Who-to-retweet recommendation and popularity bias

Task 4: recommend accounts to retweet. Train on days 0–20, evaluate Recall@20 on the authors each user newly retweets in days 20–30. Popularity baseline vs. LightGCN-style propagation (no training: score = 2-step propagation on the normalised user–author bipartite graph, the "K = 2, untrained" special case of L11). We also measure how concentrated the recommendations are (share of all recommendations that go to the top 1% of authors).

Predict first. Which wins on recall, and which concentrates attention more on already-popular accounts?

_cut = 20 * 24
_past = world.retweets[world.retweets.time <= _cut]
_fut = world.retweets[world.retweets.time > _cut]
R = sp.csr_matrix((np.ones(len(_past)), (_past.retweeter, _past.author)), shape=(n, n))
R.data[:] = 1.0
_seen = {(int(a), int(b)) for a, b in zip(_past.retweeter, _past.author, strict=True)}
_truth = {}
for a, b in zip(_fut.retweeter, _fut.author, strict=True):
    if (int(a), int(b)) not in _seen and a != b:
        _truth.setdefault(int(a), set()).add(int(b))
_du, _da = np.asarray(R.sum(1)).ravel(), np.asarray(R.sum(0)).ravel()
_Rn = sp.diags(1 / np.sqrt(np.maximum(_du, 1))) @ R @ sp.diags(1 / np.sqrt(np.maximum(_da, 1)))
_S_lgcn = (_Rn @ _Rn.T @ _Rn).toarray()  # user -> author scores after user-author-user-author propagation
_pop = np.tile(_da, (n, 1)).astype(float)

def _evaluate(S, K=20):
    S = S.copy()
    S[R.nonzero()] = -np.inf  # do not recommend what they already retweet
    np.fill_diagonal(S, -np.inf)
    users_ = [u for u in _truth if _du[u] > 0]
    top = np.argsort(-S[users_], axis=1)[:, :K]
    rec = np.mean([len(set(t) & _truth[u]) / len(_truth[u]) for u, t in zip(users_, top, strict=True)])
    counts = np.bincount(top.ravel(), minlength=n)
    top1 = np.argsort(-_da)[: max(1, n // 100)]
    return rec, counts[top1].sum() / counts.sum()

rec_table = pd.DataFrame(
    [{"model": m, "Recall@20": r, "share of recs to top-1% authors": c}
     for m, (r, c) in [("popularity", _evaluate(_pop)), ("LightGCN-style propagation (K = 2, untrained)", _evaluate(_S_lgcn))]]
).round(3)
rec_table
_p, _l = rec_table.iloc[0], rec_table.iloc[1] · 7 lines
_p, _l = rec_table.iloc[0], rec_table.iloc[1]
mo.md(rf"""
Popularity: Recall@20 {_p["Recall@20"]:.3f}, with {_p["share of recs to top-1% authors"]:.0%} of recommendations going to the
top 1% of accounts. Propagation: {_l["Recall@20"]:.3f} and {_l["share of recs to top-1% authors"]:.0%}. A recommender is an
edge-rewiring operator on the social graph (L11 frontier): the concentration column is a first measure of what it does to
attention, and the one regulators ask about.
""")

L14 · Trust: conformal prediction sets for "is this a ring member?"

A moderator needs an error rate, not a score. Split conformal (L14 proof workshop): calibrate the bot classifier from §L06 (feature set 3) on half the accounts, then report, on the other half, prediction sets with a 95% coverage guarantee.

Predict first. Will empirical coverage be near 0.95? What fraction of accounts get the uninformative set {human, bot}?

_Xs = (X_feat - X_feat.mean(0)) / X_feat.std(0)
_X = np.hstack([_Xs, np.column_stack([recip, np.log1p(tri)])])
_y = is_bot.astype(int)
_rows = []
for _seed in range(20):
    _r = np.random.default_rng(_seed)
    _idx = _r.permutation(len(_y))
    _tr, _cal, _te = np.split(_idx, [len(_y) // 2, 3 * len(_y) // 4])
    _clf = LogisticRegression(max_iter=4000, class_weight="balanced").fit(_X[_tr], _y[_tr])
    _s_cal = 1 - _clf.predict_proba(_X[_cal])[np.arange(len(_cal)), _y[_cal]]
    _q = np.quantile(_s_cal, np.ceil((len(_cal) + 1) * 0.95) / len(_cal), method="higher")
    _P = _clf.predict_proba(_X[_te])
    _sets = (1 - _P) <= _q
    _rows.append({"coverage": _sets[np.arange(len(_te)), _y[_te]].mean(), "share {human, bot}": _sets.all(1).mean(),
                  "bots caught as {bot}": (_sets[_y[_te] == 1] == [False, True]).all(1).mean() if (_y[_te] == 1).any() else np.nan})
conf_table = pd.DataFrame(_rows).mean().round(3).to_frame("mean over 20 splits")
conf_table
_c = conf_table.iloc[:, 0] · 7 lines
_c = conf_table.iloc[:, 0]
mo.md(rf"""
Coverage {_c["coverage"]:.3f} (target 0.95), as the theorem promises whatever the model. {_c["share {human, bot}"]:.0%} of accounts get
the honest answer "could be either", and {_c["bots caught as {bot}"]:.0%} of ring members are flagged as bots alone. A better
model makes the sets smaller, never less valid. The caveat for social data: exchangeability fails when rings are split across
calibration and test non-randomly, or when tomorrow's bots differ from today's.
""")

L20 · Few labels: label propagation from a handful of known ring members

A platform usually knows a few confirmed accounts. Seed label propagation (L20) with 3 known ring members and 30 known humans on (a) the undirected retweet graph and (b) the mutual graph, where an edge needs retweets both ways.

Predict first. Which graph propagates the bot label better, and why?

def _lp(Ad, seeds_pos, seeds_neg, alpha=0.9, iters=50):
    d = np.asarray(Ad.sum(1)).ravel()
    Dm = sp.diags(1 / np.sqrt(np.maximum(d, 1)))
    S = Dm @ Ad @ Dm
    Y = np.zeros((n, 2))
    Y[seeds_neg, 0] = 1
    Y[seeds_pos, 1] = 1
    Fm = Y.copy()
    for _ in range(iters):
        Fm = alpha * (S @ Fm) + (1 - alpha) * Y
    return Fm[:, 1] - Fm[:, 0]

_rows = []
for _seed in range(10):
    _r = np.random.default_rng(_seed)
    _pos = _r.choice(np.flatnonzero(is_bot), 3, replace=False)
    _neg = _r.choice(np.flatnonzero(~is_bot), 30, replace=False)
    _mask = np.ones(n, bool)
    _mask[np.r_[_pos, _neg]] = False
    for _name, _Ad in [("undirected retweet graph", U), ("mutual graph", M)]:
        _rows.append({"graph": _name, "AUC": roc_auc_score(is_bot[_mask], _lp(_Ad, _pos, _neg)[_mask])})
lp_table = pd.DataFrame(_rows).groupby("graph").AUC.agg(["mean", "std"]).round(3)
lp_table
_m, _u = lp_table.loc["mutual graph", "mean"], lp_table.loc["undirected retweet graph", "m · 10 lines
_m, _u = lp_table.loc["mutual graph", "mean"], lp_table.loc["undirected retweet graph", "mean"]
_reach = (np.asarray(M.sum(1)).ravel()[is_bot] > 0).mean()
mo.md(rf"""
Mutual graph AUC {_m:.3f} vs undirected {_u:.3f}. Two forces pull in opposite directions. The mutual graph is more
**homophilous** for the bot label (an edge needs coordination both ways), but it is **sparse**: only {_reach:.0%} of ring members
have a mutual tie at all, and labels cannot flow to the others, which score like isolated humans. The undirected graph reaches
everyone through the targets they all retweet, at the price of camouflage edges. Here the {"undirected graph's reach wins" if _u > _m else "mutual graph's homophily wins"}.
Choosing the graph is choosing the homophily assumption *and* the reach: the representation decision of L01, made with the
question in mind. (A natural next step: propagate on a weighted mix, $U+\lambda M$.)
""")

L21 · Motifs: what is over-represented, against which null?

Count reciprocal dyads and triangles in the mutual graph, and compare with a degree-preserving randomisation of the directed retweet graph (double-edge swaps keep every in- and out-degree). The Z-score says how many null standard deviations the observed count is away from the null mean.

def _counts(g):
    mutual = sum(1 for u, v in g.edges() if u < v and g.has_edge(v, u))
    m = nx.Graph((u, v) for u, v in g.edges() if g.has_edge(v, u))
    tri = sum(nx.triangles(m).values()) // 3
    return mutual, tri

def _swap(g, n_swaps, seed):
    r = np.random.default_rng(seed)
    e = np.array(g.edges())
    for _ in range(n_swaps):
        i, j = r.integers(len(e), size=2)
        (a, b), (c, d) = e[i], e[j]
        if len({a, b, c, d}) == 4:
            e[i], e[j] = (a, d), (c, b)
    h = nx.DiGraph()
    h.add_edges_from(map(tuple, e))
    return h

_obs = _counts(G)
_null = np.array([_counts(_swap(G, 3 * G.number_of_edges(), s)) for s in range(5)])
_bots_only = _counts(G.subgraph(np.flatnonzero(is_bot)))
motif_table = pd.DataFrame(
    {"observed": _obs, "null mean": _null.mean(0), "null std": _null.std(0) + 1e-9,
     "inside rings": _bots_only},
    index=["reciprocal dyads", "triangles in the mutual graph"],
)
motif_table["Z"] = (motif_table.observed - motif_table["null mean"]) / motif_table["null std"]
motif_table.round(1)
_t = motif_table.loc["triangles in the mutual graph"] · 38 lines
_t = motif_table.loc["triangles in the mutual graph"]
mo.md(rf"""
Mutual triangles: observed {int(_t.observed)}, null {_t["null mean"]:.0f} ± {_t["null std"]:.0f} (Z ≈ {_t.Z:.0f}); {int(_t["inside rings"])}
of them are inside rings. The motif signal is real and mostly *is* the rings. On real data this is the first plot to make: it
tells you whether there is any coordination structure to find before you train anything.

---

## Swap in real data, and what to ask next

```python
from gl.retweet import load_retweets_csv
world = load_retweets_csv("data/my_retweets.csv")   # columns: retweeter, author, time[, hashtag]
```
Labels (`is_bot`, `community`) are unknown on real data: use §L21 (motif Z-scores) and §L20 (propagation from a few confirmed
accounts) first, then validate on any labelled subset you have.

A real retweet log ships with the repo's dataset catalogue: `gl.datasets.higgs_world()` loads the Higgs boson Twitter
activity (1–7 July 2012, De Domenico et al. 2013) as a `World`, so every section above runs on it unchanged except
the ones that need labels. The cell below loads it and computes two label-free signals of §L21. For bot-detection
labels, the catalogue (`data/README.md`) lists the benchmarks that need an application to their authors (TwiBot-22,
TwiBot-20, MGTAB) and how to get them.

**Research extensions.** Each section of this notebook has a pre-registered experiment in `experiments/` that asks
the open question behind it (hypotheses written before running, verdicts computed, refuted ones kept):

| section | experiment | question |
|---|---|---|
| L02–L04 | `experiments/rings-in-communities` | do rings look like tiny communities, or hide inside real ones? |
| L05 | `experiments/honest-splits` | how much does a random split overestimate a bot detector? (+ 5 real temporal networks) |
| L06–L07 | `experiments/rege-vs-wl` | regular equivalence (sets) vs WL (multisets) for finding rings |
| L11 | `experiments/recsys-amplification` | does the recommender amplify the rings, and by how much? |
| L14 | `experiments/conformal-shift` | do conformal guarantees hold when bots change tactics? (+ Twitch countries) |
| L19 | `experiments/design-transfer` | are Cora's best GNN designs also best for bot detection? (+ 4 real social graphs) |
| L20–L21 | `experiments/motif-signatures` | which small patterns give rings away? |

`experiments/README.md` has the findings in one table.
""")

higgs_run = mo.ui.run_button(label="Load the Higgs retweet log (downloads 4 MB once; ≈ 1 min; then cached)")
higgs_run

from gl.datasets import higgs_world
from gl.exp import adjacency

_cache = mo.notebook_dir() / "cache" / "higgs_summary.json"
mo.stop(
    not (higgs_run.value or os.environ.get("GL_RUN_EXPERIMENTS") == "1") and not _cache.exists(),
    mo.md("Press the button. A summary is cached in `notebooks/project/cache/higgs_summary.json`."),
)
if not _cache.exists():
    _w = higgs_world()
    _A, _U, _M = adjacency(_w)
    _tri = 0.5 * (_M @ _M @ _M).diagonal()
    _per_hour = _w.retweets.time.floordiv(1).value_counts().sort_index()
    _cache.parent.mkdir(exist_ok=True)
    _cache.write_text(json.dumps({
        "users": len(_w.users), "retweets": len(_w.retweets), "hours": float(_w.retweets.time.max()),
        "per_hour": {int(k): int(v) for k, v in _per_hour.items()},
        "retweet_pairs": int(_A.nnz), "mutual_pairs": int(_M.nnz // 2),
        "accounts_mutual": int((np.asarray(_M.sum(1)).ravel() > 0).sum()),
        "accounts_in_mutual_triangles": int((_tri > 0).sum()),
        "top_author_share": float(_w.retweets.author.value_counts().iloc[0] / len(_w.retweets)),
    }))
higgs = json.loads(_cache.read_text())
_ph = pd.Series({int(k): v for k, v in higgs["per_hour"].items()}).sort_index()
_fig, _ax = plt.subplots(figsize=(8, 2.4))
_ax.plot(_ph.index / 24, _ph.values, lw=0.8)
_ax.set_xlabel("days since 1 July 2012, 00:00 UTC (approx.)")
_ax.set_ylabel("retweets per hour")
_ax.set_title("Higgs boson retweets: one burst around the 4 July announcement", fontsize=9)
plt.tight_layout()
mo.vstack([
    _fig,
    mo.md(f"""
**{higgs["users"]:,} accounts, {higgs["retweets"]:,} retweets over {higgs["hours"] / 24:.1f} days**, peaking at
{_ph.max():,} retweets in one hour. The single most-retweeted account receives {higgs["top_author_share"]:.1%} of all
retweets. Of {higgs["retweet_pairs"]:,} retweeter → author pairs only {higgs["mutual_pairs"]:,} are mutual;
{higgs["accounts_mutual"]:,} accounts have a mutual tie and {higgs["accounts_in_mutual_triangles"]:,} sit in a mutual
triangle, the label-free pattern §L21 found almost only inside rings. On real data such accounts are a *shortlist* to
inspect, not a list of bots: friends, colleagues and fan accounts retweet each other too.
"""),
])