diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 5d4d49785..b424c990d 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4696,7 +4696,6 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool bool error_model_b_unmeasured = false; // b had no leverage; ISa is not reported bool error_model_b_resolved = true; // b^2 stood clear of its own standard error std::vector &samples = em_samples; - samples.clear(); // A full's counting variance at intensity I_for_b. Its own sigma carries its own Poisson fluctuation: // a full that came out low has a smaller sigma, so 1/sigma^2 weights it up and the merged mean drifts // below . Rebuild the variance at the reflection's EXPECTED intensity - var_bkg + var_per_I*, @@ -4780,6 +4779,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for (int g = 0; g < n_groups; ++g) if (cnt[g] >= 2 && sw[g] > 0.0) em_mean[g] = swI[g] / sw[g]; + samples.clear(); samples.reserve(fulls.size()); for (const auto &o : fulls) { if (!usable_merge(o) || cnt[o.group] < 2) continue; @@ -5079,8 +5079,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // and whether it is usable at all are fixed here, and neither depends on the error model that the // merge below refits. Rank parity within the reflection - see AssignHalvesOfMembers - so the halves are // balanced and identical on the host and device paths. - std::vector obs_key(fulls.size()); - std::vector merge_half(fulls.size(), 0); + std::vector &obs_key = merge_obs_key; + std::vector &merge_half = merge_half_set; + obs_key.resize(fulls.size()); + merge_half.resize(fulls.size()); bool did_gpu_halves = false; #ifdef JFJOCH_USE_CUDA // On the device where the merge is: integer work, the same halves. @@ -5092,8 +5094,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (!did_gpu_halves) { ParallelChunks(static_cast(fulls.size()), ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) { - for (int i = lo; i < hi; ++i) + for (int i = lo; i < hi; ++i) { obs_key[i] = ObservationKey(fulls[i].h, fulls[i].k, fulls[i].l, fulls[i].frame); + merge_half[i] = 0; + } }); // mf.group is already this array: it was filled with the same `usable_merge(o) ? o.group : -1` // above, on all threads, so a second serial walk over the fat records only rebuilds it. @@ -5139,9 +5143,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // model has been refitted on the reflections that survive the cut. ---- std::vector &acc = merge_acc; Result result; - std::vector merged_I; + std::vector &merged_I = merge_merged_I; size_t reject_count = 0; - std::vector rejected_obs; // per-full outlier-rejected flag (both paths) + std::vector &rejected_obs = merge_rejected; // per-full outlier-rejected flag (both paths) // ---- The Wilson test (WilsonOutliers.h): every observation against the intensity distribution of // its shell, beside the median test. It reaches what the median cannot - a reflection measured @@ -5180,15 +5184,22 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // vector doubles its way up to a few hundred megabytes on a large cell, copying everything it // already holds each time. result.merged.reserve(n_groups); - merged_I.assign(n_groups, NAN); + merged_I.resize(n_groups); + rejected_obs.resize(fulls.size()); + ParallelChunks(n_groups, ThreadsForWork(n_groups, nthreads), [&](int lo, int hi) { + std::fill(merged_I.begin() + lo, merged_I.begin() + hi, NAN); + }); + ParallelChunks(n_full, ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) { + std::fill(rejected_obs.begin() + lo, rejected_obs.begin() + hi, uint8_t{0}); + }); reject_count = 0; - rejected_obs.assign(fulls.size(), 0); if (wilson_test) { // Each observation with its OWN sigma, the error model evaluated at its own intensity: the // reflection's mean, which the merge weights use, is exactly what an artefact inflates. // Which fulls take part, in order, and then each one's record on its own - on all threads. wilson_full = KeptIndices(n_full, nthreads, [&](int i) { return mf.group[i] >= 0; }); - std::vector wobs(wilson_full.size()); + std::vector &wobs = merge_wobs; + wobs.resize(wilson_full.size()); ParallelChunks(static_cast(wobs.size()), ThreadsForWork(wobs.size(), nthreads), [&](int lo, int hi) { for (int j = lo; j < hi; ++j) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index a523ca256..70dbf3546 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -17,6 +17,7 @@ #include "../IntegrationOutcome.h" #include "ErrorModel.h" +#include "WilsonOutliers.h" #include "Merge.h" // MergedReflection, MergeStatistics #include "SphericalHarmonicSurface.h" #ifdef JFJOCH_USE_CUDA @@ -427,6 +428,15 @@ private: }; MergeAccumStaging merge_accum; + // MergeAndStats' other per-full and per-group arrays, members for the same reason: five merges a + // run, each over up to tens of millions of fulls and millions of groups, and as locals each one + // allocated, faulted in and zeroed them again on one thread. As members they are only resized, and + // reset on all threads where a pass needs them reset. + std::vector merge_merged_I; + std::vector merge_rejected, merge_half_set; + std::vector merge_obs_key; + std::vector merge_wobs; + // Per-group scatter for the strong-reflection ISa asymptote (see MergeAndStats). A member for the // same reason: 32 bytes a group, once per merge. struct GroupScatter { double sum = 0, sum_sq = 0, sum_var = 0; int n = 0; };