The global dataset

Two cross-checks fell out of this that we did not plan for, and both matter. R3 = 4,615 points against NNPDF4.0's published 4,618 — three points apart, so our independent assembly of the 82 datasets reproduces theirs essentially exactly. And R0 = 3,089 points against the ≈3,092 we have quoted for full-DIS all along — confirming that the control rung really is Phase A's dataset, so comparisons across rungs are like-for-like.

What "hadronic" means for cost. A DIS prediction contracts the FK table with the PDF once. A hadronic one contracts it with the PDF twice — two partons come in, one from each proton. The FK tables are correspondingly larger (13–122 MB versus 0.1–1.4 MB) and the convolution is ~100× more expensive per data point. Caveat: that ratio comes from a numpy emulation of the contraction, not Colibri's own convolution — the FK sizes are exact, the ratio is indicative.

What this rules out as a constraint. 868 MB of FK tables against 251 GB of RAM, and 231 GB of free disk. Storage and memory were the things I could not size in advance, and they turn out to be non-issues. Compute is the binding constraint — which is exactly why PPDF-32 (Fisher, no sampler) runs before any sampling.

On the parity check. The original gate (PPDF-10, 96 checks to ≲10⁻⁹) compares against the FPPDF reference package, which is not installed on this box and was not rebuilt. What is decisive here is cross-machine determinism: the model reproduces its Phase A verified values to 2.5×10⁻¹³ under a different OS and different library versions. If the port had drifted, that number would not hold.

# box: gadget4-twin-c (Nebius) — 64 vCPU, 251 GB RAM, no GPU (dispatch-bound work)
ssh -i ~/.ssh/nebius_fdm ubuntu@<ip>
# env built in ~4 min: jax 0.11.0 (CpuDevice), validphys, colibri
micromamba create -y -n colibri-dev -f ~/colibri/environment.yml

# stage the theory — 762 MB, 25 s
micromamba run -n colibri-dev vp-get theoryID 40_000_000
# -> fastkernel/ 759 MB (254 pineappl + 652 .dat), cfactor/ 2.6 MB, yamldb/ 528 KB

# the three probes
python ~/fk_probe.py             # 7/7 load, 5 hadronic=True
python ~/fk_total.py             # all four rungs, FK size + ndata
python ~/cross_machine_parity.py # port vs PPDF-10 values -> 2.5e-13

Rung definitions: analysis/phaseB/rungs.py · generated runcards: analysis/phaseB/runcards/ (24 = 4 rungs × 3 methods × closure/real, from make_runcards.py)

5 Aug. Two of my own estimates were wrong in the same direction — I had guessed the environment build would take about an hour (it took 4 minutes) and that hadronic FK tables would be "much larger", implying tens of GB (the whole theory is 762 MB). Both were guesses where a measurement was cheap, which is the argument for making this a gate rather than an assumption.

The genuinely useful surprise was the point-count agreement: 4,615 vs 4,618 was not something I set out to test, and it is stronger evidence that the dataset assembly is right than anything I had planned.

Why this was a gate and not a nicety. Phase A's own limitations slide says the real-data verdict is "transferred from closure with no external anchor" — coverage can only be measured where a truth is planted, so nothing in Phase A independently confirmed the fit was good on real data. This is that anchor. It also pre-empts the obvious objection to Phase B: that a rigid 52-parameter polynomial simply cannot carry global data, which would make any method comparison at global scale meaningless. It can.

On being 5.7% worse than NNPDF. That is the expected direction and a reasonable margin. They fit a 763-parameter neural network; we fit 52 polynomial coefficients in the MSHT20 functional form. More flexibility buys a better χ². The point is not to win — it is to be in the same regime, which we are.

The unplanned second result: the penalty block is load-bearing. Both penalty-free arms fail — R0 returns a non-finite χ², R3 lands at 16.35. Same model, same data, same optimiser and schedule; the only change is removing positivity and integrability. Those constraints are not cosmetic bookkeeping here, they are what keeps a 52-parameter polynomial inside the region where the likelihood is finite. No learning-rate schedule rescued it — fixed and exponentially-decayed both went to nan.

Correction, 2026-08-06. This result originally read "positivity is load-bearing." That over-attributed it: our penalty-free arm removes positivity and integrability in one switch, so it cannot separate them. A follow-up run that removes positivity only (integrability left on) gives epoch-50 loss 3.27×10⁶ vs 3.75×10⁶ with both on — so positivity does contribute on its own, but the R3 figure of 16.35 belongs to the two constraints together, not to positivity alone. The χ² numbers above are unaffected.

