Merge branch 'model-null-speed' into rc173 (model validation null: replicate threads, residual reuse, sd floor)

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C

rugnux/ModelValidation.cpp

@@ -8,6 +8,8 @@
 #include <cmath>
 #include <complex>
 #include <array>
+#include <atomic>
+#include <future>
 #include <numeric>
 #include <optional>
 #include <random>
@@ -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<MergedReflection> &
         // 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<double> 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<MergedReflection> &
                 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<int> next_replicate{0};
+        std::vector<std::future<void>> replicate_threads;
+        for (size_t t = 0; t < std::min<size_t>(std::max<size_t>(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<void> &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;
@@ -770,13 +786,14 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
         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);
@@ -794,8 +811,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
             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.

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<std::vector<double>> shifted(6, std::vector<double>(n));
@@ -244,6 +253,8 @@ private:
     Evaluator &ev_;
     std::vector<Evaluator> &columns_;
     size_t nthreads_;
+    mutable std::array<double, 6> last_q_{};
+    mutable std::vector<double> 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