Jungfraujoch/image_analysis/scale_merge/ErrorModel.h

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#pragma once

#include <cstddef>
#include <cstdint>
#include <optional>
#include <vector>

// The error model of the rotation merge: var = a * s2 + b^2 * <I>^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*<I> 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<ErrorModelSample> &pool, std::vector<ErrorModelBinned> &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<double> &bs2, const std::vector<double> &bI2, const std::vector<double> &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<ErrorModelB2Interval> BootstrapErrorModelB2(const std::vector<ErrorModelSample> &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 * <I>^2 + c^2 * ((1 - p) / p)^2 * <I>^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<StillsErrorModelSample> &pool, size_t nthreads = 1,
                                        const StillsErrorModelFit &start = {});