Surfaces: compare held-out half-set CC on Fisher's z; fit modulation, time, goniometer frame in that order

The correction-surface gate averaged the per-shell change of the held-out half-set CC over 20
equal-occupancy shells and adopted a surface on its sign. On the shells where a multiplicative error
shows, CC1/2 sits at 0.99-0.9996 and a large reduction of the error moves it in the fourth decimal;
the shells beyond the data's reach (up to 12 of 20 at CC ~ 0) each add +-0.01, so the sign was set
by noise (measured spread of the mean +-0.001-0.003 against effects of -0.0003..-0.0008). It refused
the 24x24 detector surface on four EIGER2 16M sweeps. The change is now averaged on atanh(CC) - the
candidate from the hq-bisect investigation.

The order of the three overlapping surfaces also decides what is adopted. With the gate fixed,
goniometer-frame first took a 5 keV insulin sweep from ISa 31.8 to 20.1 (the time surface then
refused; R_meas 7.8 -> 8.9%, CC to the model 0.8307 -> 0.8282) and cost a 0.1 deg-sliced thaumatin
sweep R_meas 18.2 -> 19.6%; time x detector first refused the detector surface on the 16 keV
thaumatin sweep (ISa 33 instead of 40). Modulation, time, goniometer frame avoids both.

Bare runs, --report-resolution at XDS's range, against rc173 (0252880e0): ISa / lowest-shell R_meas
thaumatin 16 keV 32.9/.030 -> 39.8/.028, lysozyme 90 deg sweeps 10.1/.059 -> 14.1/.048 and
10.1/.060 -> 13.8/.049, cytochrome c 19.2 -> 20.4 and 17.3 -> 19.3, thaumatin 3.8 keV 13.3 -> 14.8;
CC to a deposited model (standard-protein models on 14 in-house sets, the deposition on 13 open sets)
within +-0.0045, R_model_shell_scaled within +-0.002 except where d_min moved. 46 sets, one regression:
a half-image lysozyme test set, CC1/2 0.878 -> 0.862.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C

docs/ACKNOWLEDGEMENT.md