# four fits in parallel on 64 vCPU — penalties on/off × R0/R3
for job in R0_gd_nopos R0_gd_pos R3_gd_nopos R3_gd_pos; do
  micromamba run -n colibri-dev msht20_exe phaseB_${job}_real.yaml -o out_${job} &
done; wait

# the learning-rate schedule that made 60k epochs stable.
# A FIXED lr=2e-3 descends to loss 3,876 by epoch 26k and then diverges to nan before 36k.
scheduler:
  name: exponential_decay
  params: {init_value: 2e-3, transition_steps: 4000, decay_rate: 0.75}
max_epochs: 60000

# NB: run_gradient_descent cannot be a standalone action — it demands `initial_parameters`
# directly. Use run_hessian_fit and harvest the chi2 its GD phase logs BEFORE the
# eigendecomposition, which stalls on a slow XLA compile (the PPDF-20 wall).
grep "Minimum chi2 found" R3_gd_pos.log   # -> 5707.6675364220355

Runcards: analysis/phaseB/runcards/phaseB_R*_gd_{pos,nopos}_real.yaml, generated by make_runcards.py · logs on box ~/phaseB/R*_gd_*.log

5 Aug. Getting here cost six sequential schema failures — integrability wanted a mapping not a list, the initialiser wanted bounds, uniform initialisation produced a non-finite χ² from every restart, hessian_settings turned out to be mandatory, run_gradient_descent was the wrong action, and then the optimiser diverged at 60k epochs. Each was read out of the Colibri source rather than guessed, and all six are written up in analysis/phaseB/README.md so they cost an hour once instead of an hour each time.

The genuinely useful surprise was the penalty finding. I had it backwards for most of the afternoon — I read the 260:1 penalty-to-χ² ratio as evidence the penalties were breaking the fit, when in fact they were the only thing holding it together. The control that was supposed to validate the port ended up re-characterising Phase A's own convergence.

The question this answers. Phase A established that the three uncertainty methods diverge because a large number of the 52 directions are unconstrained by DIS data (12, exactly re-measured — we long quoted 14). The obvious follow-up — does more data fix it? — was left open, because the Tier-2 sampling run came back under-converged at ESS 30 and no claim was made (PPDF-26). This card answers it by measurement.

The answer: more data helps, then stops. Drell–Yan/W/Z removes three directions. Jets and direct photon remove one more. Top-quark data — 23 datasets, the largest single block added — removes none. The full NNPDF4.0-like global set still leaves 8 of 52 directions unconstrained.

Against the pre-registered brackets (≤8 cured, ≥12 structural) R3 = 8 lands exactly on the CURED threshold. We are calling that a boundary result, not a pass: landing precisely on a threshold is the least informative place to land, eight directions are still unconstrained by every dataset we have, and a one-decade change of numerical step size would have moved it. The bracket was never the finding. The saturation is — and that is what survived exact re-measurement.

Which data, not how much. Per dataset added: DY/W/Z removes 0.086 flat directions, jets and photon 0.200 — about 2.3× more effective, not the 7× we first reported — and top 0.000. That is physically coherent: jets and direct photon are the gluon-sensitive observables, and much of the degeneracy lives in the gluon. Once the gluon directions are pinned, adding more gluon-sensitive data does nothing, and the remaining flat directions are elsewhere — the strange sector Phase A localised.

An unplanned validation. In the original forward-difference run, R0 = 14 and R1 = 12 reproduced PPDF-19/PPDF-20 exactly — on a different box, jax 0.10.2→0.11.0, and with the dataset assembled independently from NNPDF's own global runcard rather than from Phase A's card. Condition numbers (1.2–1.6×10¹¹) match Phase A's ≈10¹¹ too. Phase A's counts were forward-differenced by the same code, so this agreement validates the machinery's determinism, not the counts themselves — the exact central-difference values are 12 and 9. It was not designed as a reproduction test; it is still stronger evidence for the Fisher machinery than anything that was.

# one rung at a time — never parallel; concurrent Hessians took the box down once
for rung in R0 R1 R2 R3; do
  ( ulimit -v $((150*1024*1024))          # cap the job, not the box
    micromamba run -n colibri-dev python fisher_fd.py runcards/fisher_${rung}.yaml $rung )
done

