diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index 097512a26..c9568d994 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -313,8 +313,10 @@ are: as TLSMD cuts them (by dynamic programming over the least-squares fit of an isotropic TLS model, a quadratic B(r), to the refined B factors), each cut where it gives more groups than the one before; and a U per atom, at the U weight balanced as $w_{B,0}$ is (on the U shaken by 20% of its B) times $\{1/4, 1\}$. The cuts between one group - per chain and the finest run only where the finest fits the free set better (lower $-\ln L$) than one per chain; - on 49 refinements of 33 structures neither this nor leaving out the U weight $\times 4$ changed a choice taken. A + per chain and the finest run only where the finest fits the free set better (lower $-\ln L$) than one per chain, + and the U at its balanced weight only where the U at a quarter of it lowers R-free or $-\ln L$ of the free + reflections from where the weights ended; on 49 refinements of 33 structures neither the cuts' rule nor leaving + out the U weight $\times 4$ changed a choice taken. A trial displaces the one before it (the first: where the weights ended) only where its geometry is within bounds and both its R-free (beyond the tie) and the $-\ln L$ of the free reflections are lower: the free set tests the parameters a model adds, so the choice is made on the evidence of each run, not at a resolution (Merritt 2012). @@ -322,6 +324,12 @@ are: `REFINE_ANISO_ATOMS`, `REFINE_TLS_GROUPS`). It is decided after the solvent build (§15.10), on the model with its solvent: a water a round put in refines a U of its own in the per-atom trial as every refinable atom does, and no TLS group holds a water. +5. **The model as given** is the first candidate in order of trust: where the refinement does not displace it by the + rule every trial is held to (its geometry within bounds and R-free lower beyond the tie), the model is written as + it came, scaled to the data, and the solvent build's edits are undone (`REFINE_AS_GIVEN`). What comes out is never + worse on the free set than what went in. This holds only where the free set is free to the model as given: a model + placed as given on a free set drawn by the run (`MODEL_FREE_SET`; in `--mode refine`, a file without free flags) has + most likely been refined against those reflections, and its R-free as given does not compete. Where the model has copies (§15.5), the NCS restraints are one more weight: every rung of the first ladder is run without them, and then with them at $\sigma = 0.1$ Å at the B factors' weight the rungs without chose (NCS moves the diff --git a/docs/RUGNUX_MODEL_AND_REFERENCE.md b/docs/RUGNUX_MODEL_AND_REFERENCE.md index 7aca9af63..8c3da7c74 100644 --- a/docs/RUGNUX_MODEL_AND_REFERENCE.md +++ b/docs/RUGNUX_MODEL_AND_REFERENCE.md @@ -339,7 +339,13 @@ columns are found, not named: amplitudes with their sigmas (F/SIGF, FP/SIGFP, FO map coefficients such as FC, FWT, DELFWT are passed over) where the file has them, else intensities (IMEAN/SIGIMEAN, I/SIGI, ...) made into amplitudes by French-Wilson, else the mean of the anomalous pairs. The free set is read from a FreeR_flag-like column (flag 0 is the test set; where 0 marks the majority, as phenix and CNS write it, the -minority value is); a file without one gets rugnux's own 5 % set. The refinement is the one run after a merge on a placed model, and +minority value is; in a structure-factor mmCIF, `_refln.status` 'f' first); a file without one gets rugnux's own 5 % set. +Where a structure-factor mmCIF's flags are classes (`_refln.pdbx_r_free_flag` 0-19), the class the file marks is not +always the one the depositors held out, so the model as given is scored on each class: the mark is kept where the model +fits that class worst or the test-set count the model file records (REMARK 3, `_refine.ls_number_reflns_R_free`) is +nearest its count; otherwise the class the model fits worst is the test set where that count is nearest it, or where +it is above the marked class by more than twice the scatter of the other classes' R (`REFINE_FREE_SET_SOURCE= +STATUS|REMARK3_COUNT|WORST_FIT`, `REFINE_FREE_SET_CLASS`). The refinement is the one run after a merge on a placed model, and writes the same files - `_refined_model.cif`/`.pdb`, `_refined_maps.mtz`, `_blobs.json` - with the `REFINE_` keys in `_report.txt`, beside `REFINE_FRENCH_WILSON=` and `REFINE_FREE_FLAGS= FILE|GENERATED` saying what had to be made of the input. diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 956d1aca6..79ddbfb6b 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1235,6 +1235,7 @@ Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md | `REFINE_RIGID_BODIES` | The chains moved as rigid bodies, where that step was kept (0: not kept) | | `REFINE_R_WORK_PLACEMENT`, `REFINE_R_FREE_PLACEMENT` | After the placement stage (rigid bodies, then one B per residue, each kept only where it lowers R-free) | | `REFINE_R_WORK`, `REFINE_R_FREE` | After the restrained stage: the refined model, the one written as `_refined_model.cif` | +| `REFINE_AS_GIVEN` | `TRUE` where no refinement lowered R-free beyond the tie of the model as given, which is then the model written (as given, scaled to the data) | | `REFINE_RMSD_BONDS`, `REFINE_RMSD_ANGLES` | The refined model's rmsd from the dictionary's bonds (Å) and angles (°); the weight is chosen within 0.02 Å and 2° | | `REFINE_RMSD_METAL_BONDS` | The rmsd of the metal-donor bonds made from the model's LINK records (those within the metal's inner sphere), apart from the covalent bonds. Present where there are any | | `REFINE_RMSD_BONDS_START`, `REFINE_RMSD_ANGLES_START` | *dev*. The same for the model as given | @@ -1307,4 +1308,7 @@ resolution shells. `rugnux --mode refine --refine-score-only` writes this score In `--mode refine` the input's keys come first: `REFINE_REFLECTIONS` and `REFINE_FREE_REFLECTIONS` (those read, after any `--refine-d-min`), `REFINE_FRENCH_WILSON` (`TRUE` where the file had intensities only and the amplitudes were made by French-Wilson), `REFINE_FREE_FLAGS` (`FILE`, or `GENERATED` where the file had no free set and -rugnux drew its own) and `REFINE_D_MIN_LIMIT` (the `--refine-d-min` given). +rugnux drew its own), `REFINE_FREE_SET_SOURCE` and `REFINE_FREE_SET_CLASS` (where a structure-factor mmCIF's +flags are classes: the class taken as the test set, and whether the file's own mark (`STATUS`), the model file's +test-set count with the worst fit (`REMARK3_COUNT`) or the worst fit alone (`WORST_FIT`) chose it) and +`REFINE_D_MIN_LIMIT` (the `--refine-d-min` given). diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index dc18038e8..59b52ed94 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -2994,9 +2994,12 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec copy_g = {{}}; copy_rg = {{}}; } - glog.Info("Molecular replacement: {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", + // a placement that clashes too much to be one is not refined, and has no R-factors + glog.Info("Molecular replacement: {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}% {}", sg.xhm(), rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * f.clash, - score.r_work, score.r_free); + std::isfinite(score.r_free) + ? fmt::format("R-work {:.3f} R-free {:.3f}", score.r_work, score.r_free) + : std::string("not refined")); } // The likelihood of the placement kept, judged against the same orientation elsewhere - only where // its translation function's Z-score does not already stand out. @@ -3166,9 +3169,13 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec trial.refined = LikelihoodTranslationZ(lsetup(), *lknown, known_llg, trial.g, centroid(trial.copy), model.error_A, nthreads); } - glog.Info("Molecular replacement: {:<12} copy {} TF Z {:5.1f} LLG gain {:7.1f}, refined {:7.1f} Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", - sg.xhm(), n, tf2[trial.i].z, trial.gain, trial.refined.gain, - trial.refined.z, 100 * trial.clash, trial.score.r_work, trial.score.r_free); + // the likelihood's own Z only where the translation function's does not already stand out + glog.Info("Molecular replacement: {:<12} copy {} TF Z {:5.1f} LLG gain {:7.1f}, {} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", + sg.xhm(), n, tf2[trial.i].z, trial.gain, + std::isfinite(trial.refined.z) + ? fmt::format("refined {:7.1f} Z {:5.1f}", trial.refined.gain, trial.refined.z) + : std::string("likelihood Z not needed"), + 100 * trial.clash, trial.score.r_work, trial.score.r_free); // The first, in that order, that stands out, packs and lowers R-free. if ((tf2[trial.i].z >= TF_Z_PLACED || trial.refined.z >= TF_Z_PLACED) && trial.clash <= CLASH_MAX && trial.score.r_free <= cand.r_free - COPY_R_FREE_GAIN) { @@ -3176,8 +3183,12 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec break; } } - if (!best) + if (!best) { + glog.Info("Molecular replacement: {:<12} copy {}: none of {} trials stands out, packs and lowers " + "R-free by {:.3f}; {} {} kept", sg.xhm(), n, trials.size(), COPY_R_FREE_GAIN, + cand.copies, cand.copies == 1 ? "copy" : "copies"); break; + } copies.push_back(best->copy); copy_rotation.push_back(best->rotation); copy_position.push_back(centroid(best->copy)); diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 9fdf78de0..2a76ce97d 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -2209,6 +2209,27 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector given_atoms; + std::vector> given_residues; + const SolventPool given_pool = pool; + if (settings.restrained && !settings.score_only) + for (const gemmi::Chain &ch : ctx.Model().chains) + for (const gemmi::Residue &res : ch.residues) { + given_residues.emplace_back(res.name, res.entity_type); + given_atoms.insert(given_atoms.end(), res.atoms.begin(), res.atoms.end()); + } + // The waters and ions the build relabelled and kept (SolventEdit::site their old name), for the metal records the // refined model is written with. std::vector relabels; @@ -2385,6 +2406,11 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector hyps = AdpHypotheses(ctx, topo, input_tls); const double w = state.trial.weight, adp = state.trial.adp_weight > 0 ? state.trial.adp_weight : adp_w0; // the U weight balanced on the model with every refinable atom given a U @@ -2956,13 +2984,19 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector= 0 && outcomes[k].ran && outcomes[k].trial.within_bounds && outcomes[k].trial.r_free < state.trial.r_free; }; - const bool tls = lowers(coarsest) || lowers(finest), u_full = lowers(u_quarter); + // the U at a quarter fitting the free reflections better than where the weights ended asks for the + // balanced weight as well: the choice below is made on -ln L free as much as on R-free + auto fits_better = [&](int k) { + return k >= 0 && outcomes[k].ran && outcomes[k].trial.within_bounds && + outcomes[k].trial.nll_free < state.trial.nll_free; + }; + const bool tls = lowers(coarsest) || lowers(finest), u_full = lowers(u_quarter) || fits_better(u_quarter); const bool cuts = tls && finest > coarsest && outcomes[finest].ran && outcomes[coarsest].ran && outcomes[finest].trial.nll_free < outcomes[coarsest].trial.nll_free; std::vector needed(trials.size(), 0); @@ -2992,6 +3026,28 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vectorevaluations; ctx.evaluations -= ahead_evaluations; diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index 6205c02ca..191e0f6fd 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -78,6 +78,10 @@ struct RefinementSettings { // The data's own intensity statistics read twinned - the caller's verdict, from the merge (the L-test): the law // is then taken from 2 standard errors over its null instead of 3 (AdoptTwinLaw). bool twin_data_indicated = false; + // The model as given a candidate of the final choice (the restrained stage, 5.): only where its R-free is a free + // statistic - false where the model has most likely been refined against the free reflections too (a deposited + // model of this crystal scored on a free set drawn here). + bool as_given_competes = true; bool gpu = true; size_t nthreads = 1; size_t gpu_resident_bytes = 0; // held on the card by what runs beside @@ -159,6 +163,7 @@ struct RefinementResult { double r_work_placement = 1, r_free_placement = 1; // After the restrained stage (NaN where it did not run). bool restrained = false; + bool as_given = false; // the model as given kept: no trial displaced it (the restrained stage ran) bool cancelled = false; // cut short by RefinementSettings::cancel double r_work = NAN, r_free = NAN; double weight = 0; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 686203368..06bc30f0f 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -569,6 +569,7 @@ namespace { if (rf.restrained) { Add(s, KeyReal("REFINE_R_WORK", rf.r_work, "{:.4f}")); Add(s, KeyReal("REFINE_R_FREE", rf.r_free, "{:.4f}")); + Add(s, KeyBool("REFINE_AS_GIVEN", rf.as_given)); Add(s, KeyReal("REFINE_RMSD_BONDS", rf.geometry.rmsd_bonds, "{:.4f}")); if (rf.geometry.metal_bonds > 0) Add(s, KeyReal("REFINE_RMSD_METAL_BONDS", rf.geometry.rmsd_metal, "{:.4f}")); @@ -3438,6 +3439,10 @@ ReportDocument BuildRefineReportDocument(const std::string &output_prefix, const Add(s, KeyInt("REFINE_FREE_REFLECTIONS", r.free_reflections)); Add(s, KeyBool("REFINE_FRENCH_WILSON", r.french_wilson)); Add(s, KeyEnum("REFINE_FREE_FLAGS", r.free_flags_generated ? "GENERATED" : "FILE", {"FILE", "GENERATED"})); + if (!r.free_set_source.empty()) { + Add(s, KeyEnum("REFINE_FREE_SET_SOURCE", r.free_set_source, {"STATUS", "REMARK3_COUNT", "WORST_FIT"})); + Add(s, KeyInt("REFINE_FREE_SET_CLASS", r.free_class)); + } if (r.d_min_limit > 0) Add(s, KeyReal("REFINE_D_MIN_LIMIT", r.d_min_limit, "{:.3f}")); AddRefinementKeys(s, rf, output_prefix); diff --git a/rugnux/RugnuxRefine.cpp b/rugnux/RugnuxRefine.cpp index fd8e51d90..fe217da2b 100644 --- a/rugnux/RugnuxRefine.cpp +++ b/rugnux/RugnuxRefine.cpp @@ -3,13 +3,134 @@ #include "RugnuxRefine.h" +#include +#include #include +#include + +#include "gemmi/symmetry.hpp" #include "../common/UnitCell.h" #include "../image_analysis/scale_merge/FrenchWilson.h" #include "../image_analysis/scale_merge/RfreeFlags.h" #include "../image_analysis/structure_refinement/ModelValidation.h" // SameIntensities +namespace { + +// The test set of a structure-factor mmCIF whose free flags are classes (0-19): the file marks one (_refln.status +// 'f', or flag 0), but that mark is the archive's reading of the flags and is not always the class the depositors' +// refinement held out. The model as given tells: it fits worst the class it was not refined against, and the +// test-set count its file records (REMARK 3 FREE R VALUE TEST SET COUNT, _refine.ls_number_reflns_R_free) is that +// class's count within the refinement's resolution range. The count alone does not single a class out (of 20 classes +// of about the same size another often has it too, or the marked one is off by a reflection or two), so the mark is +// kept ("STATUS") where the model fits it worst or the count is nearest it. Else the class the model fits worst is +// taken where the count is nearest it ("REMARK3_COUNT"), or where it is above the marked one by more than twice the +// scatter of the R of the other classes, so that the marked class is fitted like a working one ("WORST_FIT"; the +// model file gives no count, or the count points nowhere). `class_of`: the class of each reflection of the file; `counts`: per class, its reflections within that range. +void ChooseFreeClass(std::vector &refl, const std::map &class_of, + const std::map &counts, double d_hi, double d_lo, const gemmi::Structure &model, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg, double wavelength_A, const RefinementSettings &settings, + Logger &logger, RefineOnlyResult &out) { + std::vector cls(refl.size(), -1); + for (size_t i = 0; i < refl.size(); i++) + if (const auto c = class_of.find({{refl[i].h, refl[i].k, refl[i].l}}); c != class_of.end()) + cls[i] = c->second; + // the class the file marks: that of most reflections flagged free + std::map flagged; + for (size_t i = 0; i < refl.size(); i++) + if (refl[i].rfree_flag && cls[i] >= 0) + flagged[cls[i]]++; + int marked = -1; + for (const auto &[c, n] : flagged) + if (marked < 0 || n > flagged[marked]) + marked = c; + + // the model as given, scored against these flags: R per class + RefinementSettings score = settings; + score.score_only = true; + score.twin = false; + score.solvent = false; + logger.Info("Refinement input: the free flags are {} classes; the model as given scored on each", counts.size()); + const RefinementResult s = RefineModel(model, refl, cell, sg, wavelength_A, logger, score); + if (!s.ok) + return; + const gemmi::GroupOps gops = sg.operations(); + const gemmi::ReciprocalAsu asu(&sg); + auto in_asu = [&](gemmi::Miller h) { return asu.is_in(h) ? h : asu.to_asu(h, gops).first; }; + std::map fmodel; + for (size_t m = 0; m < s.score_start.hkl.size(); m++) + fmodel[in_asu(s.score_start.hkl[m])] = + s.score_start.k_total[m] * std::abs(s.score_start.fcalc[m] + s.score_start.fmask[m]); + std::map> sums; // class: sum |Fo - |Fmodel||, sum Fo + for (size_t i = 0; i < refl.size(); i++) { + const auto f = fmodel.find(in_asu({{refl[i].h, refl[i].k, refl[i].l}})); + if (cls[i] < 0 || f == fmodel.end() || refl[i].d < d_hi || refl[i].d > d_lo) + continue; + sums[cls[i]][0] += std::fabs(refl[i].F - f->second); + sums[cls[i]][1] += refl[i].F; + } + int worst = -1; + std::map r; + for (const auto &[c, v] : sums) { + r[c] = v[1] > 0 ? v[0] / v[1] : 0; + if (worst < 0 || r[c] > r[worst]) + worst = c; + } + + // the classes whose count is nearest the model file's + const int n_test = model.meta.refinement.empty() ? -1 : model.meta.refinement[0].rfree_set_count; + std::vector nearest; + if (n_test > 0) { + int best = -1; + for (const auto &[c, n] : counts) { + const int dist = std::abs(n - n_test); + if (best < 0 || dist < best) { + best = dist; + nearest.clear(); + } + if (dist == best) + nearest.push_back(c); + } + } + auto is_nearest = [&](int c) { return std::find(nearest.begin(), nearest.end(), c) != nearest.end(); }; + + // the scatter of the classes but the worst: how much a class the model was refined against varies + double mean = 0, var = 0; + for (const auto &[c, v] : r) + if (c != worst) + mean += v / static_cast(r.size() - 1); + for (const auto &[c, v] : r) + if (c != worst) + var += (v - mean) * (v - mean) / static_cast(r.size() - 2); + const double scatter = std::sqrt(var); + const bool fitted_like_work = r.count(marked) > 0 && r[worst] - r[marked] > 2 * scatter; + + int chosen = marked; + out.free_set_source = "STATUS"; + if (worst != marked && !is_nearest(marked) && (is_nearest(worst) || fitted_like_work)) { + chosen = worst; + out.free_set_source = is_nearest(worst) ? "REMARK3_COUNT" : "WORST_FIT"; + } + out.free_class = chosen; + + std::string by_class, near; + for (const auto &[c, v] : r) + by_class += fmt::format("{}{} {:.3f}", by_class.empty() ? "" : ", ", c, v); + for (int c : nearest) + near += fmt::format("{}{}", near.empty() ? "" : ", ", c); + logger.Info("Refinement input: R as given by free-flag class: {}; the file marks class {}, the model fits class {} " + "worst ({:.1f} times the scatter {:.4f} of the others above the marked one){} -> class {} is the test " + "set ({})", by_class, marked, worst, scatter > 0 ? (r[worst] - r[marked]) / scatter : 0.0, scatter, + n_test > 0 ? fmt::format(", the model file's test-set count {} is nearest class {}", n_test, near) + : std::string(", the model file gives no test-set count"), + chosen, out.free_set_source); + if (chosen != marked) + for (size_t i = 0; i < refl.size(); i++) + refl[i].rfree_flag = cls[i] == chosen; +} + +} // namespace + RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &model, double d_min_limit, const RefinementSettings &settings, Logger &logger) { RefineOnlyResult out; @@ -31,6 +152,26 @@ RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &m out.space_group = sg->xhm(); out.wavelength_A = data.wavelength_A > 0 ? data.wavelength_A : 1.0; std::vector &refl = data.reflections; + // free-flag classes, by reflection, and the reflections of each within the resolution range the model file + // records its refinement to, counted before any limit here + std::map class_of; + std::map class_counts; + double d_hi = 0, d_lo = INFINITY; + if (!data.free_column.empty() && data.free_class.size() == refl.size()) { + if (!model.meta.refinement.empty()) { + const gemmi::RefinementInfo &ri = model.meta.refinement[0]; + if (std::isfinite(ri.resolution_high) && ri.resolution_high > 0) + d_hi = ri.resolution_high - 1e-3; + if (std::isfinite(ri.resolution_low) && ri.resolution_low > 0) + d_lo = ri.resolution_low + 1e-3; + } + for (size_t i = 0; i < refl.size(); i++) + if (data.free_class[i] >= 0) { + class_of[{{refl[i].h, refl[i].k, refl[i].l}}] = data.free_class[i]; + if (refl[i].d >= d_hi && refl[i].d <= d_lo) + class_counts[data.free_class[i]]++; + } + } if (d_min_limit > 0) std::erase_if(refl, [&](const MergedReflection &r) { return r.d < d_min_limit; }); if (!data.amplitudes) { @@ -65,7 +206,17 @@ RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &m if (!st.models.empty() && st.models[0].chains.size() != chains) logger.Info("The model's strict-NCS operators generate the rest of the asymmetric unit: {} chains -> {}", chains, st.models[0].chains.size()); - out.refinement = RefineModel(st, refl, data.cell, *sg, out.wavelength_A, logger, settings); + if (class_counts.size() > 2) { + ChooseFreeClass(refl, class_of, class_counts, d_hi, d_lo, st, data.cell, *sg, out.wavelength_A, settings, + logger, out); + out.free_reflections = 0; + for (const MergedReflection &r : refl) + out.free_reflections += r.rfree_flag; + } + // the model has most likely been refined against every reflection of a free set drawn here + RefinementSettings refine = settings; + refine.as_given_competes = !out.free_flags_generated; + out.refinement = RefineModel(st, refl, data.cell, *sg, out.wavelength_A, logger, refine); } catch (const std::exception &e) { out.refinement = RefinementResult(); out.refinement.failure_reason = std::string("refinement failed: ") + e.what(); diff --git a/rugnux/RugnuxRefine.h b/rugnux/RugnuxRefine.h index 78e2bf2af..64336d304 100644 --- a/rugnux/RugnuxRefine.h +++ b/rugnux/RugnuxRefine.h @@ -20,6 +20,11 @@ struct RefineOnlyResult { double wavelength_A = 0; // the file's, or 1 A where it states none bool french_wilson = false; // amplitudes made here from the file's intensities bool free_flags_generated = false; // the file has no free set: one was drawn here (AssignRfreeFlags) + // Where the file's free flags are classes (0-19), the class taken as the test set and what chose it (ChooseFreeClass + // in RugnuxRefine.cpp): "STATUS" the file's own mark, "REMARK3_COUNT" the class whose count the model file + // records and that the model fits worst, "WORST_FIT" the class the model fits worst. -1 / empty: not classes. + int free_class = -1; + std::string free_set_source; int reflections = 0; // amplitudes refined against, after the resolution limit int free_reflections = 0; double d_min_limit = 0; // the limit asked for (0: none) diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 9577ec095..4d031458d 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -5424,6 +5424,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // refinement takes a law at a smaller lead over its null (AdoptTwinLaw). settings.twin_data_indicated = TwinningVerdictKey(result) == "INDICATED" && result.twinning.twin_fraction_source == TwinFractionSource::LTest; + // a model placed as given on a free set drawn here has seen those reflections (MODEL_FREE_SET): its R-free + // as given is no free statistic to keep it by + settings.as_given_competes = !(result.model_placement == "RIGID_BODY" && result.free_flags_source == "DRAWN"); const PlacedModelOutput p = *placed_model_out; // a copy: a rival taken below replaces *placed_model_out // In the group RefinableIn chose: the model's own where it has the rotations of the group written - // the screws are then the model's, as its fit already took them - so a model whose screws the merge diff --git a/tests/RefineOnlyTest.cpp b/tests/RefineOnlyTest.cpp index 37a62d6a9..1377319aa 100644 --- a/tests/RefineOnlyTest.cpp +++ b/tests/RefineOnlyTest.cpp @@ -221,6 +221,106 @@ TEST_CASE("RunRefineOnly scores a model against intensities, with amplitudes and CHECK(text.find("REFINE_FRENCH_WILSON= TRUE") != std::string::npos); } +TEST_CASE("ParseRugnuxReflections keeps a structure-factor mmCIF's free-flag classes beside its status", "[refine_only]") { + // status marks class 5 free; the classes are h % 3 (three of them) + std::string cif = "data_r\n_cell.length_a 31\n_cell.length_b 37\n_cell.length_c 43\n_cell.angle_alpha 90\n" + "_cell.angle_beta 90\n_cell.angle_gamma 90\n_symmetry.space_group_name_H-M 'P 21 21 21'\n" + "loop_\n_refln.index_h\n_refln.index_k\n_refln.index_l\n_refln.status\n" + "_refln.pdbx_r_free_flag\n_refln.F_meas_au\n_refln.F_meas_sigma_au\n"; + for (int h = 1; h <= 6; h++) + for (int k = 1; k <= 3; k++) + cif += fmt::format("{} {} 1 {} {} {} 1\n", h, k, h % 3 == 2 ? 'f' : 'o', h % 3 == 2 ? 5 : h % 3, 10 * h); + const RugnuxReflections r = ParseRugnuxReflections(cif, "synthetic-sf.cif"); + REQUIRE(r.reflections.size() == 18); + REQUIRE(r.free_class.size() == 18); + CHECK(r.n_free == 6); + for (size_t i = 0; i < r.reflections.size(); i++) { + CHECK(r.reflections[i].rfree_flag == (r.reflections[i].h % 3 == 2)); + CHECK(r.free_class[i] == (r.reflections[i].h % 3 == 2 ? 5 : r.reflections[i].h % 3)); + } +} + +namespace { + +// The synthetic input with its reflections in 20 free-flag classes, class 5 marked free, and the intensities of +// class `held_out` disturbed by 30 % at random: the class a model was not refined against, which it fits worst. +// `count`: the reflections of that class. +void ClassedRefineInput(gemmi::Structure &st, RugnuxReflections &data, int held_out, int &count) { + SyntheticRefineInput(st, data); + std::mt19937 rng(11); + std::uniform_real_distribution noise(0.7f, 1.3f); + data.free_column = "FreeR_flag"; + count = 0; + for (MergedReflection &m : data.reflections) { + const int c = (m.h * 7 + m.k * 13 + m.l * 17 + 1000) % 20; + data.free_class.push_back(c); + m.rfree_flag = c == 5; + data.n_free += m.rfree_flag; + if (c == held_out) { + m.I *= noise(rng) * noise(rng); + count++; + } + } +} + +} // namespace + +TEST_CASE("RunRefineOnly takes the free-flag class the model fits worst, by the model file's test-set count", + "[refine_only]") { + gemmi::Structure st; + RugnuxReflections data; + int count = 0; + ClassedRefineInput(st, data, 3, count); + st.meta.refinement.emplace_back(); + st.meta.refinement[0].rfree_set_count = count; + RefinementSettings settings; + settings.score_only = true; + settings.gpu = false; + settings.nthreads = 2; + Logger logger("RefineOnlyTest"); + const RefineOnlyResult r = RunRefineOnly(data, st, 0, settings, logger); + REQUIRE(r.refinement.ok); + CHECK(r.free_class == 3); + CHECK(r.free_set_source == "REMARK3_COUNT"); + CHECK(r.free_reflections == count); + const std::string text = RenderReportText(BuildRefineReportDocument("out", "x-sf.cif", "model.cif", r), true); + CHECK(text.find("REFINE_FREE_SET_SOURCE= REMARK3_COUNT") != std::string::npos); + CHECK(text.find("REFINE_FREE_SET_CLASS= 3") != std::string::npos); +} + +TEST_CASE("RunRefineOnly takes the worst-fit free-flag class where the model file gives no count", "[refine_only]") { + gemmi::Structure st; + RugnuxReflections data; + int count = 0; + ClassedRefineInput(st, data, 3, count); + RefinementSettings settings; + settings.score_only = true; + settings.gpu = false; + settings.nthreads = 2; + Logger logger("RefineOnlyTest"); + const RefineOnlyResult r = RunRefineOnly(data, st, 0, settings, logger); + REQUIRE(r.refinement.ok); + CHECK(r.free_class == 3); + CHECK(r.free_set_source == "WORST_FIT"); +} + +TEST_CASE("RunRefineOnly keeps the marked free-flag class where the model fits it worst", "[refine_only]") { + gemmi::Structure st; + RugnuxReflections data; + int count = 0; + ClassedRefineInput(st, data, 5, count); + RefinementSettings settings; + settings.score_only = true; + settings.gpu = false; + settings.nthreads = 2; + Logger logger("RefineOnlyTest"); + const RefineOnlyResult r = RunRefineOnly(data, st, 0, settings, logger); + REQUIRE(r.refinement.ok); + CHECK(r.free_class == 5); + CHECK(r.free_set_source == "STATUS"); + CHECK(r.free_reflections == count); +} + TEST_CASE("RunRefineOnly scores a model with strict-NCS operators as the expanded model", "[refine_only]") { gemmi::Structure st; RugnuxReflections data; diff --git a/tools/battery/model_check.py b/tools/battery/model_check.py index 686ae4eb2..58efc4444 100644 --- a/tools/battery/model_check.py +++ b/tools/battery/model_check.py @@ -454,12 +454,59 @@ def french_wilson(workdir, hkl, i, sigi, flags, cell, sg): return out.make_miller_array(), col("F"), col("SIGF"), col("FreeR_flag").astype(int) +def free_class(hkl, f, cls, marked, xyz, cell, sg): + """The test set of deposited structure factors whose free flags are classes (0-19), as rugnux + --mode refine takes it (RugnuxRefine.cpp, ChooseFreeClass): the class the file marks (status 'f', + else flag 0) is not always the one the depositors held out. The deposited model fits the class it + was not refined against worst, and the test-set count its file records is that class's count, + both within the refinement's resolution range. The mark is kept where the model fits it worst or the + count is nearest it (STATUS); else the worst-fit class is taken where the count is nearest it + (REMARK3_COUNT) or where it is above the marked one by more than twice the scatter of the other + classes' R (WORST_FIT). The fit here: |Fcalc| without bulk solvent, an overall scale and B, on + the reflections finer than 6 A (the finer half on data coarser than that). Returns (class, source).""" + st = read_model(xyz) + ref = st.meta.refinement[0] if st.meta.refinement else None + d = cell.calculate_d_array(hkl) + d_hi = ref.resolution_high if ref is not None and ref.resolution_high > 0 else 0 + d_lo = ref.resolution_low if ref is not None and ref.resolution_low > 0 else np.inf + in_range = (d >= d_hi - 1e-3) & (d <= d_lo + 1e-3) + classes = np.unique(cls[cls >= 0]) + counts = {c: int(np.sum(in_range & (cls == c))) for c in classes} + n_test = ref.rfree_set_count if ref is not None else -1 + nearest = [] + if n_test > 0: + dist = {c: abs(n - n_test) for c, n in counts.items()} + nearest = [c for c in classes if dist[c] == min(dist.values())] + # the model as given on each class + ghkl, gf, _, gcls = to_group(hkl, f, np.ones(len(f)), cls, cell, sg, np.arange(len(f))) + fc_hkl, fc = fcalc(st, max(float(np.min(d)), 1.0)) + _, i, j = np.intersect1d(hkl_key(ghkl), hkl_key(fc_hkl), return_indices=True) + fo, fcm, c_of = gf[i], fc[j], gcls[i] + dg = cell.calculate_d_array(ghkl[i]) + s2 = 1.0 / dg ** 2 + sel = (s2 > min(1 / 36.0, np.median(s2))) & (fo > 0) & (fcm > 0) & (c_of >= 0) & \ + (dg >= d_hi - 1e-3) & (dg <= d_lo + 1e-3) + slope, lnk = np.polyfit(s2[sel], np.log(fo[sel] / fcm[sel]), 1) + resid = np.abs(fo - np.exp(lnk + slope * s2) * fcm) + r = {c: np.sum(resid[sel & (c_of == c)]) / np.sum(fo[sel & (c_of == c)]) for c in classes} + worst = max(r, key=r.get) + scatter = np.std([x for c, x in r.items() if c != worst], ddof=1) + if worst == marked or marked in nearest: + return marked, "STATUS" + if worst in nearest: + return worst, "REMARK3_COUNT" + if r.get(marked, np.inf) < r[worst] - 2 * scatter: + return worst, "WORST_FIT" + return marked, "STATUS" + + def depositor_data(pdb_id, cache_dir, workdir, cell_m, sg_m): """The first reflection block of the entry's -sf.cif as an MTZ (F, SIGF, FreeR_flag) in the model's cell and group, and what the amplitudes were made from. Amplitudes: F_meas, else the mean of F(+)/F(-); where only intensities were deposited (I, else the mean of I(+)/I(-)) they - are converted by French-Wilson. Free flags from _refln.status ('f' free), else from a numeric - _refln.pdbx_r_free_flag (0 when it takes more than two values - the CCP4 convention - else the + are converted by French-Wilson. Where a numeric _refln.pdbx_r_free_flag takes more than two + values, the free set is the class free_class decides; else free flags from _refln.status ('f' + free), else from a numeric _refln.pdbx_r_free_flag (0 when it takes more than two values - the CCP4 convention - else the minority value). Where neither marks a free set (e.g. a SHELXL entry whose free set was not deposited), a synthetic one is made by synthetic_free and noted. Returns None, reason when there is nothing usable, else the data and a note (None, or what was synthesised).""" @@ -483,7 +530,19 @@ def depositor_data(pdb_id, cache_dir, workdir, cell_m, sg_m): if "status" in labels else np.array([]) note = None v = b.make_int_array("pdbx_r_free_flag", -1) if "pdbx_r_free_flag" in labels else None - if len(status) == len(hkl) and np.any(status == "f"): + classes = np.unique(v[v >= 0]) if v is not None else np.array([]) + if len(classes) > 2 and (len(status) != len(hkl) or np.any(status == "f")): + # free flags in classes: the test set by evidence (free_class) + marked = np.bincount(v[(status == "f") & (v >= 0)]).argmax() if len(status) == len(hkl) else 0 + use = (v >= 0) & ((status == "f") | (status == "o") if len(status) == len(hkl) else True) + ok = use & np.isfinite(f) & (f > 0) + amp = f[ok] if amplitudes else np.sqrt(f[ok]) + free, source = free_class(hkl[ok], amp, v[ok], marked, + deposition(pdb_id, cache_dir)["xyz"], cell_m, sg_m) + flags = np.where(v == free, 0, 1) + if free != marked: + note = f"free set: flag class {free} ({source}), not the class {marked} the file marks" + elif len(status) == len(hkl) and np.any(status == "f"): use = (status == "f") | (status == "o") flags = np.where(status == "f", 0, 1) elif v is not None and len(np.unique(v[v >= 0])) >= 2: @@ -493,14 +552,17 @@ def depositor_data(pdb_id, cache_dir, workdir, cell_m, sg_m): flags = np.where(v == free, 0, 1) else: use = np.ones(len(hkl), bool) if len(status) != len(hkl) else (status == "o") - flags = np.ones(len(hkl), int) + flags = None note = f"deposited structure factors carry no free set: a synthetic {SYNTHETIC_FREE:.0%} one is used" use &= np.isfinite(f) & np.isfinite(s) & (s > 0) if not np.allclose(b.cell.parameters, cell_m.parameters, rtol=0.02, atol=0.5): return None, "deposited structure factors in another cell" + synthetic = flags is None + if synthetic: + flags = np.ones(len(hkl), int) hkl, f, s, flags = to_group(hkl[use], f[use], s[use], flags[use], cell_m, sg_m, np.arange(use.sum())) - if note: + if synthetic: flags = synthetic_free(hkl) if not amplitudes: hkl, f, s, flags = french_wilson(workdir, hkl, f, s, flags, cell_m, sg_m) diff --git a/writer/RugnuxReflectionFile.cpp b/writer/RugnuxReflectionFile.cpp index 5fd297502..a6f8a1aaf 100644 --- a/writer/RugnuxReflectionFile.cpp +++ b/writer/RugnuxReflectionFile.cpp @@ -107,6 +107,8 @@ std::vector SfCifSpec(const gemmi::cif::Loop &loop) { continue; spec.push_back(l); } + // the numeric classes as well, whichever marks the free set (RugnuxReflections::free_class) + spec.emplace_back("pdbx_r_free_flag FreeR_class I 0"); return spec; } @@ -143,6 +145,7 @@ RugnuxReflections RugnuxReflectionsFromMtz(const gemmi::Mtz &mtz, const std::str } const gemmi::Mtz::Column *free = FindFree(mtz); + const gemmi::Mtz::Column *free_class = mtz.column_with_label("FreeR_class", nullptr, 'I'); std::vector raw_free; const size_t stride = mtz.columns.size(); for (size_t r = 0; r < static_cast(mtz.nreflections); r++) { @@ -184,6 +187,10 @@ RugnuxReflections RugnuxReflectionsFromMtz(const gemmi::Mtz &mtz, const std::str const float f = row[free->idx]; raw_free.push_back(std::isfinite(f) ? static_cast(std::lround(f)) : -1); } + if (free_class != nullptr) { + const float f = row[free_class->idx]; + out.free_class.push_back(std::isfinite(f) ? static_cast(std::lround(f)) : -1); + } } if (out.reflections.empty()) throw std::runtime_error(name + " holds no reflection with a value in " + out.value_column); diff --git a/writer/RugnuxReflectionFile.h b/writer/RugnuxReflectionFile.h index 08a2de429..7d080db88 100644 --- a/writer/RugnuxReflectionFile.h +++ b/writer/RugnuxReflectionFile.h @@ -29,6 +29,10 @@ struct RugnuxReflections { std::string value_column, sigma_column; std::string free_column; // empty: the file has no free-set flags int n_free = 0; + // A structure-factor mmCIF's numeric free-set classes (_refln.pdbx_r_free_flag) per reflection, -1 where it + // gives none; empty for an MTZ or where the file has no such column. free_column marks one class of them + // (_refln.status 'f', or the flag by the convention above), which need not be the depositors' test set. + std::vector free_class; }; // The columns of `mtz` (an MTZ, or an mmCIF converted to one), found as RugnuxReflections describes: