mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

RotationScaleMerge: refit the error model about the merge's own mean

Browse Source 
The error model was fitted once, on deviations from a mean that weighs every full by its COUNTING
variance, and the outlier test's median took the same weights. Where equivalents disagree beyond
counting statistics, the low-count observations dominate that centre: on a strongly absorbing crystal
(YAG, Ia-3d, equivalents spread over a factor 100 after scaling) the mean of (0 4 0) sat at 21k among
observations from 11k to 1.9M, a ran into its bound (100), b came out at 800% internally, and the
six-sigma test about the biased median removed 63% of the observations - the strong ones.

Now, after the first fit, the model is refitted about the model-weighted mean (each full weighted by
the variance the fitted model gives it at the reflection's mean) until a and b settle, and the
rejection median takes the same weights. Where counting statistics are right nothing moves. Following
Blessing (1997) J. Appl. Cryst. 30, 421-426. GPU path: the refitted means are uploaded (SetEmMean).

Measured (SHELXL R1(>4sigma) against the COD model, sm-a's harness; rc174 scaling):
  YAG 0.556 -> 0.127 (XDS 0.083 merged), rejected 8055 -> 53, normalised deviations calibrated
  (median |z| 0.62-0.69 in every intensity decile); aspirin 20 keV 0.0964 -> 0.0958; citric acid
  0.161 -> 0.159; HEPES 0.0903 -> 0.0899; L-cystine 25 keV 0.1425 -> 0.1456;
  aspirin 25 keV 0.094 -> 0.106 (fixed-model R1 0.107 -> 0.173): its strong equivalents split into two
  frame-dependent populations from the per-frame partial scaling (sm-a's dq-smallmol), which the old
  under-sized sigmas happened to cut; with that scaling fixed (f69339ce6 + pooling, --no-scale-partials)
  this change is neutral to better on every small molecule (aspirin 20 .0618 -> .0586, aspirin 25
  .0456 -> .0454, citric .1093 -> .1006, HEPES .0703 -> .0697, YAG .649 -> .222; SHELXL GooF ~1.1).
  => ship together with the scaling fix.
Proteins (GPU full runs): CC1/2 and R_meas unchanged to 0.002; ISa myob 9.06 -> 8.20, thau 52.5 -> 47.5,
cytc 25.8 -> 25.4, lyso 29.4 -> 29.4 (still above XDS's 5.2 / 44.5 / 31.8 / 28.3 except cytc).
CPU build gives the same statistics as the GPU build on aspirin 20 keV and myob.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB
This commit is contained in:
leonarski_f
2026-10-04 00:48:35 +02:00
co-authored by Claude Opus 5.5
parent 359e749074
commit 1baf926066
3 changed files with 118 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+107 -4
View File
@@ -4448,6 +4448,104 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
resid * resid / factor, o.d});
}
}
fit_error_model(samples);
// ---- The error model refitted about the mean the merge takes. ----
// The fit above is made about a mean that weighs each full by its COUNTING variance. Where the
// equivalents of a reflection agree to within counting statistics that is the merge's own centre,
// and the model it gives is the model the merge applies. Where they do not - a systematic error
// the counting variance does not know about - it is the wrong centre, and in the worst way: the
// full with the fewest counts has the smallest variance, so the mean, and every deviation the
// model is fitted on, is pulled toward the observations that read lowest. Measured on a strongly
// absorbing crystal whose equivalents spread over a factor 100: the mean of (0 4 0) sat at 21k
// among observations from 11k to 1.9M, the deviations of the rest read as hundreds of counting
// sigmas, a ran into its bound and the six-sigma test about the median taken with the same
// weights deleted 63% of the observations - the strong ones (SHELXL R1 0.50 on that merge).
// So refit about the mean the MERGE takes - each full weighted by the variance the fitted model
// gives it, evaluated at the reflection's mean - until the model stops moving; the outlier
// test's median below takes the same weights. Where counting statistics were right the first fit
// is the answer and nothing moves; where b dominates, the weights even out, as the model says.
// Following Blessing (1997) J. Appl. Cryst. 30, 421-426: the centre an outlier is judged from
// is weighted as the merge weighs.
if (error_model_active) {
// Fulls in group order, so each pass over a group's members is one contiguous walk and the
// groups split across threads without every thread reading every full.
std::vector<int32_t> gstart(n_groups + 1, 0);
for (int i = 0; i < n_full; ++i)
if (mf.group[i] >= 0 && cnt[mf.group[i]] >= 2) ++gstart[mf.group[i] + 1];
for (int g = 0; g < n_groups; ++g) gstart[g + 1] += gstart[g];
std::vector<int32_t> member(gstart[n_groups]);
{
std::vector<int32_t> fill(gstart.begin(), gstart.end() - 1);
for (int i = 0; i < n_full; ++i)
if (mf.group[i] >= 0 && cnt[mf.group[i]] >= 2) member[fill[mf.group[i]]++] = i;
}
std::vector<double> next_mean;
std::vector<Sample> next;
constexpr int REFIT_ROUNDS = 10;
for (int round = 0; round < REFIT_ROUNDS; ++round) {
const double a = error_model_a, b = error_model_b;
const auto model_var = [&](int i, double mean) {
const double sc = static_cast<double>(mf.sigma[i]) * mf.corr[i];
return a * counting_variance(fulls[i], mean, sc * sc) + (b * mean) * (b * mean);
};
std::vector<std::vector<Sample>> part(ThreadsForWork(member.size(), nthreads));
const int nt = static_cast<int>(part.size());
next_mean.assign(n_groups, NAN);
// Each group's mean and samples come from its own members in index order, so the result
// does not depend on how the groups were split.
ParallelFor(nt, nt, [&](int t) {
const int g0 = static_cast<int>(static_cast<int64_t>(n_groups) * t / nt);
const int g1 = static_cast<int>(static_cast<int64_t>(n_groups) * (t + 1) / nt);
for (int g = g0; g < g1; ++g) {
if (gstart[g + 1] - gstart[g] < 2 || !std::isfinite(em_mean[g])) continue;
double sw = 0.0, swI = 0.0, swh[2] = {0.0, 0.0}, swIh[2] = {0.0, 0.0};
int nh[2] = {0, 0};
for (int q = gstart[g]; q < gstart[g + 1]; ++q) {
const int i = member[q];
const double v = model_var(i, em_mean[g]);
if (!(v > 0.0)) continue;
const double I_corr = static_cast<double>(mf.I[i]) * mf.corr[i];
sw += 1.0 / v; swI += I_corr / v;
if (merge_friedel && group_has_hands[g]) {
swh[obs_hand[i]] += 1.0 / v; swIh[obs_hand[i]] += I_corr / v; nh[obs_hand[i]]++;
}
}
if (!(sw > 0.0)) continue;
const double mean = swI / sw;
next_mean[g] = mean;
// The samples about the hand's own mean where it has two of its own, as above.
for (int q = gstart[g]; q < gstart[g + 1]; ++q) {
const int i = member[q];
const double v = model_var(i, em_mean[g]);
if (!(v > 0.0)) continue;
const int hh = obs_hand[i];
const bool on_hand = merge_friedel && group_has_hands[g] && nh[hh] >= 2 && swh[hh] > 0.0;
const double centre = on_hand ? swIh[hh] / swh[hh] : mean;
const double factor = 1.0 - (1.0 / v) / (on_hand ? swh[hh] : sw);
if (factor < 0.05) continue;
const double sc = static_cast<double>(mf.sigma[i]) * mf.corr[i];
const double resid = static_cast<double>(mf.I[i]) * mf.corr[i] - centre;
part[t].push_back({counting_variance(fulls[i], mean, sc * sc), centre * centre,
resid * resid / factor, mf.d[i]});
}
}
});
for (int g = 0; g < n_groups; ++g)
if (!std::isfinite(next_mean[g])) next_mean[g] = em_mean[g];
next.clear();
for (auto &v : part) next.insert(next.end(), v.begin(), v.end());
em_mean.swap(next_mean);
samples.swap(next);
fit_error_model(samples);
if (std::fabs(error_model_a - a) <= 1e-3 * a
&& std::fabs(error_model_b - b) <= 1e-3 * std::max(b, 1e-6))
break;
}
#ifdef JFJOCH_USE_CUDA
if (use_gpu_merge) gpu_->SetEmMean(em_mean.data());
#endif
}
// Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on
// the GPU; cheap here, cnt >= 3 filter from the em pass). Fed to the merge accumulate.
@@ -4510,15 +4608,21 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
// Inverse-variance WEIGHTED median: a frame the crystal barely diffracted on is scaled up by
// 1/G together with its sigma, so a plain median lets two such observations outvote one
// well-measured one - and the test below then rejects the well-measured one against its
// own small sigma.
std::vector<std::pair<float, float>> iv(start[n_sets]); // (I*corr, 1/(sigma*corr)^2)
// own small sigma. The variance is the model's, at the reflection's mean - the weight the
// merge gives the observation (see the centre refit above).
std::vector<std::pair<float, float>> iv(start[n_sets]); // (I*corr, 1/model variance)
{
std::vector<int32_t> fill(start.begin(), start.end() - 1);
for (int i = 0; i < n_full; ++i) {
const int g = mf.group[i];
if (g < 0) continue;
const float sc = mf.sigma[i] * mf.corr[i];
const std::pair<float, float> v{mf.I[i] * mf.corr[i], sc > 0.0f ? 1.0f / (sc * sc) : 0.0f};
const double mean = std::isfinite(em_mean[g]) ? em_mean[g] : static_cast<double>(mf.I[i]) * mf.corr[i];
const double mv = error_model_active
? error_model_a * counting_variance(fulls[i], mean, static_cast<double>(sc) * sc)
+ (error_model_b * mean) * (error_model_b * mean)
: static_cast<double>(sc) * sc;
const std::pair<float, float> v{mf.I[i] * mf.corr[i], mv > 0.0 ? static_cast<float>(1.0 / mv) : 0.0f};
if (cnt[g] >= 3) iv[fill[g]++] = v;
if (!pair_needed.empty() && pair_needed[pair_of_group[g]])
iv[fill[n_groups + pair_of_group[g]]++] = v;
@@ -4545,7 +4649,6 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
reject_median[g] = cnt[g] >= 3 ? set_median[g]
: !pair_needed.empty() ? set_median[n_groups + pair_of_group[g]] : NAN;
}
fit_error_model(samples);
}
// The full's sigma under the error model, with the variance evaluated at intensity I_for_b.
image_analysis/scale_merge/RotationScaleMergeGPU.cu
+8
View File
@@ -1110,6 +1110,14 @@ void RotationScaleMergeGPU::SetFrameCellOk(const uint8_t *frame_cell_ok) {
// The per-group inv-var mean (em_mean) + the per-full leverage-corrected error-model samples over the
// resident+scaled fulls. Stashes the filter context for the later MergeAccum/MergeRmeas calls.
void RotationScaleMergeGPU::SetEmMean(const double *em_mean) {
DeviceGuard guard(impl_->device, impl_->available);
auto &d = *impl_;
if (d.n_groups > 0)
CopyAndWait(d.m_em_mean.get(), em_mean, size_t(d.n_groups) * sizeof(double), cudaMemcpyHostToDevice,
impl_->s(), "ul em_mean");
}
void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partiality,
const uint8_t *hand, const uint8_t *has_hands,
double *em_mean_out, int32_t *cnt_out, double *s2_out,
image_analysis/scale_merge/RotationScaleMergeGPU.h
+3
View File
@@ -103,6 +103,9 @@ public:
// half-set weights multiplied by it. Requires MergeEmSamples first (em_mean resident).
// reject_var_add (n_groups) widens the pooled cut by the shell's own measured Bijvoet
// variance; null leaves the plain n-sigma test.
// Replace the per-group means MergeEmSamples left on the device (n_groups values): the merge's
// model sigmas are evaluated at them.
void SetEmMean(const double *em_mean);
void MergeAccum(double error_model_a, double error_model_b, bool error_model_active,
bool reject_outliers, double reject_nsigma, const float *reject_median,
const float *reject_var_add,