diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 0138440fb..84bbb1d45 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; @@ -770,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); @@ -794,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. 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