diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7885d767c..b1371ee4e 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4478,104 +4478,6 @@ 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. @@ -4638,21 +4540,15 @@ 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. 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) + // own small sigma. + std::vector> iv(start[n_sets]); // (I*corr, 1/(sigma*corr)^2) { 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 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}; + const std::pair v{mf.I[i] * mf.corr[i], sc > 0.0f ? 1.0f / (sc * sc) : 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; @@ -4679,6 +4575,7 @@ 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 257886bc3..0ecf08707 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -1126,14 +1126,6 @@ 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 349343b1a..98aa1b18a 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -103,9 +103,6 @@ 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,