mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 64 Wiki Activity

Add global (joint) scaling path: jfjoch_process --global-scale

Browse Source 
Refine all per-image scales, per-image mosaicities and the shared per-HKL
true intensities (plus the global wedge) together in one Ceres problem
(GlobalScale.{h,cpp}, GlobalScaleCeres), as an alternative to the existing
alternating merge-then-scale-each-image loop. The result feeds back into each
reflection's image_scale_corr exactly as ScaleOnTheFly does, so the rot3d
combine and MergeOnTheFly run unchanged and every metric stays comparable.

On HEWL crystal 2 the joint fit converges to the same place as the alternating
loop (anomalous S-peak vs XDS 0.53x -> 0.54x, ISa 9.4 both), confirming the
alternating path already reaches the joint optimum. Kept as a validated
alternative and the home for future global corrections.

A shared quadratic absorption surface in detector position was prototyped here
and dropped: it fit large non-physical coefficients (radially degenerate with
the per-HKL resolution structure), lowered the scaling-fit residual but raised
the error-model b (ISa 9.4 -> 7.3) and did not improve anomalous accuracy.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-25 20:43:04 +02:00
co-authored by Claude Opus 4.8
parent 2041ce2bd5
commit 50ebf1694e
6 changed files with 632 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
process/JFJochProcess.cpp
+74 -5
View File
@@ -27,6 +27,7 @@
#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/GlobalScale.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/scale_merge/Combine3D.h"
#include "../image_analysis/WriteReflections.h"
@@ -54,6 +55,69 @@ namespace {
ret.assign(unique_ordinals.begin(), unique_ordinals.end());
return ret;
}
// Global (joint) scaling: refine all per-image scales/mosaicities and the shared per-HKL
// intensities in one Ceres problem (GlobalScaleCeres), then write the result back into each
// reflection's image_scale_corr exactly as ScaleOnTheFly does - recompute the rotation
// partiality from the refined per-image mosaicity, then image_scale_corr = rlp/(partiality*G).
// The downstream rot3d combine and MergeOnTheFly are unchanged, so every metric stays
// comparable with the alternating path; only how G and mosaicity are found differs.
void RunGlobalScaling(const DiffractionExperiment &experiment,
std::vector<IntegrationOutcome> &outcomes, Logger &logger) {
const auto &s = experiment.GetScalingSettings();
const bool rotation = experiment.GetPartialityModel() == PartialityModel::Rotation;
std::vector<std::vector<Reflection>> observations;
observations.reserve(outcomes.size());
for (const auto &o: outcomes)
observations.push_back(o.reflections);
GlobalScaleOptions opt;
opt.merge_friedel = s.GetMergeFriedel();
if (const auto sg = experiment.GetGemmiSpaceGroup())
opt.space_group = *sg;
opt.d_min_limit_A = s.GetHighResolutionLimit_A().value_or(0.0);
opt.mosaicity_init_deg = s.GetDefaultMosaicity();
// The joint problem is large (per-image G + mosaicity + all Itrue + wedge); give the solver
// enough room to converge rather than the per-image default.
opt.max_solver_time_s = 300.0;
opt.max_num_iterations = 50;
double wedge = 0.0;
if (rotation) {
opt.partiality_model = GlobalScaleOptions::PartialityModel::Rotation;
wedge = experiment.GetRotationWedgeForScaling().value_or(0.0);
opt.wedge_deg = wedge;
opt.mosaicity_min_deg = s.GetMinMosaicity();
opt.mosaicity_max_deg = s.GetMaxMosaicity();
} else {
opt.partiality_model = GlobalScaleOptions::PartialityModel::Fixed;
}
const auto result = GlobalScaleCeres(observations, opt);
const double half_wedge = wedge / 2.0;
for (size_t i = 0; i < outcomes.size(); ++i) {
const double G = result.image_scale_g[i];
const double mos = result.mosaicity_deg[i];
for (auto &r: outcomes[i].reflections) {
if (rotation && std::isfinite(r.delta_phi_deg) && std::isfinite(r.zeta) && mos > 1e-6) {
const double c1 = r.zeta / std::sqrt(2.0);
r.partiality = static_cast<float>(
(std::erf((r.delta_phi_deg + half_wedge) * c1 / mos)
- std::erf((r.delta_phi_deg - half_wedge) * c1 / mos)) / 2.0);
}
const double denom = r.partiality * G;
r.image_scale_corr = (std::isfinite(r.rlp) && std::isfinite(denom) && denom > 0.0)
? static_cast<float>(r.rlp / denom) : NAN;
}
if (std::isfinite(G))
outcomes[i].image_scale_g = static_cast<float>(G);
if (rotation && std::isfinite(mos))
outcomes[i].mosaicity_deg = static_cast<float>(mos);
}
logger.Info("Global scaling complete ({} images)", outcomes.size());
}
}
JFJochProcess::JFJochProcess(JFJochHDF5Reader &reader, DiffractionExperiment experiment,
@@ -380,11 +444,16 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) {
// ScaleOnTheFly is only for the classical, no-reference path; with a reference (or
// PixelRefine) each image is already scaled, so we merge directly.
if (config_.reference_data.empty() && !pixel_refine_path) {
logger.Info("Running scaling ...");
ScalingResult scale_result(0);
for (int i = 0; i < config_.scaling_iter; i++) {
auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false);
scale_result = indexer->ScaleAllImages(merge_result);
if (config_.global_scaling) {
logger.Info("Running global (joint) scaling ...");
RunGlobalScaling(experiment_, indexer->GetIntegrationOutcome(), logger);
} else {
logger.Info("Running scaling ...");
ScalingResult scale_result(0);
for (int i = 0; i < config_.scaling_iter; i++) {
auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false);
scale_result = indexer->ScaleAllImages(merge_result);
}
}
}