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

rugnux --model: score the twin hypothesis too, and take the maps from the merge the model prefers

Browse Source 
Where a twin hypothesis (<prefix>_twin_<SG>) was written, the model is also
validated against that merge: as merged, and as twinned by its law -
|F_twin(h)|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 from the placed, scaled
model, one overall scale refitted per a, a scanned 0..0.5 on the working set
and R-free read at it (ValidateAgainstModel's new twin_law argument; GPU
structure factors as before). The maps and the placed model are then written
from whichever merge has the lower R-free; the report gives MAPS_FROM and the
subgroup's plain and twinned R-work/R-free and fraction. The subgroup file
follows its own validation's indexing, as the answer follows its own.

The model chooses the maps only: the space group, <prefix>.mtz and R_FREE stay
the answer's. Documented in RUGNUX_REPORT.md, RUGNUX_ADVANCED.md and
CPU_DATA_ANALYSIS_DECISIONS.md (14.9); Yeates (1997) credited.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
This commit is contained in:
leonarski_f
2026-10-09 19:40:22 +02:00
co-authored by Claude Opus 5.5
parent 475df9e5cb
commit 8442d55f6f
11 changed files with 189 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/structure_refinement/ModelValidation.cpp
+61 -3
View File
@@ -423,6 +423,7 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
double wavelength_A,
const std::vector<float> &report_shell_d_min,
const ModelValidationSchedule &schedule,
const gemmi::Op *twin_law,
bool rigid_body_gpu) {
ModelValidationResult result;
result.model_path = model_path;
@@ -1228,6 +1229,61 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
result.b_sol = best.b_sol;
result.k_overall = best.k_overall;
// --- the twin target, where a twin law was given ---
// The intensity a twinned crystal records at h is (1-a) I(h) + a I(Th), so the model is compared
// against the data as twinned: |F_twin(h)|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 from the
// placed, scaled model, with one overall scale refitted at each a. The fraction is scanned in steps
// of 0.01 and chosen on the working set; R-free is read at it. Nothing else is refitted against the
// twinned target - the placement, the anisotropic B and the bulk solvent are those of the untwinned
// fit - so this is the same few-parameter comparison as R above, not a twin refinement.
// Following Yeates (1997) Methods Enzymol. 276, 344-358
if (twin_law != nullptr) {
std::unordered_map<long, double> model_i;
model_i.reserve(fmodel.v.size());
for (const auto &hv : fmodel.v)
model_i[hkl_key(hv.hkl)] = std::norm(hv.value);
struct TwinTerm { double fo, i1, i2; bool free; };
std::vector<TwinTerm> twin_terms;
for (const auto &hv : fmodel.v) {
const auto it = obs_by_hkl.find(hkl_key(hv.hkl));
if (it == obs_by_hkl.end()) continue;
const gemmi::Miller mate = asu.to_asu(twin_law->apply_to_hkl(hv.hkl), gops).first;
const auto m = model_i.find(hkl_key(mate));
if (m == model_i.end()) continue;
twin_terms.push_back({it->second.F, std::norm(hv.value), m->second, it->second.free});
}
double best_r_work = 2.0;
for (int step = 0; step <= 50 && !twin_terms.empty(); ++step) {
const double a = 0.01 * step;
double fo_fc = 0, fc_fc = 0;
for (const auto &t : twin_terms)
if (!t.free) {
const double fc = std::sqrt((1.0 - a) * t.i1 + a * t.i2);
fo_fc += t.fo * fc;
fc_fc += fc * fc;
}
const double k = fc_fc > 0 ? fo_fc / fc_fc : 1.0;
double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
for (const auto &t : twin_terms) {
const double diff = std::fabs(t.fo - k * std::sqrt((1.0 - a) * t.i1 + a * t.i2));
if (t.free) { num_f += diff; den_f += t.fo; }
else { num_w += diff; den_w += t.fo; }
}
if (den_w > 0 && den_f > 0 && num_w / den_w < best_r_work) {
best_r_work = num_w / den_w;
result.r_work_twin = best_r_work;
result.r_free_twin = num_f / den_f;
result.twin_fraction = a;
}
}
result.twin_target = std::isfinite(result.r_free_twin);
if (result.twin_target)
logger.Info("Model validation: against the data as twinned by {} - fraction {:.2f} (scanned on the "
"working set) - R-work {:.4f}, R-free {:.4f}; untwinned {:.4f} / {:.4f}",
twin_law->as_hkl().triplet('h'), result.twin_fraction, result.r_work_twin,
result.r_free_twin, result.r_work, result.r_free);
}
// --- CC(model, data) by resolution shell ---
// The correlation of the merged intensities with |F_model|^2, binned on the merge table's own
// shells so the two tables line up row for row. Nearly free: the scaled Fmodel is already here,
@@ -1655,14 +1711,16 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
size_t nthreads,
double wavelength_A,
const std::vector<float> &report_shell_d_min,
const ModelValidationSchedule &schedule) {
const ModelValidationSchedule &schedule,
const gemmi::Op *twin_law) {
#ifdef JFJOCH_USE_CUDA
// A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through:
// the validation reports why it did not finish, and a validation that did not finish decides nothing,
// so the reflection files are those of a run without a model.
try {
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group,
probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, true);
probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law,
true);
} catch (const RigidBodyGPUFailure &e) {
cuda_clear_error();
ModelValidationResult failed;
@@ -1673,7 +1731,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
}
#else
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, probe_indexing_ambiguity,
nthreads, wavelength_A, report_shell_d_min, schedule, false);
nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false);
#endif
}