RESULT-CHANGING: --mode refine: the test set of a structure-factor mmCIF whose free flags are classes, decided by the model as given

A deposited structure-factor mmCIF marks the free set with _refln.status 'f', and the reader took that mark over
_refln.pdbx_r_free_flag. Where the flags are classes (0-19), the status can mark a class the depositors' refinement
did not hold out: the deposited model then fits the marked class like working data, R-work comes out above R-free,
and refinement seems to make the model worse.

The reader now keeps the numeric classes beside the status (RugnuxReflections::free_class, an extra FreeR_class
column in the conversion). Where there are more than two classes, RunRefineOnly scores the model as given on each
class (score-only, untwinned, no solvent build), within the resolution range the model file records its refinement
to, and counts each class's reflections in that range:

- STATUS: the marked class is the one the model fits worst, or its count is the nearest to the model file's test-set
  count (REMARK 3 FREE R VALUE TEST SET COUNT / _refine.ls_number_reflns_R_free);
- REMARK3_COUNT: else the worst-fit class, where its count is the nearest;
- WORST_FIT: else the worst-fit class, where it is above the marked one by more than twice the scatter of the other
  classes' R (the marked class is fitted like a working one); otherwise the mark stands.

The count is corroboration only: over the battery's PDB cache (106 files with 20 classes) the marked class is within
a few reflections of the count on nearly all, and another class has the count exactly on about 7 of them (20 classes
of about the same size), so the count alone would move correct marks. The choice is logged with R per class and
reported as REFINE_FREE_SET_SOURCE / REFINE_FREE_SET_CLASS.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

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 - `<prefix>_refined_model.cif`/`.pdb`, `<prefix>_refined_maps.mtz`, `<prefix>_blobs.json` - with
 the `REFINE_` keys in `<prefix>_report.txt`, beside `REFINE_FRENCH_WILSON=` and `REFINE_FREE_FLAGS= FILE|GENERATED`
 saying what had to be made of the input.

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).

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);

rugnux/RugnuxRefine.cpp

@@ -3,13 +3,134 @@
 
 #include "RugnuxRefine.h"
 
+#include <algorithm>
+#include <array>
 #include <cmath>
+#include <map>
+
+#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<MergedReflection> &refl, const std::map<gemmi::Miller, int> &class_of,
+                     const std::map<int, int> &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<int> 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<int, int> 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<gemmi::Miller, float> 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<int, std::array<double, 2>> 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<int, double> 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<int> 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<double>(r.size() - 1);
+    for (const auto &[c, v] : r)
+        if (c != worst)
+            var += (v - mean) * (v - mean) / static_cast<double>(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<MergedReflection> &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<gemmi::Miller, int> class_of;
+    std::map<int, int> 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();

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)

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<float> 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;

writer/RugnuxReflectionFile.cpp

@@ -107,6 +107,8 @@ std::vector<std::string> 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<int> raw_free;
     const size_t stride = mtz.columns.size();
     for (size_t r = 0; r < static_cast<size_t>(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<int>(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<int>(std::lround(f)) : -1);
+        }
     }
     if (out.reflections.empty())
         throw std::runtime_error(name + " holds no reflection with a value in " + out.value_column);

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<int> free_class;
 };
 
 // The columns of `mtz` (an MTZ, or an mmCIF converted to one), found as RugnuxReflections describes: