diff --git a/common/ImageBuffer.cpp b/common/ImageBuffer.cpp index 98b0101f4..6d448e809 100644 --- a/common/ImageBuffer.cpp +++ b/common/ImageBuffer.cpp @@ -109,10 +109,10 @@ bool ImageBuffer::CheckIfBufferReturned(std::chrono::microseconds timeout) { return (slots_preparation == 0) && (slots_sending == 0); } -ZeroCopyReturnValue *ImageBuffer::GetImageSlot() { +ZeroCopyReturnValue *ImageBuffer::GetImageSlot(int64_t keep_free) { std::unique_lock ul(m); - if (!disable_requests && !free_slots.empty()) { + if (!disable_requests && static_cast(free_slots.size()) > keep_free) { auto i = free_slots.front(); free_slots.pop(); diff --git a/common/ImageBuffer.h b/common/ImageBuffer.h index d0a92c463..8c954d419 100644 --- a/common/ImageBuffer.h +++ b/common/ImageBuffer.h @@ -76,7 +76,8 @@ public: void StartMeasurement(size_t buffer_location_size); void StartMeasurement(const DiffractionExperiment& experiment); - ZeroCopyReturnValue *GetImageSlot(); + // Nullptr where no more than `keep_free` slots are free. + ZeroCopyReturnValue *GetImageSlot(int64_t keep_free = 0); void ReadyToSend(uint32_t location, int64_t image_number, size_t image_size, bool indexed); void ReleaseSlot(uint32_t location); diff --git a/docs/JFJOCH_BROKER.md b/docs/JFJOCH_BROKER.md index 7ea6040e2..200315425 100644 --- a/docs/JFJOCH_BROKER.md +++ b/docs/JFJOCH_BROKER.md @@ -190,10 +190,11 @@ back to the broker, so with it the route answers 409. ## Rugnux live mode -With rugnux live mode on, rotation collections from a DECTRIS detector are processed by rugnux inside -the broker while they arrive, and serial (stills) collections in the FPGA workflow are analysed on the GPU -with the geometry and mask a rugnux pre-pass decided (see [Serial collections](#serial-collections-experimental) -below). It is off by default; the configuration file switches it on with +With rugnux live mode on, collections from a DECTRIS detector are processed inside the broker while they +arrive: rotation collections by rugnux, and serial (stills) collections on the GPU with the geometry and mask a +rugnux pre-pass decided (see [Serial collections](#serial-collections-experimental) below). In the FPGA workflow +(PSI JUNGFRAU and EIGER) the mode changes nothing: collections run image by image, as with it off, and +`/rugnux/status` stays idle. It is off by default; the configuration file switches it on with `"rugnux": {"enable": true, "reserved_cpus": 4}`, where `reserved_cpus` is the number of CPU threads left to the receiver and the image pusher. @@ -251,9 +252,9 @@ estimate and background of every image in `GET /result/scan` and the master file ### Serial collections (experimental) Live processing of serial collections is **experimental**. With rugnux live mode on, a serial -(non-rotation) collection in the FPGA workflow, in the standard detector mode, is analysed as frames from a -DECTRIS detector are: the converted frames go to GPU spot finding, indexing and integration, with rugnux's -stills spot finder (adaptive threshold, single-pixel spots kept). Before any frame is analysed, the first +(non-rotation) collection from a DECTRIS detector, in the standard detector mode, is analysed frame by frame: the +frames as the stream delivers them go to GPU spot finding, indexing and integration, with rugnux's stills spot +finder (adaptive threshold, single-pixel spots kept). Before any frame is analysed, the first 500 frames are run through rugnux, in memory, as `rugnux --no-merge` would run on them. That pre-pass decides the beam centre, refines the detector distance (where the dataset settings give a unit cell) and masks the beam-stop shadow and defective pixels; every frame of the collection is then analysed with that @@ -262,23 +263,25 @@ the centre and distance used (`refined_beam_center_x/y`, `refined_detector_dista [CBOR](CBOR.md)). The broker's own distance is taken as the beamline's: unlike offline rugnux, the pre-pass does not search the whole distance range. It refines the centre and distance for a detector square to the beam, as the frames are then analysed. At a serial hit rate the distance needs the 500 frames to be refined at -all: on fewer, too few index for the geometry refinement to run, and the stated distance is kept. The image -buffer has to hold them, and the frames that arrive while the pre-pass runs, in half its slots (for a JUNGFRAU -8M, a slot is the size of an uncompressed frame); a smaller buffer starts the pre-pass on as many as it holds. +all: on fewer, too few index for the geometry refinement to run, and the stated distance is kept. The pre-pass +holds and fits on at most half the slots of the image buffer (a slot is the size of an uncompressed +frame), so a smaller buffer starts it on as many as that. -Even 500 frames give the pre-pass only a few dozen indexed ones at a serial hit rate, and none where the -first frames are empty (a fixed-target chip whose crystals sit further in). So, where a cell is given and the -pre-pass refined the geometry on fewer than 200 frames, the centre and distance are refined again on the frames -as they are analysed, from their spots and lattices, without holding their images: once 200 strong frames (20 -indexed spots or more) are in, or 5000 indexed, the geometry is fitted, as rugnux fits it on its sample, and the -analysis switches to it; a second round does the same on frames analysed at that geometry. Two fits in all, a -fixed cost. The frames analysed before the last switch are written as analysed but left out of the merge (their -reflections were predicted at the earlier geometry, and the images are no longer there to integrate again); -the switches are notes of the live status, and the end message carries the geometry the run ended at. +Where a unit cell is given and those frames are too few for the geometry fit (200 indexed frames), frames go on +being held and the pre-pass runs again on all of them each time they have doubled, until the fit has its frames, half +the image buffer is held, or the collection ends. The analysis starts when the pre-pass ends and catches up from the +frames waiting, so it starts later on a low hit rate. -Frames arriving during the pre-pass, or faster than the analysis, wait in the image buffer, in at most half of -its slots; a frame that finds no free slot is dropped and counted as for a full image buffer, so analysis -slower than the frame rate costs frames, not memory. A pre-pass that fails leaves the stated geometry and +The geometry is fixed when the pre-pass ends and does not change during the collection, so every frame is analysed +and merged at it and the radial profile has one binning throughout. Where the pre-pass could not refine it (too few +of its frames indexed, as with a small image buffer, a low hit rate or a chip whose first frames are empty, or no unit +cell given), the distance in the header is used for the whole run: this is a warning in the log and a note of the +live status, and a wrong header distance then costs resolution and indexed frames. + +Frames arriving during the pre-pass, or faster than the analysis, wait in the image buffer, in all its slots but +one per analysis thread, kept for sending; a frame that finds it full is dropped and counted as for a full image +buffer, so analysis slower than the frame rate costs frames, not memory. Where the writer lags behind, the analysis +waits for slots to come back. A pre-pass that fails leaves the stated geometry and mask in use. Where the start message cannot be sent after the pre-pass (no writer listening) or the analysis cannot be set up, the collection ends with that error, reported by `/status` and `/wait_till_done`, and the next `/start` runs as usual. @@ -296,7 +299,34 @@ JUNGFRAU 8M, against about half a minute for the same merge made from nothing. T writer in the end message, as for rotation (`.mtz`), and the run's step "Merge" carries the final statistics. -`GET /rugnux/status` reports such a run as three steps (pre-pass, GPU analysis, merge) and a `stills` block: the +With an uploaded model, the model step runs on every one of these merges, as `rugnux --mode scale --model` +places a model and `rugnux --mode refine` refines it: the model is placed by its rigid body where it fits the +data as given - its indexing ambiguity probed by R-free, and the merge put in the model's indexing - and by +molecular replacement where it does not (or is a prediction); the 2mFo-DFc and mFo-DFc maps are made, and the +placed model is refined, with the blobs of its difference map. Without a space group in the dataset settings +the frames are merged in the model's own, where the model's cell is the consensus cell of the run on the same +axes; otherwise the model is left out, with a warning in the live status saying why. The step on a provisional +merge runs on a thread of its own, on 8 CPU threads and the GPU, beside the analysis and the merge, which do not +wait for it: its result - the placed and the refined model, their maps, the scores - goes to the live status +(`model`, `/rugnux/model.cif`, `/rugnux/map.ccp4`) as soon as it is made, so the density can be looked at while +the collection goes on. A step still running when the next merge is made is cut (a refinement cut short keeps +the weight trials it finished), and each step starts from the placement of the one before: the merge is put in +the same indexing, and the model goes in where it was placed, with no null, indexing probe or search made again. +After the last frame the step on the merge of the collection runs on every thread, and the end message waits +for it until `deadline_s` after the last image; then it is cut as for rotation - the refinement first (cut +short, or left out where it has not started), and a molecular replacement that has not started. The end message +carries the merge in the model's indexing (`.mtz`), the map coefficients, the placed and the +refined model, the refined maps and the blobs, written as for rotation; one that would not fit the pushers' +buffer goes without them. The GPU memory the analysis's indexers hold for the broker's life is counted against +the card when the model step decides what goes on it. `POST /cancel` cuts it +the same way, at any time (`cut` in `/rugnux/status`): the frames that came are still analysed and merged. + +The frames wait in the image buffer compressed, and go on to the writer byte for byte, with their analysis. `efficiency` is the +images received over the images expected; a stream has no module packets, so per-image collection +efficiency, ADU histograms and pulse IDs are absent. The dark frames of the mask sequence stay outside this mode. + +`GET /rugnux/status` reports such a run as three steps (pre-pass, GPU analysis, merge; with a model, the steps of +the model on the merge of the collection after them) and a `stills` block: the pre-pass state and duration, the stated and decided beam centre and distance, the pixels masked, and the images analysed and indexed with the analysis rate. The "Rugnux live" page of the web frontend shows the same. `tools/jfjoch_stills_gpu_replay` replays a stored serial collection through the same analysis. diff --git a/frontend/src/components/RugnuxPage.tsx b/frontend/src/components/RugnuxPage.tsx index abf3d164a..2322c28af 100644 --- a/frontend/src/components/RugnuxPage.tsx +++ b/frontend/src/components/RugnuxPage.tsx @@ -119,13 +119,18 @@ function SummaryTiles({s}: { s: rugnux_live_status }) { function ModelSection({s}: { s: rugnux_live_status }) { const m = s.model!; const mrSeconds = s.steps.find(step => step.title === "Molecular replacement")?.duration_s; + // Serial data place the model on each provisional merge while collecting: the viewer reloads when the model + // step has a new result, not every time the status is polled. + const resultKey = s.stills + ? `${m.placement} ${m.r_work} ${m.r_free} ${m.refinement?.r_free} ${s.refined_model}` + : `${s.elapsed_s}`; return (
{s.placed_model && Loading viewer...}> - @@ -166,7 +171,7 @@ function RugnuxPage({measuring}: { measuring: boolean }) {
- {s.model && s.state !== 'Running' && } + {s.model && (s.state !== 'Running' || s.stills) && }
diff --git a/frontend/src/components/RugnuxSettings.tsx b/frontend/src/components/RugnuxSettings.tsx index 8fd76bbca..c3f6e8ab0 100644 --- a/frontend/src/components/RugnuxSettings.tsx +++ b/frontend/src/components/RugnuxSettings.tsx @@ -50,7 +50,7 @@ function RugnuxSettings({measuring}: { measuring: boolean }) { onUpload={() => settings.submit({enable, reserved_cpus: cpus, deadline_s: deadline})} uploadDisabled={measuring || cpusError || deadlineError || settings.pending} snackbar={settings.snackbar}> - setEnable(e.target.checked)}/>}/> &IndexAndRefine::GetIntegrationOutcome() { return integration_outcome; } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 0467e31ad..f90c3267b 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -309,7 +309,5 @@ public: // The retained outcome of one image, moved out: what is left in its place is empty. Safe beside the // analysis of other images. IntegrationOutcome TakeIntegrationOutcome(int64_t image); - // The retained outcome of one image dropped, and its cell with it: it is not in the consensus. - void DropIntegrationOutcome(int64_t image); const std::vector &GetIntegrationOutcome() const; }; diff --git a/image_analysis/indexing/FFTIndexerGPU.cu b/image_analysis/indexing/FFTIndexerGPU.cu index 234820928..d048e093f 100644 --- a/image_analysis/indexing/FFTIndexerGPU.cu +++ b/image_analysis/indexing/FFTIndexerGPU.cu @@ -263,6 +263,16 @@ FFTIndexerGPU::FFTIndexerGPU(const IndexingSettings &settings) } +size_t FFTIndexerGPU::DeviceBytes() const { + const size_t batch = std::min(nDirections, FFT_BATCH_DIRECTIONS); + const size_t out_len = histogram_size / 2 + 1; + const size_t directions = nDirections; + return batch * histogram_size * sizeof(float) + batch * out_len * sizeof(cufftComplex) + + directions * out_len * sizeof(float) + directions * sizeof(FFTResult) + + 3 * FFT_MAX_SPOTS * sizeof(float) + 3 * directions * sizeof(float) + + plan.WorkBytes() + (nDirections % batch != 0 ? plan_last.WorkBytes() : 0); +} + void FFTIndexerGPU::DirectionsChanged() { std::vector dir_x(nDirections), dir_y(nDirections), dir_z(nDirections); for (int i = 0; i < nDirections; i++) { diff --git a/image_analysis/indexing/FFTIndexerGPU.h b/image_analysis/indexing/FFTIndexerGPU.h index 99cf0c3e5..03658bacf 100644 --- a/image_analysis/indexing/FFTIndexerGPU.h +++ b/image_analysis/indexing/FFTIndexerGPU.h @@ -58,6 +58,7 @@ public: FFTIndexerGPU(const FFTIndexerGPU &i) = delete; const FFTIndexerGPU &operator=(const FFTIndexerGPU &i) = delete; ~FFTIndexerGPU() override = default; + size_t DeviceBytes() const override; }; diff --git a/image_analysis/indexing/Indexer.h b/image_analysis/indexing/Indexer.h index 903a38737..0ecb7ec57 100644 --- a/image_analysis/indexing/Indexer.h +++ b/image_analysis/indexing/Indexer.h @@ -86,6 +86,8 @@ protected: virtual std::optional RunSeverityOnly(const std::vector &coord) { return {}; } public: virtual ~Indexer() = default; + // The device memory the indexer holds for as long as it lives; 0 for one on the CPU. + virtual size_t DeviceBytes() const { return 0; } void Setup(const DiffractionExperiment& experiment); // severity_only = true runs RunSeverityOnly instead of indexing: no lattice comes back and // `executed` stays false, since nothing was an indexing attempt. Used when the seed escalation diff --git a/image_analysis/indexing/IndexerThreadPool.cpp b/image_analysis/indexing/IndexerThreadPool.cpp index ea4150583..25b761483 100644 --- a/image_analysis/indexing/IndexerThreadPool.cpp +++ b/image_analysis/indexing/IndexerThreadPool.cpp @@ -246,6 +246,10 @@ void IndexerThread::Worker(int threadid) { { std::unique_lock lock(m); + device_ = get_gpu(); + for (const auto *i : {fft_indexer.get(), ffbidx_indexer.get(), fftw_indexer.get()}) + if (i) + device_bytes_ += i->DeviceBytes(); state = TaskState::IDLE; } c_running.notify_all(); @@ -392,6 +396,16 @@ IndexerThreadPool::IndexerThreadPool(const IndexingSettings &settings, IndexerCo tasks.emplace_back(std::make_unique(std::cref(settings), i, construction, beside)); } +size_t IndexerThreadPool::DeviceBytes() const { + std::map by_card; + for (const auto &t : tasks) + by_card[t->Device()] += t->DeviceBytes(); + size_t most = 0; + for (const auto &[card, bytes] : by_card) + most = std::max(most, bytes); + return most; +} + void IndexerThreadPool::ReleaseIndexers() { for (auto &t : tasks) t->Release(); diff --git a/image_analysis/indexing/IndexerThreadPool.h b/image_analysis/indexing/IndexerThreadPool.h index 023f6971d..7e3b7fdd6 100644 --- a/image_analysis/indexing/IndexerThreadPool.h +++ b/image_analysis/indexing/IndexerThreadPool.h @@ -98,6 +98,9 @@ class IndexerThread { std::unique_ptr result = nullptr; std::unique_ptr task_input = nullptr; std::thread worker_thread; + // The card the worker is pinned to, and the device memory of the indexers it built before it was ready + int device_ = 0; + size_t device_bytes_ = 0; void Worker(int threadid); public: @@ -109,6 +112,8 @@ public: // them again on the next frame that needs one. Returns when they are gone. void Release(); void Finalize(); + [[nodiscard]] int Device() const { return device_; } + [[nodiscard]] size_t DeviceBytes() const { return device_bytes_; } }; class IndexerThreadPool { @@ -137,6 +142,10 @@ public: // 0.6-1.2 GB of it - for a stretch with nothing to index. A frame that comes after all builds its // indexer again, as on first use, and the indexer it gets computes what the dropped one would. void ReleaseIndexers(); + // The device memory the indexers built with the pool (IndexerConstruction::Preconstruct) hold, on the + // card they hold most on: what is resident beside anything else that plans its GPU memory on a card while + // the pool lives (model validation, molecular replacement, refinement). + [[nodiscard]] size_t DeviceBytes() const; }; diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index 80c03c117..d22f64de4 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -3,14 +3,9 @@ #include "JFJochReceiverFPGA.h" -#include "../image_analysis/WriteReflections.h" - #include #include "ImageMetadata.h" -#include "JFJochReceiverLite.h" -#include "../common/CUDAWrapper.h" -#include "../image_analysis/MXAnalysisWithoutFPGA.h" JFJochReceiverFPGA::JFJochReceiverFPGA(const DiffractionExperiment &in_experiment, const PixelMask &in_pixel_mask, @@ -24,9 +19,7 @@ JFJochReceiverFPGA::JFJochReceiverFPGA(const DiffractionExperiment &in_experimen ImageBuffer &in_send_buf_ctrl, ZMQPreviewSocket *in_zmq_preview_socket, ZMQMetadataSocket *in_zmq_metadata_socket, - IndexerThreadPool *indexing_thread_pool, - bool in_gpu_analysis, - RugnuxLiveStatus *in_rugnux_status) + IndexerThreadPool *indexing_thread_pool) : JFJochReceiver(in_experiment, in_send_buf_ctrl, in_image_sender, @@ -46,10 +39,7 @@ JFJochReceiverFPGA::JFJochReceiverFPGA(const DiffractionExperiment &in_experimen pedestal_nthreads((experiment.GetStorageCellNumber() > 2) ? 1 : 4), summation_nthreads(2), frame_transformation_ready((experiment.GetImageNum() > 0) ? frame_transformation_nthreads : 0), - data_acquisition_ready(ndatastreams), - gpu_analysis(in_gpu_analysis && experiment.GetImageNum() > 0 - && experiment.GetDetectorMode() == DetectorMode::Standard), - rugnux_status(in_rugnux_status) { + data_acquisition_ready(ndatastreams) { // The FPGA integration core has a fixed bin-count limit. Reject configurations it cannot // handle before starting acquisition, unless azimuthal integration is forced onto the CPU. const auto &azint_settings = experiment.GetAzimuthalIntegrationSettings(); @@ -128,37 +118,8 @@ JFJochReceiverFPGA::JFJochReceiverFPGA(const DiffractionExperiment &in_experimen } logger.Debug("Pedestal threads ready ({} threads/module*SC)", pedestal_nthreads); } else if (experiment.GetImageNum() > 0) { - if (gpu_analysis) { - // The start message and the calibration go out once the pre-pass is done. - gpu = std::make_unique( - experiment, pixel_mask, MakeStartMessage(experiment, pixel_mask, *az_int_mapping), - indexing_thread_pool, std::min(16, frame_transformation_nthreads), - // Half the image buffer for the frames waiting, the other half for those sent. - image_buffer.GetStatus().total_slots / 2, logger, - [this](const StillsGPUAnalysis::Setup &s) { - // The profiles are binned on the pre-pass's geometry. - plots.Setup(s.experiment, *s.mapping); - SendStartMessage(MakeStartMessage(s.experiment, s.mask, *s.mapping)); - SendCalibration(); - }, - [this](DataMessage &m, const AzimuthalIntegrationProfile &p) { GPUSink(m, p); }, - [this] { - // The broker's settings, with the spot finder rugnux uses on stills (the pre-pass - // decided the centre and mask with it): adaptive threshold, single-pixel spots kept. - // With the broker's default fixed threshold 3 points fewer frames index. - SpotFindingSettings s = GetSpotFindingSettings(); - s.adaptive_threshold = true; - s.min_pix_per_spot = std::nullopt; - return s; - }, - // The start message did not go out (or the analysis could not be set up): nothing can - // be written, so the collection stops, as it does on an error in a frame thread. - [this](const JFJochException &e) { Cancel(e); }, - rugnux_status); - } else { - SendStartMessage(); - SendCalibration(); - } + SendStartMessage(); + SendCalibration(); for (int i = 0; i < experiment.GetImageNum(); i++) images_to_go.Put(i); @@ -459,34 +420,6 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { experiment.GetImageMode(), CompressionAlgorithm::NO_COMPRESSION); - if (gpu_analysis) { - // Compressed here into a slot of the image buffer, where it waits for the GPU analysis; GPUSink - // writes it, with its analysis, to a second slot, which is sent. A frame that finds no slot - // free, or as many frames waiting as the analysis is to hold, is dropped, as the receiver drops - // any frame that finds the image buffer full. - ++images_collected; - uncompressed_size += experiment.GetModulesNum() * RAW_MODULE_SIZE * experiment.GetByteDepthImage(); - auto held = std::make_shared(); - held->slot = image_buffer.GetImageSlot(); - if (held->slot == nullptr) { - writer_queue_full = true; - continue; - } - const size_t size = transformation.CompressImage(static_cast(held->slot->GetImage()), - experiment.GetImageBufferLocationSize()); - held->message = std::move(message); - held->message.image = CompressedImage(held->slot->GetImage(), size, - experiment.GetXPixelsNum(), experiment.GetYPixelsNum(), - experiment.GetImageMode(), - experiment.GetCompressionAlgorithm()); - if (!gpu->Hold(held)) { - writer_queue_full = true; - continue; - } - compressed_size += size; - continue; - } - // No-op unless the CPU backend is forced; fills az_int_profile_image from the assembled image analyzer->RunAzimuthalIntegration(transformation.GetImage(), az_int_profile_image); @@ -665,22 +598,10 @@ void JFJochReceiverFPGA::FinalizeMeasurement() { logger.Debug("All processing threads done"); } - // Where the start message never went out (or the analysis could not be set up), there is no end - // message either: the devices are stopped all the same, and the error is the collection's result. - std::exception_ptr gpu_error; - if (gpu) { - try { - gpu->Finish(); - } catch (const std::exception &) { - gpu_error = std::current_exception(); - } - } - current_status.SetProgress(1.0); current_status.SetStatus(GetStatus()); - if (!gpu_error) - SendEndMessage(); + SendEndMessage(); if (experiment.GetImageNum() > 0) { for (int d = 0; d < ndatastreams; d++) @@ -699,18 +620,10 @@ void JFJochReceiverFPGA::FinalizeMeasurement() { throw JFJochException(JFJochExceptionCategory::ZeroMQ, "Send commands not finalized in " + std::to_string(SendFinalizeTimeout.count()) + " seconds"); } - if (gpu_error) { - current_status.SetProgress({}); - current_status.SetStatus(GetStatus()); - std::rethrow_exception(gpu_error); - } logger.Info("Receiving data done"); if (push_images_to_writer) writer_error = image_pusher.Finalize(); - // Last image -> all written, as the rugnux receiver reports it - if (gpu) - gpu->Written(); current_status.SetProgress({}); current_status.SetStatus(GetStatus()); @@ -790,96 +703,3 @@ void JFJochReceiverFPGA::LoadCalibrationToFPGA(uint16_t data_stream) { // Initialize data processing acquisition_device[data_stream].InitializeDataProcessing(experiment, *az_int_mapping); } - -void JFJochReceiverFPGA::GPUSink(DataMessage &message, const AzimuthalIntegrationProfile &profile) { - try { - const auto status = image_buffer.GetStatus(); - message.receiver_buf_available = status.available_slots; - message.receiver_buf_in_preparation = status.preparation_slots; - message.receiver_buf_in_sending = status.sending_slots; - if (message.image.GetCompressedSize() > 0) - message.compression_ratio = static_cast(experiment.GetPixelsNum() * experiment.GetByteDepthImage()) - / static_cast(message.image.GetCompressedSize()); - - if (message.image_collection_efficiency == 1.0f) { - saturated_pixels.Add(message.saturated_pixel_count); - error_pixels.Add(message.error_pixel_count); - if (message.roi.contains("beam")) { - roi_beam_npixel.Add(message.roi["beam"].pixels); - roi_beam_sum.Add(message.roi["beam"].sum); - } - } - - scan_result.Add(message); - - if (!serialmx_filter.ApplyFilter(message)) - ++images_skipped; - else if (auto loc = image_buffer.GetImageSlot(); loc == nullptr) - writer_queue_full = true; - else { - auto writer_buffer = (uint8_t *) loc->GetImage(); - CBORStream2Serializer serializer(writer_buffer, experiment.GetImageBufferLocationSize()); - if (JFJochReceiverLite::SerializeImageOrDrop(serializer, message, *loc, - experiment.GetImageBufferLocationSize(), logger)) { - loc->SetImageNumber(message.number); - loc->SetImageSize(serializer.GetBufferSize()); - loc->SetIndexed(message.indexing_result.value_or(false)); - loc->ReadyToSend(); - - if (zmq_preview_socket != nullptr) - zmq_preview_socket->SendImage(writer_buffer, serializer.GetBufferSize()); - if (zmq_metadata_socket != nullptr) - zmq_metadata_socket->AddDataMessage(message); - - if (push_images_to_writer) { - if (image_pusher.SendImage(*loc)) - ++images_sent; - } else - loc->release(); - UpdateMaxImageSent(message.number); - } - } - - plots.Add(message, profile); - current_status.SetProgress(GetProgress()); - current_status.SetStatus(GetStatus()); - } catch (const JFJochException &e) { - logger.ErrorException(e); - Cancel(e); - } catch (const std::exception &e) { - const JFJochException error(JFJochExceptionCategory::WrongDAQState, - fmt::format("Image {} not written - {}", message.number, e.what())); - logger.ErrorException(error); - Cancel(error); - } -} - -void JFJochReceiverFPGA::AddToEndMessage(EndMessage &message) { - if (!gpu) - return; - message.unit_cell = gpu->GetConsensusUnitCell(); - message.refined_beam_center_x = gpu->GetBeamX_pxl(); - message.refined_beam_center_y = gpu->GetBeamY_pxl(); - message.refined_detector_distance_mm = gpu->GetDetectorDistance_mm(); - // The merged data, for the writer to put beside the master file. - if (const auto &m = gpu->GetMerged(); m && gpu->GetMergeExperiment().GetUnitCell()) - message.merged = MergedMtzTable(m->merge.merged, *gpu->GetMergeExperiment().GetUnitCell(), - gpu->GetMergeExperiment(), OutputProvenance{}); - if (!message.merged) - return; - // Every pusher serialises END into a buffer of MESSAGE_SIZE_FOR_START_END, and an END that does not fit there - // fails the collection. So it is serialised here first, as JFJochReceiverRugnux does; one that does not fit - // goes without the merged data. - const auto scratch = std::make_unique_for_overwrite(MESSAGE_SIZE_FOR_START_END); - CBORStream2Serializer serializer(scratch.get(), MESSAGE_SIZE_FOR_START_END); - try { - serializer.SerializeSequenceEnd(message); - } catch (const JFJochException &e) { - const std::string text = fmt::format("the merged data do not fit in the end message ({} MiB) and are not " - "written", MESSAGE_SIZE_FOR_START_END / (1024 * 1024)); - logger.Error("GPU analysis: {} - {}", text, e.what()); - if (rugnux_status) - rugnux_status->Warning(text); - message.merged.reset(); - } -} diff --git a/receiver/JFJochReceiverFPGA.h b/receiver/JFJochReceiverFPGA.h index 167fe50af..d4bcdf3d0 100644 --- a/receiver/JFJochReceiverFPGA.h +++ b/receiver/JFJochReceiverFPGA.h @@ -40,7 +40,6 @@ #include "../image_analysis/MXAnalysisAfterFPGA.h" #include "../common/MovingAverage.h" #include "JFJochReceiverOutput.h" -#include "StillsGPUAnalysis.h" class JFJochReceiverFPGA : public JFJochReceiver { const JFCalibration *calibration; @@ -87,14 +86,6 @@ class JFJochReceiverFPGA : public JFJochReceiver { void Cancel(const JFJochException &e) override; void FinalizeMeasurement(); - // Serial data with rugnux on: the converted frames are analysed on the GPU (StillsGPUAnalysis), and - // the start message waits for its pre-pass. - const bool gpu_analysis; - RugnuxLiveStatus *const rugnux_status; - std::unique_ptr gpu; - void GPUSink(DataMessage &message, const AzimuthalIntegrationProfile &profile); - void AddToEndMessage(EndMessage &message) override; - void SendCalibration(); void SendPedestal(const std::string &prefix, const std::vector &v, int gain, int sc); void RetrievePedestal(std::vector &output) const; @@ -111,9 +102,7 @@ public: ImageBuffer &image_buffer, ZMQPreviewSocket *zmq_preview_socket, ZMQMetadataSocket *zmq_metadata_socket, - IndexerThreadPool *indexing_thread_pool, - bool gpu_analysis = false, - RugnuxLiveStatus *rugnux_status = nullptr); + IndexerThreadPool *indexing_thread_pool); ~JFJochReceiverFPGA(); JFJochReceiverFPGA(const JFJochReceiverFPGA &other) = delete; JFJochReceiverFPGA& operator=(const JFJochReceiverFPGA &other) = delete; diff --git a/receiver/JFJochReceiverRugnux.cpp b/receiver/JFJochReceiverRugnux.cpp index 6afcbd87a..6ebd1a821 100644 --- a/receiver/JFJochReceiverRugnux.cpp +++ b/receiver/JFJochReceiverRugnux.cpp @@ -4,6 +4,8 @@ #include "JFJochReceiverRugnux.h" #include +#include +#include #include "spdlog/fmt/fmt.h" @@ -12,6 +14,7 @@ #include "../frame_serialize/CBORStream2Deserializer.h" #include "../frame_serialize/CBORStream2Serializer.h" #include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../rugnux/RugnuxDefaults.h" #include "../rugnux/RugnuxLive.h" #include "../rugnux/ResultReport.h" #include "../writer/RugnuxResultBlocks.h" @@ -195,7 +198,8 @@ JFJochReceiverRugnux::JFJochReceiverRugnux(const DiffractionExperiment &in_exper JFJochReceiverPlots &in_plots, ImageBuffer &in_image_buffer, ZMQPreviewSocket *in_zmq_preview_socket, - ZMQMetadataSocket *in_zmq_metadata_socket) + ZMQMetadataSocket *in_zmq_metadata_socket, + IndexerThreadPool *in_indexing_thread_pool) : JFJochReceiver(in_experiment, in_image_buffer, in_image_pusher, in_current_status, in_plots, in_spot_finding_settings, in_logger, in_pixel_mask, in_zmq_preview_socket, in_zmq_metadata_socket, nullptr), @@ -208,16 +212,21 @@ JFJochReceiverRugnux::JFJochReceiverRugnux(const DiffractionExperiment &in_exper deadline_s(in_deadline_s), live_status(in_live_status), measurement_started(1), - forward_started(forward_nthreads) { - // GPU FFT indexers and cuFFT plans are kept for reuse for the life of the broker, not of one run as - // in `rugnux` (FFTIndexerReuse, FFTPlanReuse): every pass of every collection builds the same pools, - // and making a plan takes a process-wide lock and took over a second on a card busy with the online - // analysis. Never destroyed - at exit the CUDA runtime may already be gone. - [[maybe_unused]] static const FFTIndexerReuse *fft_indexers = new FFTIndexerReuse; - [[maybe_unused]] static const FFTPlanReuse *fft_plans = new FFTPlanReuse; - logger.Info("rugnux live: starting {} forward threads; rugnux runs on {} threads while the images arrive, " - "planning against {:.1f} GB of host memory", forward_nthreads, rugnux_nthreads, - static_cast(host_budget) / 1e9); + forward_started(forward_nthreads), + stills_run(!RugnuxRotationData(in_experiment)), + indexing_thread_pool(in_indexing_thread_pool), + stills_nthreads(std::min(16, in_forward_nthreads)) { + if (!stills_run) { + // GPU FFT indexers and cuFFT plans are kept for reuse for the life of the broker, not of one run as + // in `rugnux` (FFTIndexerReuse, FFTPlanReuse): every pass of every collection builds the same pools, + // and making a plan takes a process-wide lock and took over a second on a card busy with the online + // analysis. Never destroyed - at exit the CUDA runtime may already be gone. + [[maybe_unused]] static const FFTIndexerReuse *fft_indexers = new FFTIndexerReuse; + [[maybe_unused]] static const FFTPlanReuse *fft_plans = new FFTPlanReuse; + logger.Info("rugnux live: starting {} forward threads; rugnux runs on {} threads while the images arrive, " + "planning against {:.1f} GB of host memory", forward_nthreads, rugnux_nthreads, + static_cast(host_budget) / 1e9); + } for (int i = 0; i < forward_nthreads; i++) forward_futures.emplace_back(std::async(std::launch::async, &JFJochReceiverRugnux::ForwardThread, this, i)); measurement = std::async(std::launch::async, &JFJochReceiverRugnux::MeasurementThread, this); @@ -268,24 +277,28 @@ void JFJochReceiverRugnux::MeasurementThread() { JFJochReceiverLite::ApplyStartMessage(experiment, msg->cbor->start_message.value()); image_buffer.StartMeasurement(experiment); start_time = std::chrono::system_clock::now(); - SendStartMessage(); + if (stills_run) + StartStills(); + else { + SendStartMessage(); - // Rugnux reads the images back from the buffer, so it needs one slot per image: a slot taken - // again would hand it a later image. Without that the collection is only forwarded. - const auto slots = image_buffer.GetStatus().total_slots; - // The image time as the detector produces images, summation included - const double image_time_s = std::chrono::duration(experiment.GetImageTime()).count(); - std::optional expected_collection_s; - if (image_time_s > 0) - expected_collection_s = static_cast(experiment.GetImageNum()) * image_time_s; - live_status.Start(experiment.GetImageNum(), expected_collection_s, msg->received); - if (slots < experiment.GetImageNum()) { - skipped_warning = fmt::format("no rugnux analysis: the image buffer holds {} images, the collection " - "has {}", slots, experiment.GetImageNum()); - logger.Error("{}", skipped_warning); - live_status.Finish(RugnuxLiveState::Skipped, "the image buffer cannot hold the whole collection"); - } else - StartRugnux(rugnux_done, rugnux_out, probe_done); + // Rugnux reads the images back from the buffer, so it needs one slot per image: a slot taken + // again would hand it a later image. Without that the collection is only forwarded. + const auto slots = image_buffer.GetStatus().total_slots; + // The image time as the detector produces images, summation included + const double image_time_s = std::chrono::duration(experiment.GetImageTime()).count(); + std::optional expected_collection_s; + if (image_time_s > 0) + expected_collection_s = static_cast(experiment.GetImageNum()) * image_time_s; + live_status.Start(experiment.GetImageNum(), expected_collection_s, msg->received); + if (slots < experiment.GetImageNum()) { + skipped_warning = fmt::format("no rugnux analysis: the image buffer holds {} images, the collection " + "has {}", slots, experiment.GetImageNum()); + logger.Error("{}", skipped_warning); + live_status.Finish(RugnuxLiveState::Skipped, "the image buffer cannot hold the whole collection"); + } else + StartRugnux(rugnux_done, rugnux_out, probe_done); + } measurement_started.count_down(); } catch (const JFJochException &e) { logger.ErrorException(e); @@ -299,6 +312,11 @@ void JFJochReceiverRugnux::MeasurementThread() { throw; } + if (stills_run) { + StillsMeasurement(); + return; + } + logger.Info("Receiving started"); for (auto &f: forward_futures) f.get(); @@ -365,6 +383,212 @@ void JFJochReceiverRugnux::MeasurementThread() { current_status.SetStatus(GetStatus()); } +void JFJochReceiverRugnux::StartStills() { + const int64_t slots = image_buffer.GetStatus().total_slots; + stills_keep_free = std::min(stills_nthreads, slots / 2); + auto analysis = std::make_unique( + experiment, pixel_mask, MakeStartMessage(experiment, pixel_mask, *az_int_mapping), + indexing_thread_pool, stills_nthreads, + // What the pre-pass holds and fits on; the frames beyond it wait in the image buffer as well. + slots / 2, logger, + [this](const StillsGPUAnalysis::Setup &s) { + // The profiles are binned on the pre-pass's geometry. + plots.Setup(s.experiment, *s.mapping); + std::unique_lock ul(calibration_m); + SendStartMessage(MakeStartMessage(s.experiment, s.mask, *s.mapping)); + start_sent = true; + if (push_images_to_writer) + for (const auto &c: calibration_held) + image_pusher.SendCalibration(c.cbor->calibration.value()); + calibration_held.clear(); + }, + [this](DataMessage &m, const AzimuthalIntegrationProfile &p) { SendAnalysedImage(m, p); }, + [this] { + // The broker's settings, with the spot finder rugnux uses on stills (the pre-pass decided the + // centre and mask with it): adaptive threshold, single-pixel spots kept. With the broker's + // default fixed threshold 3 points fewer frames index. + SpotFindingSettings settings = GetSpotFindingSettings(); + settings.adaptive_threshold = true; + settings.min_pix_per_spot = std::nullopt; + return settings; + }, + // The start message did not go out (or the analysis could not be set up): nothing can be + // written, so the collection stops. + [this](const JFJochException &e) { Cancel(e); }, + &live_status, true, model, deadline_s); + std::unique_lock ul(rugnux_m); + stills = std::move(analysis); + if (cut) + stills->Cut("cancel"); +} + +void JFJochReceiverRugnux::StillsMeasurement() { + logger.Info("Receiving started"); + for (auto &f: forward_futures) + f.get(); + images_done = true; + logger.Info("Receiving finished, {} images", images_collected.load()); + + // Where the start message never went out (or the analysis could not be set up), there is no end message + // either, and the error is the collection's result. + std::exception_ptr error; + try { + stills->Finish(); + } catch (const std::exception &) { + error = std::current_exception(); + } + + current_status.SetProgress(1.0); + current_status.SetStatus(GetStatus()); + + if (!error) + SendEndMessage(); + + end_time = std::chrono::system_clock::now(); + + if (!image_buffer.CheckIfBufferReturned(SendFinalizeTimeout)) { + logger.Error("Send commands not finalized in {} seconds", SendFinalizeTimeout.count()); + throw JFJochException(JFJochExceptionCategory::ZeroMQ, "Send commands not finalized in " + + std::to_string(SendFinalizeTimeout.count()) + " seconds"); + } + if (error) { + current_status.SetProgress({}); + current_status.SetStatus(GetStatus()); + std::rethrow_exception(error); + } + if (push_images_to_writer) + writer_error = image_pusher.Finalize(); + stills->Written(); + + current_status.SetProgress({}); + current_status.SetStatus(GetStatus()); +} + +void JFJochReceiverRugnux::HoldStillsImage(const DataMessage &message) { + // The analysis keeps per-image results by number, for the images the collection declares; the stream is not ours. + if (message.number < 0 || message.number >= experiment.GetImageNum()) { + logger.Error("Dropping image {} - the collection has {} images", message.number, experiment.GetImageNum()); + return; + } + ++images_collected; + const size_t size = message.image.GetCompressedSize(); + compressed_size += size; + uncompressed_size += message.image.GetUncompressedSize(); + UpdateMaxImageReceived(message.number); + + // Copied into a slot of the image buffer, where it waits for the analysis; SendAnalysedImage writes it, with + // its analysis, to a second slot, which is sent. A frame that finds the buffer full, but for a slot kept + // free for each analysis thread to send with, is dropped, as any frame that finds the image buffer full. + auto held = std::make_shared(); + held->slot = image_buffer.GetImageSlot(stills_keep_free); + if (held->slot == nullptr) { + writer_queue_full = true; + return; + } + if (size > static_cast(experiment.GetImageBufferLocationSize())) { + logger.Error("Dropping image {} - {} B do not fit in buffer slot {} B", message.number, size, + experiment.GetImageBufferLocationSize()); + return; + } + memcpy(held->slot->GetImage(), message.image.GetCompressed(), size); + held->message = message; + held->message.image = CompressedImage(held->slot->GetImage(), size, message.image.GetWidth(), + message.image.GetHeight(), message.image.GetMode(), + message.image.GetCompressionAlgorithm(), message.image.GetChannel()); + held->message.original_number = message.number; + held->message.user_data = experiment.GetImageAppendix(); + held->message.run_number = experiment.GetRunNumber(); + held->message.run_name = experiment.GetRunName(); + held->message.receiver_aq_dev_delay = image_puller.GetCurrentFifoUtilization(); + stills->Hold(held); +} + +void JFJochReceiverRugnux::StillsCalibration(const ImagePullerOutput &msg) { + std::unique_lock ul(calibration_m); + if (!start_sent) + calibration_held.push_back(msg); + else if (push_images_to_writer) + image_pusher.SendCalibration(msg.cbor->calibration.value()); +} + +ZeroCopyReturnValue *JFJochReceiverRugnux::AnalysedImageSlot() { + // An analysed frame is not dropped because the writer lags behind the catch-up from a full buffer: the analysis + // waits for the slots being sent to come back, and the frames arriving meanwhile find the buffer full. A + // writer taking nothing for this long loses the frame. + constexpr auto WRITER_WAIT = std::chrono::seconds(2); + const auto t0 = std::chrono::steady_clock::now(); + auto loc = image_buffer.GetImageSlot(); + while (loc == nullptr && !cancelled && std::chrono::steady_clock::now() - t0 < WRITER_WAIT) { + std::this_thread::sleep_for(std::chrono::milliseconds(1)); + loc = image_buffer.GetImageSlot(); + } + return loc; +} + +void JFJochReceiverRugnux::SendAnalysedImage(DataMessage &message, const AzimuthalIntegrationProfile &profile) { + try { + const auto status = image_buffer.GetStatus(); + message.receiver_buf_available = status.available_slots; + message.receiver_buf_in_preparation = status.preparation_slots; + message.receiver_buf_in_sending = status.sending_slots; + if (message.image.GetCompressedSize() > 0) + message.compression_ratio = static_cast(experiment.GetPixelsNum() * experiment.GetByteDepthImage()) + / static_cast(message.image.GetCompressedSize()); + + // A stream that does not state it (DECTRIS) delivers whole images. + if (message.image_collection_efficiency.value_or(1.0f) == 1.0f) { + saturated_pixels.Add(message.saturated_pixel_count); + error_pixels.Add(message.error_pixel_count); + if (message.roi.contains("beam")) { + roi_beam_npixel.Add(message.roi["beam"].pixels); + roi_beam_sum.Add(message.roi["beam"].sum); + } + } + + scan_result.Add(message); + + if (!serialmx_filter.ApplyFilter(message)) + ++images_skipped; + else if (auto loc = AnalysedImageSlot(); loc == nullptr) + writer_queue_full = true; + else { + auto writer_buffer = (uint8_t *) loc->GetImage(); + CBORStream2Serializer serializer(writer_buffer, experiment.GetImageBufferLocationSize()); + if (JFJochReceiverLite::SerializeImageOrDrop(serializer, message, *loc, + experiment.GetImageBufferLocationSize(), logger)) { + loc->SetImageNumber(message.number); + loc->SetImageSize(serializer.GetBufferSize()); + loc->SetIndexed(message.indexing_result.value_or(false)); + loc->ReadyToSend(); + + if (zmq_preview_socket != nullptr) + zmq_preview_socket->SendImage(writer_buffer, serializer.GetBufferSize()); + if (zmq_metadata_socket != nullptr) + zmq_metadata_socket->AddDataMessage(message); + + if (push_images_to_writer) { + if (image_pusher.SendImage(*loc)) + ++images_sent; + } else + loc->release(); + UpdateMaxImageSent(message.number); + } + } + + plots.Add(message, profile); + current_status.SetProgress(GetProgress()); + current_status.SetStatus(GetStatus()); + } catch (const JFJochException &e) { + logger.ErrorException(e); + Cancel(e); + } catch (const std::exception &e) { + const JFJochException error(JFJochExceptionCategory::WrongDAQState, + fmt::format("Image {} not written - {}", message.number, e.what())); + logger.ErrorException(error); + Cancel(error); + } +} + void JFJochReceiverRugnux::StartRugnux(std::future &rugnux_done, std::future &rugnux_out, std::future &probe_done) { // A failure here costs the analysis only: the images go on to the writer as they would without it. @@ -516,6 +740,10 @@ void JFJochReceiverRugnux::ProbeThread() { } void JFJochReceiverRugnux::AddToEndMessage(EndMessage &message) { + if (stills) { + stills->FillEndMessage(message); + return; + } // A failure here costs the end message the analysis only: it goes out as the receiver made it. try { std::unique_lock ul(rugnux_m); @@ -605,6 +833,10 @@ void JFJochReceiverRugnux::ForwardThread(uint32_t id) { continue; if (msg->cbor->end_message) { end_message_received = true; + } else if (stills_run && msg->cbor->calibration) { + StillsCalibration(*msg); + } else if (stills_run && msg->cbor->data_message) { + HoldStillsImage(msg->cbor->data_message.value()); } else if (msg->cbor->calibration) { if (push_images_to_writer) image_pusher.SendCalibration(msg->cbor->calibration.value()); @@ -686,6 +918,10 @@ void JFJochReceiverRugnux::Cut(const std::string &reason) { if (cut || rugnux_result) return; cut = true; + if (stills) { + stills->Cut(reason); + return; + } reader.Ended(); // no image is waited for any more if (rugnux != nullptr) rugnux->Cancel(); @@ -702,6 +938,15 @@ void JFJochReceiverRugnux::StopReceiver() { } } +void JFJochReceiverRugnux::SetSpotFindingSettings(const SpotFindingSettings &in_spot_finding_settings) { + // A rotation run takes rugnux's own spot finding; the frames of a serial run are analysed with these. + if (!stills_run) + return; + std::unique_lock ul(spot_finding_settings_mutex); + DiffractionExperiment::CheckDataProcessingSettings(in_spot_finding_settings); + spot_finding_settings = in_spot_finding_settings; +} + float JFJochReceiverRugnux::GetEfficiency() const { if (experiment.GetFrameNum() == 0) return 0; diff --git a/receiver/JFJochReceiverRugnux.h b/receiver/JFJochReceiverRugnux.h index dabd7f6e5..c316012a5 100644 --- a/receiver/JFJochReceiverRugnux.h +++ b/receiver/JFJochReceiverRugnux.h @@ -11,6 +11,7 @@ #include "../reader/JFJochImageBufferReader.h" #include "../rugnux/Rugnux.h" #include "RugnuxLiveStatus.h" +#include "StillsGPUAnalysis.h" // What the run made of the model it was given, for the live status; none where no model was given. std::optional LiveModel(const ProcessResult &result); @@ -20,6 +21,11 @@ std::optional LiveModel(const ProcessResult &result); // image buffer and send it on, and one rugnux thread reads the collection back from that buffer. The // end message waits for rugnux, and carries what it found - the lattice, refined geometry, per-image // results and its merged data, map coefficients and placed model (see AddToEndMessage). +// A serial collection (no goniometer scan) is analysed frame by frame instead (StillsGPUAnalysis): each image +// the puller delivers waits in the image buffer for the pre-pass, is analysed at the geometry it decides, and +// goes to the writer with its analysis; the start message waits for the pre-pass. The merges and the model step +// run beside it, and the end message carries the last of them. Only what the puller delivers is used: a +// DataMessage with a compressed image, whatever produced it. class JFJochReceiverRugnux : public JFJochReceiver { using JFJochReceiver::Cancel; @@ -51,6 +57,23 @@ class JFJochReceiverRugnux : public JFJochReceiver { bool cut = false; std::atomic images_done{false}; // the forward threads have ended + const bool stills_run; + IndexerThreadPool *const indexing_thread_pool; + const int64_t stills_nthreads; + int64_t stills_keep_free = 0; // slots frames do not take, for the analysis threads to send with + // Calibration that came while the pre-pass held the start message back. + std::mutex calibration_m; + std::vector calibration_held; + bool start_sent = false; + std::unique_ptr stills; // under rugnux_m until the frames are in, for Cut + + void StartStills(); + void StillsMeasurement(); + void HoldStillsImage(const DataMessage &message); + void StillsCalibration(const ImagePullerOutput &msg); + ZeroCopyReturnValue *AnalysedImageSlot(); + void SendAnalysedImage(DataMessage &message, const AzimuthalIntegrationProfile &profile); + void MeasurementThread(); // Opens the reader on the collection and starts the rugnux and probe threads; a failure is the // analysis's alone (live status Failed), never the collection's. rugnux_out is ready once the run's @@ -79,7 +102,8 @@ public: JFJochReceiverPlots &plots, ImageBuffer &image_buffer, ZMQPreviewSocket *zmq_preview_socket = nullptr, - ZMQMetadataSocket *zmq_metadata_socket = nullptr); + ZMQMetadataSocket *zmq_metadata_socket = nullptr, + IndexerThreadPool *indexing_thread_pool = nullptr); ~JFJochReceiverRugnux() override; void StopReceiver() override; @@ -89,5 +113,5 @@ public: float GetProgress() const override; JFJochReceiverOutput GetFinalStatistics() const override; - void SetSpotFindingSettings(const SpotFindingSettings &spot_finding_settings) override {} + void SetSpotFindingSettings(const SpotFindingSettings &spot_finding_settings) override; }; diff --git a/receiver/JFJochReceiverService.cpp b/receiver/JFJochReceiverService.cpp index d3f8f0fac..61c619f4c 100644 --- a/receiver/JFJochReceiverService.cpp +++ b/receiver/JFJochReceiverService.cpp @@ -4,7 +4,6 @@ #include "JFJochReceiverService.h" #include "JFJochReceiverRugnux.h" #include "../common/ThreadAffinity.h" -#include "../rugnux/RugnuxDefaults.h" #include "JFJochReceiverFPGA.h" #include "JFJochReceiverLite.h" #include "../preview/JFJochJPEG.h" @@ -114,11 +113,7 @@ void JFJochReceiverService::Start(const DiffractionExperiment &experiment, image_buffer, zmq_preview_socket.get(), zmq_metadata_socket.get(), - indexer_thread_pool.get(), - // Serial data with rugnux on: the frames are - // analysed on the GPU, after a pre-pass. - rugnux_enable && !RugnuxRotationData(experiment), - &rugnux_status); + indexer_thread_pool.get()); break; case DetectorType::DECTRIS: if (puller) @@ -127,12 +122,13 @@ void JFJochReceiverService::Start(const DiffractionExperiment &experiment, image_puller = std::make_shared( experiment.GetDetectorSetup().GetDECTRISStream2Addr()); } - if (rugnux_enable && RugnuxRotationData(experiment)) { + // The dark frames of the mask sequence take Lite's mask path, with rugnux on as well. + if (rugnux_enable && experiment.GetDetectorMode() == DetectorMode::Standard) { receiver = std::make_unique( experiment, pixel_mask, *image_puller, image_pusher, logger, nthreads_local, std::max(1, static_cast(AvailableCpus()) - rugnux_reserved_cpus), rugnux_model, rugnux_deadline_s, rugnux_status, spot_finding_settings, receiver_status, plots, image_buffer, - zmq_preview_socket.get(), zmq_metadata_socket.get()); + zmq_preview_socket.get(), zmq_metadata_socket.get(), indexer_thread_pool.get()); break; } receiver = std::make_unique(experiment, diff --git a/receiver/StillsGPUAnalysis.cpp b/receiver/StillsGPUAnalysis.cpp index e5ccebe65..e005604d3 100644 --- a/receiver/StillsGPUAnalysis.cpp +++ b/receiver/StillsGPUAnalysis.cpp @@ -8,18 +8,24 @@ #include "../common/CUDAWrapper.h" #include "../common/ThreadAffinity.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/WriteReflections.h" #include "../rugnux/RugnuxDefaults.h" +#include "../frame_serialize/CBORStream2Serializer.h" +#include "../writer/RugnuxResultBlocks.h" +#include "JFJochReceiverRugnux.h" +#include "RugnuxStructure.h" StillsGPUAnalysis::StillsGPUAnalysis(const DiffractionExperiment &in_experiment, const PixelMask &in_mask, const StartMessage &in_start, IndexerThreadPool *in_indexing_thread_pool, int64_t nthreads, int64_t max_held, Logger &in_logger, Ready in_ready, Sink in_sink, Settings in_settings, Failed in_failed, - RugnuxLiveStatus *in_status, bool in_merge) + RugnuxLiveStatus *in_status, bool in_merge, + std::shared_ptr in_model, double in_deadline_s) : experiment(in_experiment), mask(in_mask), start(in_start), indexing_thread_pool(in_indexing_thread_pool), logger(in_logger), ready(std::move(in_ready)), sink(std::move(in_sink)), settings(std::move(in_settings)), setup_failed(std::move(in_failed)), status(in_status), merge_frames(in_merge), - queue(std::max(1, max_held)) { - live.images_needed = std::min(PREPASS_IMAGES, experiment.GetImageNum()); + hold_limit(std::max(1, max_held)), model(std::move(in_model)), deadline_s(in_deadline_s) { + live.images_needed = std::min({PREPASS_IMAGES, experiment.GetImageNum(), hold_limit}); live.file_beam_x_pxl = experiment.GetBeamX_pxl(); live.file_beam_y_pxl = experiment.GetBeamY_pxl(); live.file_distance_mm = experiment.GetDetectorDistance_mm(); @@ -34,7 +40,7 @@ StillsGPUAnalysis::StillsGPUAnalysis(const DiffractionExperiment &in_experiment, prepass_future = std::async(std::launch::async, &StillsGPUAnalysis::PrePassThread, this); for (int64_t i = 0; i < nthreads; i++) analysis_futures.emplace_back(std::async(std::launch::async, &StillsGPUAnalysis::AnalysisThread, this)); - logger.Info("GPU analysis: {} threads, the first {} images held for the pre-pass", nthreads, PREPASS_IMAGES); + logger.Info("GPU analysis: {} threads, at least the first {} images held for the pre-pass", nthreads, PREPASS_IMAGES); } StillsGPUAnalysis::~StillsGPUAnalysis() { @@ -44,6 +50,7 @@ StillsGPUAnalysis::~StillsGPUAnalysis() { } catch (const std::exception &e) { logger.Error("GPU analysis: {}", e.what()); } + StopModel(true); if (status && outcome) status->Finish(outcome->first, outcome->second); } @@ -65,43 +72,33 @@ void StillsGPUAnalysis::Publish() { status->Stills(live); } -bool StillsGPUAnalysis::Hold(const std::shared_ptr &image) { +void StillsGPUAnalysis::Hold(const std::shared_ptr &image) { if (status) status->ImageReceived(std::chrono::steady_clock::now()); - if (!queue.Put(image)) { - // As many frames wait as the buffer is to hold: the pre-pass, if it has not run, runs on what it has. - StartPrePass(); - return false; - } + queue.PutBlocking(image); { std::unique_lock ul(prepass_mutex); - if (!prepass_go) { - if (image->message.number < PREPASS_IMAGES) { - prepass_images[image->message.number] = image; - std::unique_lock ul_live(live_mutex); - live.images_held = static_cast(prepass_images.size()); - } - ++images_held; - // Every one of the first images is there - or, where some never come, twice as many are. - prepass_go = (prepass_images.size() == std::min(PREPASS_IMAGES, experiment.GetImageNum())) - || (images_held >= 2 * PREPASS_IMAGES); + if (!prepass_over && static_cast(prepass_images.size()) < hold_limit) { + prepass_images[image->message.number] = image; + std::unique_lock ul_live(live_mutex); + live.images_held = static_cast(prepass_images.size()); } } prepass_cv.notify_all(); if (status) Publish(); - return true; } void StillsGPUAnalysis::StartPrePass() { { std::unique_lock ul(prepass_mutex); - prepass_go = true; + prepass_last = true; } prepass_cv.notify_all(); } void StillsGPUAnalysis::Finish() { + const auto finish_start = std::chrono::steady_clock::now(); finished = true; StartPrePass(); prepass_future.get(); @@ -124,6 +121,26 @@ void StillsGPUAnalysis::Finish() { std::nullopt); } } + if (model_future.valid()) { + // The model step on the merge of the collection, until the deadline after the last image: the step on a + // provisional merge still running is cut, and this one starts from its placement. + if (merged && model_left_out.empty()) { + const auto deadline = last_image.value_or(finish_start) + std::chrono::duration_cast< + std::chrono::steady_clock::duration>(std::chrono::duration(deadline_s)); + const auto cut = SubmitModel(*merged, "the collection", true); + StopModel(false); + if (deadline_s > 0 && model_future.wait_until(deadline) == std::future_status::timeout) { + *cut = true; + logger.Warning("Model step cut at the deadline, {:.1f} s after the last image", deadline_s); + if (status) + status->Cut("deadline"); + } + model_future.get(); + } else + StopModel(true); + std::unique_lock ul(model_mutex); + model_over = true; + } if (status) { if (setup_error) status->Finish(RugnuxLiveState::Failed, "GPU analysis setup failed"); @@ -151,71 +168,97 @@ void StillsGPUAnalysis::Written() { } void StillsGPUAnalysis::PrePassThread() { - std::map> sample; - { - std::unique_lock ul(prepass_mutex); - prepass_cv.wait(ul, [this] { return prepass_go; }); - sample.swap(prepass_images); - } - const auto masked = [](const PixelMask &m) { return std::count_if(m.GetMask().begin(), m.GetMask().end(), [](uint32_t v) { return v != 0; }); }; setup.experiment = experiment; setup.mask = mask; - { - std::unique_lock ul(live_mutex); - live.prepass = RugnuxLivePrePass::Running; - } - if (status) - Publish(); - bool failed = sample.empty(); - int geometry_frames = 0; - const auto prepass_start = std::chrono::steady_clock::now(); - if (!sample.empty()) { + const int64_t all_images = experiment.GetImageNum(); + std::optional result; + std::chrono::steady_clock::time_point prepass_start = std::chrono::steady_clock::now(); + // At a low hit rate the first frames may not be enough to fit the geometry on, or may not hold a crystal: with + // a cell to fit, the pre-pass runs again on all that is held once that has doubled, until the fit has its + // frames, or it holds hold_limit frames (the last try), or the collection ended. + const int64_t most_held = std::min(all_images, hold_limit); + for (int64_t needed = std::min(PREPASS_IMAGES, most_held); ;) { + std::map> sample; + bool last; + { + std::unique_lock ul(prepass_mutex); + prepass_cv.wait(ul, [&] { + return prepass_last || static_cast(prepass_images.size()) >= needed; + }); + sample = prepass_images; + last = prepass_last || static_cast(sample.size()) >= most_held; + } + if (sample.empty()) + break; + { + std::unique_lock ul(live_mutex); + live.prepass = RugnuxLivePrePass::Running; + } + if (status) + Publish(); try { const auto t0 = std::chrono::steady_clock::now(); - const auto r = RunStillsPrePass(start, sample, static_cast(std::thread::hardware_concurrency())); - setup.experiment.BeamX_pxl(r.beam_x_pxl).BeamY_pxl(r.beam_y_pxl).DetectorDistance_mm(r.distance_mm); - setup.mask = r.mask; - geometry_frames = r.geometry_frames; - const auto &m = r.mask.GetMask(); - std::unique_lock ul(live_mutex); - live.beam_x_pxl = r.beam_x_pxl; - live.beam_y_pxl = r.beam_y_pxl; - live.distance_mm = r.distance_mm; - live.shadow_pixels = std::count_if(m.begin(), m.end(), [](uint32_t v) { - return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; - }); - live.bad_pixels = std::count_if(m.begin(), m.end(), [](uint32_t v) { - return (v & (1u << PixelMask::HotPixelBit)) != 0; - }); + auto r = RunStillsPrePass(start, sample, static_cast(std::thread::hardware_concurrency())); logger.Info("Pre-pass on {} images ({:.1f} s): beam centre ({:.2f}, {:.2f}) -> ({:.2f}, {:.2f}) px, " - "distance {:.2f} -> {:.2f} mm, masked pixels {} -> {}", r.images, + "distance {:.2f} -> {:.2f} mm on {} frames, masked pixels {} -> {}", r.images, std::chrono::duration(std::chrono::steady_clock::now() - t0).count(), experiment.GetBeamX_pxl(), experiment.GetBeamY_pxl(), r.beam_x_pxl, r.beam_y_pxl, - experiment.GetDetectorDistance_mm(), r.distance_mm, + experiment.GetDetectorDistance_mm(), r.distance_mm, r.geometry_frames, masked(mask), masked(r.mask)); + result = std::move(r); } catch (const std::exception &e) { - failed = true; - logger.Warning("Pre-pass failed, the images are analysed with the stated beam centre and mask - {}", - e.what()); + logger.Warning("Pre-pass on {} images failed - {}", sample.size(), e.what()); } + if (last || !experiment.GetUnitCell() + || (result && result->geometry_frames >= RUGNUX_REFINE_GEOMETRY_FRAMES)) + break; + needed = std::min(2 * static_cast(sample.size()), most_held); + std::unique_lock ul(live_mutex); + live.prepass = RugnuxLivePrePass::Holding; + live.images_needed = needed; } - sample.clear(); + { + // The frames are analysed from here on: their slots go back to the buffer once they are sent. + std::unique_lock ul(prepass_mutex); + prepass_over = true; + prepass_images.clear(); + } + const bool failed = !result; + int geometry_frames = 0; + if (result) { + const auto &r = *result; + setup.experiment.BeamX_pxl(r.beam_x_pxl).BeamY_pxl(r.beam_y_pxl).DetectorDistance_mm(r.distance_mm); + setup.mask = r.mask; + geometry_frames = r.geometry_frames; + const auto &m = r.mask.GetMask(); + std::unique_lock ul(live_mutex); + live.beam_x_pxl = r.beam_x_pxl; + live.beam_y_pxl = r.beam_y_pxl; + live.distance_mm = r.distance_mm; + live.shadow_pixels = std::count_if(m.begin(), m.end(), [](uint32_t v) { + return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; + }); + live.bad_pixels = std::count_if(m.begin(), m.end(), [](uint32_t v) { + return (v & (1u << PixelMask::HotPixelBit)) != 0; + }); + } else + logger.Warning("Pre-pass failed, the images are analysed with the stated beam centre and mask"); try { setup.mapping = std::make_unique(setup.experiment, setup.mask); // Retained for the merge, which takes each image's reflections as soon as it is analysed. indexer = std::make_unique(setup.experiment, indexing_thread_pool, merge_frames, true); - // The pre-pass's fit rests on the frames of its sample that indexed; where those are fewer than the - // running refinement bundles, the geometry is refined again, from where the pre-pass left it, on the - // frames as they are analysed. - if (geometry_frames < StillsRunningGeometry::STRONG_FRAMES && setup.experiment.GetUnitCell()) { - running_geometry = std::make_unique(setup.experiment, MERGE_THREADS); - frame_epoch.assign(experiment.GetImageNum(), 0); - geometry_state = GEOMETRY_PENDING; - logger.Info("The pre-pass refined the geometry on {} frames: it is refined again on {} strong frames as " - "they are analysed", geometry_frames, StillsRunningGeometry::STRONG_FRAMES); + if (geometry_frames == 0) { + const std::string text = fmt::format( + "The pre-pass could not refine the geometry (too few frames indexed in what the buffer holds, or no unit cell given): the " + "whole run is analysed at the distance in the header, {:.2f} mm, beam centre ({:.2f}, {:.2f}) px", + setup.experiment.GetDetectorDistance_mm(), setup.experiment.GetBeamX_pxl(), + setup.experiment.GetBeamY_pxl()); + logger.Warning("{}", text); + if (status) + status->Note(text); } if (merge_frames) { merge_experiment = setup.experiment; @@ -224,6 +267,9 @@ void StillsGPUAnalysis::PrePassThread() { setup.experiment.GetUnitCell(), std::nullopt, setup.experiment.IsDetectIceRings(), SpotFindingSettings().ice_ring_width_Q_recipA, logger); + // The model's thread first: the merge thread hands it its merges + if (model) + model_future = std::async(std::launch::async, &StillsGPUAnalysis::ModelThread, this); merge_future = std::async(std::launch::async, &StillsGPUAnalysis::MergeThread, this); } ready(setup); @@ -267,54 +313,32 @@ void StillsGPUAnalysis::AnalysisThread() { return; } - // A thread whose analysis cannot be built still passes its frames on, unanalysed. It is built again at - // the geometry the running refinement switches to. + // A thread whose analysis cannot be built still passes its frames on, unanalysed. std::unique_ptr analysis; - int epoch = -1; - const auto build = [&] { - epoch = geometry_epoch; - analysis.reset(); - try { - analysis = std::make_unique(setup.experiment, *setup.mapping, setup.mask, - *indexer, /*enable_fused_adaptive_gpu=*/true); - } catch (const std::exception &e) { - AnalysisError(std::string("analysis thread could not start, its frames are written unanalysed - ") - + e.what()); - } - }; + try { + analysis = std::make_unique(setup.experiment, *setup.mapping, setup.mask, + *indexer, /*enable_fused_adaptive_gpu=*/true); + } catch (const std::exception &e) { + AnalysisError(std::string("analysis thread could not start, its frames are written unanalysed - ") + + e.what()); + } while (auto held = queue.GetBlocking()) { DataMessage &message = held->message; std::optional profile; - bool switch_geometry = false; const auto start_time = std::chrono::steady_clock::now(); - { - std::unique_lock gate(gate_mutex); - gate_cv.wait(gate, [this] { return !switch_waiting; }); - } - { - std::shared_lock sl(geometry_mutex); - if (epoch != geometry_epoch) - build(); - profile.emplace(*setup.mapping); - try { - if (analysis) - analysis->Analyze(message, *profile, settings()); - } catch (const JFJochException &e) { - // The frame is still written, without its analysis. - logger.ErrorException(e); - } catch (const std::exception &e) { - AnalysisError(fmt::format("image {} written unanalysed - {}", message.number, e.what())); - } - if (geometry_state == GEOMETRY_PENDING) { - frame_epoch.at(message.number) = static_cast(epoch + 1); - switch_geometry = running_geometry->Add(message); - } + profile.emplace(*setup.mapping); + try { + if (analysis) + analysis->Analyze(message, *profile, settings()); + } catch (const JFJochException &e) { + // The frame is still written, without its analysis. + logger.ErrorException(e); + } catch (const std::exception &e) { + AnalysisError(fmt::format("image {} written unanalysed - {}", message.number, e.what())); } message.processing_time_s = std::chrono::duration(std::chrono::steady_clock::now() - start_time).count(); - if (switch_geometry) - SwitchGeometry(); try { sink(message, *profile); if (merge) @@ -337,56 +361,6 @@ void StillsGPUAnalysis::AnalysisThread() { } } -void StillsGPUAnalysis::SwitchGeometry() { - // One fit, in the analysis thread that offered the last frame it needed; the others go on meanwhile. - const auto t0 = std::chrono::steady_clock::now(); - const auto r = running_geometry->Refine(); - const bool last_round = ++geometry_round == GEOMETRY_ROUNDS; - if (!r) { - logger.Warning("Running geometry refinement, round {}: no fit on {} frames ({} strong); the frames are " - "analysed at the geometry they were", geometry_round, running_geometry->Offered(), - running_geometry->Strong()); - geometry_state = geometry_epoch > 0 ? GEOMETRY_SWITCHED : GEOMETRY_KEPT; - return; - } - const std::string text = fmt::format( - "Geometry refined on {} analysed frames ({} strong), round {}: beam centre ({:.2f}, {:.2f}) -> " - "({:.2f}, {:.2f}) px, distance {:.2f} -> {:.2f} mm ({:.1f} s)", running_geometry->Offered(), - running_geometry->Strong(), geometry_round, setup.experiment.GetBeamX_pxl(), - setup.experiment.GetBeamY_pxl(), r->beam_x_px, r->beam_y_px, setup.experiment.GetDetectorDistance_mm(), - r->distance_mm, std::chrono::duration(std::chrono::steady_clock::now() - t0).count()); - { - std::unique_lock gate(gate_mutex); - switch_waiting = true; - } - { - std::unique_lock ul(geometry_mutex); - setup.experiment.BeamX_pxl(static_cast(r->beam_x_px)).BeamY_pxl(static_cast(r->beam_y_px)) - .DetectorDistance_mm(static_cast(r->distance_mm)); - retired_mappings.push_back(std::move(setup.mapping)); - setup.mapping = std::make_unique(setup.experiment, setup.mask); - ++geometry_epoch; - if (last_round) - geometry_state = GEOMETRY_SWITCHED; - else - running_geometry = std::make_unique(setup.experiment, MERGE_THREADS); - } - { - std::unique_lock gate(gate_mutex); - switch_waiting = false; - } - gate_cv.notify_all(); - { - std::unique_lock ul(live_mutex); - live.beam_x_pxl = r->beam_x_px; - live.beam_y_pxl = r->beam_y_px; - live.distance_mm = r->distance_mm; - } - logger.Info("{}", text); - if (status) - status->Note(text); -} - std::optional StillsGPUAnalysis::GetConsensusUnitCell() const { return indexer->GetConsensusUnitCell(); } @@ -411,6 +385,10 @@ const DiffractionExperiment &StillsGPUAnalysis::GetMergeExperiment() const { return merge_experiment; } +const std::optional &StillsGPUAnalysis::GetModel() const { + return model_final; +} + void StillsGPUAnalysis::MergeThread() { constexpr int64_t PREPARE_EVERY = 200; // frames between two Prepare calls after a checkpoint const int64_t images = experiment.GetImageNum(); @@ -423,6 +401,19 @@ void StillsGPUAnalysis::MergeThread() { merge_experiment.SetUnitCell(*c); merge->Cell(*c); cell_set = true; + if (model) { + const StillsModelGroup g = StillsModelSpaceGroup(merge_experiment.GetGemmiSpaceGroup(), *c, *model); + if (!g.space_group) { + model_left_out = "the model is left out: " + g.reason; + logger.Warning("{}", model_left_out); + if (status) + status->Warning(model_left_out); + } else if (!merge_experiment.GetGemmiSpaceGroup()) { + merge_experiment.SetSpaceGroup(*g.space_group); + merge->SpaceGroup(*g.space_group); + logger.Info("Merged in the model's space group {}", g.space_group->xhm()); + } + } } return cell_set; }; @@ -432,38 +423,26 @@ void StillsGPUAnalysis::MergeThread() { merge->Add(n, std::move(o)); if (next < CHECKPOINTS.size() && merge->Added() >= std::llround(CHECKPOINTS[next] * static_cast(images))) { - if (have_cell()) - PublishMerge(merge->Checkpoint(MERGE_THREADS), - fmt::format("{:.0f}% of the images (provisional)", 100.0 * CHECKPOINTS[next])); + if (have_cell()) { + StillsMergeOutput c = merge->Checkpoint(MERGE_THREADS); + const std::string title = fmt::format("{:.0f}% of the images (provisional)", + 100.0 * CHECKPOINTS[next]); + PublishMerge(c, title); + if (model_future.valid() && model_left_out.empty()) + SubmitModel(c, title, false); + } ++next; } else if (next > 0 && merge->Added() % PREPARE_EVERY == 0) merge->Prepare(MERGE_THREADS); }; - // The frames analysed at a geometry the running refinement may still replace wait until it is decided: - // left out where analysed before its last switch, merged as analysed otherwise. - std::vector held_back; - const auto release = [&] { - for (const int64_t h: held_back) - if (frame_epoch[h] - 1 < geometry_epoch) { - indexer->DropIntegrationOutcome(h); - add(h, IntegrationOutcome{}); - } else - add(h, indexer->TakeIntegrationOutcome(h)); - held_back.clear(); - }; - for (int64_t n = merge_queue.GetBlocking(); n >= 0; n = merge_queue.GetBlocking()) { - if (!frame_epoch.empty() && frame_epoch[n] > 0) - held_back.push_back(n); - else - add(n, indexer->TakeIntegrationOutcome(n)); - if (geometry_state != GEOMETRY_PENDING && !held_back.empty()) - release(); - } - release(); + for (int64_t n = merge_queue.GetBlocking(); n >= 0; n = merge_queue.GetBlocking()) + add(n, indexer->TakeIntegrationOutcome(n)); if (have_cell()) { merged = merge->Final(static_cast(AvailableCpus())); logger.Info("Stills merge after the last image: {:.2f} s", merged->seconds); } + // The frames' reflections are not read again: what the model step and the end message need is the merge + merge.reset(); } catch (const std::exception &e) { logger.Warning("Stills merge failed: {}", e.what()); merged.reset(); @@ -482,3 +461,194 @@ void StillsGPUAnalysis::PublishMerge(const StillsMergeOutput &out, const std::st if (status) status->Note(text); } + +std::shared_ptr> StillsGPUAnalysis::SubmitModel(StillsMergeOutput &out, const std::string &title, + bool final) { + ModelJob job; + // The merge of the collection is kept for the end message + job.merged = final ? out.merge.merged : std::move(out.merge.merged); + for (const auto &sh : out.merge.statistics.shells) + job.shells.push_back(sh.d_min); + job.title = title; + job.final = final; + const auto cut = job.cut; + { + std::unique_lock ul(model_mutex); + if (model_cut) + *job.cut = true; + if (model_running_cut) + *model_running_cut = true; + if (model_next) + logger.Info("Model step on the merge of {} left out: the merge of {} is there", model_next->title, title); + model_next = std::move(job); + } + model_cv.notify_all(); + return cut; +} + +void StillsGPUAnalysis::Cut(const std::string &reason) { + if (!model) + return; + { + std::unique_lock ul(model_mutex); + if (model_cut || model_over) + return; + model_cut = true; + if (model_running_cut) + *model_running_cut = true; + if (model_next) + *model_next->cut = true; + } + logger.Info("Model step cut: {}", reason); + if (status) + status->Cut(reason); +} + +void StillsGPUAnalysis::StopModel(bool cut) { + if (!model_future.valid()) + return; + { + std::unique_lock ul(model_mutex); + model_stop = true; + if (cut) { + model_next.reset(); + if (model_running_cut) + *model_running_cut = true; + } + } + model_cv.notify_all(); + if (cut) + model_future.get(); +} + +void StillsGPUAnalysis::ModelThread() { + try { + pin_gpu(); + } catch (const std::exception &e) { + logger.Warning("Error pinning GPU {}", e.what()); + } + while (true) { + ModelJob job; + { + std::unique_lock ul(model_mutex); + model_cv.wait(ul, [this] { return model_next || model_stop; }); + if (!model_next) + return; + job = std::move(*model_next); + model_next.reset(); + model_running_cut = job.cut; + } + try { + StillsModelSettings settings; + // Beside the analysis on as many threads as a provisional merge, on every thread after the last image + settings.nthreads = job.final ? AvailableCpus() : MERGE_THREADS; + // The analysis's indexers stay on the card for the broker's life + settings.gpu_resident_bytes = indexing_thread_pool ? indexing_thread_pool->DeviceBytes() : 0; + settings.cut = job.cut.get(); + if (job.final && status) + settings.step = [this](const std::string &title) { status->Step(title); }; + const DiffractionExperiment &x = merge_experiment; + StillsModelResult r = PlaceAndRefineStillsModel(job.merged, x.GetUnitCell().value(), x.GetSpaceGroupOrP1(), + *model, x.GetWavelength_A(), + x.GetScalingSettings().GetMergeFriedel(), job.shells, + settings, model_placed ? &*model_placed : nullptr, logger); + PublishModel(r, job.title); + if (r.placed_model) + model_placed = r; + if (job.final) + model_final = std::move(r); + } catch (const std::exception &e) { + logger.Warning("Model step on the merge of {} failed: {}", job.title, e.what()); + } + } +} + +void StillsGPUAnalysis::PublishModel(const StillsModelResult &r, const std::string &title) { + std::string text = fmt::format("Model on the merge of {}: ", title); + if (r.placement == "NOT_PLACED") + text += "not placed"; + else + text += fmt::format("placed by {}{}", r.placement == "MR" ? "molecular replacement" : "its rigid body", + r.warm ? " (from the earlier placement)" : ""); + if (r.validation && r.validation->ok) + text += fmt::format(", R-work {:.3f}, R-free {:.3f}", r.validation->r_work, r.validation->r_free); + if (r.refinement && r.refinement->ok) + text += fmt::format("; refined{} R-work {:.3f}, R-free {:.3f}, {} blobs", r.refinement->cancelled + ? " (cut short)" : "", r.refinement->r_work, r.refinement->r_free, + r.refinement->blobs.size()); + for (const auto &l : r.left_out) + text += "; " + l + " left out"; + text += fmt::format(" ({:.1f} s placement, {:.1f} s refinement)", r.placement_s, r.refinement_s); + logger.Info("{}", text); + if (!status) + return; + status->Note(text); + if (const auto m = LiveModel(StillsModelReport(r))) + status->Model(*m); + if (!r.placed_model) + return; + std::optional coefficients; + if (r.maps) + coefficients = r.maps->coefficients; + std::shared_ptr refined; + if (r.refinement && r.refinement->ok && r.refinement->model) { + std::optional refined_coefficients; + if (r.refinement->maps) + refined_coefficients = r.refinement->maps->coefficients; + refined = std::make_shared(r.refinement->model, r.placed_model->cell, + r.placed_model->space_group, std::move(refined_coefficients)); + } + status->Structure(std::make_shared(r.placed_model->model, r.placed_model->cell, + r.placed_model->space_group, std::move(coefficients)), + std::move(refined)); +} + +void StillsGPUAnalysis::FillEndMessage(EndMessage &message) const { + message.unit_cell = GetConsensusUnitCell(); + message.refined_beam_center_x = GetBeamX_pxl(); + message.refined_beam_center_y = GetBeamY_pxl(); + message.refined_detector_distance_mm = GetDetectorDistance_mm(); + // The merged data, for the writer to put beside the master file: in the model's frame where a model was + // placed against them, with its maps, the model as placed and as refined, and the refined model's blobs. + if (const auto &mx = GetModel()) { + DiffractionExperiment x(GetMergeExperiment()); + x.SetSpaceGroup(mx->space_group); + message.merged = MergedMtzTable(mx->merged, mx->cell, x, OutputProvenance{}); + if (mx->maps) + message.map_coefficients = mx->maps->coefficients; + if (mx->placed_model) + message.placed_model = PlacedModelFromStructure(*mx->placed_model->model, mx->placed_model->cell, + mx->placed_model->space_group.xhm()); + if (const auto &rf = mx->refinement; rf && rf->ok && rf->model) { + message.refined_model = PlacedModelFromStructure(*rf->model, mx->cell, mx->space_group.xhm()); + if (rf->maps) + message.refined_map_coefficients = rf->maps->coefficients; + message.blobs = DensityBlobsToMessage(rf->blobs); + } + } else if (const auto &m = GetMerged(); m && GetMergeExperiment().GetUnitCell()) + message.merged = MergedMtzTable(m->merge.merged, *GetMergeExperiment().GetUnitCell(), + GetMergeExperiment(), OutputProvenance{}); + if (!message.merged) + return; + // Every pusher serialises END into a buffer of MESSAGE_SIZE_FOR_START_END, and an END that does not fit there + // fails the collection. So it is serialised here first, as JFJochReceiverRugnux does for rotation; one that does not fit + // goes without the merged data, map coefficients and models. + const auto scratch = std::make_unique_for_overwrite(MESSAGE_SIZE_FOR_START_END); + CBORStream2Serializer serializer(scratch.get(), MESSAGE_SIZE_FOR_START_END); + try { + serializer.SerializeSequenceEnd(message); + } catch (const JFJochException &e) { + const std::string text = fmt::format("the merged data, map coefficients and models do not fit in the end " + "message ({} MiB) and are not written", + MESSAGE_SIZE_FOR_START_END / (1024 * 1024)); + logger.Error("GPU analysis: {} - {}", text, e.what()); + if (status) + status->Warning(text); + message.merged.reset(); + message.map_coefficients.reset(); + message.placed_model.reset(); + message.refined_model.reset(); + message.refined_map_coefficients.reset(); + message.blobs.reset(); + } +} diff --git a/receiver/StillsGPUAnalysis.h b/receiver/StillsGPUAnalysis.h index d6350127f..00dd45f2a 100644 --- a/receiver/StillsGPUAnalysis.h +++ b/receiver/StillsGPUAnalysis.h @@ -4,7 +4,6 @@ #pragma once #include -#include #include #include #include @@ -13,7 +12,6 @@ #include #include #include -#include #include #include "StillsPrePass.h" @@ -26,15 +24,16 @@ #include "../common/ThreadSafeFIFO.h" #include "../image_analysis/IndexAndRefine.h" #include "../rugnux/StillsLiveMerge.h" -#include "../rugnux/StillsRunningGeometry.h" +#include "../rugnux/StillsModel.h" -// Serial frames analysed on the GPU as the DECTRIS ones are (MXAnalysisWithoutFPGA), after a pre-pass on +// Serial frames analysed on the GPU (MXAnalysisWithoutFPGA), after a pre-pass on // the first PREPASS_IMAGES (RunStillsPrePass) has decided the beam centre and the mask they are analysed // with. Every frame given before the pre-pass is done is held; `ready` is called once it is, with the // geometry and mask the frames are then analysed at, before any of them is analysed; `sink` gets each // frame once analysed, from one of the analysis threads. Where `status` is given, the run is reported there // (pre-pass, its decisions, the analysis). -// At most `max_held` frames wait for analysis: Hold refuses one more, so the slots held stay bounded. +// Every frame given is taken and analysed: the slots of the image buffer are what bounds them. The pre-pass holds +// and fits on at most `max_held` of them; those arriving beyond that, or while it runs, wait for the analysis. // Where the setup after the pre-pass fails (`ready` included), `failed` is called at once, the frames are // dropped unanalysed and Finish throws the error. An analysis thread that fails passes its frames on // unanalysed, and the run is reported as failed. @@ -44,18 +43,17 @@ // before Finish returns (GetMerged). The cell is the run's consensus at the first merge, as rugnux merges in // it (the collection's cell picks the frames merged); the space group the collection's, else P1. // merge = false leaves it out. -// Where the pre-pass refined the centre and distance on fewer frames than StillsRunningGeometry bundles, or -// could not refine them (too few of the frames it had indexed: a small image buffer, a serial hit rate, a chip -// whose first frames are empty), and a cell is given, they are refined again while the frames are analysed -// (StillsRunningGeometry, from the spots and lattices of the frames analysed, no image held), and the analysis -// switches to that geometry as each of GEOMETRY_ROUNDS fits is made. The frames analysed before the last -// switch are not merged (their reflections were predicted at an earlier geometry; the images are gone); where -// no fit can be made they are merged as analysed. +// With a cell given, the pre-pass runs again on all the frames held each time they have doubled, until the geometry +// fit has its frames, `max_held` frames are held or the collection ends; the analysis starts after that and +// catches up from the frames held. The geometry is fixed when the pre-pass ends: every frame is analysed and merged at it. Where the pre-pass could +// not refine it (too few of the frames it had indexed: a small image buffer, a serial hit rate, a chip whose first +// frames are empty; or no cell given), the header's distance stays for the whole run, with a warning in the log +// and the status. class StillsGPUAnalysis { public: // Enough frames for the distance to be refined at a serial hit rate: of 100 frames of a collection indexed // at 18%, 7 indexed at a distance 2.6% off and the geometry refinement did not run; of 500, 34 did and - // refined it to 0.4%. They wait in the image buffer, which has to hold them in half its slots. + // refined it to 0.4%. They wait in the image buffer. constexpr static int64_t PREPASS_IMAGES = 500; constexpr static std::array CHECKPOINTS = {0.5, 0.75}; constexpr static int MERGE_THREADS = 8; @@ -74,17 +72,19 @@ public: StillsGPUAnalysis(const DiffractionExperiment &experiment, const PixelMask &mask, const StartMessage &start, IndexerThreadPool *indexing_thread_pool, int64_t nthreads, int64_t max_held, Logger &logger, Ready ready, Sink sink, Settings settings, Failed failed = {}, - RugnuxLiveStatus *status = nullptr, bool merge = true); + RugnuxLiveStatus *status = nullptr, bool merge = true, + std::shared_ptr model = nullptr, double deadline_s = 0); ~StillsGPUAnalysis(); - // False where `max_held` frames already wait: the frame is not taken. - bool Hold(const std::shared_ptr &image); + void Hold(const std::shared_ptr &image); // No frame comes after this: the pre-pass runs on what came, if it has not, and every frame is // analysed before it returns. Throws where the setup after the pre-pass failed. void Finish(); // After Finish, once the end message is out and the writer is done: the status gets the latency from the last // image and then the run's outcome. Without it the outcome is reported when the analysis goes, with no latency. void Written(); + // POST /cancel: the model step stops where it is, as at the deadline, and the frames are still analysed and merged. + void Cut(const std::string &reason); // After Finish. [[nodiscard]] std::optional GetConsensusUnitCell() const; @@ -94,6 +94,11 @@ public: // The merge of the collection, where one was made. [[nodiscard]] const std::optional &GetMerged() const; [[nodiscard]] const DiffractionExperiment &GetMergeExperiment() const; + // The model step on the merge of the collection, where one was made. + [[nodiscard]] const std::optional &GetModel() const; + // What the analysis found, for the collection's end message: the cell, geometry, and the merged data for the + // writer to put beside the master file, in the model's frame where a model was placed against them. + void FillEndMessage(EndMessage &message) const; private: const DiffractionExperiment &experiment; @@ -121,8 +126,8 @@ private: std::mutex prepass_mutex; std::condition_variable prepass_cv; - bool prepass_go = false; - int64_t images_held = 0; + bool prepass_last = false; // what is held is all there will be + bool prepass_over = false; std::map> prepass_images; std::latch prepass_done{1}; std::exception_ptr setup_error; // set before prepass_done @@ -131,6 +136,7 @@ private: std::string analysis_error; // the first, where an analysis thread failed void AnalysisError(const std::string &text); + const int64_t hold_limit; // the frames the pre-pass holds and fits on, at most ThreadSafeFIFO> queue; std::future prepass_future; std::vector> analysis_futures; @@ -143,28 +149,38 @@ private: std::optional merged; std::future merge_future; - // The geometry refined while the frames are analysed (see above), in GEOMETRY_ROUNDS rounds, each on the - // frames analysed at the last one's geometry (the spots of a frame are indexed at the geometry it was - // analysed at, and a start 2.6% off left the first fit 0.6% off). The analysis takes geometry_mutex shared - // for each frame; a switch takes it exclusive and raises geometry_epoch, which the analysis threads see to - // rebuild their analysis. geometry_state is PENDING until the last round, then SWITCHED or KEPT. - // frame_epoch: per frame number, 1 + the epoch it was analysed at while PENDING (0: after); the merge - // leaves out those analysed before the last switch. - constexpr static int GEOMETRY_ROUNDS = 2; - enum GeometryState : int { GEOMETRY_PENDING = 0, GEOMETRY_SWITCHED = 1, GEOMETRY_KEPT = 2 }; - std::unique_ptr running_geometry; - int geometry_round = 0; - std::shared_mutex geometry_mutex; - // A switch waiting for geometry_mutex closes this gate to new frames, or the frames that keep taking the - // mutex shared would hold it off to the end of the collection. - std::mutex gate_mutex; - std::condition_variable gate_cv; - bool switch_waiting = false; - std::atomic geometry_epoch{0}; - std::atomic geometry_state{GEOMETRY_KEPT}; - std::vector frame_epoch; - std::vector> retired_mappings; // the analyses before the switch read them - void SwitchGeometry(); + // Where a model is given, it is placed against each merge, mapped and refined (PlaceAndRefineStillsModel), on a + // thread of its own: the provisional merges' beside the analysis, published to the status as they are made, and the + // collection's in Finish, which waits for it until deadline_s after the last image (0: no deadline) and then cuts it - + // a refinement running is cut short, a molecular replacement or a refinement not started is left out. A step still + // running when the next merge is made is cut the same way, and one not started is replaced by it. Each step starts + // from the placement of the last one that placed the model. Without a space group stated, the frames are merged in + // the model's where its cell is the run's consensus (StillsModelSpaceGroup); otherwise the model is left out. + const std::shared_ptr model; + const double deadline_s; + std::string model_left_out; // why the model is left out; set by the merge thread before any step + struct ModelJob { + std::vector merged; + std::vector shells; // the merge statistics' shell bounds + std::string title; + bool final = false; + std::shared_ptr> cut = std::make_shared>(false); + }; + std::mutex model_mutex; + std::condition_variable model_cv; + std::optional model_next; + std::shared_ptr> model_running_cut; // of the step running + bool model_stop = false; + bool model_cut = false; // the steps from here on are cut too + bool model_over = false; // the step on the merge of the collection has ended: nothing is left to cut + std::optional model_placed; // the last step that placed the model, for the next to start from + std::optional model_final; + std::future model_future; + // The step on `out`, the merge of `title` (its reflections moved out, but for the final one); returns its cut. + std::shared_ptr> SubmitModel(StillsMergeOutput &out, const std::string &title, bool final); + void StopModel(bool cut); + void ModelThread(); + void PublishModel(const StillsModelResult &r, const std::string &title); void StartPrePass(); void PrePassThread(); diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 72ab81c3d..748515e52 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -25,8 +25,8 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxLive.h StillsLiveMerge.cpp StillsLiveMerge.h - StillsRunningGeometry.cpp - StillsRunningGeometry.h + StillsModel.cpp + StillsModel.h RugnuxRefine.cpp RugnuxRefine.h Rugnux.h diff --git a/rugnux/StillsLiveMerge.cpp b/rugnux/StillsLiveMerge.cpp index 46d317dc0..a84f7afbe 100644 --- a/rugnux/StillsLiveMerge.cpp +++ b/rugnux/StillsLiveMerge.cpp @@ -87,6 +87,10 @@ void StillsLiveMerge::Cell(const UnitCell &cell) { experiment.SetUnitCell(cell); } +void StillsLiveMerge::SpaceGroup(const gemmi::SpaceGroup &space_group) { + experiment.SetSpaceGroup(space_group); +} + void StillsLiveMerge::Add(int64_t image, IntegrationOutcome &&outcome) { ++added; if (outcome.reflections.empty()) diff --git a/rugnux/StillsLiveMerge.h b/rugnux/StillsLiveMerge.h index cf20ba9e9..ac7ec5169 100644 --- a/rugnux/StillsLiveMerge.h +++ b/rugnux/StillsLiveMerge.h @@ -49,6 +49,8 @@ public: // The cell merged in, where the experiment had none; before the first checkpoint. void Cell(const UnitCell &cell); + // The group merged in, where the experiment had none; before the first checkpoint. + void SpaceGroup(const gemmi::SpaceGroup &space_group); // Frame `image` as analysed: without reflections where it was not indexed. void Add(int64_t image, IntegrationOutcome &&outcome); // Frames added so far diff --git a/rugnux/StillsModel.cpp b/rugnux/StillsModel.cpp new file mode 100644 index 000000000..9b3077679 --- /dev/null +++ b/rugnux/StillsModel.cpp @@ -0,0 +1,192 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "StillsModel.h" + +#include +#include + +#include "spdlog/fmt/fmt.h" + +#include "../image_analysis/scale_merge/CrystalSetting.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" + +namespace { + double SecondsSince(std::chrono::steady_clock::time_point t0) { + return std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + } + + bool Cut(const StillsModelSettings &settings) { + return settings.cut != nullptr && settings.cut->load(); + } + + void Step(const StillsModelSettings &settings, const std::string &title) { + if (settings.step) + settings.step(title); + } +} + +StillsModelGroup StillsModelSpaceGroup(const std::optional &collection, const UnitCell &consensus, + const RugnuxModel &model) { + if (collection) + return {collection, ""}; + const gemmi::Structure &st = model.structure; + const gemmi::SpaceGroup *sg = st.find_spacegroup(); + if (sg == nullptr || !st.cell.is_crystal()) + return {std::nullopt, "the collection states no space group, and the model gives no crystal's"}; + if (!CellsCorrespond(st.cell, gemmi::UnitCell(consensus))) + return {std::nullopt, fmt::format("the collection states no space group, and the model's cell {:.2f} {:.2f} " + "{:.2f} {:.1f} {:.1f} {:.1f} ({}) is not the run's {:.2f} {:.2f} {:.2f} " + "{:.1f} {:.1f} {:.1f}", st.cell.a, st.cell.b, st.cell.c, st.cell.alpha, + st.cell.beta, st.cell.gamma, sg->xhm(), consensus.a, consensus.b, + consensus.c, consensus.alpha, consensus.beta, consensus.gamma)}; + return {*sg, ""}; +} + +StillsModelResult PlaceAndRefineStillsModel(const std::vector &merged, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const RugnuxModel &model, + double wavelength_A, bool merge_friedel, + const std::vector &report_shell_d_min, + const StillsModelSettings &settings, const StillsModelResult *earlier, + Logger &logger) { + const auto t0 = std::chrono::steady_clock::now(); + StillsModelResult out; + out.merged = merged; + out.cell = cell; + out.space_group = space_group; + ModelValidationSchedule placed_elsewhere; // a placement already judged: no null + placed_elsewhere.fit_test = false; + placed_elsewhere.gpu_resident_bytes = settings.gpu_resident_bytes; + + // Warm: in the earlier step's frame, from where it placed the model + if (earlier && earlier->placed_model && earlier->validation) { + const ReindexChoice c = ChooseReindex(out.merged, cell, space_group, + [&](const std::vector &m) { + const double cc = ReferenceIntensityCC(m, earlier->merged, space_group); + return std::isfinite(cc) ? cc : -1.0; + }); + std::vector moved = c.is_identity + ? out.merged + : ReindexMergedIntoAsu(out.merged, c.op, space_group, merge_friedel); + Step(settings, "Model validation"); + ModelValidationResult v = ValidateAgainstModel(moved, cell, model, logger, &earlier->space_group, + /*probe_indexing_ambiguity=*/false, settings.nthreads, + wavelength_A, report_shell_d_min, placed_elsewhere, nullptr, + earlier->placed_model->model.get()); + if (v.ok && v.placed_model) { + logger.Info("Model: placed as on the earlier merge{} - R-work {:.3f}, R-free {:.3f}", + c.is_identity ? "" : ", the merge reindexed into its frame", v.r_work, v.r_free); + KeepModelVerdict(v, *earlier->validation); + out.merged = std::move(moved); + out.space_group = earlier->space_group; + out.placement = earlier->placement; + out.mr = earlier->mr; + out.warm = true; + out.maps = std::move(v.maps); + out.placed_model = PlacedModelOutput{v.placed_model, cell, out.space_group}; + out.validation = std::move(v); + } else + logger.Warning("Model: the earlier placement does not validate on this merge ({}) - placed again", + v.failure_reason); + } + + if (!out.warm) { + const bool predicted = PrepareSearchModel(model).predicted; + // An experimental structure as given, by its rigid body; a prediction has no placement of its own + if (!predicted) { + Step(settings, "Model validation"); + ModelValidationSchedule schedule; + schedule.gpu_resident_bytes = settings.gpu_resident_bytes; + ModelValidationResult v = ValidateAgainstModel(out.merged, cell, model, logger, &space_group, + /*probe_indexing_ambiguity=*/true, settings.nthreads, + wavelength_A, report_shell_d_min, schedule); + if (v.ok && v.model_fits) { + out.space_group = *AdoptModelFrame(v, out.merged, space_group, merge_friedel, logger); + out.placement = "RIGID_BODY"; + out.maps = std::move(v.maps); + out.placed_model = PlacedModelOutput{v.placed_model, cell, out.space_group}; + } + v.maps.reset(); + out.validation = std::move(v); + } + // Molecular replacement, where the rigid body did not place it + if (out.placement == "NOT_PLACED" && Cut(settings)) + out.left_out.emplace_back("molecular replacement"); + else if (out.placement == "NOT_PLACED") { + Step(settings, "Molecular replacement"); + MolecularReplacementResult mr = RunMolecularReplacement( + model, {{"MERGE", &out.merged, cell, space_group}}, logger, settings.nthreads, + settings.gpu_resident_bytes); + if (mr.chosen >= 0) { + const MRCandidate &c = mr.candidates[mr.chosen]; + if (mr.decision == "ENANTIOMORPH") + CheckEnantiomorphHand(mr, out.merged, space_group, logger, settings.nthreads, + settings.gpu_resident_bytes); + // The merge labelled in the group the model chose, where it chose one of those the merge cannot + // tell apart (the enantiomorph past its anomalous veto, or another screw) + if (mr.decision == "ENANTIOMORPH" || mr.decision == "SPACE_GROUP") { + out.space_group = c.space_group; + const gemmi::GroupOps ops = c.space_group.operations(); + std::erase_if(out.merged, [&](const MergedReflection &r) { + return ops.is_systematically_absent({{r.h, r.k, r.l}}); + }); + } + // Its maps and R-factors: the placement validated against the merge, nothing left to decide + ModelValidationResult v = ValidateAgainstModel(out.merged, cell, model, logger, &out.space_group, + /*probe_indexing_ambiguity=*/false, settings.nthreads, + wavelength_A, report_shell_d_min, placed_elsewhere, + nullptr, c.placed.get()); + if (v.ok && v.placed_model) { + if (out.validation) { + // the rigid body's verdict stays + v.fit_tested = out.validation->fit_tested; + v.model_fits = out.validation->model_fits; + v.null_r_work_mean = out.validation->null_r_work_mean; + v.null_r_work_sd = out.validation->null_r_work_sd; + v.r_work_sigma = out.validation->r_work_sigma; + v.null_replicates = out.validation->null_replicates; + } + out.placement = "MR"; + out.maps = std::move(v.maps); + out.placed_model = PlacedModelOutput{v.placed_model, cell, out.space_group}; + out.validation = std::move(v); + } + } + out.mr = std::move(mr); + } + } + out.placement_s = SecondsSince(t0); + + if (!settings.refine || !out.placed_model) + return out; + if (Cut(settings)) { + out.left_out.emplace_back("refinement"); + return out; + } + Step(settings, "Refinement"); + const auto t1 = std::chrono::steady_clock::now(); + RefinementSettings rs; + rs.nthreads = settings.nthreads; + rs.gpu_resident_bytes = settings.gpu_resident_bytes; + rs.cancel = settings.cut; + // The free set is drawn by the merge: a model placed as given has most likely been refined against all of it + rs.as_given_competes = out.placement != "RIGID_BODY"; + try { + out.refinement = RefineModel(*out.placed_model->model, out.merged, gemmi::UnitCell(cell), out.space_group, + wavelength_A, logger, rs); + } catch (const std::exception &e) { + out.refinement = RefinementResult(); + out.refinement->failure_reason = std::string("refinement failed: ") + e.what(); + } + out.refinement_s = SecondsSince(t1); + return out; +} + +ProcessResult StillsModelReport(const StillsModelResult &result) { + ProcessResult r; + r.model_validation = result.validation; + r.model_prior = result.mr; + r.model_placement = result.placement; + r.refinement = result.refinement; + return r; +} diff --git a/rugnux/StillsModel.h b/rugnux/StillsModel.h new file mode 100644 index 000000000..56c50ea95 --- /dev/null +++ b/rugnux/StillsModel.h @@ -0,0 +1,76 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include +#include + +#include "Rugnux.h" +#include "../image_analysis/structure_refinement/MolecularReplacement.h" +#include "../image_analysis/structure_refinement/Refinement.h" + +// The model step of a merge of serial frames made in memory (StillsLiveMerge, MergeStills): the model placed +// against the merge, its maps and its refinement, as `rugnux --mode scale --model` places it and `rugnux --mode +// refine` refines it. For the live stills of the broker, which makes it at each provisional merge and after the +// last frame. No file is read or written here. + +// The group a merge of serial frames is merged in where a model is given: the collection's, where it states +// one; else the model's own, where the model's cell describes the run's consensus cell on the same axes +// (CellsCorrespond); else none, and `reason` says why the model is left out. +struct StillsModelGroup { + std::optional space_group; + std::string reason; +}; +StillsModelGroup StillsModelSpaceGroup(const std::optional &collection, const UnitCell &consensus, + const RugnuxModel &model); + +struct StillsModelSettings { + size_t nthreads = 1; + size_t gpu_resident_bytes = 0; // held on the card meanwhile (the analysis's indexers) + bool refine = true; + // Set from any thread: the step is cut - a refinement that runs is cut short, keeping the trials that ran + // (RefinementSettings::cancel), and a molecular replacement or a refinement that has not started is left out. + const std::atomic *cut = nullptr; + // Called as each part starts: "Model validation", "Molecular replacement", "Refinement" + std::function step; +}; + +struct StillsModelResult { + // The merge in the model's frame - the indexing the model chose (AdoptModelFrame), labelled in the group + // the model was placed in (its enantiomorph, where it was taken) - and its cell. + std::vector merged; + UnitCell cell{}; + gemmi::SpaceGroup space_group; + std::string placement = "NOT_PLACED"; // RIGID_BODY | MR | NOT_PLACED + bool warm = false; // placed from an earlier step's placement + std::optional validation; // its maps moved to `maps` + std::optional mr; + // Where the model was placed: the maps of the placement and the model as placed, in the merge's frame + std::shared_ptr maps; + std::optional placed_model; + std::optional refinement; // its model, maps and blobs where it ran + std::vector left_out; // what a cut left out + double placement_s = 0, refinement_s = 0; +}; + +// The model placed against `merged` (in `space_group`, the cell `cell`): by its rigid body where it fits these +// data as given, its indexing ambiguity probed and the merge put in the model's indexing; by molecular +// replacement where it does not (or is a prediction); then refined, where placed. `earlier`, the result of a +// step on an earlier merge of the same collection, is the warm start: the merge is put in its frame (the +// indexing a later merge happened to take is decided against it, ChooseReindex), the model goes in where it was +// placed there, and the verdict, the indexing and the molecular replacement's choice are its - no null, probe +// or search is made again. Where the warm placement does not validate, the step is made cold. +StillsModelResult PlaceAndRefineStillsModel(const std::vector &merged, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const RugnuxModel &model, + double wavelength_A, bool merge_friedel, + const std::vector &report_shell_d_min, + const StillsModelSettings &settings, const StillsModelResult *earlier, + Logger &logger); + +// The step's outcome in the terms of a ProcessResult (model_validation, model_prior, model_placement and +// refinement), for what reports one (LiveModel). +ProcessResult StillsModelReport(const StillsModelResult &result); diff --git a/rugnux/StillsRunningGeometry.cpp b/rugnux/StillsRunningGeometry.cpp deleted file mode 100644 index e5be725dc..000000000 --- a/rugnux/StillsRunningGeometry.cpp +++ /dev/null @@ -1,95 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "StillsRunningGeometry.h" - -#include - -#include "../image_analysis/scale_merge/CellPopulation.h" - -namespace { - // A frame the geometry refinement can use, as rugnux takes it from its sample: the spots on lattice 0 with - // their index, the others without. - GeomRefineFrame RefineFrame(const DataMessage &msg) { - GeomRefineFrame f; - f.lattice = *msg.indexing_lattice; - f.ordinal = static_cast(msg.number); - for (const auto &s: msg.spots) { - if (s.indexed && s.lattice == 0) - f.spots.push_back(GeomRefineSpot{s.x, s.y, static_cast(s.h), static_cast(s.k), - static_cast(s.l), s.panel}); - else - f.unindexed.push_back(GeomRefineSpot{s.x, s.y, 0, 0, 0, s.panel}); - } - return f; - } -} - -StillsRunningGeometry::StillsRunningGeometry(const DiffractionExperiment &in_experiment, int in_nthreads) - : experiment(in_experiment), nthreads(in_nthreads) {} - -bool StillsRunningGeometry::Add(const DataMessage &msg) { - std::optional f; - if (msg.indexing_result.value_or(false) && msg.indexing_lattice.has_value()) - f = RefineFrame(msg); - std::unique_lock ul(m); - if (done) - return false; - ++offered; - if (f && f->spots.size() >= 6) { - ++indexed; - strong += static_cast(f->spots.size()) >= MIN_STRONG_SPOTS; - frames.push_back(std::move(*f)); - } - done = static_cast(strong) >= STRONG_FRAMES || indexed >= INDEXED_BUDGET; - return done; -} - -int64_t StillsRunningGeometry::Offered() const { - std::unique_lock ul(m); - return offered; -} - -size_t StillsRunningGeometry::Strong() const { - std::unique_lock ul(m); - return strong; -} - -std::optional StillsRunningGeometry::Refine() { - std::unique_lock ul(m); - const auto cell = experiment.GetUnitCell(); - if (!cell) - return {}; - // As RugnuxStills.cpp's sample: only the frames of the given cell's population are fitted, the population - // told from the strong frames. - const double given = CrystalLattice(*cell).CalcVolume(); - std::vector strong_cells; - for (const auto &f: frames) - if (static_cast(f.spots.size()) >= MIN_STRONG_SPOTS) - strong_cells.push_back(f.lattice.GetUnitCell()); - if (strong_cells.empty()) - return {}; - const auto [low, high] = GivenCellPopulation(strong_cells, *cell); - std::vector bundle; - for (auto &f: frames) { - const double r = f.lattice.CalcVolume() / given - 1.0; - if (static_cast(f.spots.size()) >= MIN_STRONG_SPOTS && r >= low && r <= high) - bundle.push_back(std::move(f)); - } - frames.clear(); - std::sort(bundle.begin(), bundle.end(), [](const GeomRefineFrame &a, const GeomRefineFrame &b) { - if (a.spots.size() != b.spots.size()) return a.spots.size() > b.spots.size(); - return a.ordinal < b.ordinal; - }); - if (static_cast(bundle.size()) > STRONG_FRAMES) - bundle.resize(STRONG_FRAMES); - - GeometryRefinerSettings settings; - settings.crystal_system = experiment.GetCrystalSystem(); - settings.num_threads = nthreads; - settings.cell_reg_coeff = 0.0; // cell and distance both free, as rugnux fits them - const GeometryRefinerResult r = RefineGlobalGeometry(experiment.GetDiffractionGeometry(), *cell, bundle, settings); - if (!r.ok) - return {}; - return r; -} diff --git a/rugnux/StillsRunningGeometry.h b/rugnux/StillsRunningGeometry.h deleted file mode 100644 index ad3260673..000000000 --- a/rugnux/StillsRunningGeometry.h +++ /dev/null @@ -1,47 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include - -#include "../common/DiffractionExperiment.h" -#include "../common/JFJochMessages.h" -#include "../image_analysis/geom_refinement/GeometryRefiner.h" - -// The detector centre and distance of a serial collection refined from frames as they are analysed, holding -// no image: the live path's geometry where its pre-pass could not refine it (too few of the frames the image -// buffer holds indexed). The frames' indexed spots and lattices are kept until STRONG_FRAMES strong frames -// (at least MIN_STRONG_SPOTS indexed spots) are in, or INDEXED_BUDGET frames indexed (empty frames - a chip's -// wells without a crystal - do not count); then the strongest -// are bundled and the geometry refined once, as rugnux refines it on its sample (cell and distance free, the -// detector untilted), on the frames of the given cell's population. A fixed cost: one bundle, one fit. -// Thread-safe. -class StillsRunningGeometry { -public: - constexpr static int STRONG_FRAMES = 200; - constexpr static int MIN_STRONG_SPOTS = 20; - constexpr static int64_t INDEXED_BUDGET = 5000; - - // `experiment`: the geometry the frames are analysed at, the cell given. - explicit StillsRunningGeometry(const DiffractionExperiment &experiment, int nthreads); - - // A frame as analysed. True once: for the frame after which Refine is to be called. - bool Add(const DataMessage &msg); - // The refined geometry, where the fit ran. - std::optional Refine(); - // Frames offered and frames kept, for the log - [[nodiscard]] int64_t Offered() const; - [[nodiscard]] size_t Strong() const; - -private: - const DiffractionExperiment experiment; - const int nthreads; - mutable std::mutex m; - std::vector frames; - int64_t offered = 0, indexed = 0; - size_t strong = 0; - bool done = false; -}; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 217ac9a08..559864c22 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -149,6 +149,7 @@ ADD_EXECUTABLE(jfjoch_test MovingAverageTest.cpp ImageMetadataTest.cpp JFJochReceiverLiteTest.cpp + JFJochReceiverRugnuxTest.cpp GridScanSettingsTest.cpp JFJochReceiverPlotsTest.cpp GoniometerAxisTest.cpp diff --git a/tests/JFJochReceiverProcessingTest.cpp b/tests/JFJochReceiverProcessingTest.cpp index 0b949f197..ac1e72260 100644 --- a/tests/JFJochReceiverProcessingTest.cpp +++ b/tests/JFJochReceiverProcessingTest.cpp @@ -107,149 +107,6 @@ TEST_CASE("JFJochIntegrationTest_ZMQ_crystal_spot_and_index", "[JFJochReceiver]" REQUIRE_NOTHROW(writer_future.get()); } -TEST_CASE("JFJochIntegrationTest_ZMQ_crystal_spot_and_index_GPUStills", "[JFJochReceiver]") { - Logger logger(Catch::getResultCapture().getCurrentTestName()); - - RegisterHDF5Filter(); - - const uint16_t nthreads = 4; - - DiffractionExperiment experiment(DetJF4M()); - experiment.ImagesPerTrigger(5).NumTriggers(1).UseInternalPacketGenerator(true).ImagesPerFile(2) - .FilePrefix("crystal_test_gpu_stills").JungfrauConvPhotonCnt(false).SetFileWriterFormat(FileWriterFormat::NXmxVDS).OverwriteExistingFiles(true) - .DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4) - .SetUnitCell(UnitCell{.a = 36.9, .b = 78.95, .c = 78.95, .alpha =90, .beta = 90, .gamma = 90}); - experiment.SampleTemperature_K(123.0).RingCurrent_mA(115); - - PixelMask pixel_mask(experiment); - - // Load example image - HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5"); - HDF5DataSet dataset(data, "/entry/data/data"); - HDF5DataSpace file_space(dataset); - - REQUIRE(file_space.GetDimensions()[2] == experiment.GetXPixelsNum()); - REQUIRE(file_space.GetDimensions()[1] == experiment.GetYPixelsNum()); - std::vector image_conv (file_space.GetDimensions()[1] * file_space.GetDimensions()[2]); - - std::vector start = {4,0,0}; - std::vector file_size = {1, file_space.GetDimensions()[1], file_space.GetDimensions()[2]}; - dataset.ReadVector(image_conv, start, file_size); - - std::vector image_raw_geom(experiment.GetModulesNum() * RAW_MODULE_SIZE); - ConvertedToRawGeometry(experiment, image_raw_geom.data(), image_conv.data()); - logger.Info("Loaded image"); - - // Setup acquisition device - AcquisitionDeviceGroup aq_devices; - std::unique_ptr test = std::make_unique(0, 64); - for (int m = 0; m < experiment.GetModulesNum(); m++) - test->SetInternalGeneratorFrame((uint16_t *) image_raw_geom.data() + m * RAW_MODULE_SIZE, m); - - aq_devices.Add(std::move(test)); - - ZMQStream2Pusher pusher({"ipc://*"}); - pusher.WriterNotificationSocket("ipc://*"); - - ZMQImagePuller puller(pusher.GetAddress()[0]); - StreamWriter writer(logger, puller); - auto writer_future = std::async(std::launch::async, &StreamWriter::Run, &writer); - - JFJochReceiverService service(aq_devices, logger, pusher); - service.NumThreads(nthreads); - service.Indexing(experiment.GetIndexingSettings()); - // Serial data with rugnux on: the converted frames are analysed on the GPU, after the pre-pass. - service.RugnuxLive(true, nullptr, 4, 60); - - // No progress value at the start of measurement - REQUIRE(!service.GetProgress().has_value()); - - SpotFindingSettings settings = DiffractionExperiment::DefaultDataProcessingSettings(); - settings.signal_to_noise_threshold = 2.5; - settings.photon_count_threshold = 5; - settings.min_pix_per_spot = 1; - settings.max_pix_per_spot = 200; - settings.high_resolution_limit = 2.0; - settings.low_resolution_limit = 50.0; - service.SetSpotFindingSettings(settings); - - service.Start(experiment, pixel_mask, nullptr); - auto receiver_out = service.Stop(); - - CHECK(receiver_out.efficiency == 1.0); - REQUIRE(receiver_out.status.indexing_rate); - CHECK(receiver_out.status.indexing_rate.value() == 1.0); - CHECK(receiver_out.status.images_sent == experiment.GetImageNum()); - CHECK(!receiver_out.status.cancelled); - - // Done once everything is written, with the latency from the last image to then - const auto live = service.GetRugnuxStatus(); - CHECK(live.state == RugnuxLiveState::Done); - REQUIRE(live.last_image_to_written_s); - REQUIRE(live.last_image_s); - CHECK(live.elapsed_s >= live.last_image_s.value() + live.last_image_to_written_s.value() - 1e-6); - - // No progress value at the end of measurement - REQUIRE(!service.GetProgress().has_value()); - - REQUIRE_NOTHROW(writer_future.get()); -} - -TEST_CASE("JFJochIntegrationTest_ZMQ_GPUStills_writer_down", "[JFJochReceiver]") { - // The start message goes out after the pre-pass: with no writer listening, the collection ends with - // the error, and the next one runs. - Logger logger(Catch::getResultCapture().getCurrentTestName()); - - RegisterHDF5Filter(); - - DiffractionExperiment experiment(DetJF4M()); - experiment.ImagesPerTrigger(5).NumTriggers(1).UseInternalPacketGenerator(true).ImagesPerFile(2) - .FilePrefix("crystal_test_gpu_stills_writer_down").JungfrauConvPhotonCnt(false) - .SetFileWriterFormat(FileWriterFormat::NXmxVDS).OverwriteExistingFiles(true) - .DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4) - .SetUnitCell(UnitCell{.a = 36.9, .b = 78.95, .c = 78.95, .alpha =90, .beta = 90, .gamma = 90}); - - PixelMask pixel_mask(experiment); - - HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5"); - HDF5DataSet dataset(data, "/entry/data/data"); - HDF5DataSpace file_space(dataset); - std::vector image_conv (file_space.GetDimensions()[1] * file_space.GetDimensions()[2]); - std::vector start = {4,0,0}; - std::vector file_size = {1, file_space.GetDimensions()[1], file_space.GetDimensions()[2]}; - dataset.ReadVector(image_conv, start, file_size); - std::vector image_raw_geom(experiment.GetModulesNum() * RAW_MODULE_SIZE); - ConvertedToRawGeometry(experiment, image_raw_geom.data(), image_conv.data()); - - AcquisitionDeviceGroup aq_devices; - std::unique_ptr test = std::make_unique(0, 64); - for (int m = 0; m < experiment.GetModulesNum(); m++) - test->SetInternalGeneratorFrame((uint16_t *) image_raw_geom.data() + m * RAW_MODULE_SIZE, m); - aq_devices.Add(std::move(test)); - - ZMQStream2Pusher pusher({"ipc://*"}); - pusher.WriterNotificationSocket("ipc://*"); - - JFJochReceiverService service(aq_devices, logger, pusher); - service.NumThreads(4); - service.Indexing(experiment.GetIndexingSettings()); - service.RugnuxLive(true, nullptr, 4, 60); - - service.Start(experiment, pixel_mask, nullptr); - REQUIRE_THROWS_AS(service.Stop(), JFJochException); - REQUIRE(!service.GetProgress().has_value()); - - ZMQImagePuller puller(pusher.GetAddress()[0]); - StreamWriter writer(logger, puller); - auto writer_future = std::async(std::launch::async, &StreamWriter::Run, &writer); - - service.Start(experiment, pixel_mask, nullptr); - auto receiver_out = service.Stop(); - CHECK(receiver_out.status.images_sent == experiment.GetImageNum()); - CHECK(!receiver_out.status.cancelled); - REQUIRE_NOTHROW(writer_future.get()); -} - TEST_CASE("JFJochIntegrationTest_ZMQ_crystal_spot_and_index_refinement", "[JFJochReceiver]") { Logger logger(Catch::getResultCapture().getCurrentTestName()); diff --git a/tests/JFJochReceiverRugnuxTest.cpp b/tests/JFJochReceiverRugnuxTest.cpp new file mode 100644 index 000000000..0f030e89d --- /dev/null +++ b/tests/JFJochReceiverRugnuxTest.cpp @@ -0,0 +1,547 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include + +#include "../acquisition_device/HLSSimulatedDevice.h" +#include "../compression/JFJochCompressor.h" +#include "../image_analysis/structure_refinement/RugnuxModel.h" +#include "../image_puller/TestImagePuller.h" +#include "../image_pusher/ImagePusher.h" +#include "../receiver/JFJochReceiverRugnux.h" +#include "../receiver/JFJochReceiverService.h" +#include "../writer/HDF5Objects.h" + +namespace { + // What the writer would be given, in order. + class RecordingPusher : public ImagePusher { + mutable std::mutex m; + std::vector events; + std::optional start; + std::optional end; + std::vector images; + public: + void StartDataCollection(StartMessage &message) override { + std::unique_lock ul(m); + events.emplace_back("start"); + start = message; + } + bool EndDataCollection(const EndMessage &message) override { + std::unique_lock ul(m); + events.emplace_back("end"); + end = message; + return true; + } + bool SendImage(const uint8_t *, size_t, int64_t number) override { + std::unique_lock ul(m); + if (events.empty() || events.back() != "image") + events.emplace_back("image"); + images.push_back(number); + return true; + } + bool SendCalibration(const CompressedImage &) override { + std::unique_lock ul(m); + events.emplace_back("calibration"); + return true; + } + std::string PrintSetup() const override { return "recording"; } + ImagePusherType GetType() const override { return ImagePusherType::Test; } + + std::vector Events() const { std::unique_lock ul(m); return events; } + std::vector Images() const { std::unique_lock ul(m); return images; } + std::optional Start() const { std::unique_lock ul(m); return start; } + std::optional End() const { std::unique_lock ul(m); return end; } + }; + + DiffractionExperiment StillsExperiment(int64_t images) { + DiffractionExperiment x(DetDECTRIS(2068, 2164, "Test", {})); + // As JFJochReceiverLite's test: the frames are signed int16, and the stream's start message says so. + x.ImagesPerTrigger(images).NumTriggers(1).ImagesPerFile(2).FilePrefix("rugnux_stills").JungfrauConvPhotonCnt(false) + .PixelSigned(true).DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4) + .SetUnitCell(UnitCell{.a = 36.9, .b = 78.95, .c = 78.95, .alpha = 90, .beta = 90, .gamma = 90}); + return x; + } + + // The frames of the test data, as a stream delivers them: compressed, one DataMessage each, no more. + class Stream { + std::vector> frames; + size_t width, height; + public: + const std::shared_ptr puller = std::make_shared(); + + Stream(const DiffractionExperiment &x, int64_t nframes, bool blank = false) + : width(x.GetXPixelsNum()), height(x.GetYPixelsNum()) { + HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5"); + HDF5DataSet dataset(data, "/entry/data/data"); + HDF5DataSpace file_space(dataset); + std::vector image(file_space.GetDimensions()[1] * file_space.GetDimensions()[2]); + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + for (int64_t i = 0; i < nframes; i++) { + const std::vector start = {static_cast(i % file_space.GetDimensions()[0]), 0, 0}; + const std::vector size = {1, file_space.GetDimensions()[1], file_space.GetDimensions()[2]}; + if (blank) + std::fill(image.begin(), image.end(), 0); + else + dataset.ReadVector(image, start, size); + frames.push_back(compressor.Compress(image)); + } + StartMessage start_msg; + x.FillMessage(start_msg); + puller->Put(ImagePullerOutput{.cbor = std::make_shared(start_msg)}); + } + + void Frames(int64_t first, int64_t count) { + for (int64_t i = first; i < first + count; i++) { + DataMessage message; + message.number = i; + message.image = CompressedImage(frames[i].data(), frames[i].size(), width, height, + CompressedImageMode::Int16, CompressionAlgorithm::BSHUF_LZ4); + puller->Put(ImagePullerOutput{.cbor = std::make_shared(message)}); + } + } + + void Calibration() { + static const std::vector pedestal(4, 1); + auto cbor = std::make_shared(); + cbor->calibration = CompressedImage(pedestal, 2, 2); + puller->Put(ImagePullerOutput{.cbor = cbor}); + } + + void End() { + puller->Put(ImagePullerOutput{.cbor = std::make_shared(EndMessage{})}); + } + }; + + void WaitForStatus(JFJochReceiverService &service, const std::function &what) { + const auto give_up = std::chrono::steady_clock::now() + std::chrono::minutes(5); + while (!what(service.GetRugnuxStatus())) { + REQUIRE(std::chrono::steady_clock::now() < give_up); + std::this_thread::sleep_for(std::chrono::milliseconds(10)); + } + } + + bool HasNote(const RugnuxLiveSnapshot &s, const std::string &text) { + for (const auto &step: s.steps) + for (const auto ¬e: step.notes) + if (note.find(text) != std::string::npos) + return true; + return false; + } +} + +TEST_CASE("JFJochReceiverRugnux_Stills", "[JFJochReceiverRugnux]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + const int64_t images = 12; + DiffractionExperiment x = StillsExperiment(images); + PixelMask mask(x); + + Stream stream(x, images); + // Calibration ahead of the frames waits for the start message, which waits for the pre-pass. + stream.Calibration(); + stream.Frames(0, images); + stream.End(); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + const auto out = service.Stop(); + + CHECK(out.efficiency == 1.0f); + CHECK(!out.status.cancelled); + CHECK(out.status.images_collected == images); + CHECK(out.status.images_sent == images); + REQUIRE(out.status.indexing_rate); + CHECK(out.status.indexing_rate.value() > 0.5f); + CHECK(out.scan_result.images.size() == images); + + // The writer sees the start message before the calibration and the frames, the end message last. + CHECK(pusher.Events() == std::vector{"start", "calibration", "image", "end"}); + auto sent = pusher.Images(); + std::sort(sent.begin(), sent.end()); + REQUIRE(sent.size() == images); + CHECK(sent.front() == 0); + CHECK(sent.back() == images - 1); + + // The merged data travel in the end message, with the geometry the frames were analysed at. + const auto end = pusher.End(); + REQUIRE(end); + CHECK(end->merged.has_value()); + CHECK(end->unit_cell.has_value()); + CHECK(end->refined_beam_center_x.has_value()); + CHECK(end->refined_detector_distance_mm.has_value()); + REQUIRE(pusher.Start()); + CHECK(pusher.Start()->beam_center_x == Catch::Approx(end->refined_beam_center_x.value())); + + const auto live = service.GetRugnuxStatus(); + CHECK(live.state == RugnuxLiveState::Done); + CHECK(live.cut.empty()); + CHECK(HasNote(live, "Merge of 50% of the images (provisional)")); + CHECK(HasNote(live, "Merge of 75% of the images (provisional)")); + CHECK(HasNote(live, "Merge of the collection")); + REQUIRE(live.stills); + CHECK(live.stills->prepass == RugnuxLivePrePass::Done); + CHECK(live.stills->images_analysed == images); + CHECK(live.stills->images_indexed > images / 2); + REQUIRE(live.last_image_to_written_s); + CHECK(!live.space_group.empty()); + + REQUIRE(!service.GetProgress().has_value()); +} + +TEST_CASE("JFJochReceiverRugnux_Stills_CancelDuringPrePass", "[JFJochReceiverRugnux]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + const int64_t images = 8; + DiffractionExperiment x = StillsExperiment(images); + PixelMask mask(x); + Stream stream(x, images); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + int64_t given = 0; + SECTION("while it holds frames") { + given = 3; + stream.Frames(0, given); + service.Start(x, mask, nullptr, stream.puller); + WaitForStatus(service, [&](const RugnuxLiveSnapshot &s) { + return s.stills && s.stills->images_held == given; + }); + } + SECTION("while it runs") { + given = images; + stream.Frames(0, given); + service.Start(x, mask, nullptr, stream.puller); + WaitForStatus(service, [&](const RugnuxLiveSnapshot &s) { + return s.stills && s.stills->prepass == RugnuxLivePrePass::Running; + }); + } + service.Cancel(false); + const auto out = service.Stop(); + + // The frames that came are still analysed and written, and the collection ends properly. + CHECK(out.status.cancelled); + CHECK(out.status.images_sent == given); + CHECK(pusher.Events() == std::vector{"start", "image", "end"}); + const auto live = service.GetRugnuxStatus(); + CHECK(live.state == RugnuxLiveState::Done); + REQUIRE(live.stills); + CHECK(live.stills->images_analysed == given); + REQUIRE(!service.GetProgress().has_value()); +} + +TEST_CASE("JFJochReceiverRugnux_Stills_CancelDuringRun", "[JFJochReceiverRugnux]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + const int64_t images = 40; + DiffractionExperiment x = StillsExperiment(images); + PixelMask mask(x); + Stream stream(x, images); + + // Twelve slots: the pre-pass ends on six frames, and the analysis starts with most of the run to come. + const size_t buffer_MiB = 12 * x.GetImageBufferLocationSize() / (1024 * 1024) + 1; + const int64_t given = 8; + stream.Frames(0, given); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher, buffer_MiB); + service.NumThreads(2); + service.Indexing(x.GetIndexingSettings()); + // A model only so that there is a model step for the cancel to cut. + service.RugnuxLive(true, std::make_shared(), 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + WaitForStatus(service, [](const RugnuxLiveSnapshot &s) { + return s.stills && s.stills->prepass == RugnuxLivePrePass::Done && s.stills->images_analysed > 0; + }); + service.Cancel(false); + const auto out = service.Stop(); + + CHECK(out.status.cancelled); + CHECK(out.status.images_sent == given); + CHECK(pusher.Events() == std::vector{"start", "image", "end"}); + const auto live = service.GetRugnuxStatus(); + CHECK(live.state == RugnuxLiveState::Done); + CHECK(live.cut == "cancel"); + REQUIRE(!service.GetProgress().has_value()); +} + +TEST_CASE("JFJochReceiverRugnux_Stills_ShortRun", "[JFJochReceiverRugnux]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + DiffractionExperiment x = StillsExperiment(8); + PixelMask mask(x); + Stream stream(x, 3); + stream.Frames(0, 3); + stream.End(); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + const auto out = service.Stop(); + + // The pre-pass runs on what came. + CHECK(out.efficiency == Catch::Approx(3.0 / 8.0)); + CHECK(!out.status.cancelled); + CHECK(out.status.images_sent == 3); + CHECK(pusher.Events() == std::vector{"start", "image", "end"}); + const auto live = service.GetRugnuxStatus(); + CHECK(live.state == RugnuxLiveState::Done); + REQUIRE(live.stills); + CHECK(live.stills->images_analysed == 3); +} + +TEST_CASE("JFJochReceiverRugnux_Stills_FramesBeyondTheCollection", "[JFJochReceiverRugnux]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + DiffractionExperiment x = StillsExperiment(3); + PixelMask mask(x); + Stream stream(x, 5); + stream.Frames(0, 5); + stream.End(); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + const auto out = service.Stop(); + + CHECK(out.status.images_collected == 3); + CHECK(out.status.images_sent == 3); +} + +TEST_CASE("JFJochReceiverRugnux_Stills_NoFrames", "[JFJochReceiverRugnux]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + DiffractionExperiment x = StillsExperiment(8); + PixelMask mask(x); + Stream stream(x, 0); + stream.End(); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + const auto out = service.Stop(); + + CHECK(out.status.images_sent == 0); + CHECK(pusher.Events() == std::vector{"start", "end"}); + CHECK(service.GetRugnuxStatus().state == RugnuxLiveState::Done); +} + +TEST_CASE("JFJochReceiverRugnux_Stills_NoIndexableFrames", "[JFJochReceiverRugnux]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + const int64_t images = 6; + DiffractionExperiment x = StillsExperiment(images); + PixelMask mask(x); + Stream stream(x, images, true); + stream.Frames(0, images); + stream.End(); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + const auto out = service.Stop(); + + // Written, with no merge, at the geometry in the header. + CHECK(!out.status.cancelled); + CHECK(out.status.images_sent == images); + CHECK(pusher.Events() == std::vector{"start", "image", "end"}); + REQUIRE(pusher.End()); + CHECK(!pusher.End()->merged.has_value()); + CHECK(pusher.End()->refined_beam_center_x == Catch::Approx(x.GetBeamX_pxl())); + const auto live = service.GetRugnuxStatus(); + CHECK(live.state == RugnuxLiveState::Done); + CHECK(HasNote(live, "could not refine the geometry")); + REQUIRE(live.stills); + CHECK(live.stills->images_analysed == images); + CHECK(live.stills->images_indexed == 0); +} + +TEST_CASE("JFJochReceiverRugnux_Stills_WriterDown", "[JFJochReceiverRugnux]") { + // The start message goes out after the pre-pass: where it fails, the collection ends with the error. + class NoWriter : public RecordingPusher { + public: + bool down = true; + void StartDataCollection(StartMessage &message) override { + if (down) + throw JFJochException(JFJochExceptionCategory::ZeroMQ, "no writer"); + RecordingPusher::StartDataCollection(message); + } + }; + + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + const int64_t images = 4; + DiffractionExperiment x = StillsExperiment(images); + PixelMask mask(x); + Stream stream(x, images); + stream.Frames(0, images); + stream.End(); + + NoWriter pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + REQUIRE_THROWS_AS(service.Stop(), JFJochException); + CHECK(pusher.Events().empty()); + CHECK(service.GetRugnuxStatus().state == RugnuxLiveState::Failed); + + // And the next one runs. + pusher.down = false; + Stream next(x, images); + next.Frames(0, images); + next.End(); + service.Start(x, mask, nullptr, next.puller); + CHECK(service.Stop().status.images_sent == images); + CHECK(pusher.Events() == std::vector{"start", "image", "end"}); +} + +TEST_CASE("JFJochReceiverRugnux_Rotation", "[JFJochReceiverRugnux]") { + // A scan stays with rugnux's own run: the start message goes out at once, and there are no stills results. + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + const int64_t images = 5; + DiffractionExperiment x = StillsExperiment(images); + x.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, {1, 0, 0}, {})); + PixelMask mask(x); + Stream stream(x, images); + stream.Frames(0, images); + stream.End(); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.Indexing(x.GetIndexingSettings()); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, stream.puller); + const auto out = service.Stop(); + + CHECK(out.status.images_sent == images); + CHECK(pusher.Events() == std::vector{"start", "image", "end"}); + const auto live = service.GetRugnuxStatus(); + CHECK(!live.stills); + CHECK(live.images_received == images); + CHECK(live.state != RugnuxLiveState::Idle); + CHECK(live.state != RugnuxLiveState::Running); +} + +TEST_CASE("JFJochReceiverFPGA_RugnuxOnStaysPerImage", "[JFJochReceiverRugnux]") { + // Live stills are this receiver's alone until a puller feeds it: the card's run is as with rugnux off. + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + DiffractionExperiment x(DetJF(1)); + x.ImagesPerTrigger(2).NumTriggers(1).UseInternalPacketGenerator(true).FilePrefix("rugnux_on_fpga") + .DetectorDistance_mm(75).BeamX_pxl(500).BeamY_pxl(250).IncidentEnergy_keV(12.4); + PixelMask mask(x); + + AcquisitionDeviceGroup aq_devices; + aq_devices.Add(std::make_unique(0, 64)); + + RecordingPusher pusher; + JFJochReceiverService service(aq_devices, logger, pusher); + service.NumThreads(2); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr); + const auto out = service.Stop(); + + CHECK(out.status.images_sent == 2); + CHECK(pusher.Events() == std::vector{"start", "image", "end"}); + CHECK(service.GetRugnuxStatus().state == RugnuxLiveState::Idle); +} + +TEST_CASE("JFJochReceiverRugnux_DarkMaskStaysWithLite", "[JFJochReceiverRugnux]") { + // The mask sequence of a later POST /pedestal, with rugnux on: the state machine's experiment still has the + // last /start's goniometer, or none. + Logger logger(Catch::getResultCapture().getCurrentTestName()); + RegisterHDF5Filter(); + + DiffractionExperiment x(DetDECTRIS(2068, 2164, "Test", {})); + SECTION("after a serial collection") {} + SECTION("after a rotation collection") { + x.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, {1, 0, 0}, {})); + } + DarkMaskSettings mask_settings; + mask_settings.NumberOfFrames(10).MaxCounts(1).MaxFramesWithCounts(5); + x.Mode(DetectorMode::DarkMask).ImportDarkMaskSettings(mask_settings); + PixelMask mask(x); + + auto puller = std::make_shared(); + StartMessage start_msg; + x.FillMessage(start_msg); + puller->Put(ImagePullerOutput{.cbor = std::make_shared(start_msg)}); + + std::vector image(x.GetPixelsNum(), 0); + image[2] = 2; + for (int i = 0; i < x.GetFrameNum(); i++) { + DataMessage message; + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.number = i; + puller->Put(ImagePullerOutput{.cbor = std::make_shared(message)}); + } + puller->Put(ImagePullerOutput{.cbor = std::make_shared(EndMessage{})}); + + RecordingPusher pusher; + AcquisitionDeviceGroup group; + JFJochReceiverService service(group, logger, pusher); + service.NumThreads(4); + service.RugnuxLive(true, nullptr, 4, 60); + + service.Start(x, mask, nullptr, puller); + const auto out = service.Stop(); + + REQUIRE(out.dark_mask_result.size() == x.GetPixelsNum()); + CHECK(out.dark_mask_result[2] != 0); + CHECK(out.dark_mask_result[0] == 0); + CHECK(out.status.images_sent == 0); + CHECK(!service.GetRugnuxStatus().stills); +} diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index e8986a690..2fee6e828 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -34,6 +34,7 @@ #include "../writer/WriteModel.h" #include "../writer/RugnuxModelFile.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../rugnux/StillsModel.h" namespace { // A synthetic P1 cell with two carbon atoms - enough for a reader to produce a Structure with @@ -2068,3 +2069,89 @@ TEST_CASE("ModelStructureFactorsGPU_MapMatchesCPU", "[ModelValidation][gpu]") { #endif + +// The model step of a live serial merge (StillsModel.h). A merge in the other indexing of a P4 crystal is put in the +// model's by the probe; a later merge of the same collection, again in the other indexing, is put in the earlier +// step's frame and placed from where the model was, with no null or probe made again; a cut leaves out what has not +// started. The group a merge without one is merged in is the model's where its cell is the run's. +TEST_CASE("StillsModel_PlacesWarmStartsAndCuts", "[ModelValidation][RugnuxLive]") { + Logger logger("StillsModel_PlacesWarmStartsAndCuts"); + const auto path = WriteTemp("stills_model_test.pdb", + ClusterPdb("CRYST1 34.000 34.000 38.000 90.00 90.00 90.00 P 4 4\n").c_str()); + const RugnuxModel model = ReadModel(path); + auto own = ModelReferenceIntensities(model, {}, {}, 2.5, logger); + REQUIRE(own.size() > 500); + for (size_t i = 0; i < own.size(); i++) { + own[i].F = std::sqrt(std::max(0.0f, own[i].I)); + own[i].sigma = 0.02f * own[i].I + 1.0f; + own[i].sigmaF = 0.02f * own[i].F + 0.1f; + own[i].rfree_flag = (i % 20) == 0; + } + const UnitCell cell{.a = 34, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}; + const gemmi::SpaceGroup p4 = *gemmi::find_spacegroup_by_name("P 4"); + const auto laws = ReindexAmbiguityOperators(cell, p4); + REQUIRE(laws.size() == 1); + const auto other = ReindexMergedIntoAsu(own, laws.front(), p4, true); + + // The group: the collection's, else the model's on the run's cell, else none + CHECK(StillsModelSpaceGroup(*gemmi::find_spacegroup_by_name("P 4 2 2"), cell, model).space_group->xhm() == "P 4 2 2"); + CHECK(StillsModelSpaceGroup(std::nullopt, cell, model).space_group->xhm() == "P 4"); + const UnitCell longer{.a = 34, .b = 34, .c = 45, .alpha = 90, .beta = 90, .gamma = 90}; + const auto none = StillsModelSpaceGroup(std::nullopt, longer, model); + CHECK_FALSE(none.space_group); + CHECK_FALSE(none.reason.empty()); + + StillsModelSettings settings; + settings.refine = false; + const StillsModelResult first = PlaceAndRefineStillsModel(other, cell, p4, model, 1.0, true, {}, settings, nullptr, + logger); + REQUIRE(first.validation); + CHECK(first.placement == "RIGID_BODY"); + CHECK_FALSE(first.warm); + CHECK(first.validation->indexing_op == laws.front()); + REQUIRE(first.maps); + REQUIRE(first.placed_model); + // In the model's indexing: the model's own intensities + CHECK(ReferenceIntensityCC(first.merged, own, p4) > 0.99); + CHECK(first.validation->r_work < 0.15); + + // A later merge, in the other indexing again, and refined + settings.refine = true; + const StillsModelResult second = PlaceAndRefineStillsModel(other, cell, p4, model, 1.0, true, {}, settings, &first, + logger); + REQUIRE(second.validation); + CHECK(second.warm); + CHECK(second.placement == "RIGID_BODY"); + CHECK(second.validation->model_fits); + CHECK(ReferenceIntensityCC(second.merged, own, p4) > 0.99); + REQUIRE(second.refinement); + CHECK(second.refinement->ok); + CHECK(second.refinement->r_free < 0.2); + + // Cut before anything: the refinement is left out; and for data the model does not describe, the molecular + // replacement too + std::atomic cut{true}; + settings.cut = &cut; + const StillsModelResult cut_placed = PlaceAndRefineStillsModel(other, cell, p4, model, 1.0, true, {}, settings, + &first, logger); + CHECK_FALSE(cut_placed.refinement); + CHECK(cut_placed.left_out == std::vector{"refinement"}); + const auto unrelated_path = WriteTemp("stills_model_unrelated_test.pdb", + ClusterPdb("CRYST1 34.000 34.000 38.000 90.00 90.00 90.00 P 4 4\n", 7).c_str()); + auto unrelated = ModelReferenceIntensities(ReadModel(unrelated_path), {}, {}, 2.5, logger); + for (size_t i = 0; i < unrelated.size(); i++) { + unrelated[i].F = std::sqrt(std::max(0.0f, unrelated[i].I)); + unrelated[i].sigma = 0.02f * unrelated[i].I + 1.0f; + unrelated[i].sigmaF = 0.02f * unrelated[i].F + 0.1f; + unrelated[i].rfree_flag = (i % 20) == 0; + } + const StillsModelResult cut_unplaced = PlaceAndRefineStillsModel(unrelated, cell, p4, model, 1.0, true, {}, + settings, nullptr, logger); + CHECK(cut_unplaced.placement == "NOT_PLACED"); + CHECK_FALSE(cut_unplaced.mr); + CHECK(cut_unplaced.left_out == std::vector{"molecular replacement"}); + CHECK_FALSE(cut_unplaced.placed_model); + + std::filesystem::remove(path); + std::filesystem::remove(unrelated_path); +} diff --git a/tests/StillsGPUAnalysisTest.cpp b/tests/StillsGPUAnalysisTest.cpp index ef803a2e7..d4873afff 100644 --- a/tests/StillsGPUAnalysisTest.cpp +++ b/tests/StillsGPUAnalysisTest.cpp @@ -4,8 +4,10 @@ #include #include +#include #include "../common/ImageBuffer.h" +#include "../image_analysis/structure_refinement/RugnuxModel.h" #include "../receiver/StillsGPUAnalysis.h" namespace { @@ -51,7 +53,7 @@ TEST_CASE("StillsGPUAnalysis_StartFailureReleasesThreads", "[StillsGPUAnalysis]" [] { return DiffractionExperiment::DefaultDataProcessingSettings(); }, [&](const JFJochException &) { ++failed; }, &status); for (int i = 0; i < 3; i++) - REQUIRE(gpu.Hold(Frame(buffer, i))); + gpu.Hold(Frame(buffer, i)); REQUIRE_THROWS_AS(gpu.Finish(), JFJochException); } CHECK(failed == 1); @@ -75,14 +77,11 @@ TEST_CASE("StillsGPUAnalysis_HoldIsBounded", "[StillsGPUAnalysis]") { }, [](DataMessage &, const AzimuthalIntegrationProfile &) {}, [] { return DiffractionExperiment::DefaultDataProcessingSettings(); }); - // The analysis waits for the pre-pass, which waits for its sample: nothing is taken off the queue. - for (int i = 0; i < 3; i++) - REQUIRE(gpu.Hold(Frame(buffer, i))); - { - const auto refused = Frame(buffer, 3); - CHECK(!gpu.Hold(refused)); - } - CHECK(buffer.GetStatus().preparation_slots == 3); + // Past what the pre-pass holds, frames are still taken: only the image buffer bounds them. + for (int i = 0; i < 4; i++) + gpu.Hold(Frame(buffer, i)); + CHECK(buffer.GetStatus().preparation_slots == 4); + CHECK(buffer.GetImageSlot(buffer.GetStatus().available_slots) == nullptr); REQUIRE_THROWS(gpu.Finish()); CHECK(buffer.GetStatus().preparation_slots == 0); } @@ -107,7 +106,7 @@ TEST_CASE("StillsGPUAnalysis_FailingFrameDoesNotStopAnalysis", "[StillsGPUAnalys {}, &status); // The frames carry no image, so their analysis fails too - they are passed on all the same. for (int i = 0; i < 3; i++) - REQUIRE(gpu.Hold(Frame(buffer, i))); + gpu.Hold(Frame(buffer, i)); REQUIRE_NOTHROW(gpu.Finish()); CHECK(sunk == 3); // The outcome waits for everything to be written, and comes after the latency. @@ -125,3 +124,62 @@ TEST_CASE("StillsGPUAnalysis_FailingFrameDoesNotStopAnalysis", "[StillsGPUAnalys CHECK(s.steps.back().duration_s); CHECK(buffer.GetStatus().preparation_slots == 0); } + +TEST_CASE("StillsGPUAnalysis_CancelCutsTheModelStep", "[StillsGPUAnalysis]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + const DiffractionExperiment x = StillsExperiment(); + const PixelMask mask(x); + const auto run = [&](RugnuxLiveStatus &status, std::shared_ptr model) { + return std::make_unique( + x, mask, Start(x), nullptr, 1, 100, logger, [](const StillsGPUAnalysis::Setup &) {}, + [](DataMessage &, const AzimuthalIntegrationProfile &) {}, + [] { return DiffractionExperiment::DefaultDataProcessingSettings(); }, + StillsGPUAnalysis::Failed{}, &status, true, std::move(model), 60.0); + }; + + RugnuxLiveStatus during; + auto a = run(during, std::make_shared()); + a->Cut("cancel"); + CHECK(during.Get().cut == "cancel"); + REQUIRE_NOTHROW(a->Finish()); + + // Nothing is left to cut once the model step has ended, and nothing runs without a model. + RugnuxLiveStatus after; + auto b = run(after, std::make_shared()); + b->Finish(); + b->Cut("cancel"); + CHECK(after.Get().cut.empty()); + + RugnuxLiveStatus no_model; + auto c = run(no_model, nullptr); + c->Cut("cancel"); + CHECK(no_model.Get().cut.empty()); + c->Finish(); +} + +TEST_CASE("StillsGPUAnalysis_FramesDuringThePrePassAreAnalysed", "[StillsGPUAnalysis]") { + Logger logger(Catch::getResultCapture().getCurrentTestName()); + const DiffractionExperiment x = StillsExperiment(); + const PixelMask mask(x); + ImageBuffer buffer(64 * 1024 * 1024); + buffer.StartMeasurement(x); + const int64_t frames = 3 + 5; + REQUIRE(buffer.GetStatus().total_slots >= frames); + + std::atomic sunk = 0; + std::latch offered(1); + RugnuxLiveStatus status; + StillsGPUAnalysis gpu(x, mask, Start(x), nullptr, 2, 3, logger, + // The pre-pass is not over until every frame has been offered. + [&](const StillsGPUAnalysis::Setup &) { offered.wait(); }, + [&](DataMessage &, const AzimuthalIntegrationProfile &) { ++sunk; }, + [] { return DiffractionExperiment::DefaultDataProcessingSettings(); }, + {}, &status); + CHECK(status.Get().stills->images_needed == 3); + for (int i = 0; i < frames; i++) + gpu.Hold(Frame(buffer, i)); + offered.count_down(); + REQUIRE_NOTHROW(gpu.Finish()); + CHECK(sunk == frames); + CHECK(buffer.GetStatus().preparation_slots == 0); +} diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 18abcf779..a2cc71c45 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -48,7 +48,7 @@ IF(NOT JFJOCH_PORTABLE_ONLY) # Not installed: replays a stored serial dataset through the GPU analysis of converted JUNGFRAU frames. ADD_EXECUTABLE(jfjoch_stills_gpu_replay jfjoch_stills_gpu_replay.cpp) - TARGET_LINK_LIBRARIES(jfjoch_stills_gpu_replay JFJochReceiver JFJochWriter JFJochReader) + TARGET_LINK_LIBRARIES(jfjoch_stills_gpu_replay JFJochReceiver JFJochWriter JFJochReader JFJochModelWriter JFJochRugnuxFiles) ADD_EXECUTABLE(jfjoch_pcie_status jfjoch_pcie_status.cpp) TARGET_LINK_LIBRARIES(jfjoch_pcie_status JFJochDevice JFJochCommon) diff --git a/tools/jfjoch_stills_gpu_replay.cpp b/tools/jfjoch_stills_gpu_replay.cpp index 722da95d9..208710c23 100644 --- a/tools/jfjoch_stills_gpu_replay.cpp +++ b/tools/jfjoch_stills_gpu_replay.cpp @@ -1,10 +1,13 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -// Replays the frames of a stored serial collection through the broker's GPU analysis of converted -// JUNGFRAU frames (StillsGPUAnalysis): the pre-pass on the first images, then every frame analysed with the +// Replays the frames of a stored serial collection through the broker's live analysis of serial +// frames (StillsGPUAnalysis): the pre-pass on the first images, then every frame analysed with the // beam centre and mask it decided. Reports the pre-pass, the indexed fraction and the sustained rate, and -// writes the mask the frames were analysed with (_mask.u32, uint32, row by row). +// writes the mask the frames were analysed with (_mask.u32, uint32, row by row). With a model (-m), the +// model step of the merges as the broker makes it, and the files the end message would carry: .mtz (the +// merge, in the model's frame), _maps.mtz, _model.cif, _refined_model.cif, +// _refined_maps.mtz. #include #include @@ -14,6 +17,8 @@ #include #include #include +#include +#include #include "../common/ImageBuffer.h" #include "../common/Logger.h" @@ -21,6 +26,10 @@ #include "../reader/JFJochHDF5Reader.h" #include "../receiver/StillsGPUAnalysis.h" #include "../writer/HDF5Objects.h" +#include "../writer/RugnuxModelFile.h" +#include "../writer/WriteModel.h" +#include "../writer/WriteRugnuxFiles.h" +#include "../image_analysis/WriteReflections.h" void print_usage() { std::cout << "Usage: jfjoch_stills_gpu_replay [options] " << std::endl; @@ -37,6 +46,8 @@ void print_usage() { std::cout << " -S space group of the collection (the merge's)" << std::endl; std::cout << " -M no merge (StillsGPUAnalysis merge = false)" << std::endl; std::cout << " -e the first images given as empty frames (no diffraction)" << std::endl; + std::cout << " -m model (PDB/mmCIF) placed and refined against the merges" << std::endl; + std::cout << " -D deadline of the model step after the last image (default 60, 0 = none)" << std::endl; } int main(int argc, char **argv) { @@ -56,9 +67,11 @@ int main(int argc, char **argv) { std::optional space_group; bool merge = true; int64_t empty_frames = 0; + std::shared_ptr model; + double deadline_s = 60; int opt; - while ((opt = getopt(argc, argv, "o:s:n:r:N:b:x:y:C:cFS:Me:")) != -1) { + while ((opt = getopt(argc, argv, "o:s:n:r:N:b:x:y:C:cFS:Me:m:D:")) != -1) { switch (opt) { case 'o': prefix = optarg; break; case 's': first = atoll(optarg); break; @@ -73,6 +86,8 @@ int main(int argc, char **argv) { case 'S': space_group = optarg; break; case 'M': merge = false; break; case 'e': empty_frames = atoll(optarg); break; + case 'm': model = std::make_shared(ReadRugnuxModelFile(optarg)); break; + case 'D': deadline_s = atof(optarg); break; case 'C': { UnitCell c{}; if (sscanf(optarg, "%f,%f,%f,%f,%f,%f", &c.a, &c.b, &c.c, &c.alpha, &c.beta, &c.gamma) != 6) { @@ -136,15 +151,20 @@ int main(int argc, char **argv) { ImageBuffer buffer(buffer_MiB * 1024 * 1024); buffer.StartMeasurement(x); - std::atomic analysed = 0, indexed = 0, spots = 0; + std::atomic analysed = 0, indexed = 0, spots = 0, given = 0; + std::atomic caught_up = false; + double buffer_full_s = 0.0; std::atomic processing_s = 0.0; std::chrono::steady_clock::time_point t_ready; const auto t0 = std::chrono::steady_clock::now(); { // As the receiver: half the buffer for the frames waiting, the other half for those sent. + std::mutex run_profile_mutex; + std::unique_ptr run_profile; StillsGPUAnalysis gpu(x, mask, start, &pool, nthreads, buffer.GetStatus().total_slots / 2, logger, [&](const StillsGPUAnalysis::Setup &s) { t_ready = std::chrono::steady_clock::now(); + run_profile = std::make_unique(*s.mapping); const auto &m = s.mask.GetMask(); std::ofstream f(prefix + "_mask.u32", std::ios::binary); f.write(reinterpret_cast(m.data()), m.size() * sizeof(uint32_t)); @@ -152,14 +172,23 @@ int main(int argc, char **argv) { std::chrono::duration(t_ready - t0).count(), s.experiment.GetBeamX_pxl(), s.experiment.GetBeamY_pxl()); }, - [&](DataMessage &m, const AzimuthalIntegrationProfile &) { - ++analysed; + [&](DataMessage &m, const AzimuthalIntegrationProfile &profile) { + { + // As the receiver's plots do: the frames' profiles add up over the run. + std::unique_lock ul(run_profile_mutex); + *run_profile += profile; + } + // The analysis has caught up with the frames given. + if (++analysed + nthreads >= given && !caught_up.exchange(true)) + logger.Info("Analysis caught up {:.1f} s after it started, {} frames analysed", + std::chrono::duration(std::chrono::steady_clock::now() + - t_ready).count(), analysed.load()); spots += static_cast(m.spots.size()); processing_s += m.processing_time_s.value_or(0.0f); if (m.indexing_result.value_or(false)) ++indexed; }, - [&] { return spot_finding; }, {}, nullptr, merge); + [&] { return spot_finding; }, {}, nullptr, merge, model, deadline_s); JFJochReaderRawImage raw; std::vector empty; @@ -175,19 +204,24 @@ int main(int argc, char **argv) { } auto held = std::make_shared(); // Unlike the receiver, the replay waits for a slot rather than lose the frame. + const auto t_slot = std::chrono::steady_clock::now(); while ((held->slot = buffer.GetImageSlot()) == nullptr) std::this_thread::sleep_for(std::chrono::milliseconds(1)); + buffer_full_s += std::chrono::duration(std::chrono::steady_clock::now() - t_slot).count(); memcpy(held->slot->GetImage(), raw.image.GetCompressed(), raw.image.GetCompressedSize()); held->message.number = i; held->message.image = CompressedImage(held->slot->GetImage(), raw.image.GetCompressedSize(), raw.image.GetWidth(), raw.image.GetHeight(), raw.image.GetMode(), raw.image.GetCompressionAlgorithm()); held->message.image_collection_efficiency = 1.0f; - while (!gpu.Hold(held)) - std::this_thread::sleep_for(std::chrono::milliseconds(1)); + gpu.Hold(held); + ++given; } - logger.Info("Gave {} images in {:.1f} s", n, - std::chrono::duration(std::chrono::steady_clock::now() - t0).count()); + const auto t_given = std::chrono::steady_clock::now(); + // The replay waits for a slot where a receiver would drop the frame. + logger.Info("Buffer full for {:.1f} s: at -r {} Hz about {:.0f} frames would be dropped", buffer_full_s, rate, + buffer_full_s * rate); + logger.Info("Gave {} images in {:.1f} s", n, std::chrono::duration(t_given - t0).count()); gpu.Finish(); if (const auto c = gpu.GetConsensusUnitCell()) logger.Info("Consensus cell {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}", c->a, c->b, c->c, c->alpha, c->beta, @@ -201,6 +235,22 @@ int main(int argc, char **argv) { os << m->merge.statistics; logger.Info("Merging statistics:\n{}", os.str()); } + logger.Info("Last image given -> merge and model step done: {:.2f} s", std::chrono::duration( + std::chrono::steady_clock::now() - t_given).count()); + if (const auto &mx = gpu.GetModel()) { + DiffractionExperiment mx_x(gpu.GetMergeExperiment()); + mx_x.SetSpaceGroup(mx->space_group); + WriteMtz(MergedMtzTable(mx->merged, mx->cell, mx_x, OutputProvenance{}), prefix + ".mtz"); + if (mx->maps) + WriteMtz(mx->maps->coefficients, prefix + "_maps.mtz"); + if (mx->placed_model) + WritePlacedModel(*mx->placed_model->model, mx->cell, mx->space_group, prefix, logger); + if (const auto &rf = mx->refinement; rf && rf->ok && rf->model) { + WritePlacedModel(*rf->model, mx->cell, mx->space_group, prefix + "_refined", logger); + if (rf->maps) + WriteMtz(rf->maps->coefficients, prefix + "_refined_maps.mtz"); + } + } } const auto t1 = std::chrono::steady_clock::now(); const double after_ready = std::chrono::duration(t1 - t_ready).count();