diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 2af4a9d6..cedc68a8 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -637,11 +637,14 @@ PostRefineResult PostRefineRotationGeometry(const std::vectorh, pp->k, pp->l, sys)), new ceres::CauchyLoss(0.02), beam, dist, const_cast(det_rot), rot_vec, p0, p1, p2); - p.SetParameterBlockConstant(const_cast(det_rot)); - p.SetParameterBlockConstant(rot_vec); - p.SetParameterBlockConstant(p0); p.SetParameterBlockConstant(p1); p.SetParameterBlockConstant(p2); } if (p.NumResidualBlocks() == 0) { beam_out[0] = beam_x0; beam_out[1] = beam_y0; dist_out = dist0; return false; } + // Everything but the beam and the distance is held at its step-A value. Once per problem: + // the block only has to exist, and repeating it per observation was up to MAX_OBS times + // five calls for the same five blocks. + p.SetParameterBlockConstant(const_cast(det_rot)); + p.SetParameterBlockConstant(rot_vec); + p.SetParameterBlockConstant(p0); p.SetParameterBlockConstant(p1); p.SetParameterBlockConstant(p2); p.SetParameterLowerBound(dist, 0, dist0 * 0.95); p.SetParameterUpperBound(dist, 0, dist0 * 1.05); for (int j = 0; j < 2; ++j) { p.SetParameterLowerBound(beam, j, beam[j] - 15.0); p.SetParameterUpperBound(beam, j, beam[j] + 15.0); } diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 42454be5..100822ab 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -231,23 +231,18 @@ void RotationScaleMerge::Ingest() { size_t total = 0; for (const auto &o : partials_out) total += o.reflections.size(); - partials.clear(); - partials.resize(total); - frame_start.assign(n_frames, 0); - frame_count.assign(n_frames, 0); frame_cell_ok.assign(n_frames, 1); g_partial.assign(n_frames, 1.0); - // The per-frame CSR up front. Every reflection of a frame becomes one observation, so where a - // frame's slice of `partials` starts is known before anything is copied - which is what lets the - // conversion below hand a frame to each thread instead of appending, frame after frame, on one. - // (Hence resize rather than reserve above: the slices are written, not pushed.) + // Where each frame's slice of the key array starts. Every reflection of a frame becomes one key, so + // this is known before anything is read - which is what lets the sweep below hand a frame to each + // thread instead of appending, frame after frame, on one. + std::vector src_start(n_frames, 0); { int32_t at = 0; for (int o = 0; o < n_frames; ++o) { - frame_start[o] = at; - frame_count[o] = static_cast(partials_out[o].reflections.size()); - at += frame_count[o]; + src_start[o] = at; + at += static_cast(partials_out[o].reflections.size()); } } @@ -259,6 +254,15 @@ void RotationScaleMerge::Ingest() { // which is what keeps the mean bit-exact. std::vector mean_bkg(n_frames, NAN); + // The h range comes out of the same sweep - the sort below buckets by h and needs to know how many + // buckets that is - and so does the narrow key of every observation. The fat Obs is NOT built here: + // the resolution range test below throws most of it away on a crystal that integrates far past what + // it merges, and eighty bytes written for an observation nothing can ever merge is eighty bytes of + // pure waste. Only the survivors are built, from the source reflections, once the runs are known. + constexpr int INT_LO = std::numeric_limits::min(); + constexpr int INT_HI = std::numeric_limits::max(); + std::vector unsorted(total); + std::vector h_lo(n_frames, INT_HI), h_hi(n_frames, INT_LO); ParallelFor(n_frames, nthreads, [&](int o) { if (reference_cell) { const auto cell = partials_out[o].latt.GetUnitCell(); @@ -266,68 +270,33 @@ void RotationScaleMerge::Ingest() { } double bkg_sum = 0.0; int bkg_n = 0; - int32_t at = frame_start[o]; + int lmin = INT_HI, lmax = INT_LO; + int32_t at = src_start[o]; for (const auto &r : partials_out[o].reflections) { if (std::isfinite(r.bkg) && r.bkg > 0.0f) { bkg_sum += r.bkg; ++bkg_n; } - Obs &obs = partials[at++]; - obs.h = r.h; obs.k = r.k; obs.l = r.l; - obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; obs.rlp = r.rlp; - obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg; - obs.px = r.predicted_x; obs.py = r.predicted_y; - obs.image_number = r.image_number; - obs.frame = o; - obs.on_ice = r.on_ice_ring ? 1 : 0; - obs.corr = r.image_scale_corr; - obs.group = -1; + unsorted[at] = SortKey{r.h, r.k, r.l, r.image_number, at, r.d}; + lmin = std::min(lmin, r.h); + lmax = std::max(lmax, r.h); + ++at; } if (bkg_n > 0) mean_bkg[o] = bkg_sum / bkg_n; + h_lo[o] = lmin; h_hi[o] = lmax; }); MeasureIncidentFlux(mean_bkg); - // Per-obs AcceptReflection finiteness (immutable) - lets ComputeAsuGroups stamp the ASU-group id per - // obs from a flat 1-byte array instead of re-reading the fat Obs struct for every space group. - finite_ok.resize(partials.size()); - // The h range comes out of the same sweep - the sort below buckets by h and needs to know how - // many buckets that is, and this pass already reads every observation. int h_min = 0, h_max = 0; { - const int n_obs = static_cast(partials.size()); - const int nt = static_cast(std::clamp(nthreads, 1, std::max(1, n_obs))); - constexpr int INT_LO = std::numeric_limits::min(); - constexpr int INT_HI = std::numeric_limits::max(); - std::vector lo_of(nt, INT_HI), hi_of(nt, INT_LO); - const int chunk = (n_obs + nt - 1) / nt; - ParallelChunks(nt, nthreads, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - int lmin = INT_HI, lmax = INT_LO; - const int lo = t * chunk, hi = std::min(n_obs, lo + chunk); - for (int i = lo; i < hi; ++i) { - auto &o = partials[i]; - // Divide the incident flux out of rlp here: this pass already has the observation - // in hand, and rlp is the numerator of corr, so the fits below see only what is - // left. Before the finiteness test, which is where it has always been. - o.rlp = static_cast(o.rlp / frame_flux[o.frame]); - finite_ok[i] = (std::isfinite(o.I) && std::isfinite(o.rlp) && o.rlp != 0.0f - && std::isfinite(o.sigma) && o.sigma > 0.0f) ? 1 : 0; - lmin = std::min(lmin, o.h); - lmax = std::max(lmax, o.h); - } - lo_of[t] = lmin; hi_of[t] = lmax; - } - }); int gmin = INT_HI, gmax = INT_LO; - for (int t = 0; t < nt; ++t) { - gmin = std::min(gmin, lo_of[t]); - gmax = std::max(gmax, hi_of[t]); + for (int o = 0; o < n_frames; ++o) { + gmin = std::min(gmin, h_lo[o]); + gmax = std::max(gmax, h_hi[o]); } if (gmin <= gmax) { h_min = gmin; h_max = gmax; } // otherwise there are no observations } // Sort ONCE by (raw h,k,l, image_number) and split into raw-hkl runs. This is the one expensive sort; // both the 3D combine (event split) and the per-space-group ASU grouping reuse this order. - // Sorting an index array whose comparator dereferences the 80-byte Obs is a cache miss per - // comparison over a multi-GB array, so sort a packed copy of the key fields instead. // // The observation's own index is the last key, which makes the order TOTAL. Two observations can // genuinely share (h,k,l) and image_number: the predictor emits both intersections of a @@ -337,7 +306,6 @@ void RotationScaleMerge::Ingest() { // member of an event, and its per-event sums are floating point. Ordering them by arrival makes // the result independent of the sort algorithm, which is what lets the order be reproduced by a // faster one. - struct SortKey { int32_t h, k, l; float image_number; int32_t idx; }; const auto key_less = [](const SortKey &a, const SortKey &b) { if (a.h != b.h) return a.h < b.h; if (a.k != b.k) return a.k < b.k; @@ -345,16 +313,16 @@ void RotationScaleMerge::Ingest() { if (a.image_number != b.image_number) return a.image_number < b.image_number; return a.idx < b.idx; }; - perm.resize(partials.size()); + // Two key arrays for a moment: the sweep above wrote them in index order, and where a key belongs + // in the bucketed order is not known until the whole histogram is. Both together are less than half + // of what the observation array they replace used to occupy on its own. + std::vector keys(total); { // Bucket by h first, then sort each bucket. h is the comparator's leading key, so the sorted // array is exactly the buckets laid end to end - and because the order is total (see above) // the permutation is the one a single std::sort produces, whichever algorithm gets there. // That is what makes this checkable: the output is byte-identical, not merely equivalent. - // - // The keys are built straight into their bucket slot, so this replaces the build pass rather - // than adding to it and the 20-bytes-per-observation array is never duplicated. - const int n = static_cast(partials.size()); + const int n = static_cast(total); const int nt = static_cast(std::clamp(nthreads, 1, std::max(1, n))); const int chunk = (n + nt - 1) / nt; const int H = h_max - h_min + 1; @@ -363,7 +331,7 @@ void RotationScaleMerge::Ingest() { ParallelChunks(nt, nthreads, [&](int tlo, int thi) { for (int t = tlo; t < thi; ++t) { const int lo = t * chunk, hi = std::min(n, lo + chunk); - for (int i = lo; i < hi; ++i) hist[t][partials[i].h - h_min]++; + for (int i = lo; i < hi; ++i) hist[t][unsorted[i].h - h_min]++; } }); @@ -381,18 +349,15 @@ void RotationScaleMerge::Ingest() { } bstart[H] = acc; - std::vector keys(partials.size()); ParallelChunks(nt, nthreads, [&](int tlo, int thi) { for (int t = tlo; t < thi; ++t) { std::vector fill = hist[t]; const int lo = t * chunk, hi = std::min(n, lo + chunk); - for (int i = lo; i < hi; ++i) { - const auto &o = partials[i]; - keys[fill[o.h - h_min]++] = - SortKey{o.h, o.k, o.l, o.image_number, static_cast(i)}; - } + for (int i = lo; i < hi; ++i) keys[fill[unsorted[i].h - h_min]++] = unsorted[i]; } }); + unsorted.clear(); + unsorted.shrink_to_fit(); // Largest bucket first: the tail of this is one bucket, so it should be the big one. std::vector order(H); @@ -405,10 +370,6 @@ void RotationScaleMerge::Ingest() { std::sort(keys.begin() + bstart[b], keys.begin() + bstart[b + 1], key_less); }); - ParallelChunks(n, nthreads, [&](int lo, int hi) { - for (int i = lo; i < hi; ++i) perm[i] = keys[i].idx; - }); - // The runs never cross an h boundary, so each bucket owns its own and they can be counted, // scanned and written without touching each other. Sizing the arrays from the count also // removes the repeated growth the push_backs paid for. @@ -440,7 +401,7 @@ void RotationScaleMerge::Ingest() { while (j < bstart[b + 1]) { if (keys[j].h != k0.h || keys[j].k != k0.k || keys[j].l != k0.l) break; if (!std::isfinite(d)) { - const float dj = partials[keys[j].idx].d; // resolution: only until one is usable + const float dj = keys[j].d; // resolution: only until one is usable if (std::isfinite(dj) && dj > 0.0f) d = dj; } ++j; @@ -454,7 +415,7 @@ void RotationScaleMerge::Ingest() { } }); } - DropOutOfRangeObservations(); + BuildInRangeObservations(keys); rawrun_group.assign(rawrun_start.size(), -1); logger.Info("RotationScaleMerge: ingested {} partial observations from {} frames ({} distinct hkl)", partials.size(), n_frames, rawrun_start.size()); @@ -497,84 +458,95 @@ void RotationScaleMerge::Ingest() { #endif } -// Compact `partials` (and the perm / raw-hkl tables that index it) down to the observations whose -// resolution is inside the requested range. See the header for why this cannot change any result: a raw -// hkl outside the range is group -1 in every pass, and dropping WHOLE runs on the same per-hkl resolution -// ComputeAsuGroups tests leaves rawrun_d, the surviving observations and their order untouched. +// Build `partials` (and the per-frame CSR, the finiteness mask and `perm`) from the source reflections, +// skipping the observations whose resolution is outside the requested range. See the header for why this +// cannot change any result: a raw hkl outside the range is group -1 in every pass, and skipping WHOLE +// runs on the same per-hkl resolution ComputeAsuGroups tests leaves rawrun_d, the surviving observations +// and their order untouched. // // It has to come after the incident flux has been measured (the flux meter is the frame's mean background // over ALL its reflections) and after the sort (which is what defines a raw-hkl run in the first place). -void RotationScaleMerge::DropOutOfRangeObservations() { - if (!d_min_limit && !d_max_limit) - return; - +void RotationScaleMerge::BuildInRangeObservations(const std::vector &keys) { const int n_run = static_cast(rawrun_start.size()); - const int n_obs = static_cast(partials.size()); + const int n_obs = static_cast(keys.size()); + const bool limited = d_min_limit.has_value() || d_max_limit.has_value(); - std::vector keep_run(n_run, 0); - std::vector keep(n_obs, 0); - ParallelChunks(n_run, nthreads, [&](int lo, int hi) { - for (int r = lo; r < hi; ++r) { - const float d = rawrun_d[r]; - if (!std::isfinite(d) || d <= 0.0f) continue; - if (d_min_limit && d < *d_min_limit) continue; - if (d_max_limit && d > *d_max_limit) continue; - keep_run[r] = 1; - const int p0 = rawrun_start[r], p1 = p0 + rawrun_count[r]; - for (int p = p0; p < p1; ++p) keep[perm[p]] = 1; // distinct runs own disjoint observations - } - }); - - // New index of every kept observation. Compacting in index order keeps each frame's observations - // contiguous and in the order they had, so the per-frame ranges still describe them and every - // per-frame sum keeps its sequence of roundings. - const int nt = static_cast(std::clamp(nthreads, 1, std::max(1, n_obs))); - const int chunk = (n_obs + nt - 1) / nt; - std::vector base(nt + 1, 0); - ParallelChunks(nt, nthreads, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - int32_t c = 0; - for (int i = t * chunk; i < std::min(n_obs, (t + 1) * chunk); ++i) c += keep[i]; - base[t + 1] = c; - } - }); - for (int t = 0; t < nt; ++t) base[t + 1] += base[t]; - const int32_t n_keep = base[nt]; - if (n_keep == n_obs) - return; - - std::vector new_idx(n_obs); - std::vector kept_obs(n_keep); - std::vector kept_finite(n_keep); - ParallelChunks(nt, nthreads, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - int32_t at = base[t]; - for (int i = t * chunk; i < std::min(n_obs, (t + 1) * chunk); ++i) { - if (!keep[i]) continue; - new_idx[i] = at; - kept_obs[at] = partials[i]; - kept_finite[at] = finite_ok[i]; - ++at; + std::vector keep_run(n_run, 1); + std::vector keep(n_obs, 1); + if (limited) { + keep.assign(n_obs, 0); + ParallelChunks(n_run, nthreads, [&](int lo, int hi) { + for (int r = lo; r < hi; ++r) { + const float d = rawrun_d[r]; + const bool in_range = std::isfinite(d) && d > 0.0f + && !(d_min_limit && d < *d_min_limit) && !(d_max_limit && d > *d_max_limit); + keep_run[r] = in_range ? 1 : 0; + if (!in_range) continue; + const int p0 = rawrun_start[r], p1 = p0 + rawrun_count[r]; + for (int p = p0; p < p1; ++p) keep[keys[p].idx] = 1; // distinct runs own disjoint observations } - } - }); - partials.swap(kept_obs); - finite_ok.swap(kept_finite); - kept_obs.clear(); kept_obs.shrink_to_fit(); - - frame_start.assign(n_frames, 0); - frame_count.assign(n_frames, 0); - for (int i = 0; i < n_keep; ) { - const int f = partials[i].frame; - int j = i; - while (j < n_keep && partials[j].frame == f) ++j; - frame_start[f] = i; - frame_count[f] = j - i; - i = j; + }); } - // The surviving runs, in the same order, each owning the same perm entries remapped to the new - // observation indices. + // Where each frame's reflections sit in the key numbering, and how many of them survive. Building in + // index order keeps each frame's observations contiguous and in the order they had, so the per-frame + // ranges still describe them and every per-frame sum keeps its sequence of roundings. + std::vector src_start(n_frames, 0); + { + int32_t at = 0; + for (int o = 0; o < n_frames; ++o) { + src_start[o] = at; + at += static_cast(partials_out[o].reflections.size()); + } + } + frame_count.assign(n_frames, 0); + ParallelFor(n_frames, nthreads, [&](int o) { + int32_t c = 0; + const int32_t n = static_cast(partials_out[o].reflections.size()); + for (int32_t j = 0; j < n; ++j) c += keep[src_start[o] + j]; + frame_count[o] = c; + }); + frame_start.assign(n_frames, 0); + int32_t n_keep = 0; + for (int o = 0; o < n_frames; ++o) { + frame_start[o] = n_keep; + n_keep += frame_count[o]; + } + + // The observations themselves, one frame per thread. finite_ok is the per-obs AcceptReflection + // finiteness (immutable) - it lets ComputeAsuGroups stamp the ASU-group id per obs from a flat + // 1-byte array instead of re-reading the fat Obs struct for every space group. + partials.clear(); + partials.resize(n_keep); + finite_ok.assign(n_keep, 0); + std::vector new_idx(n_obs); // only written where keep[i] + ParallelFor(n_frames, nthreads, [&](int o) { + int32_t src = src_start[o]; + int32_t at = frame_start[o]; + for (const auto &r : partials_out[o].reflections) { + if (!keep[src]) { ++src; continue; } + new_idx[src++] = at; + Obs &obs = partials[at]; + obs.h = r.h; obs.k = r.k; obs.l = r.l; + obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; + // Divide the incident flux out of rlp here: rlp is the numerator of corr, so the fits + // below see only what is left. + obs.rlp = static_cast(r.rlp / frame_flux[o]); + obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg; + obs.px = r.predicted_x; obs.py = r.predicted_y; + obs.image_number = r.image_number; + obs.frame = o; + obs.on_ice = r.on_ice_ring ? 1 : 0; + obs.corr = r.image_scale_corr; + obs.group = -1; + finite_ok[at] = (std::isfinite(obs.I) && std::isfinite(obs.rlp) && obs.rlp != 0.0f + && std::isfinite(obs.sigma) && obs.sigma > 0.0f) ? 1 : 0; + ++at; + } + }); + + // The surviving runs, in the same order, each owning the same sorted keys remapped to the + // observation indices they were just built at. int n_keep_run = 0; std::vector dst_run(n_run, -1), dst_start(n_run, 0); { @@ -582,8 +554,8 @@ void RotationScaleMerge::DropOutOfRangeObservations() { for (int r = 0; r < n_run; ++r) if (keep_run[r]) { dst_run[r] = n_keep_run++; dst_start[r] = at; at += rawrun_count[r]; } } - std::vector new_perm(n_keep), ns(n_keep_run), nc(n_keep_run), - nh(n_keep_run), nk(n_keep_run), nl(n_keep_run); + perm.resize(n_keep); + std::vector ns(n_keep_run), nc(n_keep_run), nh(n_keep_run), nk(n_keep_run), nl(n_keep_run); std::vector nd(n_keep_run); ParallelChunks(n_run, nthreads, [&](int lo, int hi) { for (int r = lo; r < hi; ++r) { @@ -593,16 +565,16 @@ void RotationScaleMerge::DropOutOfRangeObservations() { ns[rr] = at; nc[rr] = rawrun_count[r]; nh[rr] = rawrun_h[r]; nk[rr] = rawrun_k[r]; nl[rr] = rawrun_l[r]; nd[rr] = rawrun_d[r]; for (int32_t p = 0; p < rawrun_count[r]; ++p) - new_perm[at + p] = new_idx[perm[rawrun_start[r] + p]]; + perm[at + p] = new_idx[keys[rawrun_start[r] + p].idx]; } }); - perm.swap(new_perm); rawrun_start.swap(ns); rawrun_count.swap(nc); rawrun_h.swap(nh); rawrun_k.swap(nk); rawrun_l.swap(nl); rawrun_d.swap(nd); - logger.Info("RotationScaleMerge: dropped {} of {} observations ({} of {} distinct hkl) outside the " - "requested resolution range - nothing downstream could have merged them", - n_obs - n_keep, n_obs, n_run - n_keep_run, n_run); + if (n_keep < n_obs) + logger.Info("RotationScaleMerge: dropped {} of {} observations ({} of {} distinct hkl) outside the " + "requested resolution range - nothing downstream could have merged them", + n_obs - n_keep, n_obs, n_run - n_keep_run, n_run); } namespace { @@ -2544,7 +2516,7 @@ namespace { } RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for_search, - bool fulls_resident) { + bool fulls_resident, bool full_stats) { // A full is usable for the merge / error model if it passes AddImage's filters (with the current // ice context). group >= 0 already encodes "not absent and passes AcceptReflection". auto usable_merge = [&](const Obs &o) { @@ -3020,9 +2992,11 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool em.a, em.b, em.isa, result.isa_asymptotic, error_model_chi2); } - AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); - // French-Wilson (F, and F(+)/F(-) from the anomalous split) is deferred until after the anomalous - // accumulator below has attached I(+)/I(-), so the two hands get their amplitudes in one pass. + // R-free flags and French-Wilson amplitudes are fields of the WRITTEN reflection, so a merge that + // is not written needs neither. French-Wilson is deferred until after the anomalous accumulator + // below has attached I(+)/I(-), so the two hands get their amplitudes in one pass. + if (full_stats) + AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); if (reject_count > 0) logger.Info("Merge outlier rejection: dropped {} observations", reject_count); @@ -3069,71 +3043,77 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } } - // R_meas: re-walk the fulls (Mask = cell only; no ice / error-model), accumulate - // |I_i - | per reflection. - struct RmeasObs { double sum_abs_dev = 0, sum_I = 0; int n = 0, shell = -1; }; - std::vector rmeas(n_groups); - bool did_gpu_rmeas = false; -#ifdef JFJOCH_USE_CUDA - if (use_gpu_merge) { - // Per-group R_meas + usable count on the GPU; the shell is assigned per group (its fulls share d). - std::vector rabsdev(n_groups), rsumI(n_groups); - std::vector rn(n_groups), rnusable(n_groups); - gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rn.data(), rnusable.data()); - for (int g = 0; g < n_groups; ++g) { - if (rnusable[g] == 0) continue; - const auto shell = shells.GetShell(acc[g].d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - // Count the MERGED population, not the R_meas one. The R_meas re-walk deliberately - // ignores the ice flag on a search pass, so its count includes observations that never - // entered `unique` - which inflates the reported multiplicity of whatever shell they - // land in. acc[g].nh is what actually went into this group's mean, and it is zero for a - // group the merge dropped entirely. (acc[g].d is NaN for such a group, so GetShell above - // already declines it; this is the same statement made where it counts.) - sa[*shell].total_obs += static_cast(acc[g].nh[0] + acc[g].nh[1]); - if (std::isfinite(merged_I[g]) && rn[g] > 0) { - auto &r = rmeas[g]; - r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.n = rn[g]; r.shell = *shell; - } - } - did_gpu_rmeas = true; - } -#endif - if (!did_gpu_rmeas) - for (const auto &o : fulls) { - if (o.group < 0) continue; - if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge -> also out of R_meas - if (!frame_cell_ok[o.frame]) continue; - if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue; - if (o.partiality < min_partiality) continue; - const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; - if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; - const auto shell = shells.GetShell(o.d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - // Only what the merge kept counts towards multiplicity - see the GPU branch above. The - // R_meas accumulation below keeps its own, wider filter. - if (usable_merge(o)) - sa[*shell].total_obs++; - if (std::isfinite(merged_I[o.group])) { - auto &r = rmeas[o.group]; - r.sum_abs_dev += std::fabs(static_cast(I_corr) - merged_I[o.group]); - r.sum_I += I_corr; r.n++; r.shell = *shell; - } - } std::vector rmeas_num(n_shells, 0.0), rmeas_den(n_shells, 0.0); double rmeas_num_all = 0.0, rmeas_den_all = 0.0; - for (const auto &r : rmeas) { - if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue; - const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); - rmeas_num[r.shell] += factor * r.sum_abs_dev; rmeas_den[r.shell] += r.sum_I; - rmeas_num_all += factor * r.sum_abs_dev; rmeas_den_all += r.sum_I; + // R_meas and the observation counts behind the reported multiplicity. Report-only, and the last + // full pass over the observations there is - so, like the correction surfaces, it is not made on a + // merge whose numbers nobody reads. + if (full_stats) { + // R_meas: re-walk the fulls (Mask = cell only; no ice / error-model), accumulate + // |I_i - | per reflection. + struct RmeasObs { double sum_abs_dev = 0, sum_I = 0; int n = 0, shell = -1; }; + std::vector rmeas(n_groups); + bool did_gpu_rmeas = false; +#ifdef JFJOCH_USE_CUDA + if (use_gpu_merge) { + // Per-group R_meas + usable count on the GPU; the shell is assigned per group (its fulls share d). + std::vector rabsdev(n_groups), rsumI(n_groups); + std::vector rn(n_groups), rnusable(n_groups); + gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rn.data(), rnusable.data()); + for (int g = 0; g < n_groups; ++g) { + if (rnusable[g] == 0) continue; + const auto shell = shells.GetShell(acc[g].d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + // Count the MERGED population, not the R_meas one. The R_meas re-walk deliberately + // ignores the ice flag on a search pass, so its count includes observations that never + // entered `unique` - which inflates the reported multiplicity of whatever shell they + // land in. acc[g].nh is what actually went into this group's mean, and it is zero for a + // group the merge dropped entirely. (acc[g].d is NaN for such a group, so GetShell above + // already declines it; this is the same statement made where it counts.) + sa[*shell].total_obs += static_cast(acc[g].nh[0] + acc[g].nh[1]); + if (std::isfinite(merged_I[g]) && rn[g] > 0) { + auto &r = rmeas[g]; + r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.n = rn[g]; r.shell = *shell; + } + } + did_gpu_rmeas = true; + } +#endif + if (!did_gpu_rmeas) + for (const auto &o : fulls) { + if (o.group < 0) continue; + if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge -> also out of R_meas + if (!frame_cell_ok[o.frame]) continue; + if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue; + if (o.partiality < min_partiality) continue; + const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; + if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; + const auto shell = shells.GetShell(o.d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + // Only what the merge kept counts towards multiplicity - see the GPU branch above. The + // R_meas accumulation below keeps its own, wider filter. + if (usable_merge(o)) + sa[*shell].total_obs++; + if (std::isfinite(merged_I[o.group])) { + auto &r = rmeas[o.group]; + r.sum_abs_dev += std::fabs(static_cast(I_corr) - merged_I[o.group]); + r.sum_I += I_corr; r.n++; r.shell = *shell; + } + } + for (const auto &r : rmeas) { + if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue; + const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); + rmeas_num[r.shell] += factor * r.sum_abs_dev; rmeas_den[r.shell] += r.sum_I; + rmeas_num_all += factor * r.sum_abs_dev; rmeas_den_all += r.sum_I; + } } // ---- Anomalous split (always, even when the merge is Friedel-averaged): for each acentric // reflection, the inverse-variance I(+)/I(-) from the SAME scaled fulls (so it never touches the // Friedel-merged IMEAN, scaling or error model above). Lets I(+)/I(-) and F(+)/F(-) be written by // default without scaling anomalously; a reflection with only one mate, or a centric, is left - // without the split. Skipped for the P1 search pass, which has no use for it. ---- + // without the split. Skipped for the P1 search pass, which has no use for it, and for a merge that + // is not written. ---- struct AnomExport { float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; }; std::unordered_map anom_export; // SigAno = <|I(+)-I(-)|>/ accumulated on the standard report shells (num=sum|dI|, @@ -3141,7 +3121,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool std::vector sig_num(n_shells, 0.0), sig_den(n_shells, 0.0); double sig_num_all = 0.0, sig_den_all = 0.0; size_t sig_n = 0; - if (!for_search) { + if (!for_search && full_stats) { const int sg_num = x.GetSpaceGroupNumber().value_or(1); const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, sg_num); const gemmi::GroupOps gops = gemmi::find_spacegroup_by_number(sg_num)->operations(); @@ -3233,15 +3213,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } // French-Wilson amplitudes for IMEAN and (now that they are attached) each Bijvoet hand. - FrenchWilsonOptions fw_opts; - fw_opts.num_threads = static_cast(nthreads); - ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1), fw_opts); + if (full_stats) { + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = static_cast(nthreads); + ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1), fw_opts); + } logger.Info("Merge complete ({} unique reflections)", result.merged.size()); return result; } -RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { +RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_stats) { const int sg_number = x.GetSpaceGroupNumber().value_or(1); HKLKeyGenerator keygen(merge_friedel, sg_number); @@ -3545,26 +3527,29 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { // diffracted-beam direction), each an alternating multiplicative fit of the fulls' corr against // the merged reference. Cheap host loops over the downloaded fulls; applied only on the final // in-symmetry merge, never the P1 space-group search pass (corrections there add risk and can - // perturb the symmetry determination). --- + // perturb the symmetry determination) and never on a merge that is not an output: a surface + // is a correction to the merged INTENSITIES, and a pass whose intensities are thrown away has + // nowhere to put one - the fitted gain lives on the fulls, which the next pass rebuilds. --- + const bool refine_surfaces = !for_search && full_stats; [[maybe_unused]] const bool corrections = - !for_search && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0); + refine_surfaces && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0); // Radiation-damage monitor: measure the per-batch relative-B on the scaled fulls BEFORE any correction // (report-only; captures the full damage signature, not a residual). Skipped on the P1 search pass. - if (!for_search) + if (refine_surfaces) MeasureRadiationDamageB(n_groups); // Sweep-quality diagnostic, on the per-frame scale the partial scaling just fitted (the flux is // already out of it) and the per-frame CC computed above. Report-only; drops nothing. - if (!for_search) + if (refine_surfaces) MeasureSweepQuality(partial_scaled, cc, cc_n); - if (!for_search && refine_decay_b) + if (refine_surfaces && refine_decay_b) RefineDecay(n_groups); - if (!for_search && relative_b_deg > 0.0) + if (refine_surfaces && relative_b_deg > 0.0) RefineRelativeB(n_groups); // per-batch relative-B on top of the single decay slope - if (!for_search && absorption_iter > 0) + if (refine_surfaces && absorption_iter > 0) RefineAbsorption(absorption_iter, n_groups); - if (!for_search && modulation_iter > 0) + if (refine_surfaces && modulation_iter > 0) RefineModulation(modulation_iter, n_groups); - if (!for_search && absorption_iter > 0) + if (refine_surfaces && absorption_iter > 0) RefineAbsorptionTime(absorption_iter, n_groups); // last: the static surfaces get first claim #ifdef JFJOCH_USE_CUDA // The corrections mutate the host fulls' corr; when the merge runs on the resident (GPU) fulls, push @@ -3580,6 +3565,6 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { #endif // --- 5. Error model + merge + statistics. --- - auto r = MergeAndStats(n_groups, for_search, combined_on_gpu && scaled_fulls_on_gpu); + auto r = MergeAndStats(n_groups, for_search, combined_on_gpu && scaled_fulls_on_gpu, full_stats); return r; } diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 5255af16..35fadbbb 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -71,7 +71,14 @@ public: // for the space group currently set on the experiment, reusing the ingested buffers. // for_search: the de-novo P1 pass whose merged intensities feed the space-group search - ice-ring // reflections are dropped from the merge and the error model (kept otherwise, for completeness). - Result Run(bool for_search); + // full_stats: these merged intensities are an OUTPUT. False on the rotation two-pass geometry + // pre-pass, whose merge exists only to choose the space group and post-refine the geometry and + // whose reflections are never written: the correction surfaces, the report-only diagnostics, the + // R_meas re-walk, the anomalous split, the R-free flags and the French-Wilson amplitudes are then + // all skipped, because computing them fills in fields nothing reads. What the pre-pass IS read for + // - the merged intensities themselves, the error model, and the completeness / CC1/2 the second + // pass is judged against - is computed either way. + Result Run(bool for_search, bool full_stats); // Override the high-resolution cut for the next Run() - used to gate the de-novo P1 search pass at // >= 1 without cutting the final in-symmetry merge. Reset to the manual limit afterwards. @@ -99,6 +106,17 @@ private: int32_t group; // dense ASU-group id for the current space group; <0 = never mergeable }; + // The narrow per-observation record the ingest sort orders: the raw hkl the runs are cut on, the + // frame position that breaks a tie inside one, the observation's own index (which makes the order + // total - see the .cpp), and the resolution the range test reads. Twenty-four bytes against the + // Obs's eighty, and it is all the ingest needs before it knows which observations survive. + struct SortKey { + int32_t h, k, l; + float image_number; + int32_t idx; + float d; + }; + const DiffractionExperiment &x; std::vector &partials_out; // written back at the end of scaling std::optional reference_cell; @@ -224,14 +242,17 @@ private: // background is a usable flux meter and what it costs when it is not. void MeasureIncidentFlux(const std::vector &mean_bkg); - // Drop the observations whose resolution can never be in range: --scaling-high/low-resolution are the - // coarsest limits any Run() uses (the space-group search only ever RAISES d_min), and an out-of-range - // raw hkl gets group -1 in every pass, which keeps it out of the scaling reference, the per-frame fit, - // the combine, the merge and the error model alike. On a crystal that integrates to the detector - // corner and merges well short of it that is most of the array. Whole raw-hkl RUNS are dropped, on the - // same per-hkl resolution ComputeAsuGroups tests, so what survives - and the order of every sum formed - // over it - is exactly what it was. No-op when no manual limit was given. - void DropOutOfRangeObservations(); + // Build the flat `partials` array (and the per-frame CSR, the finiteness mask and `perm`) from the + // source reflections, skipping the observations whose resolution can never be in range: + // --scaling-high/low-resolution are the coarsest limits any Run() uses (the space-group search only + // ever RAISES d_min), and an out-of-range raw hkl gets group -1 in every pass, which keeps it out of + // the scaling reference, the per-frame fit, the combine, the merge and the error model alike. On a + // crystal that integrates to the detector corner and merges well short of it that is most of the + // array, and an eighty-byte record built for it is eighty bytes written and then thrown away. Whole + // raw-hkl RUNS are skipped, on the same per-hkl resolution ComputeAsuGroups tests, so what survives - + // and the order of every sum formed over it - is exactly what it would have been had the whole array + // been built and then filtered. Everything is built when no manual limit was given. + void BuildInRangeObservations(const std::vector &keys); // Compute the dense ASU-group id for the current space group by grouping the (pre-sorted) raw-hkl // runs by their ASU key - one gemmi ASU reduction per distinct raw hkl, not per observation. Fills @@ -347,5 +368,6 @@ private: // Error model + merge + statistics over the fulls (the last stage). n_groups is the fulls group count. // fulls_resident: the (scaled) fulls + their group CSR are still on the GPU, so the em-stats / samples // / merge-accumulate / R_meas reductions run there (only per-group + samples come back). - Result MergeAndStats(int n_groups, bool for_search, bool fulls_resident); + // full_stats: see Run(). + Result MergeAndStats(int n_groups, bool for_search, bool fulls_resident, bool full_stats); }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 6255c59a..a8465297 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2090,7 +2090,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b auto scale_and_merge = [&](const std::string &label, bool for_search) -> ScaleMergeResult { if (rsm) { phase("Scale/combine/merge (" + label + ")"); - auto r = rsm->Run(for_search); + // The geometry pre-pass merges only to choose the space group and to give the quality + // guard something to judge the second pass against; its reflections are never written, + // so the parts of the merge that only fill in an output file are skipped there. + auto r = rsm->Run(for_search, /*full_stats=*/!geometry_prepass); result.error_model_isa = r.isa; result.error_model_isa_asymptotic = r.isa_asymptotic; result.error_model_a = r.error_model_a; @@ -2564,11 +2567,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto twin_sg_number = experiment_.GetSpaceGroupNumber(); const gemmi::SpaceGroup *twin_sg = twin_sg_number ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; - result.twinning = AnalyzeTwinning(sm.merged, twin_sg); + // Not on the geometry pre-pass: the analysis goes into that pass's statistics text and its + // written reflections, and neither survives the run. The promotion flag below is a different + // thing - it is what the SEARCH did, the second pass reads it, and it is set either way. + if (!geometry_prepass) { + result.twinning = AnalyzeTwinning(sm.merged, twin_sg); + stats_text << TwinningAnalysisToText(result.twinning) << "\n"; + } // Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the // search itself made, so the text cannot claim "no twin law exists" on its own say-so. result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group; - stats_text << TwinningAnalysisToText(result.twinning) << "\n"; // Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the // crystal was mounted, which the user can change), warned about when the two nearly coincide. @@ -2639,8 +2647,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } // Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not - // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log. - { + // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log, none + // of which the geometry pre-pass produces. + if (!geometry_prepass) { const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged); sm.statistics.wilson_b = wilson.b; sm.statistics.wilson_b_correlation = wilson.correlation; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index fedc5e09..fae3e0a6 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1424,7 +1424,7 @@ static int RunRugnux(int argc, char **argv) { RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(), scaling_iter, nthreads, logger); rsm.Ingest(); - auto r = rsm.Run(false); + auto r = rsm.Run(false, /*full_stats=*/true); merged_reflections = std::move(r.merged); merged_statistics = std::move(r.statistics); error_model_isa = r.isa;