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 AdpHypotheses(const Context &ctx, const RefinementTopology &topo, const std::vector &input) { constexpr size_t MIN_TLS_ATOMS = 50; @@ -1912,7 +1913,7 @@ std::vector AdpHypotheses(const Context &ctx, const RefinementTop const std::vector pos = ctx.Positions(); const std::vector b = ctx.B(); size_t last = 0; - for (int k : {1, 2, 4, 8, 16}) { + for (int k : {1, 4, 8, 16}) { std::vector> groups; for (const std::vector> &residues : chains) { std::vector starts = TlsSegments(residues, pos, b, k, MIN_TLS_RESIDUES); @@ -2384,6 +2385,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector= 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 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= 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 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 trials; std::vector 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(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(trials.size()); + else { if (coarsest < 0) coarsest = static_cast(trials.size()); finest = static_cast(trials.size()); @@ -2866,6 +2896,10 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector 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 &these) { std::vector t; @@ -2875,27 +2909,38 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector 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 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 needed(trials.size(), 0); + for (int k = 0; k < static_cast(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++) {