diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 91e575e14..ec42523be 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -8,6 +8,7 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxCalibration.cpp RugnuxPreScan.cpp RugnuxPasses.cpp + RugnuxImagePass.cpp RugnuxDefaults.cpp RugnuxDefaults.h Rugnux.h diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 9fdfe60dc..546231262 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1112,6 +1112,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto start_time = std::chrono::steady_clock::now(); + // The geometry pre-pass integrates only the first prepass_fraction of the sweep (pass 2 reads it all). + const int prepass_end = start_image + + static_cast(std::lround((end_image - start_image) * config_.prepass_fraction)); + const int prepass_images = (prepass_end - start_image + config_.stride - 1) / config_.stride; + const int loop_end = geometry_prepass ? prepass_end : end_image; + // Declared here rather than where they are first filled, so that the phases below can share them // through PipelineLocals: the bytes the image loop read, the end message (also written to the file), // and the unmerged MTZ when it is built beside the P1 cross-check merge (see there) with the @@ -1123,10 +1129,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::future unmerged_mtz; PipelineLocals p{ .logger = logger, .result = result, .observer = observer, - .geometry_prepass = geometry_prepass, .calibration = calibration, + .geometry_prepass = geometry_prepass, .full = full, .calibration = calibration, .calibration_spots = calibration_spots, .per_image_analysis = per_image_analysis, .write_files = write_files, .start_image = start_image, .end_image = end_image, - .images_to_process = images_to_process, .dataset = dataset, .mapping = mapping, .plots = plots, + .images_to_process = images_to_process, .prepass_images = prepass_images, .loop_end = loop_end, + .dataset = dataset, .mapping = mapping, .plots = plots, .writer = writer, .writer_queue = writer_queue, .indexer_pool = indexer_pool, .indexer = indexer, .indexing_bound_A = indexing_bound_A, .census_future = census_future, .superseded = superseded, .ambiguity_reference = ambiguity_reference, .start_time = start_time, @@ -3325,481 +3332,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b indexer->ProbeSupercell(frames); } - // Main per-image loop, spread over N worker threads pulling from a shared counter. HDF5 reads - // are serialized by the global hdf5_mutex; the analysis runs in parallel. - std::atomic next_ordinal = 0; - std::atomic finished_count = 0; - // Images are written at their own ordinal, so a frame that fails leaves a HOLE rather than - // shifting everything after it up. The extent of the written datasets is therefore the highest - // ordinal reached, not the number that succeeded - counting successes makes the file one image - // short for every failure, and the image it loses is the LAST one, not the one that failed. - std::atomic max_written_ordinal = -1; - auto note_written = [&max_written_ordinal](int ordinal) { - int prev = max_written_ordinal.load(std::memory_order_relaxed); - while (prev < ordinal - && !max_written_ordinal.compare_exchange_weak(prev, ordinal, std::memory_order_relaxed)) { - } - }; - - // What the integrator did to the reflections it was handed, summed over the workers as each one - // finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is - // done once per worker at the end rather than per image. - std::atomic bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0, - bragg_fallback = 0; - - // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and - // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of - // a fixed budget instead of everything it found - a powder ring is over-determined either way. - constexpr size_t CALIBRATION_MAX_SPOTS = 20000; - const size_t calibration_spots_per_image = - std::max(1, CALIBRATION_MAX_SPOTS / static_cast(images_to_process)); - std::mutex calibration_spot_m; - std::vector calibration_spot_list; - - // Higher-order beam contamination, counted on the spots the lattice did not take. Each worker - // keeps its own tally and adds it in when it is done, so the total is a sum of counts and does - // not depend on how the images were shared out. - std::mutex harmonic_m; - HarmonicEvidence harmonic_evidence; - - // Where the indexed spots sit against their predicted reflections, measured on the geometry - // pre-pass for the canonical pass's footprint (SpotFootprint.h). - const bool measure_offsets = geometry_prepass && prescan_footprint_.has_value(); - std::mutex footprint_m; - std::vector footprint_offsets; - - auto azint_worker = [&]() { - std::vector decompression_buffer; - std::shared_ptr img; - ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_); - ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum()); - AzIntEngineCPU azint(mapping); - AzimuthalIntegrationProfile profile(mapping); - - while (!cancelled_) { - const int ordinal = next_ordinal.fetch_add(1); - const int image_idx = start_image + ordinal * config_.stride; - if (image_idx >= end_image) break; - - // Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of - // untouched pages the decoder faults in one by one - except while the writer thread - // still holds the previous one's pixels. - if (!img || img.use_count() > 1) - img = std::make_shared(); - bool read = false; - try { - read = reader_.ReadRawImage(image_idx, *img); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Error("Failed to load image {}: {}", image_idx, e.what()); - continue; - } - if (!read) continue; - - DataMessage msg{}; - msg.image = img->image; - msg.number = ordinal; - msg.original_number = image_idx; - if (dataset->efficiency.size() > image_idx) - msg.image_collection_efficiency = dataset->efficiency[image_idx]; - total_uncompressed_bytes += msg.image.GetUncompressedSize(); - - const auto t0 = std::chrono::steady_clock::now(); - try { - const uint8_t *image_ptr = msg.image.GetUncompressedPtr(decompression_buffer); - preprocessor.Analyze(buffer, image_ptr, msg.image.GetMode()); - azint.Run(buffer, profile); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Error("Error integrating image {}: {}", image_idx, e.what()); - continue; - } - msg.azint_time_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); - msg.processing_time_s = msg.azint_time_s; - msg.az_int_profile = profile.GetResult(); - msg.az_int_profile_count = profile.GetPixelCount(); - msg.az_int_profile_std = profile.GetStd(); - msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings()); - msg.ice_ring_score = profile.GetIceRingScore(mapping.Settings(), - config_.spot_finding.ice_ring_width_Q_recipA); - msg.run_number = experiment_.GetRunNumber(); - msg.run_name = experiment_.GetRunName(); - - plots.Add(msg, profile); - if (writer_queue) writer_queue->Post(msg, img); - note_written(ordinal); - if (observer) observer->OnImageProcessed(msg); - const int done = finished_count.fetch_add(1) + 1; - if (observer) observer->OnProgress(done, images_to_process); - } - }; - - // The geometry pre-pass integrates only the first prepass_fraction of the sweep (pass 2 reads it all). - const int prepass_end = start_image - + static_cast(std::lround((end_image - start_image) * config_.prepass_fraction)); - const int prepass_images = (prepass_end - start_image + config_.stride - 1) / config_.stride; - const int loop_end = geometry_prepass ? prepass_end : end_image; - - auto full_worker = [&]() { - pin_gpu(); // round-robin per worker thread; must precede engine construction - // Before the analysis engine, which page-locks these bytes for its uploads and unregisters - // them when it is destroyed: it must not outlive the buffer it was given. - std::shared_ptr img; - MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer, - /*enable_fused_adaptive_gpu=*/true); - AzimuthalIntegrationProfile profile(mapping); - HarmonicEvidence worker_harmonic; - std::vector worker_offsets; - - while (!cancelled_) { - const int ordinal = next_ordinal.fetch_add(1); - const int image_idx = start_image + ordinal * config_.stride; - if (image_idx >= loop_end) break; - - // Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of - // untouched pages the decoder faults in one by one - except while the writer thread - // still holds the previous one's pixels. - if (!img || img.use_count() > 1) - img = std::make_shared(); - bool read = false; - try { - read = reader_.ReadRawImage(image_idx, *img); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Error("Failed to load image {}: {}", image_idx, e.what()); - continue; - } - if (!read) continue; - - DataMessage msg{}; - msg.image = img->image; - msg.number = ordinal; - msg.original_number = image_idx; - if (dataset->efficiency.size() > image_idx) - msg.image_collection_efficiency = dataset->efficiency[image_idx]; - total_uncompressed_bytes += msg.image.GetUncompressedSize(); - - const auto t0 = std::chrono::steady_clock::now(); - try { - analysis.Analyze(msg, profile, config_.spot_finding); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Error("Error analyzing image {}: {}", image_idx, e.what()); - continue; - } - msg.processing_time_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); - msg.run_number = experiment_.GetRunNumber(); - msg.run_name = experiment_.GetRunName(); - - if (calibration) { - std::lock_guard lock(calibration_spot_m); - for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i) - calibration_spot_list.push_back(msg.spots[i]); - } - - AddHarmonicEvidence(msg.spots, experiment_.GetDiffractionGeometry(), worker_harmonic); - if (measure_offsets) - MeasureFootprintOffsets(msg.spots, msg.reflections, experiment_.GetBeamX_pxl(), - experiment_.GetBeamY_pxl(), worker_offsets); - - plots.Add(msg, profile); - if (writer_queue) writer_queue->Post(msg, img); - note_written(ordinal); - if (observer) observer->OnImageProcessed(msg); - const int done = finished_count.fetch_add(1) + 1; - if (observer) observer->OnProgress(done, images_to_process); - } - - { - std::lock_guard lock(harmonic_m); - harmonic_evidence.Add(worker_harmonic); - } - { - std::lock_guard lock(footprint_m); - footprint_offsets.insert(footprint_offsets.end(), worker_offsets.begin(), worker_offsets.end()); - } - - const auto bragg = analysis.BraggCounts(); - bragg_predicted += bragg.predicted; - bragg_starved += bragg.bkg_starved; - bragg_starved_nb += bragg.bkg_starved_by_neighbour; - bragg_fallback += bragg.profile_fallback; - }; - - if (observer) - observer->OnPhase("Processing images"); - Step(fmt::format("{}{}, {} images", - !per_image_analysis ? "Azimuthal integration" - : calibration ? "Image pass" - : experiment_.IsRotationIndexing() ? "Integration" : "Indexing and integration", - geometry_prepass ? " (geometry pre-pass)" - : full && prepass_detector_geometry_ ? " (refined geometry)" : "", - geometry_prepass ? prepass_images : images_to_process)); - - std::function worker = per_image_analysis ? std::function(full_worker) - : std::function(azint_worker); - // How many workers the loop actually wants. Every one of them submits its own kernels to a card, - // and a card runs out of room to accept them long before it runs out of work: measured on two - // GPUs, the loop's own time falls to sixteen workers and then rises again, so forty-eight is - // slower than sixteen on a 16 Mpx set. Sixteen is where that measurement put the minimum, and a - // later one on a single card agrees - eight leaves the loop waiting on the queue rather than the - // card. The cap is per card because that is what the queue depth belongs to, and it only applies - // when -N was left at its default - an explicit -N is a deliberate instruction and is obeyed, - // which is what a previous attempt at this got wrong. - // Only the full-analysis worker drives a card; the azimuthal one preprocesses and integrates on - // the CPU and wants every thread it can have, so the cap must not reach it. - const int gpus = get_gpu_count(); - const int image_workers = (per_image_analysis && config_.nthreads_auto && gpus > 0) - ? std::min(config_.nthreads, std::max(16, 16 * gpus)) - : config_.nthreads; - if (image_workers < config_.nthreads) - logger.Info("Per-image loop: {} of {} threads ({} GPUs) - past what a card can take queue " - "work from, more workers make it slower, not faster; pass -N to override", - image_workers, config_.nthreads, gpus); - std::vector > futures; - futures.reserve(image_workers); - const auto image_loop_start = std::chrono::steady_clock::now(); - for (int i = 0; i < image_workers; ++i) - futures.push_back(std::async(std::launch::async, worker)); - for (auto &f: futures) - f.get(); - // Everything the workers posted has to be on disk before the end message is written and the file - // finalized, and a write that failed - out of space, above all - has to surface here rather than - // leave a truncated file behind. - if (writer_queue) - writer_queue->Finish(); - // Wall time of the per-image loop alone. The per-stage means below cover only what runs inside it, - // so without this there is nothing to compare them against and no way to see how much of a run is - // spent outside it - on the first-pass indexing and on scaling/merging. - result.image_loop_threads = image_workers; - result.image_loop_time_s = std::chrono::duration( - std::chrono::steady_clock::now() - image_loop_start).count(); - - result.cancelled = cancelled_; - result.images_processed = finished_count.load(); - Note(fmt::format("{} images, {:.0f} per second", result.images_processed, - result.images_processed / std::max(1e-6, result.image_loop_time_s))); - - // A pre-pass that stopped at loop_end is a sweep of the frames it read, and that is all the merge and - // the post-refinement are handed. Left in, the empty tail is a stretch of frames every per-frame - // smoothing in the merge extrapolates over: measured, the fulls' log-scale curve, carried straight - // across 900 empty frames, fell further every round until the scale there was zero, and the - // cross-validation of its smoothness then took every scale in the sweep to NaN. - if (geometry_prepass) - indexer->GetIntegrationOutcome().resize(prepass_images); - - // The footprint of the reflection rather than of the spot: the pre-scan's widths plus where the - // spots sit against their predictions (SpotFootprint.h). Like the widths it acts only where the - // reflection outgrows the r1 disk, and like them it is held back for the canonical pass. - if (measure_offsets && !cancelled_) { - const SpotFootprint fp = FootprintWithOffsets(*prescan_footprint_, footprint_offsets); - const BraggIntegrationSettings canonical = bragg_adaptive_.value_or(experiment_.GetBraggIntegrationSettings()); - bool outgrows = false; - for (size_t b = 0; b < fp.sigma_rad.size(); ++b) - outgrows |= BRAGG_FOOTPRINT_NSIGMA * std::max(fp.sigma_rad[b], fp.sigma_tan[b]) > canonical.GetR1(); - std::string table; - for (size_t b = 0; b < fp.sigma_rad.size(); ++b) - table += fmt::format(" {:.0f}:{:.1f}/{:.1f}", (b + 0.5f) * fp.bin_px, fp.sigma_rad[b], fp.sigma_tan[b]); - logger.Info("Reflection footprint: the spot widths plus the offsets of {} indexed spots from their " - "predicted reflections, sigma radial/tangential [px] by distance from the beam [px]:{}{}", - footprint_offsets.size(), table, - outgrows ? "" : fmt::format(" - every reflection fits the r1={:.1f} disk", canonical.GetR1())); - if (outgrows) { - if (!bragg_before_adaptive_) - bragg_before_adaptive_ = experiment_.GetBraggIntegrationSettings(); - BraggIntegrationSettings bis = canonical; - bis.Footprint(fp); - bragg_adaptive_ = bis; - } - } - - if (full && indexer && experiment_.IsRotationIndexing() && !cancelled_) { - const SupercellProbe probe = indexer->GetSupercellProbe(); - if (probe[0][0].n > 0 && probe[0][1].n > 0) { - std::string table; - for (int sh = 0; sh < 3; sh++) { - table += sh == 0 ? "\n 20-5 A:" : sh == 1 ? "\n 5-3 A:" : "\n 20-3 A:"; - const auto sums = [&](int c) { - SupercellProbeClass x = sh < 2 ? probe[c][sh] : probe[c][0]; - if (sh == 2) x += probe[c][1]; - return x; - }; - const SupercellProbeFit f0 = FitSupercellProbe(sums(0)); - table += fmt::format(" 000 {:.2f} a {:.1f} b {:.1f}

