Local Models and Interoperability

1Tutorial: Local Models and Tool Interoperability¶

This tutorial is about the half of the AI toolchain a statistician can actually own. A hosted endpoint is a moving target: it may be updated between your submission and the referee report, and the data you send it leave your control. A model you run yourself has neither problem, at the price of three subjects you could previously ignore --- memory, precision and throughput --- and one you could not previously have: the plumbing that lets an agent reach your tools at all. The first half covers the runtimes, in increasing order of how much machinery they impose; the second covers interoperability, which is what turns a model on your machine into something the rest of your workflow can call. Section that section develops the memory and quantization arithmetic in detail; the treatment here is operational. Capabilities move quickly, so the durable content is the set of distinctions --- single-stream versus batching, weights-file format, whether the runtime speaks a common API --- rather than any snapshot of features.

1.1Ollama: The Managed Local Runtime¶

• [Ollama] Local model manager. Pulls open-weight models, manages the on-disk store, and serves them behind a local HTTP API with sensible defaults. Fits: first contact with local models, single-user interactive work, and scripts that want one endpoint on a laptop. Watch: a bare model tag is not a pinned version; the same tag can resolve to a different quantization or a re-uploaded file, so record the digest and the quantization, not the friendly name.

The reason to start here is that it removes every decision except which model to run, and the reason to leave is that some of those decisions are yours to make. A typical session pulls a model and then talks to it over the local API.

ollama pull <model>:<tag>
ollama list                 # names, digests and sizes of what you have
ollama show <model>:<tag>   # architecture, parameter count, quantization
ollama run <model>:<tag>    # interactive; Ctrl-D exits
curl http://localhost:11434/api/generate \
  -d '{"model":"<model>:<tag>","prompt":"...","stream":false,"options":{"seed":0}}'

For a paper, the output of ollama show is the part that matters: it names the quantization scheme and the architecture, which are the two facts a reader needs and the friendly tag hides.

1.2llama.cpp: The Portable Single-Stream Engine¶

• [llama.cpp] Portable C/C++ inference engine. Runs quantized models in the GGUF format on CPU, consumer GPUs and Apple silicon, with a small dependency footprint and an OpenAI-compatible server mode. Fits: machines without a data-centre GPU, air-gapped or governance-restricted environments, and embedding a model inside another program. Watch: it optimises latency for one stream at a time; throughput on many documents is not what it is built for, and different quantization variants of the same model are genuinely different models.

The relevant structural fact is the file format. GGUF is a single-file container holding weights, tokenizer and metadata, which is why a GGUF model is portable and why it is not interchangeable with the multi-file layouts the research runtimes use. Conversion is one direction and one script; plan for it rather than discovering it.

python convert_hf_to_gguf.py /path/to/hf-model --outfile model-f16.gguf
llama-quantize model-f16.gguf model-q4.gguf Q4_K_M
llama-server -m model-q4.gguf --ctx-size 8192 --port 8080

1.3vLLM: The Batching Server¶

• [vLLM] High-throughput serving engine. Serves open-weight models to many concurrent requests using paged attention and continuous batching to keep the accelerator busy Kwon et al., 2023. Fits: scoring or annotating a corpus --- thousands of documents through one model --- and shared group or cluster deployments. Watch: throughput comes from batching, and batching changes reduction order, so a given prompt need not produce a bitwise identical continuation across runs even at temperature zero.

Why a separate engine exists at all is an arithmetic fact rather than a software preference, and it is worth computing. A single decode step reads the entire weight matrix from memory and does very little work with it; batching amortises that read across sequences until the arithmetic, rather than the memory traffic, becomes the constraint.

# Roofline arithmetic for one decode step. All numbers are hardware you declare.
BW_GBS = 900.0        # memory bandwidth of your accelerator, GB/s
TFLOPS = 120.0        # dense throughput at the precision you are running
W_GB = 16.0           # weight footprint you computed above

print(f"{'batch':>6} {'bytes moved':>12} {'FLOPs':>10} {'ms/step (bw)':>13} "
      f"{'ms/step (fl)':>13} {'bound by':>9} {'tok/s':>8}")
