diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 198c6366c..6e3d0b8e7 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -22,13 +22,12 @@ class ScalingSettings { bool merge_friedel = true; std::optional high_resolution_limit_A; - // Low-resolution limit for scaling and merging. The beam-stop mask does not replace this: the stop - // suppresses air scatter well beyond the shadow it masks, so in the pixels the mask leaves open the - // background is measurably depressed - 30-44% of the field value on a real geometry, recovering to - // 95% of it only around d ~ 52 A. A background ring inside that zone over-estimates the background - // and the intensity comes out negative (74% negative past 50 A there, 2% just inside). The zone is - // set by the stop's angular size, so it tracks wavelength rather than detector distance, which is - // what makes a fixed d limit reasonable. 50 A is also XDS's own default when the input leaves it out. + // Low-resolution limit of the scale FIT. Reflections coarser than it are still merged and written - + // they are what maps and refinement want of the lowest resolutions - but are kept out of every fit + // of the scale model: on a large cell a few dozen of them, measured against a background the beam + // stop depresses, bent the model enough to move the merged 10-50 A intensities by tens of percent. + // Which coarse reflections are measured at all is the beam-stop mask's decision (a reflection whose + // signal disk reaches the shadow is not integrated, MXAnalysisWithoutFPGA). 50 A is XDS's default. std::optional low_resolution_limit_A = 50.0; std::optional wedge_for_scaling; std::optional forced_mosaicity; // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index ae4bdfd32..938599b1b 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -369,8 +369,9 @@ An internal-consistency term can inflate uncertainties when multiple observation ### 10.5 Merging statistics -The shells are **nine bins of equal width in $1/d^2$**, laid between the **declared** low-resolution -limit (`--scaling-low-resolution`, or the whole sphere where it is switched off) and the +The shells are **nine bins of equal width in $1/d^2$**, laid between the whole sphere ($1/d^2 = 0$; the +low-resolution limit `--scaling-low-resolution` bounds only the scale fit, so the coarsest reflections are +merged too, and those the beam stop hides count as missing) or the declared range of `--report-resolution`, and the highest-resolution reflection the merge actually kept — XDS's rule and XDS's count, so at the same resolution limits the two programs' tables have the same shell boundaries and can be read row for row. `--resolution-shells` changes the count; the binning rule does not change with it. diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 637637fb2..cc4f8fc27 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -441,7 +441,7 @@ Scaling and merging: | `--capture-uncertainty ` | rot3d: systematic sigma on under-captured fulls, ~num·(1−captured_fraction)·I (default: 1.0 for rotation, 0 otherwise) | | `--min-captured-fraction ` | rot3d: drop a combined full whose rocking curve was captured below this fraction — edge-of-sweep truncated fulls (default: 0.7 for rotation, 0 otherwise; 0 = off) | | `--scaling-high-resolution ` | High-resolution limit for scaling, Å — manual override (default: no limit; disables the automatic cutoff below) | -| `--scaling-low-resolution ` | Low-resolution limit for scaling and merging, Å (default: 50, which is also XDS's own default; 0 removes the limit). The beam stop suppresses air scatter well beyond the shadow it casts, so the background stays depressed across pixels the beam-stop mask leaves open — measured at 30–44% of the field value, recovering only around 50 Å. A background ring in that zone over-estimates the background, so reflections coarser than the limit mostly integrate negative. Beam-stop masking therefore does not make this redundant. A large cell (a few hundred Å) does have real reflections coarser than 50 Å, but on the geometries measured here they fall in that zone and are not usable as integrated; raise the limit only if the low-resolution shell statistics justify it | +| `--scaling-low-resolution ` | Low-resolution limit of the scale fit, Å (default: 50, which is also XDS's own default; 0 removes the limit). Reflections coarser than it are merged and written but kept out of every fit of the scale model: on a large cell a few dozen of them, measured on a background the beam stop depresses, bent the model enough to move the merged 10–50 Å intensities by tens of percent. Which of them are measured at all is the beam-stop mask's decision - a reflection whose signal disk reaches the shadow is not integrated (`--detect-beam-stop`) | | `--resolution-cutoff ` | Automatic high-resolution cutoff for the written reflections and reported shells: `cc-logistic` \| `off` (default: `cc-logistic`; ignored when `--scaling-high-resolution` is set) | | `--resolution-cc-target ` | CC1/2 target defining the `cc-logistic` fall-off (default: 0.30) | | `--resolution-shells ` | Number of resolution shells in the reported statistics table (default: 9). The bins are equal steps in 1/d² between the lowest- and highest-resolution reflection merged, which is XDS's rule, and 9 is XDS's count — so at the same resolution limits the two tables have the same shells and can be read row for row | diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 16bfa0fcd..b07eb3593 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -62,7 +62,6 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) scaling_settings(x.GetScalingSettings()), indexing_settings(x.GetIndexingSettings()), high_resolution_limit(scaling_settings.GetHighResolutionLimit_A()), - low_resolution_limit(scaling_settings.GetLowResolutionLimit_A()), // A min-image-CC of 0 (the default) means "no limit": leave the optional // empty so the per-image CC cut is inactive. Otherwise a 0.0 threshold // would silently drop every image with a non-positive per-image CC. @@ -89,7 +88,7 @@ size_t MergeOnTheFly::Contributions(const IntegrationOutcome &outcome, std::vect if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) continue; - if (!AcceptReflection(r, high_resolution_limit, low_resolution_limit)) + if (!AcceptReflection(r, high_resolution_limit, std::nullopt)) continue; if (exclude_ice_rings && r.on_ice_ring) continue; @@ -267,7 +266,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc continue; if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) continue; - if (!AcceptReflection(r, high_resolution_limit, low_resolution_limit)) + if (!AcceptReflection(r, high_resolution_limit, std::nullopt)) continue; if (exclude_ice_rings && r.on_ice_ring) continue; @@ -641,7 +640,9 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m auto d_min_limit_A = d_min_override.has_value() ? d_min_override : scaling_settings.GetHighResolutionLimit_A(); - auto d_max_limit_A = scaling_settings.GetLowResolutionLimit_A(); + // The low-resolution limit bounds the scale fit, not the merge (ScalingSettings), so the statistics + // run to the coarsest reflection unless a range is declared. + std::optional d_max_limit_A; if (declared_range) { d_min_limit_A = declared_range->d_min; d_max_limit_A = std::isfinite(declared_range->d_max) diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 6c8ca78a7..3fe75bea5 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -288,7 +288,6 @@ class MergeOnTheFly { std::optional reference_cell; std::optional high_resolution_limit; - std::optional low_resolution_limit; std::optional image_cc_limit; // Apply image_cc_limit in Mask(). One flag for the whole engine, not a per-call argument, so the // merge, the error model and MergeStats can never disagree about which images are in. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 4761d0954..a0e2b527d 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -350,7 +350,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, const auto s = x.GetScalingSettings(); min_partiality = s.GetMinPartiality(); d_min_limit = s.GetHighResolutionLimit_A(); - d_max_limit = s.GetLowResolutionLimit_A(); + d_max_scale = s.GetLowResolutionLimit_A(); merge_friedel = s.GetMergeFriedel(); capture_uncertainty_coeff = s.GetCaptureUncertaintyCoeff(); min_captured_fraction = s.GetMinCapturedFraction(); @@ -606,7 +606,7 @@ void RotationScaleMerge::IngestHost() { // genuinely share (h,k,l) and image_number: the predictor emits both intersections of a // reflection's rotation circle with the Ewald sphere, and near the blind region the two are close // enough in angle that both land on the same frame. Without a final key their relative order was - // whatever the sort happened to produce, and the combine reads it - it takes on_ice from the first + // whatever the sort happened to produce, and the combine reads it - it takes unscaled from the first // 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. @@ -755,7 +755,7 @@ void RotationScaleMerge::IngestOnDevice() { std::vector mean_bkg(n_frames); gpu_->IngestMeanBackground(mean_bkg.data()); MeasureIncidentFlux(mean_bkg); - const int n = gpu_->IngestSort(d_min_limit, d_max_limit); + const int n = gpu_->IngestSort(d_min_limit, std::nullopt); n_partials_obs = n; RotationScaleMergeGPU::IngestLayout layout; gpu_->IngestGetLayout(layout); @@ -777,7 +777,8 @@ void RotationScaleMerge::IngestOnDevice() { return; auto pos = std::make_shared(); pos->px.resize(n); pos->py.resize(n); - gpu_->IngestObservations(refl.data(), count.data(), nthreads, frame_flux.data(), pos->px.data(), pos->py.data()); + gpu_->IngestObservations(refl.data(), count.data(), nthreads, frame_flux.data(), + d_max_scale ? static_cast(*d_max_scale) : INFINITY, pos->px.data(), pos->py.data()); obs_position = std::move(pos); const SmoothedGeometry sg = SmoothedFrameGeometry(); @@ -997,7 +998,7 @@ void RotationScaleMerge::IngestDevice(const RotationScaleMerge *observations_fro fbuf[3][i] = o.partiality; fbuf[4][i] = r.zeta; fbuf[5][i] = o.corr; fbuf[6][i] = r.bkg; fbuf[7][i] = r.var_bkg; fbuf[8][i] = r.image_number; fbuf[9][i] = r.d; - flags[i] = static_cast((r.on_ice_ring ? 1 : 0) | (ClippedFlags(r) << 1)); // SetObsFlags + flags[i] = static_cast(Unscaled(r) | (ClippedFlags(r) << 1)); // SetObsFlags } }); constexpr F fields[NFLOAT] = {F::I, F::Sigma, F::PrescalingCorr, F::Partiality, @@ -1042,7 +1043,7 @@ void RotationScaleMerge::IngestDevice(const RotationScaleMerge *observations_fro void RotationScaleMerge::BuildInRangeObservations(std::unique_ptr &keys, int n_obs, const SmoothedGeometry &sg) { const int n_run = static_cast(rawrun_start.size()); - const bool limited = d_min_limit.has_value() || d_max_limit.has_value(); + const bool limited = d_min_limit.has_value(); std::vector keep_run(n_run, 1); // Filled once, with the value the branch below actually wants: a resolution-limited run used to @@ -1053,7 +1054,7 @@ void RotationScaleMerge::BuildInRangeObservations(std::unique_ptr &ke 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); + && !(d_min_limit && d < *d_min_limit); keep_run[r] = in_range ? 1 : 0; if (!in_range) continue; const int p0 = rawrun_start[r], p1 = p0 + rawrun_count[r]; @@ -1555,7 +1556,6 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) { const float d = rawrun_d[r]; // resolution is a per-raw-hkl property (all its partials share d) 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; key[r] = keygen(rawrun_h[r], rawrun_k[r], rawrun_l[r]).pack(); eligible[r] = 1; } @@ -3978,11 +3978,11 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int const SphericalHarmonicBasis *sh_basis) { // An observation on an ice ring is excluded from scaling, as the run's log says: its background is // the ring's, uneven round the azimuth, so it would teach a surface the ring rather than the - // detector. The `Is > 0` filter the fit once had hid this omission by dropping the over-subtracted + // detector. So is one coarser than the scale fit's low-resolution limit (Unscaled). The `Is > 0` filter the fit once had hid this omission by dropping the over-subtracted // half of those observations. auto usable = [&](const Obs &o) { return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && o.partiality >= min_partiality - && !o.on_ice; + && !o.unscaled; }; // The observations any of the passes below can use, in fulls order, with the four fields those // passes read copied out beside them. Around twenty passes follow, and each one used to reach back @@ -4713,7 +4713,7 @@ void RotationScaleMerge::Combine() { float peak_frame = partials[ev[i]].r->image_number; float peak_px = partials[ev[i]].r->predicted_x, peak_py = partials[ev[i]].r->predicted_y; float peak_partiality = -1.0f; - const bool on_ice = partials[ev[i]].r->on_ice_ring; + const uint8_t unscaled = Unscaled(*partials[ev[i]].r); uint8_t flags = 0; for (size_t m = i; m < kk; ++m) flags |= ClippedFlags(*partials[ev[m]].r); // A saturated pixel anywhere in the rocking curve means its brightest part was not measured: @@ -4797,7 +4797,7 @@ void RotationScaleMerge::Combine() { full.image_number = peak_frame; full.px = peak_px; full.py = peak_py; full.frame = peak_outcome; - full.on_ice = on_ice ? 1 : 0; + full.unscaled = unscaled; full.clipped = clipped ? 1 : 0; full.group = group; // the raw hkl's ASU group for the current space group (<0 = absent) out.push_back(full); @@ -4933,7 +4933,7 @@ void RotationScaleMerge::ComputePerFrameCC(const std::vector &partial_gr const PartialObs &o = partials[i]; const Reflection &r = *o.r; const int g = partial_group[i]; - if (r.on_ice_ring) continue; + if (Unscaled(r)) continue; if (g < 0) continue; if (o.partiality < min_partiality) continue; if (!std::isfinite(r.I) || !std::isfinite(o.corr) || o.corr <= 0.0f) continue; @@ -5028,7 +5028,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // A full is usable for the merge / error model if it passes AddImage's filters (with the current // ice context). auto usable_merge = [&](const Obs &o) { - return UsableFull(o) && !(for_search && o.on_ice); + return UsableFull(o) && !(for_search && o.unscaled); }; // The em-stats / samples / merge-accumulate / R_meas reductions run on the resident, scaled fulls @@ -5731,7 +5731,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // MergeOnTheFly does on the stills path): the intensity is still merged and written, // the flag only lets a consumer that must not read ring signal as crystal signal - // the CC1/2 resolution-cutoff fit - leave it out. - if (o.on_ice) a.on_ice = true; + if (o.unscaled & UNSCALED_ICE) a.on_ice = true; if (!std::isfinite(a.d) && std::isfinite(o.d) && o.d > 0.0f) a.d = o.d; } } @@ -6551,8 +6551,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // infinite d_max - when no limit was declared. Reflections the beam stop, a mask or the limit ate // then count as missing. The high end stays the finest d reached. See MergeStats in Merge.cpp. const float d_min_pad = sd_min * 0.999f; - const float d_max_declared = d_max_limit ? static_cast(*d_max_limit) - : std::numeric_limits::infinity(); + const float d_max_declared = std::numeric_limits::infinity(); // The innermost shell is labelled by the bound it reaches, except when that bound is infinite // and the coarsest reflection measured is the only finite thing to say. const float d_max_label = std::isfinite(d_max_declared) ? d_max_declared : sd_max; @@ -6798,7 +6797,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool frame_scaled_scratch.assign(n_frames, 1); partial_loop.converged = true; } else if (!scaled_on_gpu) { - const PartialLoopKey key{Group().xhm(), merge_friedel, d_min_limit, d_max_limit, + const PartialLoopKey key{Group().xhm(), merge_friedel, d_min_limit, d_max_scale, min_partiality, smooth_window, scaling_iter}; const auto memo = std::find_if(partial_loop_memos.begin(), partial_loop_memos.end(), [&](const PartialLoopMemo &m) { return m.key == key; }); @@ -6872,7 +6871,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // model needs zeta > 0. auto fitted_obs = [&](int i) { const PartialObs &o = partials[i]; - return partial_group[i] >= 0 && !o.r->on_ice_ring && o.partiality >= min_partiality + return partial_group[i] >= 0 && !Unscaled(*o.r) && o.partiality >= min_partiality && std::isfinite(o.r->zeta) && o.r->zeta > 0.0f; }; const int n_fit_frames = static_cast(frame_start.size()); @@ -7171,7 +7170,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool FullsStaging &st = fulls_staging; st.Resize(nf); gpu_->GetFulls(st.run.data(), st.I.data(), st.sigma.data(), st.d.data(), - st.image_number.data(), st.frame.data(), st.on_ice.data(), st.clipped.data(), + st.image_number.data(), st.frame.data(), st.unscaled.data(), st.clipped.data(), st.group.data()); gpu_->GetFullsPeak(st.peak.data()); const ObsPosition &pos = *obs_position; @@ -7188,7 +7187,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool o.corr = scaled_fulls_on_gpu ? st.corr[i] : 1.0f; o.image_number = st.image_number[i]; o.frame = st.frame[i]; o.px = pos.px[st.peak[i]]; o.py = pos.py[st.peak[i]]; - o.on_ice = st.on_ice[i]; o.clipped = st.clipped[i]; o.group = st.group[i]; + o.unscaled = st.unscaled[i]; o.clipped = st.clipped[i]; o.group = st.group[i]; // A full has no rocking geometry of its own and carries no background; the combine on // the device does not produce these, and a reused buffer would keep the last chain's. o.zeta = 0.0f; o.delta_phi = 0.0f; o.bkg = 0.0f; @@ -7244,7 +7243,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // The partials' filters for the Unity model, which needs no zeta, and whose coefficient is the // mean itself: the partiality is folded into the full, and prescaling_corr is 1. auto fitted_full = [&](const Obs &o) { - return o.group >= 0 && !o.on_ice && o.partiality >= min_partiality; + return o.group >= 0 && !o.unscaled && o.partiality >= min_partiality; }; std::vector fit_start(n_frames + 1, 0); ParallelFor(n_frames, nthreads, [&](int f) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 4ef64c4fb..ced81eccb 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -338,7 +338,7 @@ private: float px, py; // predicted detector position (for the absorption surface; CPU path only) float image_number; // fractional frame position (for 3D-combine contiguity) int32_t frame; // index of the outcome whose per-frame scale G applies to this obs - uint8_t on_ice; + uint8_t unscaled; // UNSCALED_ICE | UNSCALED_LOW_RES: kept out of the scale fit, still merged uint8_t clipped; // OBS_CLIPPED | OBS_OVERLOADED (Reflection::clipped / ::overloaded); a full // is clipped if any of its partials is, and an overloaded event is no full float corr; // image_scale_corr (working; updated by scaling) @@ -362,6 +362,14 @@ private: static uint8_t ClippedFlags(const Reflection &r) { return (r.clipped ? OBS_CLIPPED : 0) | (r.overloaded ? OBS_OVERLOADED : 0); } + // Why an observation stays out of the scale fit (Obs::unscaled), with the bit values the device flags + // use: it sits on an ice ring, whose background is the ring's, or it is coarser than the scale fit's + // low-resolution limit (ScalingSettings::low_resolution_limit_A). Either is still merged and written. + static constexpr uint8_t UNSCALED_ICE = 1; + static constexpr uint8_t UNSCALED_LOW_RES = 8; + uint8_t Unscaled(const Reflection &r) const { + return (r.on_ice_ring ? UNSCALED_ICE : 0) | (d_max_scale && r.d > *d_max_scale ? UNSCALED_LOW_RES : 0); + } // One leverage-corrected error-model sample per usable full: its counting variance at the group's // mean, that mean squared, its squared deviation from the mean - and the resolution it sits at, @@ -410,7 +418,9 @@ private: int n_frames = 0; double min_partiality = 0.02; std::optional d_min_limit; - std::optional d_max_limit; + // The low-resolution limit of the scale FIT. Coarser observations are merged and written, but a + // handful of them on a large cell can bend the whole scale model (Unscaled). + std::optional d_max_scale; bool merge_friedel = true; double capture_uncertainty_coeff = 0.0; double min_captured_fraction = 0.0; @@ -566,12 +576,12 @@ private: struct FullsStaging { FillVector run, peak, frame, group; FillVector I, sigma, d, image_number, corr, var_bkg, var_per_I, capture; - FillVector on_ice, clipped; + FillVector unscaled, clipped; void Resize(int n) { run.resize(n); peak.resize(n); frame.resize(n); group.resize(n); I.resize(n); sigma.resize(n); d.resize(n); image_number.resize(n); corr.resize(n); var_bkg.resize(n); var_per_I.resize(n); capture.resize(n); - on_ice.resize(n); clipped.resize(n); + unscaled.resize(n); clipped.resize(n); } }; FullsStaging fulls_staging; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index e14faab8b..909ceb8bf 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -61,11 +61,14 @@ namespace { struct alignas(16) ObsRecord { float I, sigma, partiality, prescaling_corr; }; // floats 0, 1, 2, 3 // A partial's flags, one byte: on an ice ring, then RotationScaleMerge::OBS_CLIPPED and - // OBS_OVERLOADED one bit up. A full's on_ice byte is 0 or OBS_ON_ICE_FLAG, so the kernels shared by - // the two read either the same way. + // OBS_OVERLOADED one bit up, then coarser than the scale fit's low-resolution limit. The two that keep + // an observation out of the scale fit are OBS_UNSCALED, and a full's unscaled byte holds those of its + // partials, so the kernels shared by the two read either the same way. constexpr uint8_t OBS_ON_ICE_FLAG = 1; constexpr uint8_t OBS_CLIPPED_FLAG = 2; constexpr uint8_t OBS_OVERLOADED_FLAG = 4; + constexpr uint8_t OBS_LOW_RES_FLAG = 8; + constexpr uint8_t OBS_UNSCALED = OBS_ON_ICE_FLAG | OBS_LOW_RES_FLAG; // An observation's frame. The fulls carry theirs (frame); the partials lie in frame order, so theirs // is the last frame starting at or before them in the frame CSR (start, n_frames entries) - an empty @@ -131,10 +134,10 @@ namespace { // where the observation is not one the reference is built from (the host RunArm says why, where it fills ScalingFit). __device__ __forceinline__ bool ScaleCoeff(int i, double min_partiality, const int32_t *__restrict__ group, const ObsRecord &o, const float *__restrict__ zeta, - const uint8_t *__restrict__ on_ice, + const uint8_t *__restrict__ unscaled, const double *__restrict__ group_mean, float &coeff) { const int g = group[i]; - if (!((g >= 0) && !(on_ice[i] & OBS_ON_ICE_FLAG) && isfinite(zeta[i]) && zeta[i] > 0.0f && o.partiality >= min_partiality)) + if (!((g >= 0) && !(unscaled[i] & OBS_UNSCALED) && isfinite(zeta[i]) && zeta[i] > 0.0f && o.partiality >= min_partiality)) return false; const double mean = group_mean[g]; if (!isfinite(mean)) return false; @@ -155,7 +158,7 @@ namespace { const int32_t *__restrict__ frame_count, const int32_t *__restrict__ perm, double min_partiality, const int32_t *__restrict__ group, const ObsRecord *__restrict__ obs, - const float *__restrict__ zeta, const uint8_t *__restrict__ on_ice, + const float *__restrict__ zeta, const uint8_t *__restrict__ unscaled, const double *__restrict__ group_mean, long min_obs, double *__restrict__ g, uint8_t *__restrict__ scaled, double *__restrict__ info) { @@ -173,7 +176,7 @@ namespace { const int a = perm ? perm[i] : i; const ObsRecord o = obs[a]; float coeff_f; - if (!ScaleCoeff(a, min_partiality, group, o, zeta, on_ice, group_mean, coeff_f)) continue; + if (!ScaleCoeff(a, min_partiality, group, o, zeta, unscaled, group_mean, coeff_f)) continue; ++cnt_local; const double coeff = coeff_f; const double w = SafeInvD(o.sigma, 1.0); @@ -214,7 +217,7 @@ namespace { double sx = 0, sy = 0, sx2 = 0, sy2 = 0, sxy = 0; long nl = 0; for (int i = lo + threadIdx.x; i < hi; i += blockDim.x) { - if (flags[i] & OBS_ON_ICE_FLAG) continue; + if (flags[i] & OBS_UNSCALED) continue; const int g = group[i]; if (g < 0) continue; const ObsRecord o = obs[i]; @@ -329,7 +332,7 @@ namespace { const int32_t *rr_offset; // emit pass: per-run base offset into the fulls arrays int32_t *f_run, *f_peak, *f_frame, *f_group; float *f_I, *f_sigma, *f_d, *f_img, *f_var_bkg, *f_var_per_I, *f_capture; - uint8_t *f_on_ice, *f_clipped; + uint8_t *f_unscaled, *f_clipped; }; // One thread's raw-hkl run: split its usable partials (already in image_number order within the run) @@ -381,7 +384,7 @@ namespace { return double(p.obs[i].I) + n_bkg * (double(p.bkg[i]) - pooled_bkg); }; - // Pass B: seed F (inverse-variance mean of pooled_I*corr), plus peak / d / on_ice / partiality. + // Pass B: seed F (inverse-variance mean of pooled_I*corr), plus peak / d / unscaled / partiality. double sum_w = 0.0, sum_wI = 0.0, sum_partiality = 0.0; float d = NAN; const int first = p.perm[ev_start]; @@ -389,7 +392,7 @@ namespace { float peak_frame = p.image_number[first]; int peak_obs = first; float peak_partiality = -1.0f; - const bool on_ice = p.flags[first] & OBS_ON_ICE_FLAG; + const uint8_t unscaled = p.flags[first] & OBS_UNSCALED; uint8_t flags = 0; for (int m = ev_start; m <= ev_end; ++m) flags |= p.flags[p.perm[m]]; // An overloaded event is not a measurement of the reflection; see the host combine. @@ -466,7 +469,7 @@ namespace { p.f_img[o] = peak_frame; p.f_peak[o] = peak_obs; p.f_frame[o] = peak_outcome; - p.f_on_ice[o] = on_ice ? 1 : 0; + p.f_unscaled[o] = unscaled; p.f_clipped[o] = clipped ? 1 : 0; p.f_group[o] = group; } @@ -586,13 +589,13 @@ namespace { // Per full: 1 if it gives a term (all, even frame, odd frame) - the host's `usable(o) && cell >= 0`. // A full's partiality is 1 (the combine sets it so on both paths). - __global__ void SurfaceKeepKernel(int n, const int32_t *group, const float *corr, const uint8_t *on_ice, + __global__ void SurfaceKeepKernel(int n, const int32_t *group, const float *corr, const uint8_t *unscaled, const int32_t *frame, const int32_t *cell, double min_partiality, int32_t *keep, int32_t *keep_even, int32_t *keep_odd) { const int i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= n) return; const float c = corr[i]; - const bool k = group[i] >= 0 && c > 0.0f && isfinite(c) && 1.0f >= min_partiality && !on_ice[i] + const bool k = group[i] >= 0 && c > 0.0f && isfinite(c) && 1.0f >= min_partiality && !unscaled[i] && cell[i] >= 0; const bool odd = frame[i] & 1; keep[i] = k; keep_even[i] = k && !odd; keep_odd[i] = k && odd; @@ -699,7 +702,7 @@ namespace { const float *I, *sigma, *corr, *d, *reject_median, *var_bkg, *var_per_I, *capture; const float *reject_var_add; // per group: the shell's measured share of dI^2/4 const int32_t *group, *frame; - const uint8_t *on_ice, *frame_cell_ok, *half; + const uint8_t *unscaled, *frame_cell_ok, *half; const uint8_t *hand, *has_hands; // per-full Bijvoet hand; per-group "two hands to split" const int32_t *gperm, *gstart, *gcount; const double *em_mean, *merged_I, *cc_factor; @@ -727,7 +730,7 @@ namespace { if (!p.frame_cell_ok[p.frame[i]]) return false; const float c = p.corr[i]; if (!(c > 0.0f) || !isfinite(c)) return false; - if (p.for_search && p.on_ice[i]) return false; + if (p.for_search && p.unscaled[i]) return false; if (FULL_PARTIALITY < p.min_partiality) return false; const float I_corr = p.I[i] * c, sigma_corr = p.sigma[i] * c; return isfinite(I_corr) && isfinite(sigma_corr) && sigma_corr > 0.0f; @@ -982,7 +985,7 @@ namespace { swI += wI; sw += w; if (half) { swIh1 += wI; swh1 += w; swht1 += wt; ++nh1; } else { swIh0 += wI; swh0 += w; swht0 += wt; ++nh0; } - if (p.on_ice[i]) onice = 1; // see the host merge loop + if (p.unscaled[i] & OBS_ON_ICE_FLAG) onice = 1; // see the host merge loop if (!isfinite(dd) && isfinite(p.d[i]) && p.d[i] > 0.0f) dd = p.d[i]; } p.a_swI[g] = swI; p.a_sw[g] = sw; p.a_swIh0[g] = swIh0; p.a_swIh1[g] = swIh1; @@ -1295,7 +1298,7 @@ struct RotationScaleMergeGPU::Impl { // run: the raw-hkl run a full is of; peak: the partial its position is taken from CudaDevicePtr f_run, f_peak, f_frame, f_group; CudaDevicePtr f_I, f_sigma, f_d, f_img, f_var_bkg, f_var_per_I, f_capture; - CudaDevicePtr f_on_ice, f_clipped; + CudaDevicePtr f_unscaled, f_clipped; // scale-fulls (Unity model, kept resident): all-ones partiality/zeta and records with partiality and // prescaling_corr 1, so the shared scaling kernels yield coeff=mean, plus the working corr and the // fulls frame/group CSRs (built on the host from the small f_frame/f_group key arrays, over the @@ -1426,7 +1429,7 @@ void RotationScaleMergeGPU::Impl::ReleaseFulls() { n_fulls = 0; f_run = {}; f_peak = {}; f_frame = {}; f_group = {}; f_I = {}; f_sigma = {}; f_d = {}; f_img = {}; - f_var_bkg = {}; f_var_per_I = {}; f_capture = {}; f_on_ice = {}; f_clipped = {}; + f_var_bkg = {}; f_var_per_I = {}; f_capture = {}; f_unscaled = {}; f_clipped = {}; f_corr = {}; f_zeta = {}; f_obs = {}; f_frame_perm = {}; f_frame_start = {}; f_frame_count = {}; f_gperm = {}; f_gstart = {}; f_gcount = {}; @@ -1647,7 +1650,7 @@ void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partialit p.for_search = d.merge_for_search; p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); - p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); + p.unscaled = d.f_unscaled.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.em_mean = d.m_em_mean.get(); p.sw = sw.get(); p.swI = swI.get(); p.em_mean_out = d.m_em_mean.get(); p.cnt = cnt.get(); @@ -1730,7 +1733,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ p.reject_var_add = var_add.get(); p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); - p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); + p.unscaled = d.f_unscaled.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.half = d_half.get(); p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); p.capture = d.f_capture.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); @@ -1825,7 +1828,7 @@ void RotationScaleMergeGPU::MergeBijvoetExcess(double error_model_a, double erro p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); - p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); + p.unscaled = d.f_unscaled.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.hand = d.m_hand.get(); p.has_hands = d.m_has_hands.get(); p.error_model_a = error_model_a; p.error_model_b = error_model_b; p.error_model_active = 1; @@ -1857,7 +1860,7 @@ void RotationScaleMergeGPU::MergeGroupScatter(double *sum, double *sum_sq, doubl p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); - p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); + p.unscaled = d.f_unscaled.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.em_mean = d.m_em_mean.get(); p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); p.capture = d.f_capture.get(); @@ -1884,7 +1887,7 @@ void RotationScaleMergeGPU::MergeHalves(uint64_t *obs_key, uint8_t *half) { p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); - p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); + p.unscaled = d.f_unscaled.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); const int obs_blocks = std::min(65535, (nf + BLK - 1) / BLK); @@ -2041,7 +2044,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti d.f_I = d.Alloc(nf); d.f_sigma = d.Alloc(nf); d.f_d = d.Alloc(nf); d.f_img = d.Alloc(nf); d.f_var_bkg = d.Alloc(nf); d.f_var_per_I = d.Alloc(nf); d.f_capture = d.Alloc(nf); - d.f_on_ice = d.Alloc(nf); d.f_clipped = d.Alloc(nf); + d.f_unscaled = d.Alloc(nf); d.f_clipped = d.Alloc(nf); d.f_corr = d.Alloc(nf); d.f_zeta = d.Alloc(nf); d.f_obs = d.Alloc(nf); @@ -2052,7 +2055,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.f_frame = d.f_frame.get(); p.f_group = d.f_group.get(); p.f_I = d.f_I.get(); p.f_sigma = d.f_sigma.get(); p.f_d = d.f_d.get(); p.f_img = d.f_img.get(); p.f_var_bkg = d.f_var_bkg.get(); p.f_var_per_I = d.f_var_per_I.get(); p.f_capture = d.f_capture.get(); - p.f_on_ice = d.f_on_ice.get(); p.f_clipped = d.f_clipped.get(); + p.f_unscaled = d.f_unscaled.get(); p.f_clipped = d.f_clipped.get(); if (d.n_fulls > 0) { CombineKernel<<s()>>>(p); CudaCheck(cudaGetLastError(), "combine emit launch"); @@ -2062,7 +2065,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti } void RotationScaleMergeGPU::GetFulls(int32_t *run, float *I, float *sigma, float *d, - float *image_number, int32_t *frame, uint8_t *on_ice, + float *image_number, int32_t *frame, uint8_t *unscaled, uint8_t *clipped, int32_t *group) const { DeviceGuard guard(impl_->device, impl_->available); const auto &dd = *impl_; @@ -2075,7 +2078,7 @@ void RotationScaleMergeGPU::GetFulls(int32_t *run, float *I, float *sigma, float dl(group, dd.f_group.get(), n * sizeof(int32_t)); dl(I, dd.f_I.get(), n * sizeof(float)); dl(sigma, dd.f_sigma.get(), n * sizeof(float)); dl(d, dd.f_d.get(), n * sizeof(float)); dl(image_number, dd.f_img.get(), n * sizeof(float)); - dl(on_ice, dd.f_on_ice.get(), n * sizeof(uint8_t)); + dl(unscaled, dd.f_unscaled.get(), n * sizeof(uint8_t)); dl(clipped, dd.f_clipped.get(), n * sizeof(uint8_t)); } @@ -2111,7 +2114,7 @@ void RotationScaleMergeGPU::FitFullsScale(double min_partiality) { d.f_obs.get(), d.f_corr.get(), d.group_mean.get()); CudaCheck(cudaGetLastError(), "ReduceGroupMeansKernel launch"); FitPerFrameGKernel<<s()>>>(d.n_frames, d.f_frame_start.get(), d.f_frame_count.get(), - d.f_frame_perm.get(), min_partiality, d.f_group.get(), d.f_obs.get(), d.f_zeta.get(), d.f_on_ice.get(), + d.f_frame_perm.get(), min_partiality, d.f_group.get(), d.f_obs.get(), d.f_zeta.get(), d.f_unscaled.get(), d.group_mean.get(), 1L, d.g.get(), d.scaled.get(), d.info.get()); CudaCheck(cudaGetLastError(), "scale fulls launch"); CudaCheck(cudaStreamSynchronize(impl_->s()), "scale fulls sync"); @@ -2250,7 +2253,7 @@ int RotationScaleMergeGPU::SurfaceBuildTerms(const int32_t *cell, const float *c keep = d.Alloc(std::max(1, nf)); even = d.Alloc(std::max(1, nf)); odd = d.Alloc(std::max(1, nf)); if (nf > 0) { - SurfaceKeepKernel<<<(nf + BLK - 1) / BLK, BLK, 0, st>>>(nf, d.f_group.get(), d_corr.get(), d.f_on_ice.get(), + SurfaceKeepKernel<<<(nf + BLK - 1) / BLK, BLK, 0, st>>>(nf, d.f_group.get(), d_corr.get(), d.f_unscaled.get(), d.f_frame.get(), d_cell.get(), min_partiality, keep.get(), even.get(), odd.get()); CudaCheck(cudaGetLastError(), "surface keep launch"); } @@ -2511,6 +2514,7 @@ namespace { ObsRecord *obs; float *zeta, *corr, *bkg, *var_bkg, *image_number, *d, *dphi; uint8_t *flags; + float d_max_scale; // coarser than this is kept out of the scale fit (OBS_LOW_RES_FLAG) const int32_t *new_idx; // the cut's place of each reflection (an exclusive prefix over the kept) int n, n_keep; // reflections and observations }; @@ -2538,7 +2542,7 @@ namespace { a.d[i] = r.d; a.dphi[i] = r.delta_phi_deg; a.flags[i] = (r.on_ice_ring ? OBS_ON_ICE_FLAG : 0) | (r.clipped ? OBS_CLIPPED_FLAG : 0) - | (r.overloaded ? OBS_OVERLOADED_FLAG : 0); + | (r.overloaded ? OBS_OVERLOADED_FLAG : 0) | (r.d > a.d_max_scale ? OBS_LOW_RES_FLAG : 0); } @@ -3117,14 +3121,14 @@ void RotationScaleMergeGPU::IngestGetLayout(IngestLayout &out) const { } void RotationScaleMergeGPU::IngestObservations(const Reflection *const *refl, const int32_t *count, size_t nthreads, - const double *frame_flux, float *px, float *py) { + const double *frame_flux, float d_max_scale, float *px, float *py) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; const int n_keep = d.in_n_keep; d.AllocObservations(n_keep); d.in_dphi = d.Alloc(std::max(1, n_keep)); const IngestArrays a{d.obs.get(), d.zeta.get(), d.corr.get(), d.bkg.get(), d.var_bkg.get(), d.image_number.get(), - d.d_obs.get(), d.in_dphi.get(), d.flags.get(), d.in_new_idx.get(), + d.d_obs.get(), d.in_dphi.get(), d.flags.get(), d_max_scale, d.in_new_idx.get(), d.in_src_start[d.n_frames], n_keep}; d.StageReflections(refl, count, nthreads, [&](int cnt, int base, const Reflection *src) { IngestExtractKernel<<<(cnt + BLK - 1) / BLK, BLK, 0, d.s()>>>(cnt, base, src, a); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index b479fe174..d28fe3476 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -88,8 +88,9 @@ public: // its order, each with its deterministic correction divided by its frame's flux (length n_frames). // px/py (length IngestSort's count) receive each observation's predicted position, which is not // kept on the card (GetFullsPeak). + // d_max_scale: observations coarser than this are flagged to stay out of the scale fit (INFINITY: none). void IngestObservations(const Reflection *const *refl, const int32_t *count, size_t nthreads, - const double *frame_flux, float *px, float *py); + const double *frame_flux, float d_max_scale, float *px, float *py); // The ingest's geometry smoothing, on the observations IngestSort left (RotationScaleMerge::IngestOnDevice): // each delta_phi re-solved against its frame's smoothed lattice where the frame has one - recip holds // a*, b*, c* per frame (9 floats), ok which frames, and the goniometer frame as SmoothedGeometry @@ -240,7 +241,7 @@ public: // of. The working corr is downloaded separately by GetFullsCorr (it is only meaningful after the // fulls scaling; otherwise the caller sets it). void GetFulls(int32_t *run, float *I, float *sigma, float *d, - float *image_number, int32_t *frame, uint8_t *on_ice, uint8_t *clipped, + float *image_number, int32_t *frame, uint8_t *unscaled, uint8_t *clipped, int32_t *group) const; // Download the observation each full's predicted detector position is that of (its peak partial), diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index be30055f0..2a3bb54b5 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -222,7 +222,7 @@ void print_usage() { std::cout << " --no-expected-variance-merge stills: disable the default expected-variance merge weighting (which rebuilds each weak observation's signal variance at the reflection mean to de-bias the inverse-variance merge); restores observed-sigma weighting" << std::endl; std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; std::cout << " --scaling-high-resolution High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl; - std::cout << " --scaling-low-resolution Low resolution limit for scaling/merging, in A (default: 50, also XDS's own default; 0 = no limit). The beam stop suppresses air scatter well beyond the shadow the beam-stop mask covers, so reflections coarser than this are measured on a depressed background and mostly integrate negative; masking does not reach them" << std::endl; + std::cout << " --scaling-low-resolution Low resolution limit of the scale fit, in A (default: 50, also XDS's own default; 0 = no limit). Coarser reflections are merged and written, but kept out of the scale model, which a few of them on a large cell can bend" << std::endl; std::cout << " --resolution-cutoff Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl; std::cout << " --resolution-cc-target CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl; std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 9, as XDS)" << std::endl; diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index c2c930911..d8fffe5f0 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -162,19 +162,25 @@ namespace { return out; } + // The statistics over a declared low-resolution bound (--report-resolution), or over the whole sphere + // without one: the scaling's own low-resolution limit bounds the scale fit, not the merge. MergeStatistics StatsWithLowLimit(const gemmi::SpaceGroup &sg, const std::optional &cell, const std::vector &merged, std::optional low_limit, bool merge_friedel) { DiffractionExperiment x; x.SetSpaceGroup(sg); ScalingSettings s = x.GetScalingSettings(); - s.LowResolutionLimit_A(low_limit); s.MergeFriedel(merge_friedel); x.ImportScalingSettings(s); MergeOnTheFly merge(x); merge.ReferenceCell(cell); - return merge.MergeStats(merged, {}); + const auto own = merge.MergeStats(merged, {}); + if (!low_limit) + return own; + // Down to the data's own d_min, padded as MergeStats pads it. + return merge.MergeStats(merged, {}, {}, std::nullopt, + ReportResolutionRange{.d_min = own.overall.d_min * 0.999, .d_max = *low_limit}); } double Completeness(const MergeStatisticsShell &s) { diff --git a/tests/RotationIngestGPUTest.cpp b/tests/RotationIngestGPUTest.cpp index 7b8588480..f0d4f973a 100644 --- a/tests/RotationIngestGPUTest.cpp +++ b/tests/RotationIngestGPUTest.cpp @@ -161,7 +161,7 @@ TEST_CASE("RotationIngestGPU_SortAndCutMatchHost", "[RotationScale][gpu]") { } std::vector px(n_keep), py(n_keep); - gpu.IngestObservations(ptr.data(), count.data(), 4, flux.data(), px.data(), py.data()); + gpu.IngestObservations(ptr.data(), count.data(), 4, flux.data(), INFINITY, px.data(), py.data()); for (int i = 0; i < n; ++i) if (keep[i]) { CHECK(px[new_idx[i]] == all[i]->predicted_x);