// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only #pragma once #include #include #include #include // The error model of the rotation merge: var = a * s2 + b^2 * ^2, fitted so that the scatter of // symmetry equivalents matches it - in each bin of counting I/sigma, the median of dev2 / var is the // median of a chi^2 with one degree of freedom. // // One sample per observation: s2 is its counting variance as the merge rebuilds it at the reflection's // mean intensity, I2 that mean squared, dev2 its leverage-corrected squared deviation from the mean. // // Two choices, each against a measured failure of the fit this replaces, which took the medians of s2, // I2 and dev2 separately over bins of I2: // * Bins are equal counts in COUNTING I/sigma (I2 / s2), not in I2. b is identified by how far the // bins reach into the regime where b* dominates the counting term, and that is what this ranks // on; ranking on I2 mixes weak high-resolution reflections into the strong bins. // * The median is taken of the NORMALISED deviation, each sample divided by the variance the previous // iteration gave it. A median of dev2 over observations whose variances differ is not 0.455 times // their mean variance, and the ratio of three medians read a too low and b too high: on most of 28 // sets the core of the normalised deviations then scattered at up to 1.8 times its stated variance // in the weak and middle bins and at 0.1-0.7 times it in the strongest. // A median, not a mean: on heavy-tailed data (split and powder-contaminated crystals) a mean calibrates // the sigmas on the tails, the merge's outlier test widens with them, and CC1/2 fell 0.80 -> 0.71 on one // such set where this fit raises it to 0.83. struct ErrorModelSample { double s2, I2, dev2; float d; int32_t group = 0; // the unique reflection it belongs to; only the bootstrap below reads it }; struct ErrorModelFit { bool active = false; // the pool was large enough to fit double a = 1.0, b2 = 0.0; bool b_measured = true; // false: no bin reached where b could be seen; b is held at 0 bool b_resolved = true; // b^2 stood clear of two of its standard errors, with a inside its bounds // and b measured bool a_at_bound = false; // a stopped on its clamp, so the b^2 beside it is not the solution of any fit }; // Scratch: the partitioned copy of the pool (reused between calls to save the allocation), with the rank // key and the current normalised deviation. struct ErrorModelBinned { double snr2, s2, I2, dev2, z2; }; // The fit's equal-count bins of counting I/sigma; the strongest bin's median (I/sigma)^2 must exceed // ERROR_MODEL_B_LEVER_MIN for b to be fitted (see FitErrorModel); and the median of a chi^2 with one // degree of freedom, which the median normalised deviation of a bin is compared with. constexpr int ERROR_MODEL_BINS = 16; constexpr double ERROR_MODEL_B_LEVER_MIN = 4.0; constexpr double ERROR_MODEL_CHI2_1_MEDIAN = 0.454936; ErrorModelFit FitErrorModel(const std::vector &pool, std::vector &scratch, size_t nthreads); // One update of the fit from its per-bin terms - bs2 and bI2 the sums of s2 / var and I2 / var under // the current (a, b^2), bd2 the bin's median normalised deviation expressed as the sum it stands for: // the new (a, b^2) into a, b2 and fit. False when the fit stops, settled or no longer determined. bool ErrorModelUpdate(const std::vector &bs2, const std::vector &bI2, const std::vector &bd2, double &a, double &b2, ErrorModelFit &fit); // The 95% interval of b^2 (b2 being the pool's own fit) by a bootstrap over the reflections: each of the // replicates refits the pool with every group drawn a Poisson number of times, its samples going or // staying together, since symmetry mates are deviations from one shared mean. A pool over max_samples is // resampled thinned to about that size, which bounds the cost. The draws are a hash of (replicate, // group), so the interval does not depend on the threads or on the order of the pool's groups. // Report-only: nothing fitted or merged reads it. Empty if no replicate gave a fit. struct ErrorModelB2Interval { double low, high; }; std::optional BootstrapErrorModelB2(const std::vector &pool, double b2, int replicates, size_t max_samples, size_t nthreads); // The error model of the stills merge: the rotation model plus a partiality term, // var = a * s2 + b^2 * ^2 + c^2 * ((1 - p) / p)^2 * ^2, // p the partiality the observation is divided by. A still records one sample of its rocking curve, and the // partiality model's error, carried through the division by p, grows like 1/p in the tail; without the term // the observations below p ~ 0.5 cannot be merged at all, with it they add to the merge. Measured on four // serial and rotation-as-serial sets, a term growing like 1/p fits the scatter better than the one a // Gaussian rocking curve propagates (|eps| sigma_eps / sigma^2, which grows only like sqrt(-ln p)) and gives // the stronger anomalous map; the curve's real tail is heavier than its Gaussian core. // Fitted as FitErrorModel fits the rotation model - equal-count bins of counting I/sigma, the median // normalised deviation of each, every bin relative, iterated - with each bin split in four by // ((1 - p) / p)^2, so that c is told apart from b by observations of one strength and different partiality. struct StillsErrorModelSample { double s2, I2, P2, dev2; // P2 = ((1 - p) / p)^2 * I2 }; struct StillsErrorModelFit { bool active = false; double a = 1.0, b2 = 0.0, c2 = 0.0; bool b_resolved = false; // as ErrorModelFit::b_resolved }; // start: where the iteration begins - identity by default, or the fit of a previous round on a pool that // has moved little since. StillsErrorModelFit FitStillsErrorModel(const std::vector &pool, size_t nthreads = 1, const StillsErrorModelFit &start = {});