for b in (1, 4, 16, 64, 256):
    bytes_moved = W_GB                     # weights stream once per step
    gflops = 2 * (W_GB * 1e9 / 2) * b / 1e9   # 2 FLOPs per fp16 weight per seq
    t_bw = bytes_moved / BW_GBS * 1000
    t_fl = gflops / (TFLOPS * 1000) * 1000
    t = max(t_bw, t_fl)
    print(f"{b:6d} {bytes_moved:10.1f} GB {gflops:8.0f} G {t_bw:13.2f} "
          f"{t_fl:13.2f} {'memory' if t_bw > t_fl else 'compute':>9} "
          f"{b / t * 1000:8.0f}")

 batch  bytes moved      FLOPs  ms/step (bw)  ms/step (fl)  bound by    tok/s
     1       16.0 GB       16 G         17.78          0.13    memory       56
     4       16.0 GB       64 G         17.78          0.53    memory      225
    16       16.0 GB      256 G         17.78          2.13    memory      900
    64       16.0 GB     1024 G         17.78          8.53    memory     3600
   256       16.0 GB     4096 G         17.78         34.13   compute     7500

Between batch one and batch sixty-four the throughput rises by a factor of sixty-four at no cost in latency, because the step time is set by streaming the weights and that happens once regardless. This is why a corpus-annotation job belongs on a batching server and an interactive session does not care.

The second half of the memory question is the cache, which grows with context and with batch and is the usual reason a model that loads still fails.

GIB = 2 ** 30

def weights_gib(n_params_b, bits):
    return n_params_b * 1e9 * bits / 8 / GIB

def kv_gib(layers, kv_heads, head_dim, ctx, batch, bits=16):
    return 2 * layers * kv_heads * head_dim * ctx * batch * bits / 8 / GIB

# a mid-sized open-weight decoder; substitute your own config from config.json
CFG = dict(layers=32, kv_heads=8, head_dim=128, params_b=8.0)
VRAM = 24.0

print(f"{'bits':>5} {'weights':>9} {'ctx':>7} {'batch':>6} {'KV':>7} "
      f"{'total':>7}  fits in {VRAM:.0f} GiB?")
for bits in (16, 8, 4):
    w = weights_gib(CFG["params_b"], bits)
    for ctx, batch in ((4096, 1), (32768, 1), (4096, 16), (32768, 16)):
        kv = kv_gib(CFG["layers"], CFG["kv_heads"], CFG["head_dim"], ctx, batch)
        tot = (w + kv) * 1.15                      # activations + fragmentation
        print(f"{bits:5d} {w:9.2f} {ctx:7d} {batch:6d} {kv:7.2f} {tot:7.2f}  "
              f"{'yes' if tot < VRAM else 'NO'}")

 bits   weights     ctx  batch      KV   total  fits in 24 GiB?
   16     14.90    4096      1    0.50   17.71  yes
   16     14.90   32768      1    4.00   21.74  yes
   16     14.90    4096     16    8.00   26.34  NO
   16     14.90   32768     16   64.00   90.74  NO
    8      7.45    4096      1    0.50    9.14  yes
    8      7.45   32768      1    4.00   13.17  yes
    8      7.45    4096     16    8.00   17.77  yes
    8      7.45   32768     16   64.00   82.17  NO
    4      3.73    4096      1    0.50    4.86  yes
    4      3.73   32768      1    4.00    8.88  yes
    4      3.73    4096     16    8.00   13.48  yes
    4      3.73   32768     16   64.00   77.88  NO

Quantizing to four bits frees eleven gigabytes of weights, and a long-context batch of sixteen consumes sixty-four in cache --- so the quantization that made the model fit does not make the workload fit. Compute this table before provisioning anything.

1.4Hugging Face Transformers: The Research Runtime¶

• [Hugging Face Transformers] Model library and hub client. Provides reference implementations of published architectures with a common loading interface, backed by a versioned repository of weights Wolf et al., 2020. Fits: anything that needs the model internals --- log-probabilities, hidden states, attention, custom heads, fine-tuning --- and any pipeline where exact provenance of the weights matters. Watch: from_pretrained("name") follows a moving branch. Pass an explicit revision commit hash, or your “fixed” model is not fixed.