# prior-normalised so "flat" is scale-independent: lambda' < 2 means the direction
# is less constrained than the prior itself
H[:,i] = (grad(truth + eps*e_i) - grad(truth)) / eps      # FD of the autodiff gradient
Hs     = H * outer(prior_halfwidth, prior_halfwidth)      # dimensionless
nflat  = (eigvalsh(Hs) < 2).sum()

# jax.hessian is too slow at this data size — the PPDF-20 lesson, reused
XLA_PYTHON_CLIENT_PREALLOCATE=false   # and do not let JAX grab all the RAM

Script: analysis/msht20/fisher_fd.py (Phase A, reused unchanged) · runcards analysis/phaseB/runcards/fisher_R*.yaml from make_fisher_cards.py · eigenvalues ~/fisher_R*_eig.npy

5 Aug, overnight. This was the card the whole phase was designed around, and it was deliberately placed before any sampling: Fisher needs no sampler, so it cannot fail the way PPDF-26 failed. That ordering paid off — the sampling card the same night (PPDF-33) died on a config error, while this one produced a clean four-point curve.

The result I did not expect was the per-dataset efficiency. Five jet and photon datasets did the work of thirty-five Drell–Yan datasets, and twenty-three top datasets did nothing at all. That has a direct consequence for what to ask for next: the case for targeted data (ν-DIS dimuon, W+charm) over bulk data now rests on a measurement rather than an intuition.

The dominant risk in all of Phase B. Our only calibration point at larger data is PPDF-26: 73 datasets, ESS 30, 178/3000 divergences — a failure. 82 datasets with quadratic hadronic observables is strictly harder. If the Bayesian arm does not converge there is no reference method and the comparison is meaningless.

Which is why whitening is not optional here. Phase A established that a pilot-covariance mass matrix, not more samples, is what makes these posteriors converge (PPDF-28). Phase B starts there rather than discovering it again.

Reporting a failure is a valid outcome. Phase A already reports one — nested sampling escalated four times and never settled. That non-convergence is itself a transferable result.

# args: <runcard> <out> <warmup> <samples> <target_accept> <max_tree_depth> <level_1_seed>
# NOTE argv[7] is level_1_seed = the L1 NOISE DRAW. Varying it changes the DATASET, not the
# chain. Independent chains on one posterior need a fixed argv[7] and a varying
# NUTSF_CHAIN_SEED (which seeds the warmup/sampling PRNGKeys; 0 reproduces the old keys).
for CS in 1 2 3; do
  NUTSF_CENTRE=map NUTSF_MAP_EPOCHS=20000 NUTSF_CHAIN_SEED=$CS \
  python analysis/phaseB/nuts_fisher.py \
    phaseB_R0_bayes_closureL1.yaml out_v11_chain$CS  2000 6000  0.9 8  123456
done
# the gate: pooled per-parameter ESS, and true multi-chain R-hat
python analysis/phaseB/combine_chains.py "out_v11_chain*" closure_truth_phaseB_R0.json

# the retraction, from artifacts already on disk — no new sampling
python analysis/phaseB/centre_scatter_projection.py \
  analysis/phaseB/results_backup/chains  phaseB_R0_bayes_closureL1.yaml

Harness ~/phaseB/harness11.sh · the Hessian and MC arms differ from the Bayesian arm only in the inference block (site/code/msht_hessian_real.yaml vs msht_mc_real.yaml differ by one line)

Gate written down before the run. Coverage is not reported below the ESS bar — that rule is what makes PPDF-26 a non-claim rather than a bad claim.

Why five draws, again. Phase A measured the swing: coverage moves 47%↔76% between noise realisations, essentially uncorrelated with sampling quality. One draw is noise. This is a transferable warning — closure tests in the literature are commonly reported on a single draw.

Closure only. Coverage requires a known truth, so this is synthetic data by construction. The real-data statement comes separately, in PPDF-35, and the two must never be quoted as though they were the same modality.

# five draws — level_1_seed is the seed that actually varies the noise
for s in 1 2 3 4 5; do
  python3 analysis/phaseB/nuts_global.py --rung R3 --white --level1-seed $s
done
python3 analysis/phaseB/coverage.py --draws 5

Pattern: analysis/msht20/phase1/p1_analyze.py

The metric and the gate are fixed in context/PHASE_B_PREREG.md before the run, so neither can be tuned after seeing the answer.

What we can and cannot say today. PPDF-38 shows data still has purchase — removing any of four classes makes the degeneracy worse — so it is not purely structural. But whether new strangeness data drives 8 → 0 or only 8 → 7 is exactly what R4 would settle, and we cannot run it.

