From 9a1f729993eaa9c8acf08109e4625dcde597cfac Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 14:34:46 +0200 Subject: [PATCH] rugnux: do not make the first validation's maps where the forecast one is kept Where the validation in the model's setting, started beside the first on the forecast (b5980c9ba), is going to be kept, it writes the maps, the map coefficients and the placed model on the written axes, and the first validation's p_2fofc/p_fofc/p_anom.ccp4 and p_maps.mtz were overwritten unread. ModelValidationSchedule::maps_superseded is asked once the verdict, the setting and the indexing are decided, before any map is made; the first validation then makes none (the map coefficients and the MTZ are not built, the three map transforms not run, nothing written). The answer is the same rule that keeps the forecast validation afterwards (forecast_kept, now one lambda for both). The kept validation also brings back the reflections it relabelled, free set included, so relabel_output does not reindex and redraw them a second time. Where the forecast is discarded or was never started, everything runs as before, and modelpar's delete-before-rewrite stays for those cases. If the kept validation then fails, a warning says no maps were written. Log lines that disappear (first validation only, kept-forecast runs): "mean 2mFo-DFc density at atom centres = ...", "anomalous difference map from ... Bijvoet pairs; strongest density ...", "wrote p_2fofc.ccp4, ..."; in their place: "no maps made here - the validation in the model's setting writes them". (The "anomalous density ... inverted" warning would also go; it cannot be raised where a forecast is kept.) Verified on 8sa8, 8xtg, 9ea5 (forecast kept), 8tyy, 7qis, myob_x10sa (no forecast): p.mtz, p.hkl, p.cif, the three maps, p_maps.mtz and p_model.cif/pdb md5-identical to the merge commit, repeated and at -N 8; the discard path (8sa8 with a rotated model) identical files and log. Time saved is small with the maps already on the GPU: 8sa8 7.6 -> 7.4 s, 8xtg 5.1 -> 5.0 s, 9ea5 within noise. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- .../structure_refinement/ModelValidation.cpp | 312 +++++++++--------- .../structure_refinement/ModelValidation.h | 3 + rugnux/RugnuxScaleMerge.cpp | 44 ++- 3 files changed, 198 insertions(+), 161 deletions(-) diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index 1f0cb5ce2..3c5a2a611 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -1409,113 +1409,119 @@ ModelValidationResult Validate(const std::vector &merged, result.failure_reason = "superseded before its files were written"; return result; } + // Where the validation in the model's setting, started beside this one, is going to be kept, it + // writes the maps and the map coefficients on the axes the files are written on, and this one's would + // be overwritten unread: they are not made at all. + const bool maps_superseded = schedule.maps_superseded && schedule.maps_superseded(result); - // --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) --- - // The three maps - this one, the difference map and the anomalous map - share nothing but what - // they read, so the other two are made and written on threads of their own beside this one. - std::future fofc_map = std::async(std::launch::async, [&] { - write_ccp4(map_from_coefficients(mapfofc, sf_gpu), output_prefix + "_fofc.ccp4"); - }); + if (!maps_superseded) { + // --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) --- + // The three maps - this one, the difference map and the anomalous map - share nothing but what + // they read, so the other two are made and written on threads of their own beside this one. + std::future fofc_map = std::async(std::launch::async, [&] { + write_ccp4(map_from_coefficients(mapfofc, sf_gpu), output_prefix + "_fofc.ccp4"); + }); - // --- anomalous difference map, where the merge kept the Bijvoet split --- - // Coefficients F(+) - F(-) carried on the model phase turned back by 90 degrees. Its peaks sit on - // the anomalous scatterers, so reading the map at each of the model's own atoms names them, - // rather than leaving a list of coordinates for someone to look up. - // Following ANODE, Thorn & Sheldrick (2011) J. Appl. Cryst. 44, 1285-1287 - std::future anomalous_map = std::async(std::launch::async, [&] { - // Read from the merged reflections as they came in, and carry each one into the model's frame - // here: the hand each Bijvoet difference belongs to is a property of the frame the merge was - // made in, and both operators can change it. F(+) and F(-) are attached to the + index of the - // Friedel ASU of that frame, so an anomalous merge - which keeps each mate as a row of its own, - // both carrying the same pair - is read on its + rows only. Taking the - rows as well would - // give one reflection both signs of its difference, and the last row written would decide. - const gemmi::SpaceGroup *data_sg = data_space_group != nullptr ? data_space_group : sg; - const gemmi::ReciprocalAsu data_asu(data_sg); - const gemmi::GroupOps data_gops = data_sg->operations(); + // --- anomalous difference map, where the merge kept the Bijvoet split --- + // Coefficients F(+) - F(-) carried on the model phase turned back by 90 degrees. Its peaks sit on + // the anomalous scatterers, so reading the map at each of the model's own atoms names them, + // rather than leaving a list of coordinates for someone to look up. + // Following ANODE, Thorn & Sheldrick (2011) J. Appl. Cryst. 44, 1285-1287 + std::future anomalous_map = std::async(std::launch::async, [&] { + // Read from the merged reflections as they came in, and carry each one into the model's frame + // here: the hand each Bijvoet difference belongs to is a property of the frame the merge was + // made in, and both operators can change it. F(+) and F(-) are attached to the + index of the + // Friedel ASU of that frame, so an anomalous merge - which keeps each mate as a row of its own, + // both carrying the same pair - is read on its + rows only. Taking the - rows as well would + // give one reflection both signs of its difference, and the last row written would decide. + const gemmi::SpaceGroup *data_sg = data_space_group != nullptr ? data_space_group : sg; + const gemmi::ReciprocalAsu data_asu(data_sg); + const gemmi::GroupOps data_gops = data_sg->operations(); - std::unordered_map danom_by_hkl; - for (const MergedReflection &r : merged) { - if (!std::isfinite(r.F_plus) || !std::isfinite(r.F_minus)) - continue; - gemmi::Op::Miller h{{r.h, r.k, r.l}}; - if (data_gops.is_reflection_centric(h)) // a centric reflection has no anomalous difference - continue; - if (!data_asu.to_asu_sign(h, data_gops).second) - continue; - if (!(result.indexing_op == gemmi::Op::identity())) - h = result.indexing_op.apply_to_hkl(h); - const auto [hasu, plus] = asu.to_asu_sign(h, gops); - danom_by_hkl[hkl_key(hasu)] = plus ? r.F_plus - r.F_minus : r.F_minus - r.F_plus; - } + std::unordered_map danom_by_hkl; + for (const MergedReflection &r : merged) { + if (!std::isfinite(r.F_plus) || !std::isfinite(r.F_minus)) + continue; + gemmi::Op::Miller h{{r.h, r.k, r.l}}; + if (data_gops.is_reflection_centric(h)) // a centric reflection has no anomalous difference + continue; + if (!data_asu.to_asu_sign(h, data_gops).second) + continue; + if (!(result.indexing_op == gemmi::Op::identity())) + h = result.indexing_op.apply_to_hkl(h); + const auto [hasu, plus] = asu.to_asu_sign(h, gops); + danom_by_hkl[hkl_key(hasu)] = plus ? r.F_plus - r.F_minus : r.F_minus - r.F_plus; + } - gemmi::AsuData> mapanom; - mapanom.unit_cell_ = ucell; - mapanom.spacegroup_ = sg; - for (const auto &hv : fmodel.v) { - const auto it = danom_by_hkl.find(hkl_key(hv.hkl)); - if (it == danom_by_hkl.end()) - continue; - const auto phi = static_cast(std::arg(hv.value) - PI / 2); - mapanom.v.push_back({hv.hkl, it->second * std::polar(1.0f, phi)}); - } - result.anomalous_pairs = static_cast(mapanom.v.size()); + gemmi::AsuData> mapanom; + mapanom.unit_cell_ = ucell; + mapanom.spacegroup_ = sg; + for (const auto &hv : fmodel.v) { + const auto it = danom_by_hkl.find(hkl_key(hv.hkl)); + if (it == danom_by_hkl.end()) + continue; + const auto phi = static_cast(std::arg(hv.value) - PI / 2); + mapanom.v.push_back({hv.hkl, it->second * std::polar(1.0f, phi)}); + } + result.anomalous_pairs = static_cast(mapanom.v.size()); - if (!mapanom.v.empty()) { - const gemmi::Grid grid = map_from_coefficients(mapanom, sf_gpu); - const double rms = write_ccp4(grid, output_prefix + "_anom.ccp4"); - std::vector sites; - double scatterer_sum = 0; + if (!mapanom.v.empty()) { + const gemmi::Grid grid = map_from_coefficients(mapanom, sf_gpu); + const double rms = write_ccp4(grid, output_prefix + "_anom.ccp4"); + std::vector sites; + double scatterer_sum = 0; + for (gemmi::Model &m : st.models) + for (gemmi::Chain &ch : m.chains) + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) { + if (a.is_hydrogen()) // hydrogen scatters no anomalous signal + continue; + sites.push_back({fmt::format("{} {} {}{}", a.name, r.name, ch.name, + r.seqid.str()), + rms > 0 ? grid.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER) / rms + : 0.0}); + if (a.element.atomic_number() >= ANOMALOUS_SCATTERER_MIN_Z) { + scatterer_sum += sites.back().sigma; + result.anomalous_scatterers++; + } + } + if (result.anomalous_scatterers > 0) + result.anomalous_scatterer_mean_sigma = scatterer_sum / result.anomalous_scatterers; + std::sort(sites.begin(), sites.end(), + [](const auto &x, const auto &y) { return x.sigma > y.sigma; }); + + // A model and a dataset in opposite hands turn every anomalous peak into a trough, so a + // map whose deepest hole at an atom is both deep and deeper than its highest peak says + // the two disagree about the hand. That is worth reporting: it is real evidence about + // the crystal, and the alternative - reindexing until the two agree - would erase it. + if (!sites.empty()) { + const auto &deepest = sites.back(); + if (deepest.sigma < -ANOMALOUS_INVERSION_SIGMA && -deepest.sigma > sites.front().sigma) { + result.anomalous_hands_disagree = true; + result.anomalous_deepest_site = deepest.label; + result.anomalous_deepest_sigma = deepest.sigma; + } + } + + if (sites.size() > MAX_ANOMALOUS_SITES) + sites.resize(MAX_ANOMALOUS_SITES); + result.anomalous_sites = std::move(sites); + } + }); + + const gemmi::Grid grid2fofc = map_from_coefficients(map2fofc, sf_gpu); + const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4"); + { + double s = 0; int n = 0; for (gemmi::Model &m : st.models) for (gemmi::Chain &ch : m.chains) for (gemmi::Residue &r : ch.residues) - for (gemmi::Atom &a : r.atoms) { - if (a.is_hydrogen()) // hydrogen scatters no anomalous signal - continue; - sites.push_back({fmt::format("{} {} {}{}", a.name, r.name, ch.name, - r.seqid.str()), - rms > 0 ? grid.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER) / rms - : 0.0}); - if (a.element.atomic_number() >= ANOMALOUS_SCATTERER_MIN_Z) { - scatterer_sum += sites.back().sigma; - result.anomalous_scatterers++; - } - } - if (result.anomalous_scatterers > 0) - result.anomalous_scatterer_mean_sigma = scatterer_sum / result.anomalous_scatterers; - std::sort(sites.begin(), sites.end(), - [](const auto &x, const auto &y) { return x.sigma > y.sigma; }); - - // A model and a dataset in opposite hands turn every anomalous peak into a trough, so a - // map whose deepest hole at an atom is both deep and deeper than its highest peak says - // the two disagree about the hand. That is worth reporting: it is real evidence about - // the crystal, and the alternative - reindexing until the two agree - would erase it. - if (!sites.empty()) { - const auto &deepest = sites.back(); - if (deepest.sigma < -ANOMALOUS_INVERSION_SIGMA && -deepest.sigma > sites.front().sigma) { - result.anomalous_hands_disagree = true; - result.anomalous_deepest_site = deepest.label; - result.anomalous_deepest_sigma = deepest.sigma; - } - } - - if (sites.size() > MAX_ANOMALOUS_SITES) - sites.resize(MAX_ANOMALOUS_SITES); - result.anomalous_sites = std::move(sites); + for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; } + result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0; } - }); - - const gemmi::Grid grid2fofc = map_from_coefficients(map2fofc, sf_gpu); - const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4"); - { - double s = 0; int n = 0; - for (gemmi::Model &m : st.models) - for (gemmi::Chain &ch : m.chains) - for (gemmi::Residue &r : ch.residues) - for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; } - result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0; + fofc_map.get(); + anomalous_map.get(); } - fofc_map.get(); - anomalous_map.get(); // --- the enantiomorph, now that there is something to decide it on --- // Two things have to hold before the model's hand is written on these data. The model has to @@ -1536,55 +1542,57 @@ ModelValidationResult Validate(const std::vector &merged, ? "the anomalous density at the model's atoms is inverted" : "the model does not fit these data"); - // --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. --- - try { - gemmi::Mtz mtz(true); - // The group the REFLECTIONS end up in, which is the model's only where its hand was adopted - - // AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes - // identically, so the coefficients are the same numbers either way and only the label moves; - // but the two groups have different screw translations, so a reader that expands symmetry out - // of this file works in the wrong one if the label disagrees with the .mtz beside it. - mtz.spacegroup = result.adopted_model_enantiomorph - ? sg - // No data group means the caller merged in P1, and P1 is what the - // reflection files beside this one carry - not the model's group. - : (data_space_group != nullptr ? data_space_group - : gemmi::find_spacegroup_by_number(1)); - mtz.set_cell_for_all(ucell); - mtz.add_dataset("model_validation"); - mtz.datasets.back().wavelength = wavelength_A; - mtz.add_column("FP", 'F', -1, -1, false); - mtz.add_column("FC", 'F', -1, -1, false); - mtz.add_column("PHIC", 'P', -1, -1, false); - mtz.add_column("FWT", 'F', -1, -1, false); - mtz.add_column("PHWT", 'P', -1, -1, false); - mtz.add_column("DELFWT", 'F', -1, -1, false); - mtz.add_column("PHDELWT", 'P', -1, -1, false); - // The figure of merit the coefficients carry, so the weighting can be read off - and undone - - // from the file rather than having to be taken on trust. - mtz.add_column("FOM", 'W', -1, -1, false); - mtz.add_column("FREE", 'I', -1, -1, false); - std::vector data; - for (size_t i = 0; i < terms.size(); i++) { - // [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg - // returns (-180, 180] and a P column is not supposed to. - const double phi_wrapped = terms[i].phi * 180.0 / PI; - const auto phi_deg = static_cast(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped); - data.insert(data.end(), {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), - static_cast(terms[i].hkl[2]), - static_cast(terms[i].fo), static_cast(terms[i].fc), - phi_deg, - fwt[i], phi_deg, - delfwt[i], phi_deg, - static_cast(sigma_a.weight[i].m), - terms[i].free ? 0.0f : 1.0f}); + if (!maps_superseded) { + // --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. --- + try { + gemmi::Mtz mtz(true); + // The group the REFLECTIONS end up in, which is the model's only where its hand was adopted - + // AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes + // identically, so the coefficients are the same numbers either way and only the label moves; + // but the two groups have different screw translations, so a reader that expands symmetry out + // of this file works in the wrong one if the label disagrees with the .mtz beside it. + mtz.spacegroup = result.adopted_model_enantiomorph + ? sg + // No data group means the caller merged in P1, and P1 is what the + // reflection files beside this one carry - not the model's group. + : (data_space_group != nullptr ? data_space_group + : gemmi::find_spacegroup_by_number(1)); + mtz.set_cell_for_all(ucell); + mtz.add_dataset("model_validation"); + mtz.datasets.back().wavelength = wavelength_A; + mtz.add_column("FP", 'F', -1, -1, false); + mtz.add_column("FC", 'F', -1, -1, false); + mtz.add_column("PHIC", 'P', -1, -1, false); + mtz.add_column("FWT", 'F', -1, -1, false); + mtz.add_column("PHWT", 'P', -1, -1, false); + mtz.add_column("DELFWT", 'F', -1, -1, false); + mtz.add_column("PHDELWT", 'P', -1, -1, false); + // The figure of merit the coefficients carry, so the weighting can be read off - and undone - + // from the file rather than having to be taken on trust. + mtz.add_column("FOM", 'W', -1, -1, false); + mtz.add_column("FREE", 'I', -1, -1, false); + std::vector data; + for (size_t i = 0; i < terms.size(); i++) { + // [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg + // returns (-180, 180] and a P column is not supposed to. + const double phi_wrapped = terms[i].phi * 180.0 / PI; + const auto phi_deg = static_cast(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped); + data.insert(data.end(), {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), + static_cast(terms[i].hkl[2]), + static_cast(terms[i].fo), static_cast(terms[i].fc), + phi_deg, + fwt[i], phi_deg, + delfwt[i], phi_deg, + static_cast(sigma_a.weight[i].m), + terms[i].free ? 0.0f : 1.0f}); + } + mtz.nreflections = static_cast(terms.size()); + mtz.data = std::move(data); + remove_before_rewriting(output_prefix + "_maps.mtz"); + mtz.write_to_file(output_prefix + "_maps.mtz"); + } catch (const std::exception &e) { + logger.Warning("Model validation: could not write map MTZ: {}", e.what()); } - mtz.nreflections = static_cast(terms.size()); - mtz.data = std::move(data); - remove_before_rewriting(output_prefix + "_maps.mtz"); - mtz.write_to_file(output_prefix + "_maps.mtz"); - } catch (const std::exception &e) { - logger.Warning("Model validation: could not write map MTZ: {}", e.what()); } result.ok = true; @@ -1596,7 +1604,8 @@ ModelValidationResult Validate(const std::vector &merged, result.r_work, result.n_work, result.r_free, result.n_free); logger.Info("Model validation: bulk solvent k_sol={:.3f} b_sol={:.1f}, k_overall={:.3f}", result.k_sol, result.b_sol, result.k_overall); - logger.Info("Model validation: mean 2mFo-DFc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma); + if (!maps_superseded) + logger.Info("Model validation: mean 2mFo-DFc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma); logger.Info("Model validation: sigma_A weighting over {} resolution shell(s) (estimated on the {} " "free reflections): sigma_A {:.2f} at low resolution, {:.2f} at high, mean FOM {:.3f}", sigma_a.shells, sigma_a.free_reflections, sigma_a.sigma_a_lowest_shell, @@ -1615,10 +1624,13 @@ ModelValidationResult Validate(const std::vector &merged, "agree - either the model is the wrong enantiomorph for this crystal, or the " "data were indexed in the wrong hand, and reindexing would hide which", result.anomalous_deepest_site, result.anomalous_deepest_sigma); - logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4{}, {}_maps.mtz", - output_prefix, output_prefix, - result.anomalous_sites.empty() ? "" : ", " + output_prefix + "_anom.ccp4", - output_prefix); + if (maps_superseded) + logger.Info("Model validation: no maps made here - the validation in the model's setting writes them"); + else + logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4{}, {}_maps.mtz", + output_prefix, output_prefix, + result.anomalous_sites.empty() ? "" : ", " + output_prefix + "_anom.ccp4", + output_prefix); if (result.fit_tested && !result.model_fits) logger.Info("Model validation: the model was rejected, so nothing downstream moved - the " "reflection files are byte for byte the ones a run with no model would have " diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index a2f24f3cc..c89feb567 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -246,6 +246,9 @@ struct ModelValidationSchedule { std::function on_forecast; // Where valid, waited on before the first file is written: false returns without writing any. std::shared_future write_gate; + // Asked once the verdict, the setting and the indexing are decided, before any map is made: true + // where another validation will write this one's maps and map coefficients, which are then not made. + std::function maps_superseded; }; ModelValidationResult ValidateAgainstModel(const std::vector &merged, diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 3845d2821..cfc1d6a2d 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -2612,7 +2612,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // report states: every relabelling since the group was determined. gemmi::Op merge_to_written = gemmi::Op::identity(); gemmi::Op determined_to_written = gemmi::Op::identity(); - const auto relabel_output = [&](const gemmi::Op &cob, const char *source) { + // `relabelled`: the merged reflections are already in the new setting, free set included - made by the + // validation started on the forecast from the same reflections the same way - so only the rest moves. + const auto relabel_output = [&](const gemmi::Op &cob, const char *source, + std::vector *relabelled = nullptr) { result.setting_source = source; if (cob == gemmi::Op::identity()) return; @@ -2624,14 +2627,18 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { return; } const gemmi::Op hkl = HklOperator(cob); - sm.merged = ReindexMergedIntoAsu(sm.merged, hkl, *moved, - experiment_.GetScalingSettings().GetMergeFriedel()); + if (relabelled == nullptr) + sm.merged = ReindexMergedIntoAsu(sm.merged, hkl, *moved, + experiment_.GetScalingSettings().GetMergeFriedel()); + else + sm.merged = std::move(*relabelled); reseat_integrated(cob); relabel_group(cob, *moved); // The free set is a function of the written index (RfreeFlags.h), so that every dataset of // a crystal form written on the same axes carries the same one: drawn again on these axes. - AssignRfreeFlags(sm.merged, *moved, experiment_.GetScalingSettings().GetRfreeFraction(), 500, - result.consensus_cell, config_.nthreads); + if (relabelled == nullptr) + AssignRfreeFlags(sm.merged, *moved, experiment_.GetScalingSettings().GetRfreeFraction(), 500, + result.consensus_cell, config_.nthreads); merge_to_written = merge_to_written.combine(hkl); determined_to_written = determined_to_written.combine(hkl); result.setting_operator = IndexTriplet(determined_to_written); @@ -3257,6 +3264,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // only known at the end, from the anomalous map. struct Ahead { gemmi::Op change_of_basis_op, indexing_op; + // The reflections it validates against, relabelled as the files would be: where it is kept, + // they are the written ones, and the relabelling below is not made a second time. + std::vector merged; Logger log = Logger::Buffered(); std::future run; std::promise write; // destroyed before `run`, so a run left waiting on it is released @@ -3291,12 +3301,22 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } ModelValidationSchedule schedule; schedule.write_gate = gate; - return ValidateAgainstModel(merged, cell, config_.model_path, config_.output_prefix, ahead.log, + ahead.merged = std::move(merged); + return ValidateAgainstModel(ahead.merged, cell, config_.model_path, config_.output_prefix, ahead.log, in, /*probe_indexing_ambiguity=*/false, static_cast(config_.nthreads), wavelength, report_shell_d_min, schedule); }); }; + // Whether the validation started on the forecast is the one kept: the decision is the one forecast. + // Asked by the first validation before it makes its maps, which are then left to that one, and + // again below with the verdict in hand - the same answer, as nothing it reads changes between. + const auto forecast_kept = [&](const ModelValidationResult &v) { + return ahead.run.valid() && config_.reference_data.empty() && v.model_fits + && !(v.change_of_basis_op == gemmi::Op::identity()) && !v.adopted_model_enantiomorph + && v.change_of_basis_op == ahead.change_of_basis_op && v.indexing_op == ahead.indexing_op; + }; + first_schedule.maps_superseded = forecast_kept; // The ledger's P1 merge (started above) is taken before anything below acts on what the // validation decides: relabelling the outcomes and the group is what it must not see happen. // Every such relabelling reaches the P1 merge afterwards, through merge_to_written. @@ -3368,9 +3388,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // second one, facing a model already in the data's setting, has nothing left to decide. const bool remake = config_.reference_data.empty() && validation.model_fits && !(validation.change_of_basis_op == gemmi::Op::identity()); - const bool forecast_held = remake && ahead.run.valid() && !validation.adopted_model_enantiomorph - && validation.change_of_basis_op == ahead.change_of_basis_op - && validation.indexing_op == ahead.indexing_op; + const bool forecast_held = forecast_kept(validation); if (ahead.run.valid() && !forecast_held) { ahead.write.set_value(false); ahead.run.wait(); @@ -3378,14 +3396,15 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (remake) { gemmi::Op to_model = validation.change_of_basis_op; to_model.tran = {0, 0, 0}; - relabel_output(to_model.inverse(), "MODEL"); - const auto sg_now = experiment_.GetGemmiSpaceGroup(); ModelValidationResult remade; if (forecast_held) { ahead.write.set_value(true); remade = ahead.run.get(); + relabel_output(to_model.inverse(), "MODEL", &ahead.merged); ahead.log.ReplayInto(logger); } else { + relabel_output(to_model.inverse(), "MODEL"); + const auto sg_now = experiment_.GetGemmiSpaceGroup(); remade = ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, config_.output_prefix, logger, sg_now ? &*sg_now : nullptr, @@ -3397,6 +3416,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { KeepModelVerdict(remade, validation); validation = std::move(remade); result.model_validation = validation; + } else if (forecast_held) { + logger.Warning("Model validation: the validation in the model's setting failed ({}), and the " + "first left its maps to it - no maps were written", remade.failure_reason); } } // The coordinates that go with the maps just written, in the frame the reflection files