This is the runtime a statistician usually wants, because it is the only one of the four that hands you the quantities you would want to do inference with. A scoring loop that returns per-token log-probabilities is a few lines, and those log-probabilities are the raw material for calibration analysis, for likelihood-ratio comparisons between prompts, and for the uncertainty summaries covered in \Sthat section. The library also exposes the tokenizer as a first-class object, which matters more than it sounds.

A token count is a measurement made by an instrument, and swapping the instrument changes the number without changing the text. The example trains byte-pair encoding on a small statistical corpus and shows the same sentence measured by two vocabularies.

from collections import Counter

CORPUS = ("the estimate of the treatment effect is unbiased "
          "the estimator of the effect is consistent and unbiased "
          "estimation of treatment effects requires the estimator ") * 12

def train_bpe(text, n_merges):
    words = [tuple(w) + ("_",) for w in text.split()]
    merges = []
    for _ in range(n_merges):
        pairs = Counter()
        for w in words:
            for i in range(len(w) - 1):
                pairs[w[i], w[i + 1]] += 1
        if not pairs:
            break
        best = max(pairs.items(), key=lambda kv: (kv[1], kv[0]))[0]
        merges.append(best)
        words = [merge_word(w, best) for w in words]
    return merges

def merge_word(w, pair):
    out, i = [], 0
    while i < len(w):
        if i + 1 < len(w) and (w[i], w[i + 1]) == pair:
            out.append(w[i] + w[i + 1]); i += 2
        else:
            out.append(w[i]); i += 1
    return tuple(out)

def encode(text, merges):
    toks = []
    for w in text.split():
        p = tuple(w) + ("_",)
        for m in merges:
            p = merge_word(p, m)
        toks.extend(p)
    return toks

TEXT = "the estimator of the treatment effect is unbiased"
print(f"{'merges':>7} {'vocab':>6} {'tokens':>7}  segmentation of the test sentence")
for n in (0, 20, 60, 200):
    m = train_bpe(CORPUS, n)
    t = encode(TEXT, m)
    print(f"{n:7d} {len(set(encode(CORPUS, m))):6d} {len(t):7d}  "
          + " ".join(t[:9]) + (" ..." if len(t) > 9 else ""))

a, b = train_bpe(CORPUS, 20), train_bpe(CORPUS, 200)
na, nb = len(encode(TEXT, a)), len(encode(TEXT, b))
print(f"\nsame text, two tokenizers: {na} vs {nb} tokens, ratio {na / nb:.2f}")

 merges  vocab  tokens  segmentation of the test sentence
      0     17      50  t h e _ e s t i m ...
     20     26      25  the_ estimat o r _ of_ the_ t re ...
     60     13       8  the_ estimator_ of_ the_ treatment_ effect_ is_ unbiased_
    200     13       8  the_ estimator_ of_ the_ treatment_ effect_ is_ unbiased_

same text, two tokenizers: 25 vs 8 tokens, ratio 3.12

A factor of three in the same sentence, from the instrument alone. Any cost projection, context budget or throughput figure is therefore conditional on a tokenizer, and comparing two models by token count without saying whose tokenizer is a category error.

The reproducibility story has a second wrinkle, and it is the one that surprises people who fix a seed and expect determinism. Floating-point addition is not associative, so changing the order in which partial sums are accumulated --- which is exactly what changing batch size or kernel does --- perturbs the logits. The perturbation is tiny; it matters only when the top two logits are close.

import numpy as np

rng = np.random.default_rng(0)
V, D, STEPS = 4000, 1024, 4000
W = rng.normal(0, 0.02, (V, D)).astype(np.float32)

def logits(h, chunk):
    """Identical dot products, different summation order -- what batching alters."""
    out = np.zeros(V, dtype=np.float32)
    for s in range(0, D, chunk):
        out += W[:, s:s + chunk] @ h[s:s + chunk]
    return out

gaps, errs = [], []
for _ in range(STEPS):
    h = rng.normal(0, 1, D).astype(np.float32)
    a, b = logits(h, 1024), logits(h, 64)
    errs.append(np.abs(a - b).max())
    top = np.partition(a, -2)[-2:]
    gaps.append(abs(top[1] - top[0]))