A three-step error worth recording, because the lesson generalises. Establishing this took three attempts. (1) We grepped fastkernel/ for WCHARM, found nothing, and concluded unavailable — correct. (2) We then found W+charm names in yamldb/ and "corrected" that to available — wrong: yamldb declares what a dataset would need, not what is present. (3) We loaded the commondata (37 points) and said all five load — wrong again: commondata and FK tables are separate resources. The load test that matters is the FK test, which is precisely why PPDF-30's port gate tested FK convolution rather than data availability.

Runcard: analysis/phaseB/runcards/fisher_R4.yaml — a copy of fisher_R3.yaml with the five W+charm datasets spliced into dataset_inputs and everything else byte-identical, so the added strangeness data is the only variable. Measurement would be analysis/phaseB/fisher_eps_scan.py … central, the same method as PPDF-32 and PPDF-38.

Datasets required: ATLAS_WCHARM_7TEV_WP-YL, ATLAS_WCHARM_7TEV_WM-YL, CMS_WCHARM_13TEV_WPWM-TOT-UNNORM, CMS_WCHARM_7TEV_WPWM-RATIO, CMS_WCHARM_7TEV_WPWM-TOT. Failure record: ~/phaseB/eps_central_R4.log. Full audit trail in context/PHASE_B_LEDGER.md.

Where this experiment came from. PPDF-32 showed the degeneracy saturates: the last 23 datasets removed nothing. The obvious follow-up — would a different kind of data help? — is blocked, because W+charm's FK tables are absent from theory 40000000 (7 of its 131 dataset mappings lack tables, and W+charm accounts for five of them). Rather than wait on grid generation, this asks the inverse question, which needs no new data at all: what does each class we already have uniquely contribute?

The null control is the point of the design. One of the five variants had a known answer before it ran. If removing top-quark data had returned anything other than zero, the scan would have contradicted PPDF-32 and one of the two would have been wrong. It returned exactly zero.

Cards: analysis/phaseB/runcards/fisher_R3_no{NuDIS,LepDIS,DY,Jets,Top}.yaml, each a copy of fisher_R3.yaml with one class spliced out of dataset_inputs and everything else — theory, cuts, planted truth, prior, penalties — byte-identical. Measurement: analysis/phaseB/fisher_eps_scan.py … central. Figure: analysis/phaseB/make_figure_38.py.

Artifacts: ~/phaseB/loo.log, ~/phaseB/loo_no*.log. Full audit trail in context/PHASE_B_LEDGER.md.

A methodological trap worth naming. The first attempt at this card reported a converged minimum that had not moved at all: every backtracked Newton step was rejected. That is impossible for a genuine descent direction, which is what exposed a scaling bug — the code used diag(1/width) where diag(width) belongs, making the step ~14× too large. It did not crash; it printed “converged in 1s” and a plausible spectrum.

And |∇χ²| is the wrong convergence measure here. It is expressed in raw parameter units across 52 parameters spanning orders of magnitude, so its size means nothing. The right criterion is the Newton decrement — the χ² a full Newton step would recover. Here that is 0.0104 out of 3,539.8, one part in 340,000. By that measure the point was at the minimum all along, and an apparently alarming |∇χ²| = 3.6×10⁴ was a scaling artefact.

Script: analysis/phaseB/hessian_at_minimum.py — Newton with backtracking to 10⁻⁴, prior-normalised as H * outer(width, width) exactly as PPDF-32, then a central-difference Hessian at the located minimum, scanned over three step sizes.

Artifacts: ~/phaseB/ham_fixed.log, min_R0_real_fixed.json, ham_eig_1e-0{3,4,5}.npy. Full audit trail including the scaling bug in context/PHASE_B_LEDGER.md.

Why this card matters. Phase A had no external anchor on real data — every result was internal to our own machinery. A published global fit is the one number a reader outside the project can check us against. Landing at 1.058× of NNPDF4.0 with a fourteenth of the parameters is the evidence that the MSHT20 polynomial parametrisation is not the limiting factor in anything else Phase B reports.

Card: analysis/phaseB/runcards/phaseB_R3_gd_pos_real.yaml (R3 = all 82 datasets, theory 40000000, positivity + integrability on). Executed by msht20_exe under harness v7 step 5; the best-fit χ² is logged before the Hessian eigendecomposition, which is the step that exhausts memory.

Artifacts: ~/phaseB/r3_gd_v7.log, ~/phaseB/out_R3_gd_v7/. Full audit trail in context/PHASE_B_LEDGER.md.