{:.3f} n {};", f0.mean_i_over_sigma, f0.a, - f0.b, f0.mean_p, sums(0).n); - for (int c = 1; c < 8; c++) { - const SupercellProbeFit f = FitSupercellProbe(sums(c)); - table += fmt::format(" {}{}{}: {:.2f} occ {:.2f}% rock {:.2f}+-{:.2f}% flat {:.2f}%;", - (c >> 2) & 1, (c >> 1) & 1, c & 1, f.mean_i_over_sigma, - 100.0 * f.mean_i / f0.mean_i, 100.0 * f.b / f0.b, 100.0 * f.b_se / f0.b, - 100.0 * f.a / f0.a); - } - } - logger.Info("Supercell probe (the eight parity classes of the doubled primitive cell; , " - "mean intensity, and the fit I = a + b p against the lattice's own):{}", table); - - // The class with the most intensity over 20-3 A goes into the report, with the cell it would - // double the lattice to: the vectors t of the primitive lattice with (h,k,l).t even - twice - // one axis the class is odd on, and each other axis plus that one where the class is odd - // on it too - Niggli-reduced. - const auto band = [&](int c) { SupercellProbeClass x = probe[c][0]; x += probe[c][1]; return x; }; - const SupercellProbeFit f0 = FitSupercellProbe(band(0)); - const auto prim = indexer->GetSupercellProbePrimitive(); - int best = 0; - for (int c = 1; c < 8; c++) - if (best == 0 || FitSupercellProbe(band(c)).mean_i > FitSupercellProbe(band(best)).mean_i) - best = c; - if (prim && f0.mean_i > 0.0 && f0.b != 0.0) { - const SupercellProbeFit f = FitSupercellProbe(band(best)); - const int p[3] = {(best >> 2) & 1, (best >> 1) & 1, best & 1}; - const Coord v[3] = {prim->Vec0(), prim->Vec1(), prim->Vec2()}; - const int i = p[0] ? 0 : p[1] ? 1 : 2; - Coord t[3]; - for (int j = 0; j < 3; j++) - t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast(p[j]); - ProcessResult::SupercellClass sc; - sc.h = p[0]; sc.k = p[1]; sc.l = p[2]; - sc.occupancy_pct = 100.0 * f.mean_i / f0.mean_i; - sc.rock_pct = 100.0 * f.b / f0.b; - sc.rock_se_pct = 100.0 * f.b_se / f0.b; - sc.i_over_sigma = f.mean_i_over_sigma; - sc.doubled_cell = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell(); - result.supercell = sc; - } - } - } - - result.harmonic = HarmonicFromEvidence(harmonic_evidence); - result.powder = powder_; - result.powder_excluded_from_indexing = !config_.spot_finding.measured_ring_q_recipA.empty(); - if (result.powder_excluded_from_indexing) - result.powder_indexing_d_min_A = config_.spot_finding.high_resolution_limit; - if (result.harmonic.detected) - logger.Warning("Higher-order beam contamination: {:.0f}% of the {} unindexed spots this run " - "could not put on the lattice are on it at one third the wavelength (the same " - "test off a harmonic finds {:.1f}%). They sit at {:.1f} A apparent resolution", - result.harmonic.fraction, result.harmonic.tested, - result.harmonic.null_fraction, result.harmonic.mean_d_A); - - // What the stencil cost, as the integrator measured it: the reflections dropped whole for a - // background ring left below six clean pixels, and the rings the NEIGHBOURS' predicted regions - // would starve where the detector alone would not - counting the tails of rocking curves the ring - // is not kept clear of, so it measures the pattern's density. The second is what a wider r1 can - // change, and it is what the two-pass guard reads. - if (bragg_predicted > 0) { - const auto frac = [&](uint64_t n) { - return 100.0 * static_cast(n) / static_cast(bragg_predicted); - }; - bkg_starved_fraction_ = static_cast(bragg_starved_nb) / static_cast(bragg_predicted); - logger.Info("Integration at r1={:.1f} r2={:.1f} r3={:.2f}: of {} predicted reflections, " - "{:.3f}% lost their background ring; the neighbouring predictions crowd {:.3f}% of the " - "rings; {:.3f}% profile-fit fallbacks", - experiment_.GetBraggIntegrationSettings().GetR1(), - experiment_.GetBraggIntegrationSettings().GetR2(), - experiment_.GetBraggIntegrationSettings().GetR3(), - bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb), - frac(bragg_fallback)); - } - - // Every image failing is a total failure, not a run that produced nothing: it used to be reported - // only as per-image log lines while the process still exited 0 with no output file. - if (!cancelled_ && images_to_process > 0 && result.images_processed == 0) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "No image could be analyzed - see the per-image errors above"); - - result.mean_processing_time = plots.GetMeanProcessingTime(); - result.indexing_rate = plots.GetIndexingRate(); - - // Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in - // plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and - // writes the PONI file; nothing here is stored in the _process.h5. - if (calibration && !cancelled_) - FitCalibration(p, calibration_spot_list); - - // End message (also written to the file). - // The EXTENT of the per-image datasets (see note_written), not the success count: those two - // differ by exactly the number of frames that failed, and the writer sizes its virtual datasets - // from this. - end_msg.max_image_number = static_cast(max_written_ordinal.load() + 1); - end_msg.images_collected_count = result.images_processed; - end_msg.images_sent_to_write_count = result.images_processed; - end_msg.end_date = time_UTC(std::chrono::system_clock::now()); - end_msg.run_number = experiment_.GetRunNumber(); - end_msg.run_name = experiment_.GetRunName(); - end_msg.bkg_estimate = plots.GetBkgEstimate(); - end_msg.spindle_blind_fraction = plots.GetSpindleBlindFraction(); - result.spot_resolution_estimate_A = plots.GetResolutionEstimate(); - end_msg.ice_ring_score = plots.GetIceRingScoreArray(); - end_msg.ice_ring_score_mean = plots.GetIceRingScore(); - end_msg.az_int_result["dataset"] = plots.GetAzIntProfile(); - end_msg.indexing_rate = result.indexing_rate; - - if (full && !cancelled_) { - if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) { - end_msg.rotation_lattice = rot->lattice; - // The refined geometry, not the nominal StartMessage values. Two things about these five - // numbers that are not obvious from their names. The beam centre is the PONI - the foot of - // the perpendicular from the sample - so with a non-zero tilt it is NOT where the direct - // beam lands; that point is D*tan(rot)/pixel away and is GetDirectBeam_pxl(). And the tilt - // is not written here only for the record: IndexAndRefine puts it on outcome.experiment, - // which is the geometry Bragg prediction and integration then run at. It is fitted with the - // DISTANCE FROZEN and nothing checks the covariance, so the split between PONI and tilt is - // not reproducible between passes even though their sum - the direct beam - is. Neither - // number is meaningful on its own: report or feed back the whole set, or none of it. - end_msg.refined_beam_center_x = rot->geom.GetBeamX_pxl(); - end_msg.refined_beam_center_y = rot->geom.GetBeamY_pxl(); - end_msg.refined_poni_rot1 = rot->geom.GetPoniRot1_rad(); - end_msg.refined_poni_rot2 = rot->geom.GetPoniRot2_rad(); - end_msg.refined_poni_rot3 = rot->geom.GetPoniRot3_rad(); - // rot3 is not carried: a rotation about the beam is an exact null of this experiment - // (measured, smallest eigenvalue -5e-16), so the fit never moves it and reporting it would - // suggest it had been determined. - result.refined_detector_tilt_deg = std::array{ - rot->geom.GetPoniRot1_rad() * 180.0 / PI, - rot->geom.GetPoniRot2_rad() * 180.0 / PI}; - if (rot->tilt_walk && rot->tilt_walk->refused) - result.refused_detector_tilt_deg = std::array{ - rot->tilt_walk->tilt_rad[0] * 180.0 / PI, rot->tilt_walk->tilt_rad[1] * 180.0 / PI}; - if (rot->axis) - end_msg.refined_rotation_axis = rot->axis->GetAxis(); - end_msg.rotation_lattice_type = LatticeMessage{ - .centering = rot->search_result.centering, - .niggli_class = rot->search_result.niggli_class, - .crystal_system = rot->search_result.system - }; - result.rotation_lattice_found = true; - result.rotation_lattice_type = end_msg.rotation_lattice_type; - // What the equality this class asserts is worth on a REFINED cell. The indexer refines the - // same candidate a second time with nothing held (its pseudo-symmetry guard) and hands that - // cell over whenever the constrained one won; it comes back in the primitive reduced - // basis, and LengthEqualityDeparture puts it in the basis where a and b are the pair this - // class equates. Measured here, decided in RunAllPasses, which is where there is a pass to - // compare it against. - if (rot->unconstrained && ClassImposesLengthEquality(rot->search_result.system)) { - const double split = LengthEqualityDeparture(rot->search_result, - rot->unconstrained->lattice); - result.lattice_length_split = split; - logger.Info("Rotation lattice: {}-centred {} asserts a = b; the free refinement of " - "the same candidate departs from that by {:.2f} %", - rot->search_result.centering, - gemmi::crystal_system_str(rot->search_result.system), 100.0 * split); - } - if (rot->unconstrained && (rot->search_result.system == gemmi::CrystalSystem::Monoclinic - || rot->search_result.system == gemmi::CrystalSystem::Orthorhombic)) { - const double dep = AngleEqualityDeparture_deg(rot->search_result, rot->unconstrained->lattice); - result.lattice_angle_departure_deg = dep; - logger.Info("Rotation lattice: {}-centred {} holds its angles at 90 deg; the free refinement " - "of the same candidate departs from that by {:.2f} deg", - rot->search_result.centering, - gemmi::crystal_system_str(rot->search_result.system), dep); - } - // The consensus cell below is this lattice's conventional cell, so this is the centring - // it is written in - the pair the volume comparisons in Run() need. - result.consensus_centering = rot->search_result.centering; - } - result.consensus_cell = indexer->GetConsensusUnitCell(); - end_msg.unit_cell = result.consensus_cell; - } + ImagePass(p); // Scaling and merging (full analysis only). if (full && !cancelled_ && (result.indexing_rate.value_or(0.0f) > 0 || result.rotation_lattice_found) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 8d629dfb5..815e6b540 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -845,6 +845,11 @@ class Rugnux { // geometry_prepass == true stops after post-refining the detector geometry from the first integration // and applies it (+ the refined lattice) to experiment_, so the next pass re-integrates with it. ProcessResult RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass); + // The phases of RunPipeline that live in files of their own; each reads and writes RunPipeline's + // locals through `p` (RugnuxPipeline.h). + // The per-image loop over the sweep (azimuthal integration, or the full per-image analysis), what is + // read off it once it ends, and the end message. + void ImagePass(PipelineLocals &p); // Every pass of the run, returning the canonical one. Run() wraps this to time the whole thing. ProcessResult RunAllPasses(RugnuxObserver *observer); diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp new file mode 100644 index 000000000..cf698fec3 --- /dev/null +++ b/rugnux/RugnuxImagePass.cpp @@ -0,0 +1,579 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "Rugnux.h" +#include "../image_analysis/structure_refinement/ModelValidation.h" +#include "../writer/WriteModel.h" +#include "../image_analysis/indexing/SpindleCuspLoss.h" +#include "../image_analysis/bragg_integration/SpotWidth.h" +#include "../image_analysis/bragg_integration/SpotFootprint.h" +#include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA +#include "../image_analysis/hot_pixels/HotPixels.h" +#include "../image_analysis/SensorAbsorption.h" +#include "DiagnosticOutput.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../reader/JFJochHDF5Reader.h" +#include "../common/JFJochMath.h" +#include "../common/ParallelFor.h" +#include "../common/Logger.h" +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/AzimuthalIntegrationProfile.h" +#include "../common/CUDAWrapper.h" +#include "../common/JFJochException.h" +#include "../common/time_utc.h" +#include "../writer/FileWriter.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/beam_stop/ShadowFinder.h" +#include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" +#include "../image_analysis/indexing/HarmonicContamination.h" +#include "../image_analysis/geom_refinement/BeamCenterFFT.h" +#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" +#include "../image_analysis/geom_refinement/GeometryRefiner.h" +#include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" +#include "../image_analysis/spot_finding/SpotUtils.h" +#include "gemmi/twin.hpp" +#include "../image_analysis/azint/AzIntEngineCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderGPU.h" +#endif +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/RfreeFlags.h" +#include "../image_analysis/scale_merge/RotationScaleMerge.h" +#include "../image_analysis/scale_merge/ResolutionCutoff.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../image_analysis/scale_merge/ScalingResult.h" +#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/geom_refinement/PostRefine.h" +#include "../image_analysis/lattice_search/LatticeSearch.h" +#include "../image_analysis/lattice_search/LePageLattice.h" +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" +#include "../image_analysis/scale_merge/FrenchWilson.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" +#include "../image_analysis/scale_merge/TranslationalNCS.h" +#include "../image_analysis/scale_merge/HKLKey.h" +#include "../image_analysis/scale_merge/ScaleOnTheFly.h" +#include "../image_analysis/scale_merge/StillsPartialityRefine.h" +#include "../image_analysis/WriteReflections.h" +#include "../image_analysis/bragg_integration/CalcISigma.h" +#include "../common/Definitions.h" +#include "../common/CorrelationCoefficient.h" +#include +#include +#include "RugnuxPipeline.h" + +using namespace rugnux_internal; + +void Rugnux::ImagePass(PipelineLocals &p) { + auto &logger = p.logger; + auto &result = p.result; + auto &observer = p.observer; + auto &geometry_prepass = p.geometry_prepass; + auto &full = p.full; + auto &calibration = p.calibration; + auto &per_image_analysis = p.per_image_analysis; + auto &start_image = p.start_image; + auto &end_image = p.end_image; + auto &images_to_process = p.images_to_process; + auto &dataset = p.dataset; + auto &mapping = p.mapping; + auto &plots = p.plots; + auto &writer_queue = p.writer_queue; + auto &indexer = p.indexer; + auto &total_uncompressed_bytes = p.total_uncompressed_bytes; + auto &end_msg = p.end_msg; + auto &prepass_images = p.prepass_images; + auto &loop_end = p.loop_end; + // Main per-image loop, spread over N worker threads pulling from a shared counter. HDF5 reads + // are serialized by the global hdf5_mutex; the analysis runs in parallel. + std::atomic next_ordinal = 0; + std::atomic finished_count = 0; + // Images are written at their own ordinal, so a frame that fails leaves a HOLE rather than + // shifting everything after it up. The extent of the written datasets is therefore the highest + // ordinal reached, not the number that succeeded - counting successes makes the file one image + // short for every failure, and the image it loses is the LAST one, not the one that failed. + std::atomic max_written_ordinal = -1; + auto note_written = [&max_written_ordinal](int ordinal) { + int prev = max_written_ordinal.load(std::memory_order_relaxed); + while (prev < ordinal + && !max_written_ordinal.compare_exchange_weak(prev, ordinal, std::memory_order_relaxed)) { + } + }; + + // What the integrator did to the reflections it was handed, summed over the workers as each one + // finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is + // done once per worker at the end rather than per image. + std::atomic bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0, + bragg_fallback = 0; + + // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and + // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of + // a fixed budget instead of everything it found - a powder ring is over-determined either way. + constexpr size_t CALIBRATION_MAX_SPOTS = 20000; + const size_t calibration_spots_per_image = + std::max(1, CALIBRATION_MAX_SPOTS / static_cast(images_to_process)); + std::mutex calibration_spot_m; + std::vector calibration_spot_list; + + // Higher-order beam contamination, counted on the spots the lattice did not take. Each worker + // keeps its own tally and adds it in when it is done, so the total is a sum of counts and does + // not depend on how the images were shared out. + std::mutex harmonic_m; + HarmonicEvidence harmonic_evidence; + + // Where the indexed spots sit against their predicted reflections, measured on the geometry + // pre-pass for the canonical pass's footprint (SpotFootprint.h). + const bool measure_offsets = geometry_prepass && prescan_footprint_.has_value(); + std::mutex footprint_m; + std::vector footprint_offsets; + + auto azint_worker = [&]() { + std::vector decompression_buffer; + std::shared_ptr img; + ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_); + ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum()); + AzIntEngineCPU azint(mapping); + AzimuthalIntegrationProfile profile(mapping); + + while (!cancelled_) { + const int ordinal = next_ordinal.fetch_add(1); + const int image_idx = start_image + ordinal * config_.stride; + if (image_idx >= end_image) break; + + // Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of + // untouched pages the decoder faults in one by one - except while the writer thread + // still holds the previous one's pixels. + if (!img || img.use_count() > 1) + img = std::make_shared(); + bool read = false; + try { + read = reader_.ReadRawImage(image_idx, *img); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Error("Failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (!read) continue; + + DataMessage msg{}; + msg.image = img->image; + msg.number = ordinal; + msg.original_number = image_idx; + if (dataset->efficiency.size() > image_idx) + msg.image_collection_efficiency = dataset->efficiency[image_idx]; + total_uncompressed_bytes += msg.image.GetUncompressedSize(); + + const auto t0 = std::chrono::steady_clock::now(); + try { + const uint8_t *image_ptr = msg.image.GetUncompressedPtr(decompression_buffer); + preprocessor.Analyze(buffer, image_ptr, msg.image.GetMode()); + azint.Run(buffer, profile); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Error("Error integrating image {}: {}", image_idx, e.what()); + continue; + } + msg.azint_time_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + msg.processing_time_s = msg.azint_time_s; + msg.az_int_profile = profile.GetResult(); + msg.az_int_profile_count = profile.GetPixelCount(); + msg.az_int_profile_std = profile.GetStd(); + msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings()); + msg.ice_ring_score = profile.GetIceRingScore(mapping.Settings(), + config_.spot_finding.ice_ring_width_Q_recipA); + msg.run_number = experiment_.GetRunNumber(); + msg.run_name = experiment_.GetRunName(); + + plots.Add(msg, profile); + if (writer_queue) writer_queue->Post(msg, img); + note_written(ordinal); + if (observer) observer->OnImageProcessed(msg); + const int done = finished_count.fetch_add(1) + 1; + if (observer) observer->OnProgress(done, images_to_process); + } + }; + + auto full_worker = [&]() { + pin_gpu(); // round-robin per worker thread; must precede engine construction + // Before the analysis engine, which page-locks these bytes for its uploads and unregisters + // them when it is destroyed: it must not outlive the buffer it was given. + std::shared_ptr img; + MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer, + /*enable_fused_adaptive_gpu=*/true); + AzimuthalIntegrationProfile profile(mapping); + HarmonicEvidence worker_harmonic; + std::vector worker_offsets; + + while (!cancelled_) { + const int ordinal = next_ordinal.fetch_add(1); + const int image_idx = start_image + ordinal * config_.stride; + if (image_idx >= loop_end) break; + + // Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of + // untouched pages the decoder faults in one by one - except while the writer thread + // still holds the previous one's pixels. + if (!img || img.use_count() > 1) + img = std::make_shared(); + bool read = false; + try { + read = reader_.ReadRawImage(image_idx, *img); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Error("Failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (!read) continue; + + DataMessage msg{}; + msg.image = img->image; + msg.number = ordinal; + msg.original_number = image_idx; + if (dataset->efficiency.size() > image_idx) + msg.image_collection_efficiency = dataset->efficiency[image_idx]; + total_uncompressed_bytes += msg.image.GetUncompressedSize(); + + const auto t0 = std::chrono::steady_clock::now(); + try { + analysis.Analyze(msg, profile, config_.spot_finding); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Error("Error analyzing image {}: {}", image_idx, e.what()); + continue; + } + msg.processing_time_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + msg.run_number = experiment_.GetRunNumber(); + msg.run_name = experiment_.GetRunName(); + + if (calibration) { + std::lock_guard lock(calibration_spot_m); + for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i) + calibration_spot_list.push_back(msg.spots[i]); + } + + AddHarmonicEvidence(msg.spots, experiment_.GetDiffractionGeometry(), worker_harmonic); + if (measure_offsets) + MeasureFootprintOffsets(msg.spots, msg.reflections, experiment_.GetBeamX_pxl(), + experiment_.GetBeamY_pxl(), worker_offsets); + + plots.Add(msg, profile); + if (writer_queue) writer_queue->Post(msg, img); + note_written(ordinal); + if (observer) observer->OnImageProcessed(msg); + const int done = finished_count.fetch_add(1) + 1; + if (observer) observer->OnProgress(done, images_to_process); + } + + { + std::lock_guard lock(harmonic_m); + harmonic_evidence.Add(worker_harmonic); + } + { + std::lock_guard lock(footprint_m); + footprint_offsets.insert(footprint_offsets.end(), worker_offsets.begin(), worker_offsets.end()); + } + + const auto bragg = analysis.BraggCounts(); + bragg_predicted += bragg.predicted; + bragg_starved += bragg.bkg_starved; + bragg_starved_nb += bragg.bkg_starved_by_neighbour; + bragg_fallback += bragg.profile_fallback; + }; + + if (observer) + observer->OnPhase("Processing images"); + Step(fmt::format("{}{}, {} images", + !per_image_analysis ? "Azimuthal integration" + : calibration ? "Image pass" + : experiment_.IsRotationIndexing() ? "Integration" : "Indexing and integration", + geometry_prepass ? " (geometry pre-pass)" + : full && prepass_detector_geometry_ ? " (refined geometry)" : "", + geometry_prepass ? prepass_images : images_to_process)); + + std::function worker = per_image_analysis ? std::function(full_worker) + : std::function(azint_worker); + // How many workers the loop actually wants. Every one of them submits its own kernels to a card, + // and a card runs out of room to accept them long before it runs out of work: measured on two + // GPUs, the loop's own time falls to sixteen workers and then rises again, so forty-eight is + // slower than sixteen on a 16 Mpx set. Sixteen is where that measurement put the minimum, and a + // later one on a single card agrees - eight leaves the loop waiting on the queue rather than the + // card. The cap is per card because that is what the queue depth belongs to, and it only applies + // when -N was left at its default - an explicit -N is a deliberate instruction and is obeyed, + // which is what a previous attempt at this got wrong. + // Only the full-analysis worker drives a card; the azimuthal one preprocesses and integrates on + // the CPU and wants every thread it can have, so the cap must not reach it. + const int gpus = get_gpu_count(); + const int image_workers = (per_image_analysis && config_.nthreads_auto && gpus > 0) + ? std::min(config_.nthreads, std::max(16, 16 * gpus)) + : config_.nthreads; + if (image_workers < config_.nthreads) + logger.Info("Per-image loop: {} of {} threads ({} GPUs) - past what a card can take queue " + "work from, more workers make it slower, not faster; pass -N to override", + image_workers, config_.nthreads, gpus); + std::vector > futures; + futures.reserve(image_workers); + const auto image_loop_start = std::chrono::steady_clock::now(); + for (int i = 0; i < image_workers; ++i) + futures.push_back(std::async(std::launch::async, worker)); + for (auto &f: futures) + f.get(); + // Everything the workers posted has to be on disk before the end message is written and the file + // finalized, and a write that failed - out of space, above all - has to surface here rather than + // leave a truncated file behind. + if (writer_queue) + writer_queue->Finish(); + // Wall time of the per-image loop alone. The per-stage means below cover only what runs inside it, + // so without this there is nothing to compare them against and no way to see how much of a run is + // spent outside it - on the first-pass indexing and on scaling/merging. + result.image_loop_threads = image_workers; + result.image_loop_time_s = std::chrono::duration( + std::chrono::steady_clock::now() - image_loop_start).count(); + + result.cancelled = cancelled_; + result.images_processed = finished_count.load(); + Note(fmt::format("{} images, {:.0f} per second", result.images_processed, + result.images_processed / std::max(1e-6, result.image_loop_time_s))); + + // A pre-pass that stopped at loop_end is a sweep of the frames it read, and that is all the merge and + // the post-refinement are handed. Left in, the empty tail is a stretch of frames every per-frame + // smoothing in the merge extrapolates over: measured, the fulls' log-scale curve, carried straight + // across 900 empty frames, fell further every round until the scale there was zero, and the + // cross-validation of its smoothness then took every scale in the sweep to NaN. + if (geometry_prepass) + indexer->GetIntegrationOutcome().resize(prepass_images); + + // The footprint of the reflection rather than of the spot: the pre-scan's widths plus where the + // spots sit against their predictions (SpotFootprint.h). Like the widths it acts only where the + // reflection outgrows the r1 disk, and like them it is held back for the canonical pass. + if (measure_offsets && !cancelled_) { + const SpotFootprint fp = FootprintWithOffsets(*prescan_footprint_, footprint_offsets); + const BraggIntegrationSettings canonical = bragg_adaptive_.value_or(experiment_.GetBraggIntegrationSettings()); + bool outgrows = false; + for (size_t b = 0; b < fp.sigma_rad.size(); ++b) + outgrows |= BRAGG_FOOTPRINT_NSIGMA * std::max(fp.sigma_rad[b], fp.sigma_tan[b]) > canonical.GetR1(); + std::string table; + for (size_t b = 0; b < fp.sigma_rad.size(); ++b) + table += fmt::format(" {:.0f}:{:.1f}/{:.1f}", (b + 0.5f) * fp.bin_px, fp.sigma_rad[b], fp.sigma_tan[b]); + logger.Info("Reflection footprint: the spot widths plus the offsets of {} indexed spots from their " + "predicted reflections, sigma radial/tangential [px] by distance from the beam [px]:{}{}", + footprint_offsets.size(), table, + outgrows ? "" : fmt::format(" - every reflection fits the r1={:.1f} disk", canonical.GetR1())); + if (outgrows) { + if (!bragg_before_adaptive_) + bragg_before_adaptive_ = experiment_.GetBraggIntegrationSettings(); + BraggIntegrationSettings bis = canonical; + bis.Footprint(fp); + bragg_adaptive_ = bis; + } + } + + if (full && indexer && experiment_.IsRotationIndexing() && !cancelled_) { + const SupercellProbe probe = indexer->GetSupercellProbe(); + if (probe[0][0].n > 0 && probe[0][1].n > 0) { + std::string table; + for (int sh = 0; sh < 3; sh++) { + table += sh == 0 ? "\n 20-5 A:" : sh == 1 ? "\n 5-3 A:" : "\n 20-3 A:"; + const auto sums = [&](int c) { + SupercellProbeClass x = sh < 2 ? probe[c][sh] : probe[c][0]; + if (sh == 2) x += probe[c][1]; + return x; + }; + const SupercellProbeFit f0 = FitSupercellProbe(sums(0)); + table += fmt::format(" 000 {:.2f} a {:.1f} b {:.1f}

{:.3f} n {};", f0.mean_i_over_sigma, f0.a, + f0.b, f0.mean_p, sums(0).n); + for (int c = 1; c < 8; c++) { + const SupercellProbeFit f = FitSupercellProbe(sums(c)); + table += fmt::format(" {}{}{}: {:.2f} occ {:.2f}% rock {:.2f}+-{:.2f}% flat {:.2f}%;", + (c >> 2) & 1, (c >> 1) & 1, c & 1, f.mean_i_over_sigma, + 100.0 * f.mean_i / f0.mean_i, 100.0 * f.b / f0.b, 100.0 * f.b_se / f0.b, + 100.0 * f.a / f0.a); + } + } + logger.Info("Supercell probe (the eight parity classes of the doubled primitive cell; , " + "mean intensity, and the fit I = a + b p against the lattice's own):{}", table); + + // The class with the most intensity over 20-3 A goes into the report, with the cell it would + // double the lattice to: the vectors t of the primitive lattice with (h,k,l).t even - twice + // one axis the class is odd on, and each other axis plus that one where the class is odd + // on it too - Niggli-reduced. + const auto band = [&](int c) { SupercellProbeClass x = probe[c][0]; x += probe[c][1]; return x; }; + const SupercellProbeFit f0 = FitSupercellProbe(band(0)); + const auto prim = indexer->GetSupercellProbePrimitive(); + int best = 0; + for (int c = 1; c < 8; c++) + if (best == 0 || FitSupercellProbe(band(c)).mean_i > FitSupercellProbe(band(best)).mean_i) + best = c; + if (prim && f0.mean_i > 0.0 && f0.b != 0.0) { + const SupercellProbeFit f = FitSupercellProbe(band(best)); + const int p[3] = {(best >> 2) & 1, (best >> 1) & 1, best & 1}; + const Coord v[3] = {prim->Vec0(), prim->Vec1(), prim->Vec2()}; + const int i = p[0] ? 0 : p[1] ? 1 : 2; + Coord t[3]; + for (int j = 0; j < 3; j++) + t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast(p[j]); + ProcessResult::SupercellClass sc; + sc.h = p[0]; sc.k = p[1]; sc.l = p[2]; + sc.occupancy_pct = 100.0 * f.mean_i / f0.mean_i; + sc.rock_pct = 100.0 * f.b / f0.b; + sc.rock_se_pct = 100.0 * f.b_se / f0.b; + sc.i_over_sigma = f.mean_i_over_sigma; + sc.doubled_cell = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell(); + result.supercell = sc; + } + } + } + + result.harmonic = HarmonicFromEvidence(harmonic_evidence); + result.powder = powder_; + result.powder_excluded_from_indexing = !config_.spot_finding.measured_ring_q_recipA.empty(); + if (result.powder_excluded_from_indexing) + result.powder_indexing_d_min_A = config_.spot_finding.high_resolution_limit; + if (result.harmonic.detected) + logger.Warning("Higher-order beam contamination: {:.0f}% of the {} unindexed spots this run " + "could not put on the lattice are on it at one third the wavelength (the same " + "test off a harmonic finds {:.1f}%). They sit at {:.1f} A apparent resolution", + result.harmonic.fraction, result.harmonic.tested, + result.harmonic.null_fraction, result.harmonic.mean_d_A); + + // What the stencil cost, as the integrator measured it: the reflections dropped whole for a + // background ring left below six clean pixels, and the rings the NEIGHBOURS' predicted regions + // would starve where the detector alone would not - counting the tails of rocking curves the ring + // is not kept clear of, so it measures the pattern's density. The second is what a wider r1 can + // change, and it is what the two-pass guard reads. + if (bragg_predicted > 0) { + const auto frac = [&](uint64_t n) { + return 100.0 * static_cast(n) / static_cast(bragg_predicted); + }; + bkg_starved_fraction_ = static_cast(bragg_starved_nb) / static_cast(bragg_predicted); + logger.Info("Integration at r1={:.1f} r2={:.1f} r3={:.2f}: of {} predicted reflections, " + "{:.3f}% lost their background ring; the neighbouring predictions crowd {:.3f}% of the " + "rings; {:.3f}% profile-fit fallbacks", + experiment_.GetBraggIntegrationSettings().GetR1(), + experiment_.GetBraggIntegrationSettings().GetR2(), + experiment_.GetBraggIntegrationSettings().GetR3(), + bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb), + frac(bragg_fallback)); + } + + // Every image failing is a total failure, not a run that produced nothing: it used to be reported + // only as per-image log lines while the process still exited 0 with no output file. + if (!cancelled_ && images_to_process > 0 && result.images_processed == 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No image could be analyzed - see the per-image errors above"); + + result.mean_processing_time = plots.GetMeanProcessingTime(); + result.indexing_rate = plots.GetIndexingRate(); + + // Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in + // plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and + // writes the PONI file; nothing here is stored in the _process.h5. + if (calibration && !cancelled_) + FitCalibration(p, calibration_spot_list); + + // End message (also written to the file). + // The EXTENT of the per-image datasets (see note_written), not the success count: those two + // differ by exactly the number of frames that failed, and the writer sizes its virtual datasets + // from this. + end_msg.max_image_number = static_cast(max_written_ordinal.load() + 1); + end_msg.images_collected_count = result.images_processed; + end_msg.images_sent_to_write_count = result.images_processed; + end_msg.end_date = time_UTC(std::chrono::system_clock::now()); + end_msg.run_number = experiment_.GetRunNumber(); + end_msg.run_name = experiment_.GetRunName(); + end_msg.bkg_estimate = plots.GetBkgEstimate(); + end_msg.spindle_blind_fraction = plots.GetSpindleBlindFraction(); + result.spot_resolution_estimate_A = plots.GetResolutionEstimate(); + end_msg.ice_ring_score = plots.GetIceRingScoreArray(); + end_msg.ice_ring_score_mean = plots.GetIceRingScore(); + end_msg.az_int_result["dataset"] = plots.GetAzIntProfile(); + end_msg.indexing_rate = result.indexing_rate; + + if (full && !cancelled_) { + if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) { + end_msg.rotation_lattice = rot->lattice; + // The refined geometry, not the nominal StartMessage values. Two things about these five + // numbers that are not obvious from their names. The beam centre is the PONI - the foot of + // the perpendicular from the sample - so with a non-zero tilt it is NOT where the direct + // beam lands; that point is D*tan(rot)/pixel away and is GetDirectBeam_pxl(). And the tilt + // is not written here only for the record: IndexAndRefine puts it on outcome.experiment, + // which is the geometry Bragg prediction and integration then run at. It is fitted with the + // DISTANCE FROZEN and nothing checks the covariance, so the split between PONI and tilt is + // not reproducible between passes even though their sum - the direct beam - is. Neither + // number is meaningful on its own: report or feed back the whole set, or none of it. + end_msg.refined_beam_center_x = rot->geom.GetBeamX_pxl(); + end_msg.refined_beam_center_y = rot->geom.GetBeamY_pxl(); + end_msg.refined_poni_rot1 = rot->geom.GetPoniRot1_rad(); + end_msg.refined_poni_rot2 = rot->geom.GetPoniRot2_rad(); + end_msg.refined_poni_rot3 = rot->geom.GetPoniRot3_rad(); + // rot3 is not carried: a rotation about the beam is an exact null of this experiment + // (measured, smallest eigenvalue -5e-16), so the fit never moves it and reporting it would + // suggest it had been determined. + result.refined_detector_tilt_deg = std::array{ + rot->geom.GetPoniRot1_rad() * 180.0 / PI, + rot->geom.GetPoniRot2_rad() * 180.0 / PI}; + if (rot->tilt_walk && rot->tilt_walk->refused) + result.refused_detector_tilt_deg = std::array{ + rot->tilt_walk->tilt_rad[0] * 180.0 / PI, rot->tilt_walk->tilt_rad[1] * 180.0 / PI}; + if (rot->axis) + end_msg.refined_rotation_axis = rot->axis->GetAxis(); + end_msg.rotation_lattice_type = LatticeMessage{ + .centering = rot->search_result.centering, + .niggli_class = rot->search_result.niggli_class, + .crystal_system = rot->search_result.system + }; + result.rotation_lattice_found = true; + result.rotation_lattice_type = end_msg.rotation_lattice_type; + // What the equality this class asserts is worth on a REFINED cell. The indexer refines the + // same candidate a second time with nothing held (its pseudo-symmetry guard) and hands that + // cell over whenever the constrained one won; it comes back in the primitive reduced + // basis, and LengthEqualityDeparture puts it in the basis where a and b are the pair this + // class equates. Measured here, decided in RunAllPasses, which is where there is a pass to + // compare it against. + if (rot->unconstrained && ClassImposesLengthEquality(rot->search_result.system)) { + const double split = LengthEqualityDeparture(rot->search_result, + rot->unconstrained->lattice); + result.lattice_length_split = split; + logger.Info("Rotation lattice: {}-centred {} asserts a = b; the free refinement of " + "the same candidate departs from that by {:.2f} %", + rot->search_result.centering, + gemmi::crystal_system_str(rot->search_result.system), 100.0 * split); + } + if (rot->unconstrained && (rot->search_result.system == gemmi::CrystalSystem::Monoclinic + || rot->search_result.system == gemmi::CrystalSystem::Orthorhombic)) { + const double dep = AngleEqualityDeparture_deg(rot->search_result, rot->unconstrained->lattice); + result.lattice_angle_departure_deg = dep; + logger.Info("Rotation lattice: {}-centred {} holds its angles at 90 deg; the free refinement " + "of the same candidate departs from that by {:.2f} deg", + rot->search_result.centering, + gemmi::crystal_system_str(rot->search_result.system), dep); + } + // The consensus cell below is this lattice's conventional cell, so this is the centring + // it is written in - the pair the volume comparisons in Run() need. + result.consensus_centering = rot->search_result.centering; + } + result.consensus_cell = indexer->GetConsensusUnitCell(); + end_msg.unit_cell = result.consensus_cell; + } +} diff --git a/rugnux/RugnuxPipeline.h b/rugnux/RugnuxPipeline.h index 681abec3c..ed000a002 100644 --- a/rugnux/RugnuxPipeline.h +++ b/rugnux/RugnuxPipeline.h @@ -140,6 +140,7 @@ struct PipelineLocals { ProcessResult &result; RugnuxObserver *observer; const bool geometry_prepass; + const bool full; const bool calibration; const bool calibration_spots; const bool per_image_analysis; @@ -147,6 +148,8 @@ struct PipelineLocals { const int start_image; const int end_image; const int images_to_process; + const int prepass_images; + const int loop_end; const std::shared_ptr &dataset; AzimuthalIntegrationMapping &mapping; JFJochReceiverPlots &plots;