gaps, errs = np.sort(np.array(gaps)), np.array(errs)
eps = np.median(errs)

# no flip was observed directly, so estimate the rate from the small-gap density
f0 = 20 / (STEPS * gaps[19])          # 20 smallest gaps give the density at 0
p = f0 * eps
print(f"decode steps simulated              {STEPS}")
print(f"median numeric discrepancy eps      {eps:.2e}")
print(f"median top-1 / top-2 gap            {np.median(gaps):.4f}")
print(f"smallest gap observed               {gaps[0]:.2e}")
print(f"direct flips observed               {int((gaps < errs).sum())}")
print(f"density of the gap at 0             {f0:.3f}")
print(f"implied per-token flip probability  {p:.2e}")
for n in (1_000, 100_000, 10_000_000):
    print(f"  P(>=1 differing token in {n:>10,}) = {1 - (1 - p) ** n:.4f}")

decode steps simulated              4000
median numeric discrepancy eps      6.85e-07
median top-1 / top-2 gap            0.1196
smallest gap observed               1.31e-04
direct flips observed               0
density of the gap at 0             5.898
implied per-token flip probability  4.04e-06
  P(>=1 differing token in      1,000) = 0.0040
  P(>=1 differing token in    100,000) = 0.3332
  P(>=1 differing token in 10,000,000) = 1.0000

Per token the flip probability is negligible; over a corpus-scale run it is a near certainty. This is a rare-event calculation, not a bug, and it sets the honest claim: greedy decoding with a pinned model gives distributional reproducibility, not bitwise reproducibility, and a methods section should say so. Note also the estimation move --- zero events observed, so the rate is recovered from the density of the gap near zero rather than by counting.

Quantization interacts with this in a way that a likelihood-style check will under-report. Treating the precision loss as noise injected into the logits, the summary statistic moves far less than the decisions do.

import numpy as np

rng = np.random.default_rng(2)
N, K = 4000, 5                       # documents, label classes

z = rng.normal(0, 1.4, (N, K))       # full-precision logits for a labelling task
truth = z.argmax(1)

def perturb(z, sd):
    """Quantization acts like measurement error injected into the logits."""
    return z + rng.normal(0, sd, z.shape)

def nll(z, y):
    m = z.max(1, keepdims=True)
    logp = z - m - np.log(np.exp(z - m).sum(1, keepdims=True))
    return -logp[np.arange(len(y)), y].mean()

print(f"{'noise sd':>9} {'mean NLL':>9} {'NLL ratio':>10} "
      f"{'label agreement':>16} {'flipped':>8}")
base = nll(z, truth)
for sd in (0.0, 0.02, 0.05, 0.10, 0.20):
    zq = perturb(z, sd) if sd else z
    agree = (zq.argmax(1) == truth).mean()
    print(f"{sd:9.2f} {nll(zq, truth):9.4f} {nll(zq, truth)/base:10.3f} "
          f"{agree:16.3f} {int((1 - agree) * N):8d}")

 noise sd  mean NLL  NLL ratio  label agreement  flipped
     0.00    0.6475      1.000            1.000        0
     0.02    0.6479      1.001            0.989       44
     0.05    0.6486      1.002            0.982       73
     0.10    0.6498      1.004            0.955      181
     0.20    0.6659      1.028            0.901      396

A perplexity check that moves by three percent conceals a labelling disagreement of ten. If the model is an instrument in your study, validate it on the task you are using it for; a global likelihood summary is the wrong endpoint. The mechanism behind the perturbation --- and why the rounding error is not iid --- is developed in \Sthat section.

1.5Model Context Protocol: One Interface, Many Clients¶

• [Model Context Protocol (MCP)] Tool-connection protocol. A standard interface by which an agent discovers and calls tools, reads resources and loads prompts, so a capability implemented once is reachable from many different agents. Fits: exposing a lab database, a cluster submitter or an internal API to whatever agent your group uses, without rewriting an adapter each time. Watch: a server is code you granted an agent permission to run; treat installing one as you treat installing a package, and see \Sthat section for the supply-chain surface.

