v1.0.0-rc.177 (#87)

* jfjoch_broker: Faster live Rugnux analysis, with the GPU indexers prepared at the start of a collection and kept between collections; Rugnux no longer writes `.ccp4` maps (the map coefficients are in `<prefix>_maps.mtz`).

Reviewed-on: #87
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>

image_analysis/structure_refinement/ModelValidation.cpp

@@ -10,6 +10,8 @@
 #include <cstdio>
 #include <array>
 #include <atomic>
+#include <chrono>
+#include <mutex>
 #include <functional>
 #include <memory>
 #include <future>
@@ -1697,9 +1699,8 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     result.mean_fom = sigma_a.mean_fom;
     result.sigma_a_shells = sigma_a.shells;
 
-    gemmi::AsuData<std::complex<float>> map2fofc, mapfofc;
+    gemmi::AsuData<std::complex<float>> map2fofc;
     map2fofc.unit_cell_ = ucell; map2fofc.spacegroup_ = sg;
-    mapfofc.unit_cell_ = ucell;  mapfofc.spacegroup_ = sg;
     std::vector<float> fwt(terms.size()), delfwt(terms.size());
     for (size_t i = 0; i < terms.size(); i++) {
         const double m = sigma_a.weight[i].m, D = sigma_a.weight[i].d;
@@ -1710,7 +1711,6 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
         delfwt[i] = static_cast<float>(m * terms[i].fo - D * terms[i].fc);
         const std::complex<float> ph = std::polar(1.0f, static_cast<float>(terms[i].phi));
         map2fofc.v.push_back({terms[i].hkl, fwt[i] * ph});
-        mapfofc.v.push_back({terms[i].hkl, delfwt[i] * ph});
     }
 
     // A validation started before the one it follows had decided makes no maps until that one has,
@@ -1728,13 +1728,9 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     if (!maps_superseded) {
         maps = std::make_shared<ModelMaps>();
         // --- make the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) ---
-        // The three maps - this one, the difference map and the anomalous map - share nothing but what
-        // they read, so the other two are made on threads of their own beside this one. Nothing here
-        // reads the difference map, so it is made only where it is kept.
-        std::future<void> fofc_map = std::async(std::launch::async, [&] {
-            if (schedule.keep_map_grids)
-                maps->map_fofc = map_from_coefficients(mapfofc, sf_gpu);
-        });
+        // This map and the anomalous map share nothing but what they read, so the anomalous one is made
+        // on a thread of its own beside this one. Nothing here reads the difference map: its coefficients
+        // go to the file, the map is not made.
 
         // --- anomalous difference map, where the merge kept the Bijvoet split ---
         // Coefficients F(+) - F(-) carried on the model phase turned back by 90 degrees. Its peaks sit on
@@ -1749,8 +1745,6 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
             if (!mapanom.v.empty()) {
                 gemmi::Grid<float> grid = map_from_coefficients(mapanom, sf_gpu);
                 read_anomalous_map(grid, map_rms(grid), st, result);
-                if (schedule.keep_map_grids)
-                    maps->map_anom = std::move(grid);
             }
         });
 
@@ -1764,9 +1758,6 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
                         for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; }
             result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0;
         }
-        if (schedule.keep_map_grids)
-            maps->map_2fofc = std::move(grid2fofc);
-        fofc_map.get();
         anomalous_map.get();
     }
 
@@ -1948,6 +1939,19 @@ std::vector<size_t> OtherCardsTotalMemory(int device) {
     return out;
 }
 
+namespace {
+
+// The warm-up last started (WarmUpModelValidation); invalid before the first.
+std::mutex warm_up_m;
+std::shared_future<void> warm_up;
+
+std::shared_future<void> LastWarmUp() {
+    std::lock_guard lock(warm_up_m);
+    return warm_up;
+}
+
+} // namespace
+
 ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &merged,
                                            const UnitCell &cell,
                                            const RugnuxModel &model,
@@ -1960,6 +1964,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
                                            const ModelValidationSchedule &schedule,
                                            const gemmi::Op *twin_law,
                                            const gemmi::Structure *placed) {
+    if (const auto w = LastWarmUp(); w.valid())
+        w.wait();
 #ifdef JFJOCH_USE_CUDA
     // A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through:
     // the validation reports why it did not finish, and a validation that did not finish decides nothing,
@@ -1983,6 +1989,47 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
 #endif
 }
 
+void WarmUpModelValidation(std::shared_ptr<const RugnuxModel> model, size_t nthreads, Logger logger) {
+#ifdef JFJOCH_USE_CUDA
+    if (!model || get_gpu_count() == 0)
+        return;
+    std::lock_guard lock(warm_up_m);
+    std::shared_future<void> previous = warm_up;
+    warm_up = std::async(std::launch::async, [model, nthreads, logger, previous]() mutable {
+        if (previous.valid())
+            previous.wait();
+        const auto start = std::chrono::steady_clock::now();
+        // Held for as long as it runs: work beside a run waits for it where the card is too small for both.
+        GPUWorkBeside beside;
+        Logger quiet = Logger::Buffered();
+        const double d_min = model->structure.resolution > 0.0 ? model->structure.resolution : 2.0;
+        std::vector<MergedReflection> fake = ModelReferenceIntensities(*model, std::nullopt, nullptr, d_min, quiet);
+        for (size_t i = 0; i < fake.size(); i++) {
+            fake[i].F = std::sqrt(fake[i].I);
+            fake[i].rfree_flag = i % 20 == 0;
+        }
+        const gemmi::UnitCell &c = model->structure.cell;
+        const UnitCell cell{static_cast<float>(c.a), static_cast<float>(c.b), static_cast<float>(c.c),
+                            static_cast<float>(c.alpha), static_cast<float>(c.beta), static_cast<float>(c.gamma)};
+        std::string failure;
+        try {
+            failure = Validate(fake, cell, *model, quiet, nullptr, true, nthreads, 0.0, {}, {}, nullptr, true,
+                               nullptr).failure_reason;
+        } catch (const std::exception &e) {
+            cuda_clear_error();
+            failure = e.what();
+        }
+        logger.Info("Model validation warm-up on the GPU: {:.2f} s{}",
+                    std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count(),
+                    failure.empty() ? "" : " (did not finish: " + failure + ")");
+    }).share();
+#else
+    (void) model;
+    (void) nthreads;
+    (void) logger;
+#endif
+}
+
 namespace {
 
 // The reindexing operator as it reads on Miller indices ("k,h,-l" rather than "y,x,-z").