@@ -211,6 +211,13 @@ K. Diederichs, "Linking crystallographic model and data quality" (2012), Science
 P. A. Karplus, "Better models by discarding data?" (2013), Acta Cryst. D69, 1215-1222
 [doi:10.1107/S0907444913001121](https://doi.org/10.1107/S0907444913001121).
 
+**Fisher's z-transformation** — the cross-validation of the scaling correction surfaces averages the
+change of the half-set CC1/2 over resolution shells on atanh(CC), so that the shells near CC = 1, where
+a multiplicative error shows, are not outweighed by the noise of the shells without signal. R. A.
+Fisher, "Frequency distribution of the values of the correlation coefficient in samples from an
+indefinitely large population" (1915), Biometrika 10, 507-521
+[doi:10.2307/2331838](https://doi.org/10.2307/2331838).
+
 **The frame disposition** - which stretches of a rotation sweep are kept, carried at reduced weight or
 dropped from the merge - decides every exclusion on delta-CC1/2, the change in the overall CC1/2 when a
 group of images is left out, measured in the sigma-tau form so that no random half-dataset split is

docs/CPU_DATA_ANALYSIS.md

@@ -52,6 +52,7 @@ The methods draw on, and in places reimplement, solutions from:
 - A. Hennequin, B. Couturier, V. V. Gligorov & L. Lacassagne, "SparseCCL: Connected Components Labeling and Analysis for sparse images", DASIP 2019, 65-70 (the connected-component labelling of §3.4, used via ACTS/traccc).
 - S. French & K. Wilson, "On the treatment of negative intensity observations", *Acta Cryst.* **A34** (1978), 517-525 (Bayesian amplitude estimation from intensities).
 - A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free cross-validation).
+- R. A. Fisher, "Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population", *Biometrika* **10** (1915), 507-521 (the z-transformation on which the correction surfaces' held-out half-set CC1/2 is compared).
 - M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14).
 - J. P. Wright, "Experiences with GPU decompression for bitshuffle + LZ4 data", HDF5 User Group meeting (2021), and [github.com/jonwright/bslz4decoders](https://github.com/jonwright/bslz4decoders) (device-side decoding of bitshuffle+LZ4 images, §0).
 - A. Thorn & G. M. Sheldrick, "ANODE: anomalous and heavy-atom density calculation", *J. Appl. Cryst.* **44** (2011), 1285-1287 (anomalous difference density read at the model's sites).

docs/CPU_DATA_ANALYSIS_INTEGRATION.md

@@ -370,11 +370,13 @@ After scale-fulls, four **correction surfaces** are fitted on the combined fulls
 - **Decay.** Radiation damage weakens later frames more at higher resolution — a resolution×time (Debye–Waller) systematic the resolution-flat per-image scale cannot capture. A single global relative-$B$ rate is fitted, $\ln(I_\mathrm{ref}/I_\mathrm{obs}) = 2\,(\mathrm{d}B/\mathrm{d}n)\,(n-\bar n)\,s^2$ (frame $n$, $s^2 = 1/4d^2$), and folded into the scale. It engages only when the total relative-$B$ over the run exceeds a physical floor (2 Ų); below that the decay is negligible and "correcting" it only spreads symmetry equivalents (same $s^2$, different frames). An optional **per-batch relative-$B$** (`--relative-b[=deg]`, off unless requested; 10°-of-rotation batches by default) extends the single global rate to a smooth $B(n)$ curve — the same $s^2$-weighted decay fit solved independently over short frame batches, curvature-penalized so it cannot over-fit and cross-validated like the surfaces below — for crystals whose decay is non-linear in dose. Its cross-validation splits on **ASU-group parity**, not the frame parity the surfaces below use: a per-batch parameter owns whole frames and so cannot be scored on a held-out frame, whereas splitting the symmetry equivalents tests whether a batch's $B$ generalises to reflections it was not fitted on.
 - **Absorption.** A smooth multiplicative factor over the diffracted-beam direction expressed in the goniometer (crystal) frame: each full's predicted detector position gives the lab diffracted direction, de-rotated by the spindle so a fixed crystal-frame direction is sampled at many rotation angles and its grid cell is well-determined. Negligible at hard X-rays / thin crystals; it matters at low photon energy.
 - **Modulation** (detector-plane flat-field). A smooth multiplicative factor over where each reflection lands on the detector (predicted $x,y$): symmetry-equivalents land at different positions as the crystal rotates, over-determining the surface. It absorbs detector-response and geometric systematics that inflate $R_\mathrm{meas}$.
-- **Time-dependent absorption.** The same surface as *Absorption*, but over (rotation angle × detector position) instead of the crystal-frame direction alone. The two agree while the illuminated volume stays put — the incident path then depends only on the spindle angle, which the per-image scale already takes, and the exit path is fixed in the crystal frame. Once the diffracting volume drifts through the beam the exit path becomes a function of the spindle angle as well, and nothing time-independent describes it. Fitted last, on 12 rotation bins × a 10×10 equal-occupancy detector grid, so the two time-independent surfaces get first claim on what they can explain.
+- **Time-dependent absorption.** The same surface as *Absorption*, but over (rotation angle × detector position) instead of the crystal-frame direction alone. The two agree while the illuminated volume stays put — the incident path then depends only on the spindle angle, which the per-image scale already takes, and the exit path is fixed in the crystal frame. Once the diffracting volume drifts through the beam the exit path becomes a function of the spindle angle as well, and nothing time-independent describes it. Fitted on 12 rotation bins × a 10×10 equal-occupancy detector grid.
+
+The surfaces overlap, so the order decides what is adopted: modulation first (every frame measures the static detector pattern, so its fine grid is the best determined), then time-dependent absorption, then the goniometer-frame absorption, whose cells collect directions from the whole sweep and the whole resolution range and which, fitted first, takes up part of what the other two describe.
 
 Each cell's factor is fitted under a **prior pull to 1** (a Gaussian prior of width 0.1 on its logarithm): a cell moves off 1 in proportion to the information its observations carry, so a cell with little signal stays near 1 instead of being fitted to its noise. Negative observations enter the fit as measured. On the detector-plane surfaces (modulation, time-dependent absorption) the component that is a function of resolution alone is projected out within resolution shells, since the symmetry equivalents of a reflection share one resolution and such a factor cannot be determined from them.
 
-Each surface is **cross-validated** on the half-set agreement it is meant to improve: fitted on even-numbered frames and used to merge the odd ones, fitted on the odd frames and used to merge the even ones, and kept only if the correlation between the two half-set means, taken within resolution shells, rises above the same two halves merged with no surface. Each half is corrected by a surface it did not help to fit, so a surface fitted to noise lowers the correlation, and a correlation within a shell is blind to the resolution-dependent scale the surface cannot determine.
+Each surface is **cross-validated** on the half-set agreement it is meant to improve: fitted on even-numbered frames and used to merge the odd ones, fitted on the odd frames and used to merge the even ones, and kept only if the correlation between the two half-set means, taken within resolution shells, rises above the same two halves merged with no surface. Each half is corrected by a surface it did not help to fit, so a surface fitted to noise lowers the correlation, and a correlation within a shell is blind to the resolution-dependent scale the surface cannot determine. The change is averaged over the shells on Fisher's $z = \operatorname{atanh}(CC)$, not on $CC$: the strong shells, where a multiplicative error matters, sit at $CC_{1/2} \approx 0.999$, where even a large reduction of the error moves $CC$ in the fourth decimal, and averaged on $CC$ itself the shells of pure noise beyond the reach of the data decide the sign.
 
 **Radiation-damage report (rotation, report-only).** Independently of whether any decay correction is applied, Rugnux measures and reports the relative Debye–Waller $B$ across the sweep: the per-image scale's correlation to the merge and the per-image mosaicity versus frame (dose), together with a per-batch relative-$B$ curve whose first→last change is a single headline number (measured before any decay correction, against the least-damaged early wedge). It is written to the log and to the merged mmCIF as a data-quality-vs-dose diagnostic and **never** alters the merged intensities — distinct from the decay correction above, which does fold into the scale.
 

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -3463,16 +3463,25 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
     const std::vector<double> A_even = fit_surface(0), A_odd = fit_surface(1);
     const std::vector<double> odd0 = half_means(1, ident), even0 = half_means(0, ident);
     const std::vector<double> odd1 = half_means(1, A_even), even1 = half_means(0, A_odd);
+    // Each shell's change is taken on Fisher's z = atanh(CC), not on CC itself. CC1/2 is S / (S + E),
+    // and on the shells a multiplicative error matters in - the strong ones - it sits at 0.99-0.9996,
+    // where removing a third of E moves it by 0.0003: averaged as it stands, every shell counting alike,
+    // the shells of pure noise past the reach of the data (each +-0.01) set the sign. That refused the
+    // detector surface on EIGER2 16M sets where it improved 16 of 20 shells - mean change in CC -0.0003
+    // to -0.0008 - and adopting it raised ISa 33 -> 40 and 10 -> 14 and lowered the lowest-shell R_meas
+    // 0.030 -> 0.028 and 0.059 -> 0.048. On z a change of E by a given fraction is worth the same at any
+    // CC, and a shell of noise still only adds its own +-0.01.
+    // Following Fisher (1915) Biometrika 10, 507-521
     double gain = 0.0;
     int n_cc = 0;
     for (int sh = 0; sh < nsh_cc; ++sh) {
         const double c0 = shell_cc(odd0, even0, sh), c1 = shell_cc(odd1, even1, sh);
-        if (std::isfinite(c0) && std::isfinite(c1)) { gain += c1 - c0; ++n_cc; }
+        if (std::isfinite(c0) && std::isfinite(c1)) { gain += std::atanh(c1) - std::atanh(c0); ++n_cc; }
     }
     if (n_cc > 0) gain /= n_cc;
     if (!(gain > 0.0)) {
-        logger.Info("{} correction: not cross-validated (held-out half-set CC {:+.4f} per shell, skipped)",
-                    name, gain);
+        logger.Info("{} correction: not cross-validated (held-out half-set CC {:+.4f} per shell on Fisher's z, "
+                    "skipped)", name, gain);
         return;
     }
     const std::vector<double> A = fit_surface(-1);
@@ -3481,7 +3490,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
             if (cell[i] >= 0)
                 fulls[i].corr = static_cast<float>(fulls[i].corr * A[cell[i]]);
     });
-    logger.Info("{} correction: cross-validated, held-out half-set CC {:+.4f} per shell; {} in {} round(s); "
+    logger.Info("{} correction: cross-validated, held-out half-set CC {:+.4f} per shell on Fisher's z; {} in {} round(s); "
                 "resolution gauge removed a {:.2f}x ramp; {} of {} cells on the [0.25, 4] clamp",
                 name, gain, settled ? "settled" : "NOT settled",
                 rounds_used, gauge_ramp, n_clamped, ncell);
@@ -5672,12 +5681,22 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
         RefineDecay(n_groups);
     if (refine_surfaces && relative_b_deg > 0.0)
         RefineRelativeB(n_groups);   // per-batch relative-B on top of the single decay slope
-    if (refine_surfaces && absorption_iter > 0)
-        RefineAbsorption(absorption_iter, n_groups);
+    // The three surfaces overlap, so the one fitted first gets the first claim and the others are tested
+    // on what it left - the order decides what is adopted. Detector first: every frame of the sweep
+    // measures the static detector pattern, so its fine grid is the best-determined of the three. Then
+    // the time x detector surface, which takes what changes through the sweep. The goniometer-frame
+    // surface last: its cells collect directions from the whole sweep and the whole resolution range,
+    // so fitted first it soaks up part of both others and fits it where it is not. Measured with the gate
+    // below: goniometer frame first, a 5 keV sweep then refused the time surface and fell from ISa 31.8
+    // to 20.1 (R_meas 7.8 -> 8.9%, CC to the model 0.8307 -> 0.8282), and a finely sliced sweep with a
+    // large drift went from R_meas 18.2 to 19.6%; time surface first, the detector surface was refused
+    // after it on a 16 keV sweep that it takes from ISa 33 to 40. In this order neither happens.
     if (refine_surfaces && modulation_iter > 0)
         RefineModulation(modulation_iter, n_groups);
     if (refine_surfaces && absorption_iter > 0)
-        RefineAbsorptionTime(absorption_iter, n_groups);   // last: the static surfaces get first claim
+        RefineAbsorptionTime(absorption_iter, n_groups);
+    if (refine_surfaces && absorption_iter > 0)
+        RefineAbsorption(absorption_iter, n_groups);
 
     // --- 4c. The disposition. What keeping each batch of the sweep costs the merged intensities
     //         (delta-CC1/2), which frames cost too much, and what became of every frame of the run.