Keep the vocabulary straight, because these two are constantly conflated: MCP is how an agent connects; a skill is what an agent knows. A server executes and returns a result; a skill injects a procedure. Both are treated in depth in \Sthat section.

The engineering case is the same combinatorial argument that motivates any standard interface, and the statistical case is validation: a tool call is a typed record with a schema, so malformed calls are detectable before execution.

import json

# One capability, described once, in a client-agnostic schema.
TOOL = {
    "name": "cohort_query",
    "description": "Return de-identified counts for a cohort definition.",
    "input_schema": {
        "type": "object",
        "properties": {"icd10": {"type": "string"},
                       "min_age": {"type": "integer"}},
        "required": ["icd10"],
    },
}

def validate(call, schema):
    props, req = schema["properties"], schema["required"]
    missing = [k for k in req if k not in call]
    unknown = [k for k in call if k not in props]
    T = {"string": str, "integer": int}
    bad = [k for k, v in call.items()
           if k in props and not isinstance(v, T[props[k]["type"]])]
    return missing, unknown, bad

CALLS = [{"icd10": "E11", "min_age": 40},
         {"icd10": "E11"},
         {"min_age": 40},
         {"icd10": "E11", "min_age": "forty"},
         {"icd10": "E11", "minage": 40}]
for c in CALLS:
    m, u, b = validate(c, TOOL["input_schema"])
    tag = "ok" if not (m or u or b) else f"missing={m} unknown={u} type={b}"
    print(f"{json.dumps(c):<40} {tag}")

print()
for a, t in ((3, 4), (5, 10), (10, 25)):
    print(f"{a:2d} agents x {t:2d} tools: {a*t:3d} bespoke adapters, "
          f"or {a + t:3d} with a shared protocol "
          f"(factor {a*t/(a+t):.1f})")

{"icd10": "E11", "min_age": 40}          ok
{"icd10": "E11"}                         ok
{"min_age": 40}                          missing=['icd10'] unknown=[] type=[]
{"icd10": "E11", "min_age": "forty"}     missing=[] unknown=[] type=['min_age']
{"icd10": "E11", "minage": 40}           missing=[] unknown=['minage'] type=[]

 3 agents x  4 tools:  12 bespoke adapters, or   7 with a shared protocol (factor 1.7)
 5 agents x 10 tools:  50 bespoke adapters, or  15 with a shared protocol (factor 3.3)
10 agents x 25 tools: 250 bespoke adapters, or  35 with a shared protocol (factor 7.1)

Every failure mode here is a syntax failure, and every one is caught without running anything. What the validator cannot tell you is whether E11 was the right code, which is the general shape of the guarantee: schemas buy well-formedness, never correctness. That is the same asymmetry as in \Sthat section --- cheap verifiers make agentic work viable, and they only ever verify what they can see.

1.6Instruction Files: What the Agent Knows¶

Two file conventions carry procedure rather than connection, and the distinction between them is scope, not format. A project context file --- AGENTS.md and its equivalents --- describes this repository: how to build it, how to run the tests, which conventions the group uses, where the data live. A skill is a directory containing a SKILL.md instruction file plus optional scripts and references, and it is a portable procedure usable in any project, loaded on demand when a task matches its description. Both are covered properly in \Sthat section, including the retrieval framing of a skill library and the security surface of loading third-party instruction text.

For local-model work specifically, the relevant point is that a skill is where the pinning discipline belongs. The revision hash, the quantization, the decoding parameters and the seed policy are procedure, not preference, and writing them into a skill is how they survive the next person who runs your pipeline.

1.7The Statistical Payoff: A Local Model as a Cheap Instrument¶

The reason to assemble any of this is that a local model can label far more data than you can afford to label by hand, and a labelled subsample turns those cheap labels into valid inference rather than into a hopeful assertion. Prediction-powered inference Angelopoulos et al., 2023 is the clean statement of the idea: use the model everywhere, use the audited subsample to estimate and remove its bias.

import numpy as np

rng = np.random.default_rng(4)
N, n = 20_000, 400          # unlabelled notes; hand-labelled subsample

