diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 086f71c34..be65e5678 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -239,6 +239,26 @@ namespace { return out; } + // A sum over [0, n) formed on all threads, in a split that depends on n alone: add(part, i) folds + // element i into its block's partial, and the partials are merged in block order (ParallelBlocks), so + // the sum has the same bits at every thread count - though not those of the one serial walk. + template + S FixedBlockSum(int n, size_t nthreads, const S &zero, Add add, Merge merge) { + std::vector part(ReductionBlocks(n), zero); + ParallelBlocks(n, nthreads, [&](int b, int lo, int hi) { + for (int i = lo; i < hi; ++i) add(part[b], i); + }); + S total = zero; + for (const S &p : part) merge(total, p); + return total; + } + // Element-wise sum of equal-length arrays, the merge FixedBlockSum takes for per-batch / per-shell sums. + template + void AddArrays(std::array, K> &to, const std::array, K> &from) { + for (size_t k = 0; k < K; ++k) + for (size_t j = 0; j < to[k].size(); ++j) to[k][j] += from[k][j]; + } + struct ScaleObs { double coeff, Iobs, weight; }; // Per-frame scale, linear in G: the weighted least-squares slope of the observations on their @@ -2008,16 +2028,19 @@ void RotationScaleMerge::RefineDecay(int n_groups) { sw[g] = s_w; swI[g] = s_wI; } }); - double Sw = 0, Sx = 0, Sy = 0, Sxx = 0, Sxy = 0; - for (const DecayTerm &o : t.term) { - if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) continue; + using Lin = std::array; // Sw, Sx, Sy, Sxx, Sxy + const Lin S = FixedBlockSum(static_cast(t.term.size()), nthreads, Lin{}, + [&](Lin &a, int k) { + const DecayTerm &o = t.term[k]; + if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) return; const double Iref = swI[o.group] / sw[o.group]; const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; - if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue; + if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) return; const double w = (Is / sc) * (Is / sc); const double x = (o.image_number - fcenter) * s2_of(o.d), y = std::log(Iref / Is); - Sw += w; Sx += w * x; Sy += w * y; Sxx += w * x * x; Sxy += w * x * y; - } + a[0] += w; a[1] += w * x; a[2] += w * y; a[3] += w * x * x; a[4] += w * x * y; + }, [](Lin &a, const Lin &b) { for (int j = 0; j < 5; ++j) a[j] += b[j]; }); + const double Sw = S[0], Sx = S[1], Sy = S[2], Sxx = S[3], Sxy = S[4]; const double var = (Sw > 0.0) ? Sxx - Sx * Sx / Sw : 0.0; // Clamp only to guard against a near-collinear (var ~ 0) blow-up; generous enough to reach very // strong damage (slope = 2 dB/dframe, so +-1.0 admits total relative-B up to ~n_frames/2 A^2). @@ -2045,14 +2068,16 @@ void RotationScaleMerge::RefineDecay(int n_groups) { sw[g] = s_w; swI[g] = s_wI; } }); - double num = 0.0, den = 0.0; - for (const DecayTerm &o : t.term) { - if ((o.frame & 1) != parity || sw[o.group] <= 0.0) continue; + using ND = std::array; + const ND nd = FixedBlockSum(static_cast(t.term.size()), nthreads, ND{}, [&](ND &a, int k) { + const DecayTerm &o = t.term[k]; + if ((o.frame & 1) != parity || sw[o.group] <= 0.0) return; const double Iref = swI[o.group] / sw[o.group]; const double Is = static_cast(o.I) * o.corr * decay_factor(o, slope); - if (!std::isfinite(Iref) || Iref <= 0.0) continue; - num += std::abs(Is - Iref); den += Iref; - } + if (!std::isfinite(Iref) || Iref <= 0.0) return; + a[0] += std::abs(Is - Iref); a[1] += Iref; + }, [](ND &a, const ND &b) { a[0] += b[0]; a[1] += b[1]; }); + const double num = nd[0], den = nd[1]; return den > 0.0 ? num / den : 0.0; }; @@ -2158,16 +2183,20 @@ std::vector RotationScaleMerge::FitRelativeBCurve(const DecayTerms &t, i sw[g] = s_w; swI[g] = s_wI; } }); - std::vector num(n_batch, 0.0), den(n_batch, 0.0); - for (const DecayTerm &o : t.term) { - if ((gparity >= 0 && (o.group & 1) != gparity) || sw[o.group] <= 0.0) continue; + using Batches = std::array, 2>; // num, den per batch + const Batches nd = FixedBlockSum(static_cast(t.term.size()), nthreads, + Batches{std::vector(n_batch, 0.0), std::vector(n_batch, 0.0)}, + [&](Batches &a, int k) { + const DecayTerm &o = t.term[k]; + if ((gparity >= 0 && (o.group & 1) != gparity) || sw[o.group] <= 0.0) return; const double Iref = swI[o.group] / sw[o.group]; const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; - if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue; + if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) return; const double w = (Is / sc) * (Is / sc); const double s2 = s2_of(o.d), y = std::log(Iref / Is); - num[batch_of(o)] += w * s2 * y; den[batch_of(o)] += w * s2 * s2; - } + a[0][batch_of(o)] += w * s2 * y; a[1][batch_of(o)] += w * s2 * s2; + }, AddArrays<2>); + const std::vector &num = nd[0], &den = nd[1]; std::vector b = SolveCurvatureSmoothedB(num, den, RELATIVE_B_MAX); double dw = 0.0, dbw = 0.0; for (int c = 0; c < n_batch; ++c) { dw += den[c]; dbw += den[c] * b[c]; } @@ -2253,12 +2282,16 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { // spread over all of it would spend most of the grid on noise, leaving a batch with too few to fit at // all. Pooled over the whole run a shell either has signal or it has none, so walk out to the first // that has none and lay the shells inside that. - std::vector rw(MONITOR_SHELLS, 0.0), rwI(MONITOR_SHELLS, 0.0); - for (const DecayTerm &o : t.term) { + using Shells = std::array, 2>; // rw, rwI per shell + const Shells r_shell = FixedBlockSum(static_cast(t.term.size()), nthreads, + Shells{std::vector(MONITOR_SHELLS, 0.0), std::vector(MONITOR_SHELLS, 0.0)}, + [&](Shells &a, int j) { + const DecayTerm &o = t.term[j]; const double sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); const int k = shell_of(s2_of(o.d)); - rw[k] += w; rwI[k] += w * static_cast(o.I) * o.corr; - } + a[0][k] += w; a[1][k] += w * static_cast(o.I) * o.corr; + }, AddArrays<2>); + const std::vector &rw = r_shell[0], &rwI = r_shell[1]; int n_live = 0; while (n_live < MONITOR_SHELLS && rw[n_live] > 0.0 && rwI[n_live] / std::sqrt(rw[n_live]) >= MONITOR_MIN_SHELL_ISIGMA) @@ -2278,15 +2311,20 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { // mean is well determined where a single observation is not, admits negative intensities, and carries // its own I/sigma - which is what says whether the batch can be measured at all. const int n_cell = n_batch * MONITOR_SHELLS; - std::vector sw(n_cell, 0.0), swI(n_cell, 0.0), swR(n_cell, 0.0), sws2(n_cell, 0.0); - for (const DecayTerm &o : t.term) { - if (sw0[o.group] <= 0.0) continue; + using Cells = std::array, 4>; // sw, swI, swR, sws2 per (batch, shell) + const std::vector zero_cells(n_cell, 0.0); + const Cells cells = FixedBlockSum(static_cast(t.term.size()), nthreads, + Cells{zero_cells, zero_cells, zero_cells, zero_cells}, + [&](Cells &a, int j) { + const DecayTerm &o = t.term[j]; + if (sw0[o.group] <= 0.0) return; const double Iref = swI0[o.group] / sw0[o.group]; - if (!std::isfinite(Iref)) continue; + if (!std::isfinite(Iref)) return; const double s2 = s2_of(o.d), sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); const int i = batch_of(o) * MONITOR_SHELLS + shell_of(s2); - sw[i] += w; swI[i] += w * static_cast(o.I) * o.corr; swR[i] += w * Iref; sws2[i] += w * s2; - } + a[0][i] += w; a[1][i] += w * static_cast(o.I) * o.corr; a[2][i] += w * Iref; a[3][i] += w * s2; + }, AddArrays<4>); + const std::vector &sw = cells[0], &swI = cells[1], &swR = cells[2], &sws2 = cells[3]; // Per batch: slope of ln( / ) against s^2 over the shells that hold signal, weighted by // the shell's own (I/sigma)^2 since that is the precision of its log. Fitted WITH an intercept and the // normal equations centred, so that only the resolution-DEPENDENT part reaches the curve: a batch that