diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 81dabcd22..edccd2d8c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -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 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 member(gstart[n_groups]); + { + std::vector 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 next_mean; + std::vector 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(mf.sigma[i]) * mf.corr[i]; + return a * counting_variance(fulls[i], mean, sc * sc) + (b * mean) * (b * mean); + }; + std::vector> part(ThreadsForWork(member.size(), nthreads)); + const int nt = static_cast(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(static_cast(n_groups) * t / nt); + const int g1 = static_cast(static_cast(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(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(mf.sigma[i]) * mf.corr[i]; + const double resid = static_cast(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> 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> iv(start[n_sets]); // (I*corr, 1/model variance) { std::vector 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 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(mf.I[i]) * mf.corr[i]; + const double mv = error_model_active + ? error_model_a * counting_variance(fulls[i], mean, static_cast(sc) * sc) + + (error_model_b * mean) * (error_model_b * mean) + : static_cast(sc) * sc; + const std::pair v{mf.I[i] * mf.corr[i], mv > 0.0 ? static_cast(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. diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 8a0e97f85..60a2e8d76 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -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, diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 766e2f640..17b56c520 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -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,