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 e3b9ceb6f..4351c3243 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1300,4 +1300,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/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 1a85a59ba..73171fbb9 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -3393,6 +3393,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..eb75ba126 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,6 +206,13 @@ 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()); + 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; + } out.refinement = RefineModel(st, refl, data.cell, *sg, out.wavelength_A, logger, settings); } catch (const std::exception &e) { out.refinement = RefinementResult(); 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/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/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: