Files
Jungfraujoch/common/ScalingSettings.cpp
leonarski_fandClaude Opus 5.5 a84fe01712 rugnux: --sweeps joint merges several sweeps of one crystal as one dataset
Two levels, as XSCALE merges the data sets XDS made of each sweep. Each sweep is processed on its own
first (a member run: everything a run does, writing only <prefix>_s<k>_log.txt, without the model),
and its written merge hands back all of its fulls - past its resolution cut, with the reflections its
group extinguishes, and with the counting variance the merge rebuilds them from (ScaledFull var_bkg,
var_per_I, capture, d). RelateJointSweeps puts every sweep on the axes of the one with the most fulls
by ChooseReferenceSetting (axis permutations and the lattice's twin laws, scored by the merged-intensity
CC); a sweep below MIN_REFERENCE_CC on every axes is left out and named. BuildJointMergeInput lays the
sweeps' frames end to end, each full a pre-combined observation whose count, correction and variance
reproduce what its own merge gave it, and Rugnux::RunJoint runs the one pass of the pipeline on them:
no image is read, the outcomes are the indexer's, and the unchanged ScaleMergeAndSymmetry makes the
merge (RotationScaleMerge with the sweep segments, now taken from ScalingSettings::SweepSegments), the
space-group search, the cut and everything after it, into <prefix>.*. The run's log is opened first and
holds the sweeps' relation and the joint merge. Not yet: the unmerged MTZ of a joint merge, the per-sweep
report keys, and joint geometry (P1d / phase 2).

Single-sweep runs are unchanged (RunJoint's branches in RunPipeline are only taken by it; the exports
are only asked for by member runs).

Test pairs (rugnux joint against the best sweep alone, and against XSCALE, at the XDS range):
- in-house lysozyme, two 450 x 0.2 deg sweeps: P41212 (both sweeps alone too); to 1.85 A fitted against
  1.91-1.94; at 50-1.97 A R_meas 14.5% against 14.9%, I/sigma 14.1 against 9.4, ISa 24.4 against 20.8;
  XSCALE of the two XDS sweeps: R_meas 17.4%, I/sigma 9.7, ISa 16.0 / 15.2.
- 8C3E, four short CCD sweeps of a twinned P3121 crystal: the joint merge determines P3121 (the
  deposition), where the sweeps alone give P321, P6222, P321 and P622; completeness 98% against 96%.
- 5JVN, two 3600-frame sweeps: P622; to 2.32 A against 2.35; at 50-2.62 A R_meas 16.6% against 13.7%,
  I/sigma 19.7 against 17.0, ISa 11.2 against 15.3; the depositors' XSCALE: R_meas 15.5%, ISa 15.3 / 13.5.

