From 73e86d5e7df3da7df69bfa7894cd293d01bef6b8 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Wed, 28 Oct 2026 08:17:28 +0100
Subject: [PATCH] RESULT-CHANGING: Refinement: fewer weight and displacement
 trials, each family tried further only where its first members ask for it

Counted on the 368 refinements of the rc181 full battery (fee8ffa) which trial of each stage was kept,
and replayed each stage's decision without the trials left out here:

- Ladder: B refined first at the middle rung only, at the other two where it comes within two ties of
  B held there. B held was kept by 281 of 364 ladders; the rule changed the rung kept on 4 (within
  0.006 R-free at the ladder) and leaves out 412 of 1236 B-refined ladder trials.
- Ladder NCS rungs: at the rung the trials without NCS chose and the next looser one only. On 198
  ladders with copies NCS was never kept at any other rung: 232 of 594 NCS trials left out, no pick
  changed.
- Trajectory: the X-ray weight halved only in the first step (from the ladder's coarse rung). No later
  step kept the halved weight (0 of 91).
- Displacements: the 2-segment TLS cut dropped (6 choices changed, none by more than 0.001); the model's
  own TLS groups tried only where one body per chain or the finest cut lowers R-free from where the
  weights ended, the cuts between those only where also the finest fits the free set better than one
  body per chain (as before), and the U at its balanced weight only where the U at a quarter of it
  lowers R-free. Replayed: 7 of 364 choices changed (max +0.0026, mean +0.00001), 2315 -> 1527
  displacement trials.

The late trials still run ahead where the slots take them all at once (small models: no extra wave);
results do not depend on the slots or cards.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
---
 .../structure_refinement/Refinement.cpp       | 115 ++++++++++++------
 1 file changed, 80 insertions(+), 35 deletions(-)

diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp
index 1584b78dc..320a2440d 100644
--- a/image_analysis/structure_refinement/Refinement.cpp
+++ b/image_analysis/structure_refinement/Refinement.cpp
@@ -1871,8 +1871,9 @@ AdpHypothesis TlsHypothesis(const Context &ctx, std::string name, std::vector<st
 
 // The displacement models the refinement chooses among once its weights are chosen, from the model as it stands: the
 // TLS bodies of the model's own groups (`input`, where it has them); one body per chain of at least MIN_TLS_ATOMS
-// refinable atoms, and the chains cut into 2, 4, 8 and 16 TLSMD segments of at least MIN_TLS_RESIDUES residues
-// (TlsSegments) - each cut a hypothesis where it differs from the one before; and a U per atom.
+// refinable atoms, and the chains cut into 4, 8 and 16 TLSMD segments of at least MIN_TLS_RESIDUES residues
+// (TlsSegments) - each cut a hypothesis where it differs from the one before; and a U per atom. (Cut into 2, the
+// chains changed the choice on 6 of 364 refinements of the battery, none by more than 0.001 in R-free.)
 std::vector<AdpHypothesis> AdpHypotheses(const Context &ctx, const RefinementTopology &topo,
                                          const std::vector<TlsBody> &input) {
     constexpr size_t MIN_TLS_ATOMS = 50;
@@ -1912,7 +1913,7 @@ std::vector<AdpHypothesis> AdpHypotheses(const Context &ctx, const RefinementTop
         const std::vector<gemmi::Position> pos = ctx.Positions();
         const std::vector<double> b = ctx.B();
         size_t last = 0;
-        for (int k : {1, 2, 4, 8, 16}) {
+        for (int k : {1, 4, 8, 16}) {
             std::vector<std::vector<int>> groups;
             for (const std::vector<std::vector<int>> &residues : chains) {
                 std::vector<int> starts = TlsSegments(residues, pos, b, k, MIN_TLS_RESIDUES);
@@ -2384,6 +2385,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
         // Trials whose need is decided by others are run with them where the slots take them all at once, and left
         // out where it is decided they were not needed - with their evaluations, which do not count: the outcome and
         // the count are the same whatever the slots.
+        auto ahead = [&](size_t trials) { return trial_slots.size() >= trials; };
         int ahead_evaluations = 0;
         auto drop_ahead = [&](TrialOutcome &o) {
             ahead_evaluations += o.evaluations;
@@ -2417,8 +2419,12 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
         // 1. The weights from the placed model: the ladder's rungs, tighter geometry first, each with the B factors
         //    held and refined at their own balanced weight; where there are copies, each rung without NCS restraints
         //    and then with them - those with them at the B factors' weight the rungs without chose, which NCS does
-        //    not change (on 39 refinements of 21 structures with copies it never did; the X-ray weight it does) - run
-        //    ahead at every B weight where the slots take them.
+        //    not change (on 39 refinements of 21 structures with copies it never did; the X-ray weight it does), and
+        //    at the rung they chose and the next looser one (on 198 refinements of the battery with copies, NCS was
+        //    never kept at any other).
+        //    The B factors are refined first at the middle rung only, and at the others where that comes within two
+        //    ties of holding them there: on 364 refinements of the battery B held was kept by 281, and the rule
+        //    changed the rung kept on 4 (all within 0.006 in R-free).
         std::vector<TrialWeights> ladder;
         for (double rung : settings.weight_ladder)
             for (double adp : settings.adp_ladder) {
@@ -2439,20 +2445,38 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
                 if (these(ladder[i]))
                     outcomes[i] = std::move(ran[m++]);
         };
-        const bool ncs_ahead = ncs_p != nullptr && trial_slots.size() >= ladder.size();
-        run_ladder([&](const TrialWeights &l) { return l.ncs == 0 || ncs_ahead; });
+        const double mid = w0 * settings.weight_ladder[settings.weight_ladder.size() / 2];
+        auto first = [&](const TrialWeights &l) { return l.ncs == 0 && (l.adp == 0 || l.w == mid); };
+        const bool all_ahead = ahead(ladder.size());
+        run_ladder([&](const TrialWeights &l) { return first(l) || all_ahead; });
+        // B refined against B held at the middle rung (where both ran)
+        const TrialOutcome *held = nullptr, *refined = nullptr;
+        for (size_t i = 0; i < ladder.size(); i++)
+            if (ladder[i].w == mid && ladder[i].ncs == 0)
+                (ladder[i].adp == 0 ? held : refined) = &outcomes[i];
+        const bool refine_b = held == nullptr || refined == nullptr ||
+                              (refined->ran && refined->trial.within_bounds &&
+                               (!held->trial.within_bounds || refined->trial.r_free < held->trial.r_free + 2 * R_FREE_TIE));
+        auto rest_b = [&](const TrialWeights &l) { return l.ncs == 0 && !first(l); };
+        if (refine_b && !all_ahead)
+            run_ladder(rest_b);
         // the best of those without NCS restraints, those with them set aside
         std::vector<TrialOutcome> with_ncs(ladder.size());
         for (size_t i = 0; i < ladder.size(); i++)
-            if (ladder[i].ncs > 0)
+            if (ladder[i].ncs > 0 || (!refine_b && rest_b(ladder[i])))
                 std::swap(with_ncs[i], outcomes[i]);
         if (const int b = best_of(outcomes); ncs_p != nullptr && b >= 0) {
-            const double adp = outcomes[b].trial.adp_weight;
+            const double adp = outcomes[b].trial.adp_weight, w = outcomes[b].trial.weight;
+            double looser = w;
+            for (double rung : settings.weight_ladder)
+                if (w0 * rung > w && (looser == w || w0 * rung < looser))
+                    looser = w0 * rung;
+            auto these = [&](const TrialWeights &l) { return l.ncs > 0 && l.adp == adp && (l.w == w || l.w == looser); };
             for (size_t i = 0; i < ladder.size(); i++)
-                if (ladder[i].ncs > 0 && ladder[i].adp == adp)
+                if (these(ladder[i]))
                     std::swap(with_ncs[i], outcomes[i]);
-            if (!ncs_ahead)
-                run_ladder([&](const TrialWeights &l) { return l.ncs > 0 && l.adp == adp; });
+            if (!all_ahead)
+                run_ladder(these);
         }
         for (TrialOutcome &o : with_ncs)
             drop_ahead(o);
@@ -2462,10 +2486,11 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             state = std::move(outcomes[k]);
 
         // 2. The trajectory: steps of a few macro-cycles from where the last one ended, each trying the weights as they
-        //    are and, in turn, the positions' X-ray weight halved and doubled (halved: the geometry loosens as the
-        //    model improves at one weight, until it leaves the bounds) and the B factors' divided and multiplied by
-        //    four; where there are copies, also the NCS restraints changed: switched on where they are off, and in
-        //    turn tightened and loosened (to off beyond NCS_SIGMA_MAX) where they are on.
+        //    are and, in turn, the positions' X-ray weight halved and doubled (halved only in the first step, from the
+        //    ladder's coarse rung: on 364 refinements of the battery no later step kept it, of 91 that tried it) and
+        //    the B factors' divided and multiplied by four; where there are copies, also the NCS restraints changed:
+        //    switched on where they are off, and in turn tightened and loosened (to off beyond NCS_SIGMA_MAX) where
+        //    they are on.
         //    A step is taken while it lowers R-free or R-work by more than the tie and its R-free stays within the tie
         //    of the best so far; the trajectory ends at its best step (on a deposited model, whose free set it was
         //    fitted to, R-work falls while R-free rises, and it ends where it began).
@@ -2484,7 +2509,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
                     trials.push_back({w, adp / 4, ncs_s});
                     trials.push_back({w, 4 * adp, ncs_s});
                 } else {
-                    trials.push_back({w / 2, adp, ncs_s});
+                    if (step == 0)
+                        trials.push_back({w / 2, adp, ncs_s});
                     if (2 * w <= w_max)
                         trials.push_back({2 * w, adp, ncs_s});
                 }
@@ -2810,10 +2836,11 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
         //    R-free (beyond the tie) and the -ln L of the free reflections are lower: the free set is the test of the
         //    parameters it adds.
         //    Following Merritt (2012) Acta Cryst. D68, 468-477 (the displacement model the data support, decided by test)
-        //    The trials that would not change the choice are left out: the TLS cuts between one body per chain and the
-        //    finest are tried only where the finest fits the free set better than one body per chain (run ahead where
-        //    the slots take them), and the U at four times its balanced weight not at all (on 49 refinements of 33
-        //    structures, no choice taken changed).
+        //    The trials that would not change the choice are left out: the model's own TLS groups are tried only where
+        //    one body per chain or the finest cut lowers R-free, the cuts between those only where also the finest
+        //    fits the free set better, the U at its balanced weight only where the U at a quarter of it lowers R-free
+        //    (run ahead where the slots take them), and the U at four times its balanced weight not at all (on 49
+        //    refinements of 33 structures, no choice taken changed).
         if (settings.adp_hypotheses && state.ran && !ctx.Cancelled()) {
             const RefinementStage adp_stage("adp_hypotheses");
             ctx.SetPositions(state.pos);
@@ -2842,19 +2869,22 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
                 ctx.SetAniso(keep);
             }
             const double ncs_s = state.trial.ncs_sigma;
-            // every trial in order of trust; those between one TLS body per chain and the finest cut run after the
-            // others where those ask for them
+            // every trial in order of trust; the model's own TLS groups and the cuts between one TLS body per chain
+            // and the finest, and the U at its balanced weight, run after the others where those ask for them
             std::vector<TrialWeights> trials;
             std::vector<const AdpHypothesis *> which;
-            int coarsest = -1, finest = -1;
+            int own = -1, coarsest = -1, finest = -1, u_quarter = -1;
             for (const AdpHypothesis &h : hyps) {
                 if (h.kind == AdpHypothesis::Kind::Anisotropic) {
+                    u_quarter = static_cast<int>(trials.size());
                     for (double f : {0.25, 1.0}) {
                         trials.push_back({w, u_w0 * f, ncs_s});
                         which.push_back(&h);
                     }
                 } else {
-                    if (h.segments > 0) {
+                    if (h.segments == 0)
+                        own = static_cast<int>(trials.size());
+                    else {
                         if (coarsest < 0)
                             coarsest = static_cast<int>(trials.size());
                         finest = static_cast<int>(trials.size());
@@ -2866,6 +2896,10 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             std::vector<char> later(trials.size(), 0);
             for (int k = coarsest + 1; k < finest; k++)
                 later[k] = 1;
+            if (own >= 0 && coarsest >= 0)
+                later[own] = 1;
+            if (u_quarter >= 0)
+                later[u_quarter + 1] = 1;
             outcomes.assign(trials.size(), TrialOutcome());
             auto run = [&](const std::vector<char> &these) {
                 std::vector<TrialWeights> t;
@@ -2875,27 +2909,38 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
                         t.push_back(trials[k]);
                         h.push_back(which[k]);
                     }
+                if (t.empty())
+                    return;
                 std::vector<TrialOutcome> ran = RunTrials(trial_slots, t, refinable, state.pos, state.b, state.aniso,
                                                           &state.fit, settings.macro_cycles + 1, settings, ncs_p, h);
                 for (size_t k = 0, m = 0; k < trials.size(); k++)
                     if (these[k])
                         outcomes[k] = std::move(ran[m++]);
             };
-            const bool later_ahead = trial_slots.size() >= trials.size();
+            const bool later_ahead = ahead(trials.size());
             std::vector<char> now(trials.size());
             for (size_t k = 0; k < trials.size(); k++)
                 now[k] = !later[k] || later_ahead;
             run(now);
-            auto fits_better = [&](int a, int b) {
-                return outcomes[a].ran && outcomes[b].ran && outcomes[a].trial.nll_free < outcomes[b].trial.nll_free;
+            // the model's own groups where one body per chain or the finest cut lowers R-free from where the weights
+            // ended, the cuts between where also the finest fits the free set better than one body per chain, the U
+            // at its balanced weight where the U at a quarter of it lowers R-free (on 364 refinements of the battery
+            // that changed the choice on 7, none by more than 0.003 in R-free, and left out 6 in 10 of these trials)
+            auto lowers = [&](int k) {
+                return k >= 0 && outcomes[k].ran && outcomes[k].trial.within_bounds &&
+                       outcomes[k].trial.r_free < state.trial.r_free;
             };
-            if (finest > coarsest + 1 && fits_better(finest, coarsest)) {
-                if (!later_ahead)
-                    run(later);
-            } else
-                for (size_t k = 0; k < trials.size(); k++)
-                    if (later[k])
-                        drop_ahead(outcomes[k]);
+            const bool tls = lowers(coarsest) || lowers(finest), u_full = lowers(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<char> needed(trials.size(), 0);
+            for (int k = 0; k < static_cast<int>(trials.size()); k++)
+                needed[k] = later[k] && (k == u_quarter + 1 ? u_full : k == own ? tls : cuts);
+            if (!later_ahead)
+                run(needed);
+            for (size_t k = 0; k < trials.size(); k++)
+                if (later[k] && !needed[k])
+                    drop_ahead(outcomes[k]);
             int best = -1;   // where the weights ended
             auto trial = [&](int k) -> const RefinementWeightTrial & { return k < 0 ? state.trial : outcomes[k].trial; };
             for (size_t k = 0; k < outcomes.size(); k++) {