diff --git a/broker/JFJochBrokerHttp.cpp b/broker/JFJochBrokerHttp.cpp index 88cd13b6c..f0b638ba3 100644 --- a/broker/JFJochBrokerHttp.cpp +++ b/broker/JFJochBrokerHttp.cpp @@ -7,7 +7,6 @@ #include #include -#include #include #include #include @@ -17,6 +16,7 @@ #include "../preview/JFJochTIFF.h" #include "OpenAPIConvert.h" #include "../image_pusher/LiveCollection.h" +#include "../writer/RugnuxModelFile.h" #include #include "gen/model/Error_message.h" @@ -978,9 +978,9 @@ void JFJochBrokerHttp::config_rugnux_get(httplib::Response &response) { ret.setEnable(s.enable); ret.setReservedCpus(s.reserved_cpus); ret.setDeadlineS(s.deadline_s); - ret.setModelLoaded(!s.model_path.empty()); - if (!s.model_path.empty()) - ret.setModelSizeBytes(static_cast(std::filesystem::file_size(s.model_path))); + ret.setModelLoaded(s.model != nullptr); + if (s.model) + ret.setModelSizeBytes(s.model_size_bytes); ProcessOutput(ret, response); } @@ -1003,13 +1003,21 @@ void JFJochBrokerHttp::config_rugnux_model_put(const httplib::Request &request, response.set_content("Model cannot be larger than 64 MiB", "text/plain"); return; } - state_machine.SetRugnuxModel(request.body); + std::shared_ptr model; + try { + model = std::make_shared(ParseRugnuxModel(request.body, "model")); + } catch (const std::exception &e) { + response.status = 400; + response.set_content(std::string("Not an atomic model: ") + e.what(), "text/plain"); + return; + } + state_machine.SetRugnuxModel(std::move(model), static_cast(request.body.size())); logger.Info("Rugnux model uploaded ({} bytes)", request.body.size()); response.status = 200; } void JFJochBrokerHttp::config_rugnux_model_delete(httplib::Response &response) { - state_machine.SetRugnuxModel(""); + state_machine.SetRugnuxModel(nullptr, 0); logger.Info("Rugnux model removed"); response.status = 200; } diff --git a/broker/JFJochServices.cpp b/broker/JFJochServices.cpp index 8bd4ff15c..31dfaee5f 100644 --- a/broker/JFJochServices.cpp +++ b/broker/JFJochServices.cpp @@ -3,9 +3,6 @@ #include "JFJochServices.h" -#include -#include -#include #include "../common/JFJochException.h" #include "../detector_control/SLSDetectorWrapper.h" #include "../detector_control/DectrisDetectorWrapper.h" @@ -376,7 +373,7 @@ void JFJochServices::SetupIndexing(const IndexingSettings &input) { void JFJochServices::SetRugnuxLive(const RugnuxLiveSettings &input) { if (receiver) - receiver->RugnuxLive(input.enable, input.model_path, input.reserved_cpus, input.deadline_s); + receiver->RugnuxLive(input.enable, input.model, input.reserved_cpus, input.deadline_s); rugnux = input; } @@ -390,32 +387,10 @@ RugnuxLiveSnapshot JFJochServices::GetRugnuxStatus() const { return {}; } -void JFJochServices::SetRugnuxModel(const std::string &model) { +void JFJochServices::SetRugnuxModel(std::shared_ptr model, int64_t size_bytes) { RugnuxLiveSettings s = rugnux; - if (!s.model_path.empty()) - std::filesystem::remove(s.model_path); - s.model_path.clear(); - - if (!model.empty()) { - if (rugnux_model_dir.empty()) { - char dir[] = "/tmp/jfjoch_rugnux_XXXXXX"; // mkdtemp creates it 0700 - if (mkdtemp(dir) == nullptr) - throw JFJochException(JFJochExceptionCategory::FileWriteError, - "Cannot create directory for the rugnux model"); - rugnux_model_dir = dir; - } - // gemmi detects PDB/mmCIF from the content, but gzip only from the extension - const bool gzipped = model.size() >= 2 && model[0] == '\x1f' && model[1] == '\x8b'; - s.model_path = rugnux_model_dir + (gzipped ? "/model.gz" : "/model"); - { - std::ofstream f(s.model_path, std::ios::binary); - f.write(model.data(), static_cast(model.size())); - if (!f) - throw JFJochException(JFJochExceptionCategory::FileWriteError, "Cannot write the rugnux model"); - } - std::filesystem::permissions(s.model_path, - std::filesystem::perms::owner_read | std::filesystem::perms::owner_write); - } + s.model = std::move(model); + s.model_size_bytes = s.model ? size_bytes : 0; SetRugnuxLive(s); } diff --git a/broker/JFJochServices.h b/broker/JFJochServices.h index 3a732c20d..c232ff5b9 100644 --- a/broker/JFJochServices.h +++ b/broker/JFJochServices.h @@ -18,7 +18,8 @@ struct RugnuxLiveSettings { bool enable = false; int64_t reserved_cpus = 4; double deadline_s = 0; // the analysis is cut this long after the last image (0 = never) - std::string model_path; // empty = no model + std::shared_ptr model; // null = no model + int64_t model_size_bytes = 0; // as uploaded }; class JFJochServices { @@ -31,7 +32,6 @@ class JFJochServices { std::atomic cannot_stop_detector = false; RugnuxLiveSettings rugnux; - std::string rugnux_model_dir; // private directory, created at the first model upload Logger &logger; public: explicit JFJochServices(Logger &in_logger); @@ -89,7 +89,7 @@ public: void SetRugnuxLive(const RugnuxLiveSettings &input); RugnuxLiveSettings GetRugnuxLive() const; // The model file is written 0600 into a private directory; an empty model removes it - void SetRugnuxModel(const std::string &model); + void SetRugnuxModel(std::shared_ptr model, int64_t size_bytes); RugnuxLiveSnapshot GetRugnuxStatus() const; ImagePusherStatus GetImagePusherStatus() const; diff --git a/broker/JFJochStateMachine.cpp b/broker/JFJochStateMachine.cpp index 47dacea7f..97c880ded 100644 --- a/broker/JFJochStateMachine.cpp +++ b/broker/JFJochStateMachine.cpp @@ -1221,12 +1221,12 @@ void JFJochStateMachine::SetRugnuxLive(bool enable, int64_t reserved_cpus, doubl services.SetRugnuxLive(s); } -void JFJochStateMachine::SetRugnuxModel(const std::string &model) { +void JFJochStateMachine::SetRugnuxModel(std::shared_ptr model, int64_t size_bytes) { std::unique_lock ul(m); if (IsRunning()) throw WrongDAQStateException("Cannot change rugnux model during data collection"); std::unique_lock ul2(rugnux_mutex); - services.SetRugnuxModel(model); + services.SetRugnuxModel(std::move(model), size_bytes); } RugnuxLiveSettings JFJochStateMachine::GetRugnuxLive() const { diff --git a/broker/JFJochStateMachine.h b/broker/JFJochStateMachine.h index ca053658a..ed494dcaf 100644 --- a/broker/JFJochStateMachine.h +++ b/broker/JFJochStateMachine.h @@ -249,7 +249,7 @@ public: IndexingSettings GetIndexingSettings() const; void SetRugnuxLive(bool enable, int64_t reserved_cpus, double deadline_s); - void SetRugnuxModel(const std::string &model); // empty = remove the model + void SetRugnuxModel(std::shared_ptr model, int64_t size_bytes); // null = remove the model RugnuxLiveSettings GetRugnuxLive() const; RugnuxLiveSnapshot GetRugnuxStatus() const; diff --git a/image_analysis/structure_refinement/CMakeLists.txt b/image_analysis/structure_refinement/CMakeLists.txt index aacd1211b..0644f809d 100644 --- a/image_analysis/structure_refinement/CMakeLists.txt +++ b/image_analysis/structure_refinement/CMakeLists.txt @@ -5,6 +5,7 @@ ADD_LIBRARY(JFJochStructureRefinement STATIC MolecularReplacement.cpp MolecularReplacement.h $<$:MolecularReplacementGPU.cu> MolecularReplacementGPU.h ModelValidation.cpp ModelValidation.h + RugnuxModel.h RigidBodyRefine.cpp RigidBodyRefine.h $<$:RigidBodyGPU.cpp> $<$:RigidBodyGPU.cu> diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index a79dd085d..0497730c2 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -21,9 +21,7 @@ #include #include -#include // read_structure_gz #include // expand_ncs -#include // MaybeGzipped #include // IT92 x-ray form factors #include // DensityCalculator #include // get_size_for_hkl, get_f_phi_on_grid (the maps) @@ -537,7 +535,7 @@ namespace { // validation either way. ModelValidationResult Validate(const std::vector &merged, const UnitCell &cell, - const std::string &model_path, + const RugnuxModel &model, const std::string &output_prefix, Logger &logger, const gemmi::SpaceGroup *data_space_group, @@ -548,32 +546,29 @@ ModelValidationResult Validate(const std::vector &merged, const ModelValidationSchedule &schedule, const gemmi::Op *twin_law, bool rigid_body_gpu, - const gemmi::Structure *model) { + const gemmi::Structure *placed) { ModelValidationResult result; - result.model_path = model_path; + result.model_name = model.name; - // --- read the atomic model --- + // --- the atomic model --- ModelState mdl; gemmi::Structure &st = mdl.st; // the real model, the one the maps and the report describe try { - // Detect, not the default: without it GEMMI picks the format from the extension and only - // falls back to the content when it does not recognise one. A model arrives named however - // whoever produced it named it, so the file itself is the better authority. - st = model != nullptr ? *model : gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect); + st = placed != nullptr ? *placed : model.structure; expand_strict_ncs(st); } catch (const std::exception &e) { - result.failure_reason = fmt::format("cannot read model {}: {}", model_path, e.what()); + result.failure_reason = fmt::format("cannot read model {}: {}", model.name, e.what()); logger.Error("Model validation: {}", result.failure_reason); return result; } if (st.models.empty() || !st.cell.is_crystal()) { - result.failure_reason = fmt::format("model {} has no atoms or no unit cell", model_path); + result.failure_reason = fmt::format("model {} has no atoms or no unit cell", model.name); logger.Error("Model validation: {}", result.failure_reason); return result; } const gemmi::SpaceGroup *sg = st.find_spacegroup(); if (!sg) { - result.failure_reason = fmt::format("model {} has no usable space group", model_path); + result.failure_reason = fmt::format("model {} has no usable space group", model.name); logger.Error("Model validation: {}", result.failure_reason); return result; } @@ -659,7 +654,7 @@ ModelValidationResult Validate(const std::vector &merged, // and their maps and reflection files can be compared directly. Moving the model would leave // each dataset in whichever indexing its own run happened to pick. // A model handed over already placed in the data's cell (molecular replacement) is on the data's axes. - if (data_cell.is_crystal() && model == nullptr) { + if (data_cell.is_crystal() && placed == nullptr) { const gemmi::GroupOps model_gops = sg->operations(); const bool same_point_group = data_space_group != nullptr && same_rotations(data_space_group->operations(), model_gops); @@ -1921,7 +1916,7 @@ ModelValidationGPUPlan PlanModelValidationGPU(const ModelValidationGPUNeed &need ModelValidationResult ValidateAgainstModel(const std::vector &merged, const UnitCell &cell, - const std::string &model_path, + const RugnuxModel &model, const std::string &output_prefix, Logger &logger, const gemmi::SpaceGroup *data_space_group, @@ -1931,18 +1926,18 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const std::vector &report_shell_d_min, const ModelValidationSchedule &schedule, const gemmi::Op *twin_law, - const gemmi::Structure *model) { + const gemmi::Structure *placed) { #ifdef JFJOCH_USE_CUDA // A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through: // the validation reports why it did not finish, and a validation that did not finish decides nothing, // so the reflection files are those of a run without a model. try { - return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, + return Validate(merged, cell, model, output_prefix, logger, data_space_group, probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, - true, model); + true, placed); } catch (const RigidBodyGPUFailure &e) { ModelValidationResult failed; - failed.model_path = model_path; + failed.model_name = model.name; failed.failure_reason = fmt::format("the GPU failed during the validation ({})", e.what()); logger.Error("Model validation: {}", failed.failure_reason); cuda_throw_if_context_lost(); // a lost device ends the run here, not at its next user @@ -1950,8 +1945,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & return failed; } #else - return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, probe_indexing_ambiguity, - nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false, model); + return Validate(merged, cell, model, output_prefix, logger, data_space_group, probe_indexing_ambiguity, + nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false, placed); #endif } @@ -2107,7 +2102,7 @@ void KeepModelVerdict(ModelValidationResult &remade, const ModelValidationResult remade.model_space_group_number = decided.model_space_group_number; } -std::vector ModelReferenceIntensities(const std::string &model_path, +std::vector ModelReferenceIntensities(const RugnuxModel &model, const std::optional &cell, const gemmi::SpaceGroup *space_group, double d_min, @@ -2118,19 +2113,15 @@ std::vector ModelReferenceIntensities(const std::string &model return out; } - gemmi::Structure st; + gemmi::Structure st = model.structure; try { - // Detect, not the default: without it GEMMI picks the format from the extension and only - // falls back to the content when it does not recognise one. A model arrives named however - // whoever produced it named it, so the file itself is the better authority. - st = gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect); expand_strict_ncs(st); } catch (const std::exception &e) { - logger.Error("Model reference: cannot read model {}: {}", model_path, e.what()); + logger.Error("Model reference: cannot read model {}: {}", model.name, e.what()); return out; } if (st.models.empty() || !st.cell.is_crystal()) { - logger.Error("Model reference: model {} has no atoms or no unit cell", model_path); + logger.Error("Model reference: model {} has no atoms or no unit cell", model.name); return out; } @@ -2153,7 +2144,7 @@ std::vector ModelReferenceIntensities(const std::string &model st.spacegroup_hm = space_group->xhm(); const gemmi::SpaceGroup *sg = st.find_spacegroup(); if (!sg) { - logger.Error("Model reference: model {} has no usable space group", model_path); + logger.Error("Model reference: model {} has no usable space group", model.name); return out; } st.setup_cell_images(); @@ -2203,6 +2194,6 @@ std::vector ModelReferenceIntensities(const std::string &model } logger.Info("Model reference: {} intensities computed from {} to {:.2f} A, space group {}", - out.size(), model_path, d_min, sg->short_name()); + out.size(), model.name, d_min, sg->short_name()); return out; } diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index 4bc05dff7..594c3fb65 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -15,6 +15,7 @@ #include "../../common/Reflection.h" // MergedReflection #include "../../common/UnitCell.h" +#include "RugnuxModel.h" class Logger; namespace gemmi { struct Structure; } @@ -22,7 +23,7 @@ namespace gemmi { struct Structure; } // Result of validating merged intensities against an atomic model. struct ModelValidationResult { bool ok = false; - std::string model_path; // the file that was asked for, so the report can name it + std::string model_name; // RugnuxModel::name, so the report can name the model that was asked for // Why it did not run, when ok is false: the file could not be read, had no atoms, no cell or no // usable space group. A caller that asked for a model must be able to say so rather than finish // quietly with no R-free and no maps, which reads exactly like a run that was never asked. @@ -208,8 +209,7 @@ struct ModelValidationResult { std::shared_ptr placed_model; }; -// Given merged intensities and an atomic model (PDB or mmCIF, gzipped or not - the format is taken -// from the file's content), scale the model structure factors (with an +// Given merged intensities and an atomic model, scale the model structure factors (with an // optimized flat bulk-solvent contribution and an overall anisotropic B) to the WORKING observed // amplitudes, then report R-work / R-free and write 2Fo-Fc and Fo-Fc electron-density maps // (CCP4) plus an MTZ of map coefficients next to output_prefix. Nothing the reported R-free is @@ -242,8 +242,8 @@ struct ModelValidationResult { // twin_law, where given, is a twin law of the data's group (a Miller-index operator in the data's // basis): the model is then also scored against the data as twinned by it (twin_target and below). // -// model, where given, is the model itself - one already placed in the data's cell, by molecular -// replacement - and model_path only names it; no other setting is then probed for it. +// placed, where given, is the model itself - one already placed in the data's cell, by molecular +// replacement - and `model` only names it; no other setting is then probed for it. // // report_shell_d_min is the merge statistics' own shell bounds, coarse to fine, and it is what // cc_model_shells is binned on. Sharing the grid is the point: a reader has to be able to put a @@ -329,7 +329,7 @@ ModelValidationGPUPlan PlanModelValidationGPU(const ModelValidationGPUNeed &need ModelValidationResult ValidateAgainstModel(const std::vector &merged, const UnitCell &cell, - const std::string &model_path, + const RugnuxModel &model, const std::string &output_prefix, Logger &logger, const gemmi::SpaceGroup *data_space_group = nullptr, @@ -339,7 +339,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const std::vector &report_shell_d_min = {}, const ModelValidationSchedule &schedule = {}, const gemmi::Op *twin_law = nullptr, - const gemmi::Structure *model = nullptr); + const gemmi::Structure *placed = nullptr); // Reindex `merged` into the frame ValidateAgainstModel reported, so the reflection files that are // written describe the same indexing as the R-factors and the maps. Returns the space group they are @@ -378,7 +378,7 @@ ModelValidationResult ReadModelAnomalousMap(const std::vector // serial-stills run: there the ambiguity has to be broken per image, at integration time, long before // there is a merge to fit a model to. `cell` and `space_group` override the model's own, where // the run already knows them. Empty on any failure (which is logged). -std::vector ModelReferenceIntensities(const std::string &model_path, +std::vector ModelReferenceIntensities(const RugnuxModel &model, const std::optional &cell, const gemmi::SpaceGroup *space_group, double d_min, diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 47d165261..fdb9d6f15 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -18,15 +18,11 @@ #include #include #include -#include #include -#include -#include // read_structure_gz #include // remove_hydrogens, remove_alternative_conformations #include // setup_entities, remove_ligands_and_waters -#include // read_cif_gz #include // DensityCalculator #include // IT92 x-ray form factors #include // Scaling (bulk solvent + anisotropic B) @@ -892,19 +888,10 @@ void finish_search_model(SearchModel &m) { } // namespace -SearchModel PrepareSearchModel(const std::string &path) { +SearchModel PrepareSearchModel(const RugnuxModel &model_in) { SearchModel out; - bool modelcif = false; - try { - out.st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect); - if (out.st.input_format == gemmi::CoorFormat::Mmcif) { - const gemmi::cif::Document doc = gemmi::read_cif_gz(path); - modelcif = !doc.blocks.empty() && doc.blocks[0].has_mmcif_category("_ma_qa_metric"); - } - } catch (const std::exception &e) { - out.failure_reason = std::string("cannot read the model: ") + e.what(); - return out; - } + const bool modelcif = model_in.modelcif; + out.st = model_in.structure; if (out.st.models.empty()) { out.failure_reason = "the model has no atoms"; return out; @@ -978,39 +965,15 @@ SearchModel PrepareSearchModel(const std::string &path) { // --- domains from the predicted aligned error --- -std::string PaePathFor(const std::string &model_path) { - const size_t slash = model_path.find_last_of("/\\"); - const std::string dir = slash == std::string::npos ? std::string() : model_path.substr(0, slash + 1); - const std::string name = model_path.substr(dir.size()); - const size_t at = name.rfind("-model_v"); - if (at == std::string::npos) - return {}; - const size_t dot = name.find('.', at); - const std::string version = name.substr(at + 8, dot == std::string::npos ? std::string::npos : dot - at - 8); - return dir + name.substr(0, at) + "-predicted_aligned_error_v" + version + ".json"; -} - // A prediction cut into domains by its predicted aligned error before molecular replacement, as Oeffner et // al. (2022) Acta Cryst. D78, 1303-1314 do; here by average linkage. -std::vector> PaeDomains(const std::string &pae_path, const SearchModel &model) { +std::vector> PaeDomains(const std::vector> &pae, const SearchModel &model) { // Two domains a prediction does not place relative to each other with confidence are this far apart in // mean PAE (within one domain it is a few angstroms); a domain is at least this many residues. constexpr double PAE_DOMAIN_A = 10.0; constexpr int DOMAIN_MIN_RESIDUES = 40; // The clustering starts from runs of this many consecutive residues kept. constexpr int PIECE = 5; - std::vector> pae; - try { - std::ifstream in(pae_path); - if (!in) - return {}; - const nlohmann::json j = nlohmann::json::parse(in); - const nlohmann::json &m = j.is_array() ? j.at(0).at("predicted_aligned_error") : j.at("predicted_aligned_error"); - for (const auto &row : m) - pae.push_back(row.get>()); - } catch (const std::exception &) { - return {}; - } const size_t n = pae.size(); if (n != static_cast(model.residues_read) || std::any_of(pae.begin(), pae.end(), [n](const std::vector &r) { return r.size() != n; })) @@ -3323,7 +3286,7 @@ MolecularReplacementGPUPlan PlanMolecularReplacementGPU(const MolecularReplaceme return plan; } -MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &merges, +MolecularReplacementResult RunMolecularReplacement(const RugnuxModel &input, const std::vector &merges, Logger &logger, size_t nthreads, size_t gpu_resident_bytes) { const auto t0 = std::chrono::steady_clock::now(); // What the search may hold on the card: the budget model validation keeps to. What the stream-ordered pool @@ -3338,25 +3301,24 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path #else (void) gpu_resident_bytes; #endif - SearchModel model = PrepareSearchModel(model_path); + SearchModel model = PrepareSearchModel(input); if (!model.failure_reason.empty()) { MolecularReplacementResult result; - result.model_path = model_path; + result.model_name = input.name; result.residues_read = model.residues_read; result.failure_reason = model.failure_reason; return result; } model.description = "the whole model"; // Only a prediction has a predicted aligned error. - const std::string pae_path = model.predicted ? PaePathFor(model_path) : std::string(); - const std::vector> domains = pae_path.empty() ? std::vector>() - : PaeDomains(pae_path, model); + const bool pae = model.predicted && input.pae.has_value(); + const std::vector> domains = pae ? PaeDomains(*input.pae, model) : std::vector>(); const double t_model = seconds_since(t0); - if (!pae_path.empty()) - logger.Info("Molecular replacement: {} rigid domain(s) in the predicted aligned error {}", domains.size(), pae_path); + if (pae) + logger.Info("Molecular replacement: {} rigid domain(s) in the predicted aligned error", domains.size()); MolecularReplacementResult result = search_model(model, merges, logger, nthreads, gpu_budget); - result.model_path = model_path; + result.model_name = input.name; result.pae_domains = static_cast(domains.size()); result.t_model += t_model; // Where the whole model is placed nowhere and its prediction falls apart into domains whose relative @@ -3403,7 +3365,7 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path found[best].candidates[found[best].chosen].space_group.xhm()); } else if (best >= 0 && agree) { MolecularReplacementResult chosen = std::move(found[best]); - chosen.model_path = model_path; + chosen.model_name = input.name; chosen.pae_domains = result.pae_domains; chosen.t_model = result.t_model; chosen.t_transform = result.t_transform; diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index ea163349a..29105a330 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -35,6 +35,7 @@ #include "../../common/Reflection.h" // MergedReflection #include "../../common/UnitCell.h" +#include "RugnuxModel.h" class Logger; @@ -55,19 +56,16 @@ struct SearchModel { }; // A prediction where the file says so - ModelCIF's _ma_qa_metric - or where it has no cell and its B column // holds one value per residue, from 0 to 100; an experimental structure otherwise. -SearchModel PrepareSearchModel(const std::string &path); +SearchModel PrepareSearchModel(const RugnuxModel &model); -// The rigid domains of a prediction from its predicted aligned error (PAE, the AlphaFold DB's JSON: a matrix +// The rigid domains of a prediction from its predicted aligned error (PAE, the AlphaFold DB's matrix // "predicted_aligned_error" over the residues of the whole prediction): the kept residues clustered by average // linkage on the symmetrised PAE until every two clusters are a mean PAE_DOMAIN_A apart; clusters of at least -// DOMAIN_MIN_RESIDUES, largest first, each a list of positions in the prediction. Empty where the file cannot -// be read or does not match the model. -std::vector> PaeDomains(const std::string &pae_path, const SearchModel &model); +// DOMAIN_MIN_RESIDUES, largest first, each a list of positions in the prediction. Empty where the matrix does +// not match the model. +std::vector> PaeDomains(const std::vector> &pae, const SearchModel &model); // The search model of one domain: the residues of `model` at those positions, recentred. SearchModel DomainModel(const SearchModel &model, const std::vector &residues); -// The PAE file the AlphaFold DB puts beside a model: AF--F1-model_v.cif -> -// AF--F1-predicted_aligned_error_v.json in the same directory; empty where the name does not follow it. -std::string PaePathFor(const std::string &model_path); // Rotations covering SO(3) with neighbours about `step_deg` apart: the identity, then shells of // rotation angle kappa = step, 2 step, ... 180 deg, each with its axes on a golden-spiral lattice of @@ -341,7 +339,7 @@ struct MRCandidate { struct MolecularReplacementResult { bool ok = false; std::string failure_reason; - std::string model_path; + std::string model_name; // RugnuxModel::name int residues_read = 0, residues_kept = 0; double mass_da = 0; std::string search_model; // the whole model, or the PAE domain the decision was made with @@ -405,7 +403,7 @@ struct MolecularReplacementGPUPlan { // to be done rather than spilling past it. MolecularReplacementGPUPlan PlanMolecularReplacementGPU(const MolecularReplacementGPUNeed &need); -// Search `merges` (the answer's first) for the predicted model at `model_path`, in every space group with +// Search `merges` (the answer's first) for the predicted `model`, in every space group with // the rotations of each merge's group - one copy, then further copies one at a time while each new one // stands out, packs and lowers R-free - and decide. A candidate is PLACED where a copy's translation // function stands out (TF Z >= 8), the copies pack, one copy fits the cell, and its R-work beats its own @@ -413,7 +411,7 @@ MolecularReplacementGPUPlan PlanMolecularReplacementGPU(const MolecularReplaceme // every other by an R-free margin; it chooses the twin hypothesis' merge only where none of the // answer's groups can hold the model at all (Matthews, packing) and a group of the subgroup can. // `gpu_resident_bytes`: what is held on the card meanwhile by what came before (a merge engine still at work). -MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &merges, +MolecularReplacementResult RunMolecularReplacement(const RugnuxModel &model, const std::vector &merges, Logger &logger, size_t nthreads, size_t gpu_resident_bytes = 0); // The hand check of an ENANTIOMORPH decision, filled into result.hand_check: taken into the answer by diff --git a/image_analysis/structure_refinement/RugnuxModel.h b/image_analysis/structure_refinement/RugnuxModel.h new file mode 100644 index 000000000..92253b2f7 --- /dev/null +++ b/image_analysis/structure_refinement/RugnuxModel.h @@ -0,0 +1,20 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "gemmi/model.hpp" + +// An atomic model handed to a run (rugnux --model / --mr, the broker's upload), as its file held it. The +// caller parses it (ParseRugnuxModel, writer/RugnuxModelFile.h); the analysis reads no file. +struct RugnuxModel { + std::string name; // what the report calls it: the file it was given as + gemmi::Structure structure; // as read, input_format kept + bool modelcif = false; // a ModelCIF file (_ma_qa_metric): a prediction, whatever its B column holds + // A prediction's predicted aligned error over the residues of the whole prediction, where it came with one. + std::optional>> pae; +}; diff --git a/receiver/JFJochReceiverRugnux.cpp b/receiver/JFJochReceiverRugnux.cpp index 2254855ca..ec3084643 100644 --- a/receiver/JFJochReceiverRugnux.cpp +++ b/receiver/JFJochReceiverRugnux.cpp @@ -86,7 +86,7 @@ JFJochReceiverRugnux::JFJochReceiverRugnux(const DiffractionExperiment &in_exper Logger &in_logger, int64_t in_forward_nthreads, int in_rugnux_nthreads, - const std::string &in_model_path, + std::shared_ptr in_model, double in_deadline_s, RugnuxLiveStatus &in_live_status, const SpotFindingSettings &in_spot_finding_settings, @@ -102,7 +102,7 @@ JFJochReceiverRugnux::JFJochReceiverRugnux(const DiffractionExperiment &in_exper reader(in_image_buffer), forward_nthreads(JFJochReceiverLite::NumberOfDataAnalysisThreads(in_forward_nthreads, in_experiment)), rugnux_nthreads(in_rugnux_nthreads), - model_path(in_model_path), + model(std::move(in_model)), deadline_s(in_deadline_s), live_status(in_live_status), measurement_started(1), @@ -246,7 +246,7 @@ void JFJochReceiverRugnux::RugnuxThread() { parallel_detail::WorkerPool::Instance().SetNice(10); try { const auto dataset = reader.GetDataset(); - const RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, scratch_dir + "/p", model_path); + const RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, scratch_dir + "/p", model); LiveObserver observer(plots, scan_result, *az_int_mapping, live_status, logger); ProcessResult result; { @@ -288,7 +288,7 @@ void JFJochReceiverRugnux::ProbeThread() { const int64_t frames = std::lround(PROBE_DEG / std::abs(gonio->GetIncrement_deg())); if (2 * frames > static_cast(reader.GetNumberOfImages())) return; // a sweep this short is decided on soon enough by the run itself - RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, "", ""); + RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, "", nullptr); setup.config.end_image = static_cast(frames); setup.config.rotation_postrefine_geometry = false; setup.config.run_scaling = false; diff --git a/receiver/JFJochReceiverRugnux.h b/receiver/JFJochReceiverRugnux.h index 1e480e4df..22a3e9742 100644 --- a/receiver/JFJochReceiverRugnux.h +++ b/receiver/JFJochReceiverRugnux.h @@ -29,7 +29,7 @@ class JFJochReceiverRugnux : public JFJochReceiver { const int64_t forward_nthreads; const int rugnux_nthreads; - const std::string model_path; + const std::shared_ptr model; const double deadline_s; // after the last image, the analysis is cut (0 = never) std::string scratch_dir; // rugnux's own output, read back into the end message RugnuxLiveStatus &live_status; @@ -66,7 +66,7 @@ public: Logger &logger, int64_t forward_nthreads, int rugnux_nthreads, - const std::string &model_path, + std::shared_ptr model, double deadline_s, RugnuxLiveStatus &live_status, const SpotFindingSettings &spot_finding_settings, diff --git a/receiver/JFJochReceiverService.cpp b/receiver/JFJochReceiverService.cpp index ce0c633c3..36a65245d 100644 --- a/receiver/JFJochReceiverService.cpp +++ b/receiver/JFJochReceiverService.cpp @@ -20,14 +20,14 @@ JFJochReceiverService::JFJochReceiverService(AcquisitionDeviceGroup &in_aq_devic spot_finding_settings(DiffractionExperiment::DefaultDataProcessingSettings()) { } -JFJochReceiverService &JFJochReceiverService::RugnuxLive(bool enable, const std::string &model_path, +JFJochReceiverService &JFJochReceiverService::RugnuxLive(bool enable, std::shared_ptr model, int reserved_cpus, double deadline_s) { std::unique_lock ul_state(state_mutex); if (state != ReceiverState::Idle) throw JFJochException(JFJochExceptionCategory::WrongDAQState, "Cannot change rugnux settings during data collection"); rugnux_enable = enable; - rugnux_model_path = model_path; + rugnux_model = std::move(model); rugnux_reserved_cpus = reserved_cpus; rugnux_deadline_s = deadline_s; return *this; @@ -127,7 +127,7 @@ void JFJochReceiverService::Start(const DiffractionExperiment &experiment, 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_path, rugnux_deadline_s, rugnux_status, spot_finding_settings, receiver_status, plots, image_buffer, + rugnux_model, rugnux_deadline_s, rugnux_status, spot_finding_settings, receiver_status, plots, image_buffer, zmq_preview_socket.get(), zmq_metadata_socket.get()); break; } diff --git a/receiver/JFJochReceiverService.h b/receiver/JFJochReceiverService.h index 57a9985a5..1eb6d4df6 100644 --- a/receiver/JFJochReceiverService.h +++ b/receiver/JFJochReceiverService.h @@ -14,6 +14,8 @@ #include "../preview/ZMQMetadataSocket.h" #include "../preview/ZMQPreviewSocket.h" +struct RugnuxModel; + class JFJochReceiverService { std::unique_ptr receiver; AcquisitionDeviceGroup &aq_devices; @@ -26,7 +28,7 @@ class JFJochReceiverService { // Rotation collections from a DECTRIS detector analysed by rugnux while they arrive // (JFJochReceiverRugnux), instead of image by image. Off unless asked for. bool rugnux_enable = false; - std::string rugnux_model_path; + std::shared_ptr rugnux_model; int rugnux_reserved_cpus = 4; // left to the forward threads and the pusher double rugnux_deadline_s = 0; // the analysis is cut this long after the last image (0 = never) RugnuxLiveStatus rugnux_status; // of the last collection, kept after it ends @@ -53,7 +55,7 @@ public: ImagePusher &pusher, size_t send_buffer_size_MiB = 512); JFJochReceiverService& NumThreads(int64_t input); - JFJochReceiverService& RugnuxLive(bool enable, const std::string &model_path, int reserved_cpus, + JFJochReceiverService& RugnuxLive(bool enable, std::shared_ptr model, int reserved_cpus, double deadline_s); RugnuxLiveSnapshot GetRugnuxStatus() const; JFJochReceiverService& PreviewSocket(const std::string &addr, const std::optional &watermark = {}); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index bda044f24..9e56e7da6 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2329,7 +2329,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, ModelValidationResult unplaced; if (!result.model_validation) { const auto &mp = *result.model_prior; - unplaced.model_path = mp.model_path; + unplaced.model_name = mp.model_name; unplaced.failure_reason = mp.ok ? "molecular replacement did not place the model" : "the model could not be searched for: " + mp.failure_reason; } @@ -2338,7 +2338,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, const bool mr_hand = by_mr && result.model_prior->hand_check.adopted; ReportSection s; s.title = "5. MODEL VALIDATION"; - Add(s, KeyText("MODEL_FILE", mv.model_path)); + Add(s, KeyText("MODEL_FILE", mv.model_name)); if (!result.model_placement.empty()) Add(s, KeyEnum("MODEL_PLACEMENT", result.model_placement, {"RIGID_BODY", "MR", "NOT_PLACED"})); if (!mv.ok) { diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index a7c231288..23c1dfedb 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -967,12 +967,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::vector model_reference; if (!config_.reference_data.empty()) { indexer->ReferenceIntensities(config_.reference_data); - } else if (!config_.model_path.empty() && !experiment_.IsRotationIndexing() + } else if (config_.model && !experiment_.IsRotationIndexing() && experiment_.GetUnitCell().has_value() && experiment_.GetGemmiSpaceGroup().has_value()) { const auto d_min = experiment_.GetScalingSettings().GetHighResolutionLimit_A() .value_or(experiment_.GetDetectorMaxResolution_A()); model_reference = ModelReferenceIntensities( - config_.model_path, experiment_.GetUnitCell(), + *config_.model, experiment_.GetUnitCell(), &experiment_.GetSpaceGroupOrP1(), d_min, logger); if (!model_reference.empty()) { indexer->ReferenceIntensities(model_reference); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index afbdca2aa..c70acc7da 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -207,10 +207,10 @@ struct ProcessConfig { // run can write both. bool export_unmerged_partials = false; - // If set, after merging, validate the merged intensities against this atomic model (PDB or - // mmCIF, gzipped or not): compute R-free with an optimized bulk solvent and write 2Fo-Fc / - // Fo-Fc maps (--model). - std::string model_path; + // If set, after merging, validate the merged intensities against this atomic model: compute + // R-free with an optimized bulk solvent and write 2Fo-Fc / Fo-Fc maps (--model). The caller parses + // it (ParseRugnuxModel, writer/RugnuxModelFile.h). + std::shared_ptr model; // --mr: the model is searched for by molecular replacement straight away, without the rigid-body // placement first. Without it a --model that is a prediction, or an experimental structure the rigid body diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 1de2ddb02..4dd6d6dc7 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -82,8 +82,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config, if (!config.reference_column.empty()) add("--reference-column", config.reference_column); } - if (!config.model_path.empty()) - add(config.force_mr ? "--mr" : "--model", config.model_path); + if (config.model) + add(config.force_mr ? "--mr" : "--model", config.model->name); } if (config.start_image != 0) add("-s", std::to_string(config.start_image)); diff --git a/rugnux/RugnuxLive.cpp b/rugnux/RugnuxLive.cpp index afd1bc3f1..b05cf60b8 100644 --- a/rugnux/RugnuxLive.cpp +++ b/rugnux/RugnuxLive.cpp @@ -6,7 +6,7 @@ #include "RugnuxDefaults.h" RugnuxLiveSetup MakeRugnuxLiveSetup(const JFJochReaderDataset &dataset, int nthreads, - const std::string &output_prefix, const std::string &model_path) { + const std::string &output_prefix, std::shared_ptr model) { RugnuxLiveSetup s; // The experiment as rugnux_cli.cpp sets it up for a rotation sweep with no options given. @@ -48,7 +48,7 @@ RugnuxLiveSetup MakeRugnuxLiveSetup(const JFJochReaderDataset &dataset, int nthr c.beam_center_check = true; c.measure_spot_budget = true; c.run_scaling = true; - c.model_path = model_path; + c.model = std::move(model); c.live = true; c.fast_live = true; // Decided late, on the whole sweep, as it re-checks its radius there - which fast_live does not. diff --git a/rugnux/RugnuxLive.h b/rugnux/RugnuxLive.h index c4aba1297..c1e0820e8 100644 --- a/rugnux/RugnuxLive.h +++ b/rugnux/RugnuxLive.h @@ -18,4 +18,4 @@ struct RugnuxLiveSetup { }; RugnuxLiveSetup MakeRugnuxLiveSetup(const JFJochReaderDataset &dataset, int nthreads, - const std::string &output_prefix, const std::string &model_path); + const std::string &output_prefix, std::shared_ptr model); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 78ea9d081..6cb1ff152 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1571,8 +1571,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { size_t model_rotations = 0; std::optional model_group; gemmi::UnitCell model_cell; - if (!config_.model_path.empty() && write_files && !geometry_prepass) { - const SearchModel m = PrepareSearchModel(config_.model_path); + if (config_.model && write_files && !geometry_prepass) { + const SearchModel m = PrepareSearchModel(*config_.model); model_predicted = m.failure_reason.empty() && m.predicted; if (const gemmi::SpaceGroup *sg = m.st.find_spacegroup(); sg != nullptr && m.st.cell.is_crystal()) { model_rotations = sg->operations().sym_ops.size(); @@ -3073,7 +3073,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // comparison), so with one the merge in the subgroup is made wherever a metric re-ask promoted from the // indexed lattice's group, or the model is in a group with fewer rotations than the adopted one. const bool model_tests_twin = - !config_.model_path.empty() + config_.model && (twin_metric_subgroup.has_value() || (!model_predicted && model_rotations > 0 && model_rotations < result.space_group->operations().sym_ops.size())); @@ -4000,7 +4000,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { Note("cut: the model is left out"); return true; }; - const bool model_run = result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass + const bool model_run = result.consensus_cell && write_files && config_.model && !geometry_prepass && !superseded && !model_cut(); // The twin hypothesis the model is asked about: its own subgroup's merge where it has one, the data's // otherwise - and whether that merge was made only for the model, so is written only where it chooses it. @@ -4096,7 +4096,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { schedule.gpu_resident_bytes = merge_on_card; schedule.gpu_beside_bytes = beside; ahead.merged = std::move(merged); - return ValidateAgainstModel(ahead.merged, cell, config_.model_path, config_.output_prefix, ahead.log, + return ValidateAgainstModel(ahead.merged, cell, *config_.model, config_.output_prefix, ahead.log, in, /*probe_indexing_ambiguity=*/false, static_cast(config_.nthreads), wavelength, report_shell_d_min, schedule); @@ -4115,7 +4115,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // validation decides: relabelling the outcomes and the group is what it must not see happen. // Every such relabelling reaches the P1 merge afterwards, through merge_to_written. ModelValidationResult validation = - ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, + ValidateAgainstModel(sm.merged, *result.consensus_cell, *config_.model, config_.output_prefix, logger, data_sg ? &*data_sg : nullptr, /*probe_indexing_ambiguity=*/config_.reference_data.empty(), static_cast(config_.nthreads), experiment_.GetWavelength_A(), @@ -4199,7 +4199,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { ModelValidationSchedule kept_verdict; kept_verdict.fit_test = false; kept_verdict.gpu_resident_bytes = merge_on_card; - remade = ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, + remade = ValidateAgainstModel(sm.merged, *result.consensus_cell, *config_.model, config_.output_prefix, logger, sg_now ? &*sg_now : nullptr, /*probe_indexing_ambiguity=*/false, @@ -4246,7 +4246,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { twin_schedule.maps_superseded = [&](const ModelValidationResult &v) { return !maps_from(v); }; phase("Validating against model (twin hypothesis)"); auto twin_validation = ValidateAgainstModel( - mt.merged, mt.cell, config_.model_path, config_.output_prefix, logger, + mt.merged, mt.cell, *config_.model, config_.output_prefix, logger, &mt.space_group, /*probe_indexing_ambiguity=*/config_.reference_data.empty(), static_cast(config_.nthreads), experiment_.GetWavelength_A(), {}, twin_schedule, &mt.law); @@ -4328,7 +4328,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (model_twin) merges.push_back({"TWIN_SUBGROUP_MERGE", &model_twin->merged, model_twin->cell, model_twin->space_group, model_twin->law}); - MolecularReplacementResult mr = RunMolecularReplacement(config_.model_path, merges, logger, + MolecularReplacementResult mr = RunMolecularReplacement(*config_.model, merges, logger, static_cast(config_.nthreads), merge_on_card); if (mr.chosen >= 0) { const MRCandidate &c = mr.candidates[mr.chosen]; @@ -4371,7 +4371,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { placed_elsewhere.gpu_resident_bytes = merge_on_card; phase("Validating against model (molecular replacement)"); ModelValidationResult v = ValidateAgainstModel( - own ? own->merged : sm.merged, cell, config_.model_path, + own ? own->merged : sm.merged, cell, *config_.model, own ? result.model_prior_file_prefix : config_.output_prefix, logger, &data_sg, /*probe_indexing_ambiguity=*/false, static_cast(config_.nthreads), experiment_.GetWavelength_A(), report_shell_d_min, placed_elsewhere, nullptr, c.placed.get()); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e218c8a70..c77b05cb1 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -70,6 +70,7 @@ #include "RugnuxDefaults.h" #include "../image_analysis/structure_refinement/ModelValidation.h" #include "../writer/WriteModel.h" +#include "../writer/RugnuxModelFile.h" #include "ResultReport.h" #include "RugnuxConsole.h" #ifdef __APPLE__ @@ -1609,6 +1610,16 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { logger.Error("--model and --mr both name the model; give one"); exit(EXIT_FAILURE); } + // Read before anything runs: a file that holds no model ends the run here rather than after it. + std::shared_ptr model; + if (const std::string path = mr_model.empty() ? model_pdb : mr_model; !path.empty()) { + try { + model = std::make_shared(ReadRugnuxModelFile(path)); + } catch (const std::exception &e) { + logger.Error("Cannot read the model {}: {}", path, e.what()); + exit(EXIT_FAILURE); + } + } if ((indexing_algorithm == IndexingAlgorithmEnum::FFT || indexing_algorithm == IndexingAlgorithmEnum::FFBIDX) && get_gpu_count() == 0) { @@ -2300,7 +2311,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // looking exactly like one that was never given --model. std::string model_validation_failure; std::optional model_validation; - if (!output_prefix.empty() && !model_pdb.empty()) { + if (!output_prefix.empty() && model) { console->OnStep("Model validation"); const auto data_sg = experiment.GetGemmiSpaceGroup(); // With a reference MTZ the merohedral indexing was already resolved (stills per-image @@ -2310,7 +2321,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { for (const auto &sh : merged_statistics.shells) report_shell_d_min.push_back(sh.d_min); const auto validation = - ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb, + ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), *model, output_prefix, logger, data_sg ? &*data_sg : nullptr, /*probe_indexing_ambiguity=*/reference_data.empty(), @@ -3135,7 +3146,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.observation_dump_path = dump_observations; config.export_unmerged = export_unmerged; config.export_unmerged_partials = export_unmerged_partials; - config.model_path = mr_model.empty() ? model_pdb : mr_model; + config.model = model; config.force_mr = !mr_model.empty(); // When merging, the merged reflections (.mtz/.cif) are the wanted output; skip the large // _process.h5 unless explicitly requested. Without merging, the _process.h5 is the only output. diff --git a/tests/ModelScaleGPUTest.cpp b/tests/ModelScaleGPUTest.cpp index 86c167f83..90a3230b1 100644 --- a/tests/ModelScaleGPUTest.cpp +++ b/tests/ModelScaleGPUTest.cpp @@ -30,6 +30,7 @@ #include "../image_analysis/structure_refinement/ModelScaleGPU.h" #include "../image_analysis/structure_refinement/ModelValidation.h" #include "../image_analysis/structure_refinement/RigidBodyRefine.h" +#include "../writer/RugnuxModelFile.h" namespace { @@ -104,7 +105,7 @@ gemmi::Scaling ModelPoints(const char *cryst, double zone) { const std::string path = "model_scale_gpu_test.pdb"; std::ofstream(path) << pdb; Logger logger("ModelScaleGPUTest"); - const auto reference = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto reference = ModelReferenceIntensities(ReadRugnuxModelFile(path), {}, {}, 3.0, logger); gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect); std::filesystem::remove(path); st.setup_cell_images(); diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index cec04d31d..131ef6712 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -31,6 +31,7 @@ #endif #include "../image_analysis/structure_refinement/SigmaA.h" #include "../writer/WriteModel.h" +#include "../writer/RugnuxModelFile.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" namespace { @@ -107,6 +108,10 @@ ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00 f.close(); return name; } + + RugnuxModel ReadModel(const std::string &path) { + return ReadRugnuxModelFile(path); + } } // --model used to call gemmi::read_pdb unconditionally, so a deposited model handed over as mmCIF - @@ -212,8 +217,8 @@ TEST_CASE("ModelValidation_ReadsPdbAndMmcif", "[ModelValidation]") { const auto pdb = WriteTemp("model_validation_test.pdb", kPdb); const auto cif = WriteTemp("model_validation_test.cif", kMmcif); - const auto from_pdb = ModelReferenceIntensities(pdb, {}, {}, 4.0, logger); - const auto from_cif = ModelReferenceIntensities(cif, {}, {}, 4.0, logger); + const auto from_pdb = ModelReferenceIntensities(ReadModel(pdb), {}, {}, 4.0, logger); + const auto from_cif = ModelReferenceIntensities(ReadModel(cif), {}, {}, 4.0, logger); REQUIRE_FALSE(from_pdb.empty()); REQUIRE_FALSE(from_cif.empty()); @@ -222,7 +227,7 @@ TEST_CASE("ModelValidation_ReadsPdbAndMmcif", "[ModelValidation]") { // The extension is not what decides: the same mmCIF under a .pdb name still reads. const auto misnamed = WriteTemp("model_validation_test_misnamed.pdb", kMmcif); - CHECK_FALSE(ModelReferenceIntensities(misnamed, {}, {}, 4.0, logger).empty()); + CHECK_FALSE(ModelReferenceIntensities(ReadModel(misnamed), {}, {}, 4.0, logger).empty()); std::filesystem::remove(pdb); std::filesystem::remove(cif); @@ -256,8 +261,8 @@ TEST_CASE("ModelValidation_ExpandsStrictNcs", "[ModelValidation]") { const auto ncs_path = WriteTemp("model_validation_test_ncs.pdb", with_ncs); const auto full_path = WriteTemp("model_validation_test_ncs_full.pdb", expanded); - const auto from_ncs = ModelReferenceIntensities(ncs_path, {}, {}, 4.0, logger); - const auto from_full = ModelReferenceIntensities(full_path, {}, {}, 4.0, logger); + const auto from_ncs = ModelReferenceIntensities(ReadModel(ncs_path), {}, {}, 4.0, logger); + const auto from_full = ModelReferenceIntensities(ReadModel(full_path), {}, {}, 4.0, logger); REQUIRE(from_ncs.size() == from_full.size()); REQUIRE_FALSE(from_ncs.empty()); for (size_t i = 0; i < from_ncs.size(); i++) { @@ -275,22 +280,22 @@ TEST_CASE("ModelValidation_ExpandsStrictNcs", "[ModelValidation]") { TEST_CASE("ModelValidation_UnusableModelGivesAReason", "[ModelValidation]") { Logger logger("ModelValidation_UnusableModelGivesAReason"); + // A file that holds no model is refused where it is read, before anything runs. const auto empty_file = WriteTemp("model_validation_test_bogus.pdb", "not a coordinate file\n"); + CHECK_THROWS(ReadModel(empty_file)); + CHECK_THROWS(ReadModel("model_validation_test_does_not_exist.pdb")); + std::filesystem::remove(empty_file); + // A model with atoms but no crystal cell to put them in. + RugnuxModel no_cell; + no_cell.name = "no_cell.pdb"; + no_cell.structure.models.emplace_back(1); const auto result = ValidateAgainstModel({}, UnitCell{.a = 40, .b = 50, .c = 60, .alpha = 90, .beta = 90, .gamma = 90}, - empty_file, "", logger); + no_cell, "", logger); CHECK_FALSE(result.ok); CHECK_FALSE(result.failure_reason.empty()); - CHECK_THAT(result.failure_reason, Catch::Matchers::ContainsSubstring(empty_file)); - - const auto missing = ValidateAgainstModel({}, UnitCell{.a = 40, .b = 50, .c = 60, - .alpha = 90, .beta = 90, .gamma = 90}, - "model_validation_test_does_not_exist.pdb", "", logger); - CHECK_FALSE(missing.ok); - CHECK_FALSE(missing.failure_reason.empty()); - - std::filesystem::remove(empty_file); + CHECK_THAT(result.failure_reason, Catch::Matchers::ContainsSubstring("no_cell.pdb")); } // A rigid-body step is only worth taking if it can find a shift it was not told about. The check is @@ -306,7 +311,7 @@ TEST_CASE("ModelValidation_RigidBodyRecoversASmallShift", "[ModelValidation]") { // "Observed" amplitudes: the model's own structure factors, so the target's minimum is exactly // where the model started. - const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto ref = ModelReferenceIntensities(ReadModel(path), {}, {}, 3.0, logger); REQUIRE_FALSE(ref.empty()); gemmi::AsuData> fobs; fobs.unit_cell_ = st.cell; @@ -353,7 +358,7 @@ TEST_CASE("ModelValidation_RigidBodyLeavesThePolarDirectionAlone", "[ModelValida REQUIRE(sg != nullptr); st.setup_cell_images(); - const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto ref = ModelReferenceIntensities(ReadModel(path), {}, {}, 3.0, logger); REQUIRE_FALSE(ref.empty()); gemmi::AsuData> fobs; fobs.unit_cell_ = st.cell; @@ -400,7 +405,7 @@ TEST_CASE("ModelValidation_RigidBodySolventStaysPhysical", "[ModelValidation]") REQUIRE(sg != nullptr); st.setup_cell_images(); - const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto ref = ModelReferenceIntensities(ReadModel(path), {}, {}, 3.0, logger); REQUIRE_FALSE(ref.empty()); gemmi::AsuData> fobs; fobs.unit_cell_ = st.cell; @@ -487,7 +492,7 @@ TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidat // The data: the cluster's own |F|^2 in its own frame, 30 x 34 x 38, P 21 21 21. const auto data_model = WriteTemp("cob_data_test.pdb", ClusterPdb().c_str()); - auto obs = ModelReferenceIntensities(data_model, {}, {}, 2.5, logger); + auto obs = ModelReferenceIntensities(ReadModel(data_model), {}, {}, 2.5, logger); REQUIRE(obs.size() > 500); for (size_t i = 0; i < obs.size(); i++) { obs[i].F = std::sqrt(std::max(0.0f, obs[i].I)); @@ -518,7 +523,7 @@ TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidat } const auto other_setting = WriteTemp("cob_model_test.pdb", permuted.c_str()); - const auto moved = ValidateAgainstModel(obs, data_cell, other_setting, prefix, logger, data_sg, + const auto moved = ValidateAgainstModel(obs, data_cell, ReadModel(other_setting), prefix, logger, data_sg, /*probe_indexing_ambiguity=*/false, 1, 1.0); REQUIRE(moved.ok); logger.Info("Change of basis test: took {} out of {}, R-work {:.4f}", @@ -529,7 +534,7 @@ TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidat // The isomorphous case, which is what nearly every run is: the model already describes the // lattice the way the data do, the only operators that map its cell are its own symmetry, and // nothing is scored or moved. - const auto same = ValidateAgainstModel(obs, data_cell, data_model, prefix, logger, data_sg, + const auto same = ValidateAgainstModel(obs, data_cell, ReadModel(data_model), prefix, logger, data_sg, false, 1, 1.0); REQUIRE(same.ok); CHECK(same.change_of_basis_op == gemmi::Op::identity()); @@ -551,7 +556,7 @@ TEST_CASE("ModelValidation_ReindexesTheDataIntoTheModelsIndexing", "[ModelValida // Point group 4 on a tetragonal lattice (4/mmm): one alternative indexing. const auto model = WriteTemp("reidx_model_test.pdb", ClusterPdb("CRYST1 34.000 34.000 38.000 90.00 90.00 90.00 P 4 4\n").c_str()); - auto own = ModelReferenceIntensities(model, {}, {}, 2.5, logger); + auto own = ModelReferenceIntensities(ReadModel(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)); @@ -566,7 +571,7 @@ TEST_CASE("ModelValidation_ReindexesTheDataIntoTheModelsIndexing", "[ModelValida // The data as a run that picked the other indexing would have merged them. const auto obs = ReindexReflections(own, laws.front()); - const auto to_model = ValidateAgainstModel(obs, cell, model, prefix, logger, p4, + const auto to_model = ValidateAgainstModel(obs, cell, ReadModel(model), prefix, logger, p4, /*probe_indexing_ambiguity=*/true, 1, 1.0); REQUIRE(to_model.ok); CHECK(to_model.change_of_basis_op == gemmi::Op::identity()); @@ -575,7 +580,7 @@ TEST_CASE("ModelValidation_ReindexesTheDataIntoTheModelsIndexing", "[ModelValida CHECK(to_model.r_work < 0.15); // A reference already fixed the data's indexing: the data stay, and the model is moved. - const auto fixed = ValidateAgainstModel(obs, cell, model, prefix, logger, p4, false, 1, 1.0); + const auto fixed = ValidateAgainstModel(obs, cell, ReadModel(model), prefix, logger, p4, false, 1, 1.0); REQUIRE(fixed.ok); CHECK_FALSE(fixed.change_of_basis_op == gemmi::Op::identity()); CHECK(fixed.indexing_op == gemmi::Op::identity()); @@ -599,7 +604,7 @@ TEST_CASE("ModelValidation_ReindexesTheDataIntoTheModelsIndexing", "[ModelValida r.F = std::sqrt(std::max(0.0f, r.I)); } } - const auto twin = ValidateAgainstModel(twinned, cell, model, prefix, logger, p4, true, 1, 1.0); + const auto twin = ValidateAgainstModel(twinned, cell, ReadModel(model), prefix, logger, p4, true, 1, 1.0); REQUIRE(twin.ok); CHECK_FALSE(twin.indexing_decided); CHECK(twin.indexing_op == gemmi::Op::identity()); @@ -607,7 +612,7 @@ TEST_CASE("ModelValidation_ReindexesTheDataIntoTheModelsIndexing", "[ModelValida // Data already in the model's indexing: nothing moves. The null is run all the same, as on every // validation that judges its model (ModelValidationSchedule::fit_test), and the model fits. - const auto same = ValidateAgainstModel(own, cell, model, prefix, logger, p4, true, 1, 1.0); + const auto same = ValidateAgainstModel(own, cell, ReadModel(model), prefix, logger, p4, true, 1, 1.0); REQUIRE(same.ok); CHECK(same.change_of_basis_op == gemmi::Op::identity()); CHECK(same.indexing_op == gemmi::Op::identity()); @@ -736,7 +741,7 @@ TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") Logger logger("ModelValidation_CCModelFollowsTheSignalByShell"); const auto path = WriteTemp("cc_model_test.pdb", ClusterPdb().c_str()); - auto obs = ModelReferenceIntensities(path, {}, {}, 2.5, logger); + auto obs = ModelReferenceIntensities(ReadModel(path), {}, {}, 2.5, logger); REQUIRE(obs.size() > 1000); // The shells the correlation is reported on, coarse to fine, and the noise each one gets as a @@ -777,7 +782,7 @@ TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") const auto result = ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}, - path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), + ReadModel(path), prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), /*probe_indexing_ambiguity=*/false, 1, 1.0, shells); REQUIRE(result.ok); REQUIRE(result.cc_model_shells.size() == shells.size()); @@ -805,7 +810,7 @@ TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") const auto no_shells = ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}, - path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), + ReadModel(path), prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), false, 1, 1.0); CHECK(no_shells.ok); CHECK(no_shells.cc_model_shells.empty()); @@ -827,7 +832,7 @@ TEST_CASE("ModelValidation_WrongModelInAGroupOffTheDataMetric", "[ModelValidatio const auto data_path = WriteTemp("off_metric_data.pdb", ClusterPdb().c_str()); const auto wrong_path = WriteTemp("off_metric_wrong.pdb", ClusterPdb("CRYST1 40.000 40.000 40.000 90.00 90.00 90.00 I 21 3 24\n", 7).c_str()); - auto obs = ModelReferenceIntensities(data_path, {}, {}, 2.0, logger); + auto obs = ModelReferenceIntensities(ReadModel(data_path), {}, {}, 2.0, logger); REQUIRE(obs.size() > 1000); for (size_t i = 0; i < obs.size(); i++) { obs[i].sigma = std::max(1.0f, 0.05f * obs[i].I); @@ -837,7 +842,7 @@ TEST_CASE("ModelValidation_WrongModelInAGroupOffTheDataMetric", "[ModelValidatio const std::string prefix = (std::filesystem::temp_directory_path() / "off_metric_test").string(); const auto result = ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}, - wrong_path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21")); + ReadModel(wrong_path), prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21")); CHECK(result.ok); CHECK(!result.adopted_model_enantiomorph); @@ -854,7 +859,7 @@ TEST_CASE("ModelValidation_CCModelUnrelatedModelConfirmsNothing", "[ModelValidat Logger logger("ModelValidation_CCModelUnrelatedModelConfirmsNothing"); const auto data_path = WriteTemp("cc_null_data.pdb", ClusterPdb().c_str()); const auto wrong_path = WriteTemp("cc_null_wrong.pdb", ClusterPdb(kCryst, 7).c_str()); - auto obs = ModelReferenceIntensities(data_path, {}, {}, 2.0, logger); + auto obs = ModelReferenceIntensities(ReadModel(data_path), {}, {}, 2.0, logger); std::erase_if(obs, [](const MergedReflection &r) { return r.d > 8.0f; }); REQUIRE(obs.size() > 1000); for (size_t i = 0; i < obs.size(); i++) { @@ -870,7 +875,7 @@ TEST_CASE("ModelValidation_CCModelUnrelatedModelConfirmsNothing", "[ModelValidat const std::string prefix = (std::filesystem::temp_directory_path() / "cc_null_test").string(); const auto result = ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}, - wrong_path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), + ReadModel(wrong_path), prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), /*probe_indexing_ambiguity=*/false, 1, 1.0, shells); REQUIRE(result.ok); REQUIRE(result.cc_model_shells.size() == shells.size()); @@ -1003,7 +1008,7 @@ TEST_CASE("ModelValidation_RigidBodySameOnAnyNumberOfThreads", "[ModelValidation REQUIRE(sg != nullptr); st.setup_cell_images(); - const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto ref = ModelReferenceIntensities(ReadModel(path), {}, {}, 3.0, logger); REQUIRE_FALSE(ref.empty()); gemmi::AsuData> fobs; fobs.unit_cell_ = st.cell; @@ -1298,7 +1303,7 @@ namespace { gemmi::AsuData> OwnAmplitudes(const char *cryst, const gemmi::Structure &st, double d_min, Logger &logger) { const auto path = WriteTemp("rigid_body_own_amplitudes.pdb", ClusterPdb(cryst).c_str()); - const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto ref = ModelReferenceIntensities(ReadModel(path), {}, {}, 3.0, logger); std::filesystem::remove(path); gemmi::AsuData> fobs; fobs.unit_cell_ = st.cell; @@ -1327,7 +1332,7 @@ namespace { const gemmi::SpaceGroup *sg = st.find_spacegroup(); // "Observed" amplitudes: the model's own, as placed at q = 0. const auto path = WriteTemp("rigid_body_jacobian_test.pdb", ClusterPdb(cryst).c_str()); - const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto ref = ModelReferenceIntensities(ReadModel(path), {}, {}, 3.0, logger); std::filesystem::remove(path); REQUIRE_FALSE(ref.empty()); @@ -1455,7 +1460,7 @@ TEST_CASE("ModelValidation_RigidBodyReportsTheRotationItApplied", "[ModelValidat const gemmi::SpaceGroup *sg = st.find_spacegroup(); REQUIRE(sg != nullptr); st.setup_cell_images(); - const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + const auto ref = ModelReferenceIntensities(ReadModel(path), {}, {}, 3.0, logger); std::filesystem::remove(path); REQUIRE_FALSE(ref.empty()); gemmi::AsuData> fobs; diff --git a/tests/MolecularReplacementTest.cpp b/tests/MolecularReplacementTest.cpp index e03353a9f..078755225 100644 --- a/tests/MolecularReplacementTest.cpp +++ b/tests/MolecularReplacementTest.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -15,6 +16,7 @@ #include "../image_analysis/structure_refinement/MolecularReplacement.h" #include "../image_analysis/structure_refinement/RigidBodyRefine.h" #include "../common/Logger.h" +#include "../writer/RugnuxModelFile.h" #ifdef JFJOCH_USE_CUDA #include "../common/CUDAWrapper.h" #include "../image_analysis/structure_refinement/MolecularReplacementGPU.h" @@ -322,7 +324,7 @@ TEST_CASE("MolecularReplacement_ModelKind", "[MolecularReplacement]") { const float plddt[3] = {90, 85, 95}; write("MolecularReplacement_ModelKind_af.pdb", "CRYST1 1.000 1.000 1.000 90.00 90.00 90.00 P 1 1", plddt, false, false); - const SearchModel af = PrepareSearchModel("MolecularReplacement_ModelKind_af.pdb"); + const SearchModel af = PrepareSearchModel(ReadRugnuxModelFile("MolecularReplacement_ModelKind_af.pdb")); CHECK(af.failure_reason.empty()); CHECK(af.predicted); CHECK(af.residues_kept == 3); @@ -331,7 +333,7 @@ TEST_CASE("MolecularReplacement_ModelKind", "[MolecularReplacement]") { const float refined[3] = {25, 112, 40}; write("MolecularReplacement_ModelKind_xtal.pdb", "CRYST1 50.000 60.000 70.000 90.00 90.00 90.00 P 21 21 21 4", refined, true, true); - const SearchModel x = PrepareSearchModel("MolecularReplacement_ModelKind_xtal.pdb"); + const SearchModel x = PrepareSearchModel(ReadRugnuxModelFile("MolecularReplacement_ModelKind_xtal.pdb")); CHECK(x.failure_reason.empty()); CHECK(!x.predicted); CHECK(x.residues_kept == 3); @@ -375,19 +377,20 @@ TEST_CASE("MolecularReplacement_PaeDomains", "[MolecularReplacement]") { m.st.models.push_back(model); m.residues_read = 120; m.residues_kept = 110; - const std::string path = "MolecularReplacement_PaeDomains.json"; - { - std::ofstream out(path); - out << "[{\"predicted_aligned_error\":["; - for (int i = 0; i < 120; i++) { - out << (i ? "," : "") << "["; - for (int j = 0; j < 120; j++) - out << (j ? "," : "") << ((i < 55) == (j < 55) ? 2 : 25); - out << "]"; - } - out << "],\"max_predicted_aligned_error\":31.75}]"; + std::ostringstream json; + json << "[{\"predicted_aligned_error\":["; + for (int i = 0; i < 120; i++) { + json << (i ? "," : "") << "["; + for (int j = 0; j < 120; j++) + json << (j ? "," : "") << ((i < 55) == (j < 55) ? 2 : 25); + json << "]"; } - const auto domains = PaeDomains(path, m); + json << "],\"max_predicted_aligned_error\":31.75}]"; + const RugnuxModel parsed = ParseRugnuxModel( + "ATOM 1 CA ALA A 1 0.000 0.000 0.000 1.00 90.00 C\nEND\n", + "AF-P12345-F1-model_v6.pdb", json.str()); + REQUIRE(parsed.pae); + const auto domains = PaeDomains(*parsed.pae, m); REQUIRE(domains.size() == 2); CHECK(domains[0].size() == 60); CHECK(domains[0].front() == 60); @@ -395,7 +398,6 @@ TEST_CASE("MolecularReplacement_PaeDomains", "[MolecularReplacement]") { const SearchModel d = DomainModel(m, domains[1]); CHECK(d.residues_kept == 50); CHECK(d.description == "PAE domain of 50 residues (1-50)"); - std::remove(path.c_str()); } // A placement of a group carried into a subgroup's merge: in P 21 21 21 read as P 1 21 1 on the same cell each diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index d70ea191f..5d94e077d 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -383,7 +383,7 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { ModelValidationResult good; good.ok = true; - good.model_path = "model.cif"; + good.model_name = "model.cif"; good.model_space_group_number = 96; good.r_work = 0.1823; good.r_free = 0.2145; @@ -468,7 +468,7 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { mr_result.model_placement = "MR"; MolecularReplacementResult mr; mr.ok = true; - mr.model_path = "model.cif"; + mr.model_name = "model.cif"; mr.predicted_model = false; mr.decision = "ENANTIOMORPH"; MRCandidate c; @@ -568,7 +568,7 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { result.model_validation = good; ModelValidationResult failed; - failed.model_path = "broken.pdb"; + failed.model_name = "broken.pdb"; failed.failure_reason = "model broken.pdb has no atoms or no unit cell"; result.model_validation = failed; diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 74bf4c128..f6855eab4 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -274,7 +274,9 @@ TEST_CASE("RugnuxCommandLine_Full", "[process]") { CHECK(cmd.find("--model") == std::string::npos); // --mr is --model with the search forced, and is given back as itself. - config.model_path = "model.cif"; + auto model = std::make_shared(); + model->name = "model.cif"; + config.model = model; CHECK(RugnuxCommandLine(config, x, "in.h5").find("--model model.cif") != std::string::npos); config.force_mr = true; const std::string mr_cmd = RugnuxCommandLine(config, x, "in.h5"); diff --git a/viewer/windows/JFJochProcessingJobsWindow.cpp b/viewer/windows/JFJochProcessingJobsWindow.cpp index 0ba1ddbe4..b24a566e7 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.cpp +++ b/viewer/windows/JFJochProcessingJobsWindow.cpp @@ -6,6 +6,7 @@ #include "JFJochMergeStatsWindow.h" #include "../widgets/ToolbarIcons.h" #include "../../rugnux/RugnuxCommandLine.h" +#include "../../writer/RugnuxModelFile.h" #include "../../common/ThreadAffinity.h" #include @@ -396,8 +397,6 @@ ProcessConfig JFJochProcessingJobsWindow::buildConfig(const JobSpec &spec, const config.reference_space_group = inputs.reference_space_group; config.reference_path = inputs.reference_path.toStdString(); config.reference_column = inputs.reference_column.toStdString(); - // Validated at the end of the job, when there is an output to write the maps beside. - config.model_path = inputs.model_path.toStdString(); // Stills geometry-refinement two-pass. The dialog only offers it for a stills-with-cell run // (Rugnux no-ops it otherwise), so honour it whenever set; unset leaves the pass off. if (spec.refine_geometry) @@ -422,7 +421,17 @@ void JFJochProcessingJobsWindow::newJob(ProcessMode mode, CalibrationSelection c if (action == 0) return; - const ProcessConfig config = buildConfig(spec, inputs); + ProcessConfig config = buildConfig(spec, inputs); + // Validated at the end of the job, against the merge. Read here, so that a file holding no model is + // said at once rather than after the run. + if (spec.mode == ProcessMode::FullAnalysis && !inputs.model_path.isEmpty()) { + try { + config.model = std::make_shared(ReadRugnuxModelFile(inputs.model_path.toStdString())); + } catch (const std::exception &e) { + QMessageBox::warning(this, "Processing", QString("Cannot read the model: %1").arg(e.what())); + return; + } + } // The experiment carries the panel's indexing settings — including RotationIndexing set by the // rotation axis (needed so IndexAndRefine builds a rotation indexer). On top of that it gets diff --git a/writer/CMakeLists.txt b/writer/CMakeLists.txt index c96b830b7..ad4a0a3fa 100644 --- a/writer/CMakeLists.txt +++ b/writer/CMakeLists.txt @@ -36,10 +36,12 @@ ADD_LIBRARY(JFJochModelWriter STATIC WriteModel.cpp WriteModel.h) TARGET_LINK_LIBRARIES(JFJochModelWriter JFJochLogger gemmi) # rugnux's merged MTZ, map-coefficient MTZ and placed model, to and from the structured END-message -# blocks, so that jfjoch_writer writes them with gemmi alone and no analysis library. +# blocks, so that jfjoch_writer writes them with gemmi alone and no analysis library; and the atomic model +# a run is given, parsed for it. ADD_LIBRARY(JFJochRugnuxFiles STATIC WriteRugnuxFiles.cpp WriteRugnuxFiles.h - RugnuxResultBlocks.cpp RugnuxResultBlocks.h) -TARGET_LINK_LIBRARIES(JFJochRugnuxFiles JFJochModelWriter JFJochLogger gemmi) + RugnuxResultBlocks.cpp RugnuxResultBlocks.h + RugnuxModelFile.cpp RugnuxModelFile.h) +TARGET_LINK_LIBRARIES(JFJochRugnuxFiles JFJochModelWriter JFJochLogger gemmi ZLIB::ZLIB) TARGET_LINK_LIBRARIES(JFJochWriter JFJochPreview JFJochLogger JFJochHDF5Wrappers CBORStream2FrameSerialize JFJochRugnuxFiles) diff --git a/writer/RugnuxModelFile.cpp b/writer/RugnuxModelFile.cpp new file mode 100644 index 000000000..15740b333 --- /dev/null +++ b/writer/RugnuxModelFile.cpp @@ -0,0 +1,98 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RugnuxModelFile.h" + +#include +#include +#include + +#include +#include +#include + +namespace { + +// The bytes of a gzip stream inflated, every member of it. +std::string Gunzip(const std::string &bytes) { + z_stream s{}; + if (inflateInit2(&s, 15 + 32) != Z_OK) + throw std::runtime_error("cannot start zlib"); + std::string out; + char buf[1 << 16]; + s.next_in = reinterpret_cast(const_cast(bytes.data())); + s.avail_in = static_cast(bytes.size()); + int ret = Z_OK; + while (ret != Z_STREAM_END || s.avail_in > 0) { + if (ret == Z_STREAM_END) + inflateReset(&s); + s.next_out = reinterpret_cast(buf); + s.avail_out = sizeof buf; + ret = inflate(&s, Z_NO_FLUSH); + if (ret != Z_OK && ret != Z_STREAM_END) { + inflateEnd(&s); + throw std::runtime_error("the gzip stream is damaged"); + } + out.append(buf, sizeof buf - s.avail_out); + if (ret == Z_OK && s.avail_in == 0 && s.avail_out != 0) { + inflateEnd(&s); + throw std::runtime_error("the gzip stream is cut short"); + } + } + inflateEnd(&s); + return out; +} + +std::string ReadFile(const std::string &path) { + std::ifstream f(path, std::ios::binary); + if (!f) + throw std::runtime_error("cannot open " + path); + return {std::istreambuf_iterator(f), std::istreambuf_iterator()}; +} + +} // namespace + +RugnuxModel ParseRugnuxModel(const std::string &bytes, const std::string &name, const std::string &pae_json) { + RugnuxModel m; + m.name = name; + std::string text = bytes.size() >= 2 && bytes[0] == '\x1f' && bytes[1] == '\x8b' ? Gunzip(bytes) : bytes; + gemmi::cif::Document doc; + m.structure = gemmi::read_structure_from_memory(text.data(), text.size(), name, gemmi::CoorFormat::Detect, &doc); + if (m.structure.models.empty()) + throw std::runtime_error(name + " holds no atoms"); + m.modelcif = m.structure.input_format == gemmi::CoorFormat::Mmcif && !doc.blocks.empty() + && doc.blocks[0].has_mmcif_category("_ma_qa_metric"); + if (!pae_json.empty()) { + try { + const nlohmann::json j = nlohmann::json::parse(pae_json); + const nlohmann::json &pae = j.is_array() ? j.at(0).at("predicted_aligned_error") + : j.at("predicted_aligned_error"); + std::vector> rows; + for (const auto &row : pae) + rows.push_back(row.get>()); + m.pae = std::move(rows); + } catch (const std::exception &) { + } + } + return m; +} + +RugnuxModel ReadRugnuxModelFile(const std::string &path) { + const std::string pae_path = PaePathFor(path); + std::string pae_json; + if (!pae_path.empty() && std::ifstream(pae_path)) + pae_json = ReadFile(pae_path); + return ParseRugnuxModel(ReadFile(path), path, pae_json); +} + +std::string PaePathFor(const std::string &model_path) { + const size_t slash = model_path.find_last_of("/\\"); + const std::string dir = slash == std::string::npos ? std::string() : model_path.substr(0, slash + 1); + const std::string name = model_path.substr(dir.size()); + const size_t at = name.rfind("-model_v"); + if (at == std::string::npos) + return {}; + const size_t dot = name.find('.', at); + const std::string version = name.substr(at + 8, dot == std::string::npos ? std::string::npos : dot - at - 8); + return dir + name.substr(0, at) + "-predicted_aligned_error_v" + version + ".json"; +} diff --git a/writer/RugnuxModelFile.h b/writer/RugnuxModelFile.h new file mode 100644 index 000000000..5f242e6b1 --- /dev/null +++ b/writer/RugnuxModelFile.h @@ -0,0 +1,23 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../image_analysis/structure_refinement/RugnuxModel.h" + +// The atomic model of a run from the bytes of its file: PDB or mmCIF, gzipped or not, the format taken from +// the content, never from a name. `pae_json` is the AlphaFold DB's predicted-aligned-error JSON of a +// prediction, empty where there is none; one that cannot be read is left out, as one never given. `name` is +// what the report calls the model, and what gemmi names a PDB structure after. Throws where the bytes hold no +// atomic model. +RugnuxModel ParseRugnuxModel(const std::string &bytes, const std::string &name, const std::string &pae_json = {}); + +// The model file at `path`, with the predicted aligned error the AlphaFold DB puts beside it where it is there +// (PaePathFor). Named by `path`. +RugnuxModel ReadRugnuxModelFile(const std::string &path); + +// The PAE file the AlphaFold DB puts beside a model: AF--F1-model_v.cif -> +// AF--F1-predicted_aligned_error_v.json in the same directory; empty where the name does not follow it. +std::string PaePathFor(const std::string &model_path);