From ce87b98280bab5bcef1c376c3f01721a8de18702 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 27 Sep 2026 21:11:14 +0200 Subject: [PATCH 1/3] ModelValidation: run the null's replicates on threads of their own, so their inner passes go parallel The nine null replicates ran under ParallelFor, i.e. on pool workers, and a parallel pass reached from a pool worker runs inline. So the scale fit (FitModelScale's grid) and the rigid-body Jacobian of every replicate but the one on the calling thread ran serially, and the null lasted as long as its slowest serial replicate. Sampled from /proc every 0.5 s through the null of a 1.63 A P2 set (36793 atoms, 4 indexing candidates): CPU use sat at exactly 8.0 cores of 32 for most of its 220 s, falling to 1-4 at the tail; on every set measured the null averaged 6-8 busy cores. The replicates now run on std::async threads outside the pool, so their passes queue on the pool as the real fit's do. Exact: every pass splits its work the same way wherever it runs (ParallelChunks is split on the thread count, not on the caller), so p.hkl, p.mtz, p.cif, p_maps.mtz, p_model.cif are byte-identical and the model-validation section of p_report.txt unchanged on 9hnc, 8sa8, 8xtg, 9fhc, 5epe and 6oel (battery commands, --model). (6oel's report also differs in SUPERCELL_DOUBLED_CELL, 106.695/119.998 vs the equivalent 73.305/60.002 - that line flips between runs of the same code already, and 6oel never reaches the null.) Null wall clock / mean busy cores (shared 16-core/32-thread box with other jobs, so indicative): 9hnc 220 s / 7.8 -> 152 s / 15.9 8sa8 167 s / 8.2 -> 60 s / 23.7 8xtg 88 s / 7.3 -> 46 s / 20.3 9fhc 77 s / 6.7 -> 57 s / 17.7 5epe 45 s / 6.1 -> 36 s / 17.0 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- rugnux/ModelValidation.cpp | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 0138440fb..e90287ff6 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -8,6 +8,8 @@ #include #include #include +#include +#include #include #include #include @@ -29,7 +31,6 @@ #include "../common/CorrelationCoefficient.h" #include "../common/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it) #include "../common/Logger.h" -#include "../common/ParallelFor.h" // ParallelFor #include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections #include "../image_analysis/scale_merge/CrystalSetting.h" // CellMappingOperators #include "ModelFFT.h" @@ -749,7 +750,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // replicates share nothing else, so they run concurrently and the real model - which the // maps and the atom-density readout below are still taken from - is never moved. std::vector null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES); - ParallelFor(NULL_REPLICATES, nthreads, [&](int i) { + auto run_replicate = [&](int i) { ModelState rep; rep.st = mdl.st; SetModelPositions(rep.st.models[0], as_read); @@ -762,7 +763,22 @@ ModelValidationResult ValidateAgainstModel(const std::vector & identity_fit = std::move(probe.identity_fit); // replicates are placed against obs } null_r_work[i] = place_and_fit(rep, obs, std::move(identity_fit)).fit.r_work; - }); + }; + // On threads of their own, not on ParallelFor's pool. A parallel pass reached from a pool + // worker runs inline (ParallelFor.h), so there the scale fits and the rigid-body Jacobian of + // every replicate but the one on the calling thread ran serially, and the null lasted as long + // as its slowest serial replicate. From threads outside the pool those passes spread over the + // pool as the real fit's do. Each pass splits its work the same way wherever it runs, so the + // numbers are the same. + std::atomic next_replicate{0}; + std::vector> replicate_threads; + for (size_t t = 0; t < std::min(std::max(nthreads, 1), NULL_REPLICATES); t++) + replicate_threads.push_back(std::async(std::launch::async, [&] { + for (int i = next_replicate++; i < NULL_REPLICATES; i = next_replicate++) + run_replicate(i); + })); + for (std::future &f : replicate_threads) + f.get(); const auto [null_mean, null_sd] = mean_sd(null_r_work); result.fit_tested = true; result.null_replicates = NULL_REPLICATES; From dff1543bd758520f3cb854a39d82f0d896bf2cc9 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 27 Sep 2026 21:11:14 +0200 Subject: [PATCH 2/3] RigidBodyRefine: reuse the residuals of an accepted step when Ceres asks for its Jacobian Ceres evaluates some points twice - the residuals alone to decide whether to take a step, then residuals and Jacobian together at the same point - and every evaluation is a full Fcalc, a solvent mask and a scale fit. The residuals are a function of the placement alone (the zone's bulk solvent is fixed by its first evaluation), so the second computation is replaced by a copy of the first. On the null of a 1.63 A P2 set the replicates' evaluations drop 46/64/72/89/136 -> 44/59/68/82/123, and each saved evaluation is a serial one on the replicate's critical path. With the previous commit the null there goes 152 s -> 87 s (27 busy cores); 8sa8 60 -> 53 s. Exact: p.hkl, p.mtz, p.cif, p_maps.mtz, p_model.cif byte-identical to rc173 and the model-validation report unchanged on 9hnc, 8sa8, 8xtg, 9fhc, 5epe, 6oel (battery commands, --model); rigid-body rotations and translations in the log identical. [ModelValidation] tests pass. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- rugnux/RigidBodyRefine.cpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index 8b36f44d5..e60339c57 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -214,8 +214,17 @@ public: bool Evaluate(double const *const *parameters, double *residuals, double **jacobians) const override { const int n = num_residuals(); - if (!ev_.Residuals(parameters[0], residuals)) - return false; + // Ceres evaluates every step it accepts twice: the residuals alone to decide whether to take + // it, then residuals and Jacobian together at the same point. The residuals are a function of + // the placement alone, so the second time they are the first time's, copied. + if (!last_residuals_.empty() && std::equal(parameters[0], parameters[0] + 6, last_q_.begin())) { + std::copy(last_residuals_.begin(), last_residuals_.end(), residuals); + } else { + if (!ev_.Residuals(parameters[0], residuals)) + return false; + std::copy(parameters[0], parameters[0] + 6, last_q_.begin()); + last_residuals_.assign(residuals, residuals + n); + } if (jacobians != nullptr && jacobians[0] != nullptr) { const double step = ev_.JacobianStep(); std::vector> shifted(6, std::vector(n)); @@ -244,6 +253,8 @@ private: Evaluator &ev_; std::vector &columns_; size_t nthreads_; + mutable std::array last_q_{}; + mutable std::vector last_residuals_; // at last_q_; empty until an evaluation succeeds }; // The directions in which this space group's origin is free. Translating the whole cell content along From 7eaf0f21b0dc047542e1c80092759133f20847dd Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 27 Sep 2026 21:49:20 +0200 Subject: [PATCH 3/3] ModelValidation: floor the null's spread instead of zeroing the sigma below it [changes verdicts] record_fit set r_work_sigma to 0 whenever the null's sd fell under NULL_SD_FLOOR = 1e-3, and the indexing margin had the same guard at 1e-4. The floor was there so that sd == 0 and sd == 1e-7 do not land on opposite verdicts, but a zero sigma is itself a verdict: DOES NOT FIT. A large, well-measured crystal has a null that narrow for real - the replicates of a random placement agree to 0.0006-0.0010 in R-work over hundreds of thousands of reflections - so the crystal's own model, at R-work 0.17-0.29 against a null of 0.50-0.55, read "+0.00 sigma, the model DOES NOT FIT", and the decisions it gates (setting, indexing, enantiomorph label) were refused. The spread is now floored in the denominator, (mean - R) / max(sd, floor): continuous at sd -> 0, and a null with no spread because the model is too small or symmetric for a rotation to matter still reads near zero, since the real fit is then no better than its own random placements. Changed, measured on the battery commands with --model (every set in the last battery whose sigma was 0.00 for this reason; all three were REJECTED): 8sa8 +0.00 -> +328 sigma, FITS; the model's setting is adopted, R-free 0.2273 -> 0.1747 (deposited 0.1525) 9bn8 +0.00 -> +378 sigma, FITS; reindexed to the model's indexing (+137 sigma), R-free 0.5462 -> 0.1721 (deposited 0.1575) 8xtg +0.00 -> +215 sigma, FITS; the indexing it gates was already the model's ("kept"), so the reflection files are byte-identical and only the verdict changes Unchanged (files byte-identical, sigmas identical): 9hnc, 9fhc, 5epe; no margin null in the battery was under its own 1e-4 floor. [ModelValidation] tests pass. 5epe stays DOES NOT FIT for a different reason, not touched here: in F23 one of the nine random placements lies within the rigid body's reach of the twin-related orientation, walks there (held-out R-free 0.51 -> 0.18), and the null's sd becomes 0.13. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- rugnux/ModelValidation.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index e90287ff6..84bbb1d45 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -786,13 +786,14 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.null_r_work_sd = null_sd; auto record_fit = [&](const Placement &p) { - // Guarded at a floor, not at zero: sd == 0 and sd == 1e-7 are one ulp apart and land on - // opposite verdicts, and a null with no spread - a model too small or too symmetric for a - // rotation about its own centroid to move |Fcalc| - is the case that produces it. Measured - // nulls sit near 0.60 with sd ~0.009, so this is nowhere near anything genuine. + // The spread is floored, not trusted down to zero: sd == 0 and sd == 1e-7 are one ulp apart + // and would give sigmas of nothing and of millions. A null with no spread - a model too + // small or too symmetric for a rotation about its own centroid to move |Fcalc| - fits no + // better than its own random placements, so it still reads near zero here. The floor is + // not a gate: a large, well-measured crystal has a null as narrow as 0.0006-0.0010 + // (measured), and a real fit 0.3 below it is hundreds of sigma, not none. constexpr double NULL_SD_FLOOR = 1e-3; - result.r_work_sigma = - null_sd >= NULL_SD_FLOOR ? (null_mean - p.fit.r_work) / null_sd : 0.0; + result.r_work_sigma = (null_mean - p.fit.r_work) / std::max(null_sd, NULL_SD_FLOOR); result.model_fits = result.r_work_sigma >= MODEL_FIT_SIGMA; }; record_fit(real); @@ -810,8 +811,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.indexing_margin_null_mean = margin_mean; result.indexing_margin_null_sd = margin_sd; constexpr double MARGIN_SD_FLOOR = 1e-4; - result.indexing_margin_sigma = - margin_sd >= MARGIN_SD_FLOOR ? (indexing.margin - margin_mean) / margin_sd : 0.0; + result.indexing_margin_sigma = (indexing.margin - margin_mean) / std::max(margin_sd, MARGIN_SD_FLOOR); // The margin is judged against the margin a random placement produces, not against a value. // A random model also picks a winner, and on these data it picks one by a comparable lead // (measured), so the raw margin says nothing on its own.