truth = rng.random(N) < 0.22                       # true prevalence 0.22
sens, spec = 0.88, 0.94                            # the local model's behaviour
pred = np.where(truth, rng.random(N) < sens,
                rng.random(N) < 1 - spec).astype(float)

idx = rng.choice(N, n, replace=False)              # the audited subsample
y, f = truth[idx].astype(float), pred[idx]

naive = pred.mean()                                 # trust the model everywhere
classical = y.mean()                                # ignore the model
se_cl = np.sqrt(classical * (1 - classical) / n)

rect = (f - y).mean()                               # measured model bias
ppi = pred.mean() - rect
se_ppi = np.sqrt((f - y).var(ddof=1) / n + pred.var(ddof=1) / N)

print(f"true prevalence                 {truth.mean():.4f}")
print(f"model-only estimate             {naive:.4f}   (no valid interval)")
print(f"hand-labels only  n={n}        {classical:.4f} +/- {1.96*se_cl:.4f}")
print(f"prediction-powered              {ppi:.4f} +/- {1.96*se_ppi:.4f}")
print(f"\ninterval width ratio ppi/classical {se_ppi/se_cl:.3f}")
print(f"hand-labels needed to match ppi   {int(n * (se_cl/se_ppi)**2)}")

true prevalence                 0.2228
model-only estimate             0.2417   (no valid interval)
hand-labels only  n=400        0.2375 +/- 0.0417
prediction-powered              0.2192 +/- 0.0270

interval width ratio ppi/classical 0.647
hand-labels needed to match ppi   955

The model-only estimate is biased upward by two points and comes with no interval that means anything, because nothing in the pipeline measured the model’s error. The corrected estimator uses the same four hundred hand labels as the classical one and reaches the precision that would otherwise have cost nine hundred and fifty-five --- a saving of more than half the annotation budget, with validity that does not depend on the model being good. This is the shape of the argument for local models in statistical work: not that the model is accurate, but that it is cheap enough to run on everything, and that you know how to correct it.

1.8What to Record in a Paper¶

For local-model work a methods section must state: the model repository and the exact revision or digest --- not the friendly tag --- the quantization format, the runtime and its version, the hardware, the decoding parameters (temperature, top-$p$, maximum tokens) and the seed, plus the number of seeds if you report variability. It must state the tokenizer if any token count is reported. Where the model acted as an instrument, it must state the size of the hand-labelled subsample, how it was drawn, the estimated sensitivity and specificity, and the estimator used to correct for them. Where an agent reached the model through MCP, name the servers and their versions, and say which had write access. Prompts, project context files and skills belong in the supplement as files. The one-line test is whether a reader could re-run your pipeline next year and get an answer within your reported uncertainty --- which is a weaker and more honest claim than bitwise reproduction, and the only one the arithmetic supports.

1.9Exercises¶

1. Using the memory example, find the largest context length that fits a 24 GiB accelerator at batch 8 for a 13B model with $L=40$ layers, 40 key-value heads and $d_h=128$, at 4-bit weights and 16-bit cache. Then repeat with grouped-query attention reducing the key-value heads to 8, and report the ratio.

2. Take the roofline example and add a fixed per-request overhead of 5 ms. Find the batch size at which throughput is within ten percent of its asymptote, and explain why a latency-sensitive interactive workload and a corpus-annotation workload should be served by different configurations.

3. Modify the tokenizer example to compute cost rather than token count: define a rate per thousand tokens as a variable, and report the cost of processing a 500-document corpus under both vocabularies. State precisely which quantity in your answer is a property of the text and which is a property of the instrument.

4. Computational. Write the schema validator from the MCP example as a decorator that wraps a Python function, and use it to reject malformed calls before execution. Then measure something: generate 1000 calls with a ten-percent malformation rate, and report the fraction rejected, the fraction that pass validation but return wrong results, and what that says about the guarantee a schema provides.

5. Computational. Reproduce the prediction-powered example with a local model of your choice on a real annotation task with at least 200 hand labels. Report the classical interval, the corrected interval, the effective sample size gained, and the sensitivity and specificity you measured. Then show what happens to the corrected estimator when the model’s errors are differential --- when sensitivity differs between two subgroups you intend to compare.