Tests: RugnuxJointSweepsTest (a sweep indexed by the twin law put back, another crystal left out, the
segments and the pre-combined observation's count and variance).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
2026-10-28 20:17:36 +01:00

327 lines
9.8 KiB
C++

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "ScalingSettings.h"
ScalingSettings& ScalingSettings::MergeFriedel(bool input) {
merge_friedel = input;
return *this;
}
ScalingSettings& ScalingSettings::HighResolutionLimit_A(double limit) {
if (limit <= 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "High resolution limit must be positive");
high_resolution_limit_A = limit;
return *this;
}
ScalingSettings& ScalingSettings::HighResolutionLimit_A(std::optional<double> limit) {
if (limit.has_value() && limit.value() <= 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "High resolution limit must be positive");
high_resolution_limit_A = limit;
return *this;
}
ScalingSettings& ScalingSettings::LowResolutionLimit_A(std::optional<double> limit) {
if (limit.has_value() && limit.value() <= 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "Low resolution limit must be positive");
low_resolution_limit_A = limit;
return *this;
}
bool ScalingSettings::GetMergeFriedel() const {
return merge_friedel;
}
ScalingSettings &ScalingSettings::RefineRotationWedge(bool input) {
refine_wedge = input;
return *this;
}
bool ScalingSettings::GetRefineWedge() const {
return refine_wedge;
}
std::optional<double> ScalingSettings::GetHighResolutionLimit_A() const {
return high_resolution_limit_A;
}
std::optional<double> ScalingSettings::GetLowResolutionLimit_A() const {
return low_resolution_limit_A;
}
double ScalingSettings::GetMinMosaicity() const {
return 0.001;
}
double ScalingSettings::GetMaxMosaicity() const {
return 1.0;
}
double ScalingSettings::GetMinWedge() const {
return 0.001;
}
double ScalingSettings::GetMaxWedge() const {
return 10.0;
}
double ScalingSettings::GetDefaultMosaicity() const {
return 0.1;
}
ScalingSettings &ScalingSettings::RotationWedgeForScaling(std::optional<double> input) {
if (input) {
// TODO: Use fmt
if (input.value() < GetMinWedge() || input.value() > GetMaxWedge())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Wedge for scaling must be between " + std::to_string(GetMinWedge()) +
" and " + std::to_string(GetMaxWedge()));
}
wedge_for_scaling = input;
return *this;
}
std::optional<double> ScalingSettings::GetRotationWedgeForScaling() const {
return wedge_for_scaling;
}
ScalingSettings &ScalingSettings::MinPartiality(double input) {
if (!(input >= 0.0 && input <= 1.0))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Min partiality must be between 0 and 1");
min_partiality = input;
return *this;
}
double ScalingSettings::GetMinCCForImage() const {
return min_cc_for_image;
}
double ScalingSettings::GetSearchMinZeta() const {
return search_min_zeta;
}
ScalingSettings &ScalingSettings::SearchMinZeta(double input) {
if (input < 0.0 || input >= 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Search zeta limit must be in [0,1)");
search_min_zeta = input;
return *this;
}
ScalingSettings &ScalingSettings::MinCCForImage(double input) {
if (input < 0.0 || input > 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Min CC for image must be between 0 and 1");
min_cc_for_image = input;
return *this;
}
double ScalingSettings::GetOutlierRejectNsigma() const {
return outlier_reject_nsigma;
}
ScalingSettings &ScalingSettings::OutlierRejectNsigma(double input) {
outlier_reject_nsigma = input; // <= 0 disables; no upper bound (large = effectively off)
return *this;
}
ScalingSettings &ScalingSettings::ScaleFulls(bool input) {
scale_fulls = input;
return *this;
}
bool ScalingSettings::GetScaleFulls() const {
return scale_fulls;
}
ScalingSettings &ScalingSettings::AbsorptionIter(int input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Absorption iterations must be non-negative");
absorption_iter = input;
return *this;
}
int ScalingSettings::GetAbsorptionIter() const {
return absorption_iter;
}
ScalingSettings &ScalingSettings::CorrectionSurfaces(bool input) {
correction_surfaces = input;
return *this;
}
bool ScalingSettings::GetCorrectionSurfaces() const {
return correction_surfaces;
}
ScalingSettings &ScalingSettings::StillsPartialityRefine(bool input) {
stills_partiality_refine = input;
return *this;
}
bool ScalingSettings::GetStillsPartialityRefine() const {
return stills_partiality_refine;
}
ScalingSettings &ScalingSettings::ExpectedVarianceMerge(bool input) {
expected_variance_merge = input;
return *this;
}
bool ScalingSettings::GetExpectedVarianceMerge() const {
return expected_variance_merge;
}
float ScalingSettings::GetIceMinScore() const {
return ice_min_score;
}
float ScalingSettings::GetIceMinSpotRatio() const {
return ice_min_spot_ratio;
}
ScalingSettings &ScalingSettings::IceMinSpotRatio(float input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Ice spot-ratio gate must be non-negative");
ice_min_spot_ratio = input;
return *this;
}
ScalingSettings &ScalingSettings::IceMinScore(float input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Ice score gate must be non-negative");
ice_min_score = input;
return *this;
}
ScalingSettings &ScalingSettings::SmoothGDegrees(double input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Smooth-G range must be non-negative");
smooth_g_deg = input;
return *this;
}
double ScalingSettings::GetSmoothGDegrees() const {
return smooth_g_deg;
}
ScalingSettings &ScalingSettings::RelativeBDegrees(double input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Relative-B batch width must be non-negative");
relative_b_deg = input;
return *this;
}
double ScalingSettings::GetRelativeBDegrees() const {
return relative_b_deg;
}
double ScalingSettings::GetMinPartiality() const {
return min_partiality;
}
ScalingSettings &ScalingSettings::ForcedMosaicity(std::optional<double> input) {
if (input.has_value() && (input.value() < GetMinMosaicity() || input.value() > GetMaxMosaicity()))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Forced mosaicity must be between " + std::to_string(GetMinMosaicity()) +
" and " + std::to_string(GetMaxMosaicity()));
forced_mosaicity = input;
return *this;
}
std::optional<double> ScalingSettings::GetForcedMosaicity() const {
return forced_mosaicity;
}
ScalingSettings &ScalingSettings::CaptureUncertaintyCoeff(double input) {
if (input < 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Capture uncertainty coefficient must be non-negative");
capture_uncertainty_coeff = input;
return *this;
}
double ScalingSettings::GetCaptureUncertaintyCoeff() const {
return capture_uncertainty_coeff;
}
ScalingSettings &ScalingSettings::MinCapturedFraction(double input) {
if (input < 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Minimum captured fraction must be non-negative");
min_captured_fraction = input;
return *this;
}
double ScalingSettings::GetMinCapturedFraction() const {
return min_captured_fraction;
}
ScalingSettings &ScalingSettings::RfreeFraction(double input) {
if (input < 0.0 || input > 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "R-free fraction must be between 0 and 1");
rfree_fraction = input;
return *this;
}
double ScalingSettings::GetRfreeFraction() const {
return rfree_fraction;
}
ScalingSettings &ScalingSettings::ResolutionCutoff(ResolutionCutoffMethod input) {
resolution_cutoff = input;
return *this;
}
ResolutionCutoffMethod ScalingSettings::GetResolutionCutoff() const {
return resolution_cutoff;
}
ScalingSettings &ScalingSettings::ResolutionCCTarget(double input) {
if (input <= 0.0 || input >= 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Resolution CC target must be between 0 and 1");
resolution_cc_target = input;
return *this;
}
double ScalingSettings::GetResolutionCCTarget() const {
return resolution_cc_target;
}
ScalingSettings &ScalingSettings::ReportShellCount(int input) {
if (input < 1)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Number of report shells must be at least 1");
report_shell_count = input;
return *this;
}
int ScalingSettings::GetReportShellCount() const {
return report_shell_count;
}
ScalingSettings& ScalingSettings::ReportResolutionRange_A(std::optional<ReportResolutionRange> input) {
if (input && !(input->d_min > 0.0 && input->d_max > input->d_min))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Report resolution range must be positive and d_max coarser than d_min");
report_resolution_range = input;
return *this;
}
std::optional<ReportResolutionRange> ScalingSettings::GetReportResolutionRange_A() const {
return report_resolution_range;
}
ScalingSettings& ScalingSettings::SweepSegments(std::vector<SweepSegment> input) {
sweep_segments = std::move(input);
return *this;
}
const std::vector<SweepSegment> &ScalingSettings::GetSweepSegments() const {
return sweep_segments;
}