diff --git a/common/CUDAWrapper.cpp b/common/CUDAWrapper.cpp index 07a8e24c3..3c311adcd 100644 --- a/common/CUDAWrapper.cpp +++ b/common/CUDAWrapper.cpp @@ -166,6 +166,10 @@ bool get_gpu_free_memory(size_t &, size_t &) { void set_gpu(int32_t dev_id) {} +int32_t get_gpu() { + return 0; +} + bool gpu_direct_host_read() { return false; } diff --git a/common/CUDAWrapper.cu b/common/CUDAWrapper.cu index 961e706ce..028e53dc3 100644 --- a/common/CUDAWrapper.cu +++ b/common/CUDAWrapper.cu @@ -170,6 +170,13 @@ void enable_gpu_numa_binding() { gpu_numa_binding = true; } +int32_t get_gpu() { + int dev = 0; + if (get_gpu_count() > 0) + cuda_err(cudaGetDevice(&dev)); + return dev; +} + void pin_gpu(int32_t dev_id) { if (get_gpu_count() == 0) return; diff --git a/common/CUDAWrapper.h b/common/CUDAWrapper.h index a7bc51afa..537864081 100644 --- a/common/CUDAWrapper.h +++ b/common/CUDAWrapper.h @@ -43,6 +43,10 @@ bool get_gpu_free_memory(size_t &free_bytes, size_t &total_bytes); void set_gpu(int32_t dev_id); +// The calling thread's current device (set_gpu), for a thread it starts to go on with on the same card. +// 0 without CUDA or a device. +int32_t get_gpu(); + // Whether kernels on the current device read ordinary (pageable) host memory directly and at memory // speed, so a host buffer can be handed to a kernel instead of being copied to the card first: true on // a unified-memory system where the GPU walks the CPU's page tables (DGX Spark GB10, Grace Hopper GH200, diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.cpp b/image_analysis/scale_merge/AnisotropyAnalysis.cpp index d2432f797..28f840900 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.cpp +++ b/image_analysis/scale_merge/AnisotropyAnalysis.cpp @@ -1495,6 +1495,33 @@ gemmi::Mat33 AnisotropyTensorHKL(const AnisotropyResult &result, const gemmi::Un return cell.frac.mat.multiply(b).multiply(cell.frac.mat.transpose()); } +gemmi::Mat33 AnisotropyPriorHKL(const std::vector &merged, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *space_group) { + // The steps of AnalyzeAnisotropy up to the eigenvectors, and nothing else. + AnisotropyResult result; + if (merged.empty() || !cell.is_crystal()) + return AnisotropyTensorHKL(result, cell); + const auto allowed = AllowedDeviatoricBasis(cell, space_group); + if (allowed.empty()) { + for (int n = 0; n < 3; ++n) { + result.eigenvalue[n] = 0.0; + for (int j = 0; j < 3; ++j) + result.eigenvector[n][j] = n == j ? 1.0 : 0.0; + } + return AnisotropyTensorHKL(result, cell); + } + const auto refl = SelectReflections(merged, cell, space_group, N_SHELLS); + if (refl.size() < static_cast(MIN_CELLS * MIN_CELL_REFLECTIONS)) + return AnisotropyTensorHKL(result, cell); + const Eigen3 axes = Decompose(FitTensor(BuildCells(refl, cell, space_group, N_SHELLS), allowed, true).b); + for (int n = 0; n < 3; ++n) { + result.eigenvalue[n] = axes.value[n]; + for (int j = 0; j < 3; ++j) + result.eigenvector[n][j] = axes.vec[n][j]; + } + return AnisotropyTensorHKL(result, cell); +} + std::string AnisotropyToText(const AnisotropyResult &result) { std::ostringstream os; if (result.n_reflections == 0) diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.h b/image_analysis/scale_merge/AnisotropyAnalysis.h index 021eda984..a5bd8a18a 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.h +++ b/image_analysis/scale_merge/AnisotropyAnalysis.h @@ -193,3 +193,9 @@ std::string AnisotropyToText(const AnisotropyResult &result); // The fitted tensor in the Miller-index basis of `cell`: Q with s^T B s = h^T Q h, the form // FrenchWilsonOptions::anisotropy_hkl takes. Zero (the isotropic prior) when no tensor was fitted. gemmi::Mat33 AnisotropyTensorHKL(const AnisotropyResult &result, const gemmi::UnitCell &cell); + +// AnisotropyTensorHKL(AnalyzeAnisotropy(merged, ...), cell), the same bits, from the merge alone: the tensor +// is fitted to the merged intensities only, so this skips the unmerged observations and the rest of the +// analysis. For work that needs the amplitudes before the analysis is done. +gemmi::Mat33 AnisotropyPriorHKL(const std::vector &merged, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *space_group); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index ad9bc12f5..06cf4f0f8 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3473,6 +3473,142 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { cell, twin_sg, aniso_run); }); } + // fast_live, cut (the broker's deadline or a cancel): the model is left out, and the files are written + // as a run without one writes them. Asked again before the search below. + const auto model_cut = [&] { + if (!config_.fast_live || !cancelled_) + return false; + Note("cut: the model is left out"); + return true; + }; + const bool model_run = result.consensus_cell && write_output && config_.model && !geometry_prepass + && !superseded && !model_cut(); + // The merge's own shell bounds, so CC(model, data) is reported on the shells CC1/2 was. + std::vector report_shell_d_min; + for (const auto &sh : sm.statistics.shells) + report_shell_d_min.push_back(sh.d_min); + // Nothing indexes after the merge, so the image pass's indexers - an FFT indexer keeps its transform + // buffers and histograms on the card, 0.8 GB for the two of a rotation run - go before the model is + // validated or searched for, once the leftover-lattice census, the last to use them, has ended. + if (model_run) { + if (census_future.valid()) { + TimingMark census_mark(logger, "leftover-lattice census wait (model)"); + census_future.wait(); + } + p.indexer_pool.reset(); + // The model's GPU budget is 80% of the card less what the merge still holds, so from here the merge + // holds the card: a speculative probe not started yet waits for it to end (GPUMerge), and one let + // beside it is waited for. That is a probe pass of seconds: this long means it is stuck, and the + // model goes on as it would have without it. + if (p.gpu_merge) + p.gpu_merge->Holds(get_gpu_total_memory(0)); + TimingMark beside_mark(logger, "GPU work beside wait (model)"); + wait_for_gpu_work_beside(std::chrono::minutes(10)); + } + // The model's first validation reads of the merge only its amplitudes, which the analyses below make + // last, from the anisotropy tensor - a fit to the merge alone (AnisotropyPriorHKL). So it starts here, + // beside the analyses, on a copy of the merge given the same amplitudes: the same French-Wilson on the + // same tensor. Its log is held, and printed where it is taken; nothing it decides is acted on before. + std::vector model_merged; + const auto data_sg = experiment_.GetGemmiSpaceGroup(); + // With a reference MTZ the merohedral indexing was already resolved against it (rotation + // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. + // The validation in the model's setting further down depends on this one only through the + // setting and the indexing it settles, and both are forecast well before this one has finished + // (ModelFrameForecast). So it is started on that forecast, beside this one - on the next card + // where there is one - and kept only where the decision is the one forecast; its log is held, + // and its files wait, until then. On any other decision it is dropped, unwritten. The relabelling + // it runs on is the one AdoptModelFrame and relabel_output below make, made here on a copy. + // A model asserting the other enantiomorph is not forecast: whether the label is taken is + // only known at the end, from the anomalous map. + struct Ahead { + gemmi::Op change_of_basis_op, indexing_op; + // The reflections it validates against, relabelled as the files would be: where it is kept, + // they are the written ones, and the relabelling below is not made a second time. + std::vector merged; + Logger log = Logger::Buffered(); + std::future run; + std::promise write; // destroyed before `run`, so a run left waiting on it is released + } ahead; + // The merge is done with its engine - nothing below reads it but the finalist ledger's P1 merge, + // which runs on it beside the validation - so the device memory it holds goes back before the + // validation decides, from the card's total, what goes on its GPU. Where the ledger's merge keeps + // it, what it holds is counted against the card by every validation here, the same for each, so + // that one started on the forecast decides as the one it stands in for would. + if (model_run && !ledger_ahead.valid()) + rsm.reset(); + const size_t merge_on_card = rsm && rsm->OnGPU() ? rsm->DeviceBytes() : 0; + ModelValidationSchedule first_schedule; + first_schedule.gpu_resident_bytes = merge_on_card; + first_schedule.on_forecast = [&](const ModelFrameForecast &f) { + if (ahead.run.valid() || !config_.reference_data.empty() || !data_sg.has_value() + || !f.may_run_beside || f.enantiomorph_candidate + || f.change_of_basis_op == gemmi::Op::identity()) + return; + ahead.change_of_basis_op = f.change_of_basis_op; + ahead.indexing_op = f.indexing_op; + ahead.run = std::async(std::launch::async, + [&, merged = model_merged, sg = experiment_.GetSpaceGroupOrP1(), cell = *result.consensus_cell, + friedel = experiment_.GetScalingSettings().GetMergeFriedel(), + rfree_fraction = experiment_.GetScalingSettings().GetRfreeFraction(), + wavelength = experiment_.GetWavelength_A(), + gate = ahead.write.get_future().share(), beside = f.gpu_bytes]() mutable { + if (const int32_t cards = get_gpu_count(); cards > 0) + pin_gpu(1 % cards); + if (!(ahead.indexing_op == gemmi::Op::identity())) { + merged = ReindexMergedIntoAsu(merged, ahead.indexing_op, sg, friedel); + cell = CellInBasis(cell, HklOperator(ahead.indexing_op)); + } + gemmi::Op to_model = ahead.change_of_basis_op; + to_model.tran = {0, 0, 0}; + const gemmi::Op cob = to_model.inverse(); + const gemmi::SpaceGroup *in = &sg; + if (const gemmi::SpaceGroup *moved = SpaceGroupInBasis(sg, cob)) { + merged = ReindexMergedIntoAsu(merged, HklOperator(cob), *moved, friedel); + cell = CellInBasis(cell, cob); + AssignRfreeFlags(merged, *moved, rfree_fraction, 500, cell, config_.nthreads); + in = moved; + } + ModelValidationSchedule schedule; + schedule.maps_gate = gate; + schedule.fit_test = false; // the verdict is the first validation's (KeepModelVerdict) + schedule.gpu_resident_bytes = merge_on_card; + schedule.gpu_beside_bytes = beside; + ahead.merged = std::move(merged); + return ValidateAgainstModel(ahead.merged, cell, *config_.model, ahead.log, + in, /*probe_indexing_ambiguity=*/false, + static_cast(config_.nthreads), wavelength, + report_shell_d_min, schedule); + }); + }; + // Whether the validation started on the forecast is the one kept: the decision is the one forecast. + // Asked by the first validation before it makes its maps, which are then left to that one, and + // again below with the verdict in hand - the same answer, as nothing it reads changes between. + const auto forecast_kept = [&](const ModelValidationResult &v) { + return ahead.run.valid() && config_.reference_data.empty() && v.model_fits + && !(v.change_of_basis_op == gemmi::Op::identity()) && !v.adopted_model_enantiomorph + && v.change_of_basis_op == ahead.change_of_basis_op && v.indexing_op == ahead.indexing_op; + }; + first_schedule.maps_superseded = forecast_kept; + Logger model_validation_log = Logger::Buffered(); + std::future model_validation_run; // last: it is joined before what it reads goes + if (model_run && !model_predicted && !config_.force_mr) { + model_merged = sm.merged; + model_validation_run = std::async(std::launch::async, [&, device = get_gpu()] { + set_gpu(device); + TimingMark mark(model_validation_log, "model validation (beside the analyses)"); + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = config_.nthreads; + const gemmi::UnitCell cell(*result.consensus_cell); + fw_opts.anisotropy_hkl = AnisotropyPriorHKL(model_merged, cell, twin_sg); + ApplyFrenchWilson(model_merged, experiment_.GetSpaceGroupOrP1(), fw_opts); + return ValidateAgainstModel(model_merged, *result.consensus_cell, *config_.model, + model_validation_log, data_sg ? &*data_sg : nullptr, + /*probe_indexing_ambiguity=*/config_.reference_data.empty(), + static_cast(config_.nthreads), experiment_.GetWavelength_A(), + report_shell_d_min, first_schedule); + }); + } // Not on the geometry pre-pass, nor on a superseded one: the analysis goes into that pass's // statistics text and its written reflections, and neither survives the run. The promotion flag // below is a different thing - it is what the SEARCH did, the second pass reads it, and it is @@ -4068,142 +4204,26 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { std::optional model_prior_file; // The coordinates that go with the maps beside the answer's files, kept once the model is done. std::optional placed_model_out; - // fast_live, cut (the broker's deadline or a cancel): the model is left out, and the files are written - // as a run without one writes them. Asked again before the search below. - const auto model_cut = [&] { - if (!config_.fast_live || !cancelled_) - return false; - Note("cut: the model is left out"); - return true; - }; - const bool model_run = result.consensus_cell && write_output && 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. std::optional &model_twin = model_twin_merge ? model_twin_file : twin_file; std::optional &model_twin_entry = model_twin_merge ? model_twin_report : result.twin_subgroup_merge; const bool model_twin_only = model_twin_merge || twin_for_mr_only; - // The merge's own shell bounds, so CC(model, data) is reported on the shells CC1/2 was. - std::vector report_shell_d_min; - for (const auto &sh : sm.statistics.shells) - report_shell_d_min.push_back(sh.d_min); // Placed by the rigid body: the model fits against its null in the answer's merge or in the twin // hypothesis' one. Where it is not - or is a prediction, or --mr asked for the search - it is searched // for by molecular replacement below. bool placed_by_rigid_body = false; bool twin_chosen_by_model = false; - // Nothing indexes after the merge, so the image pass's indexers - an FFT indexer keeps its transform - // buffers and histograms on the card, 0.8 GB for the two of a rotation run - go before the model is - // validated or searched for, once the leftover-lattice census, the last to use them, has ended. - if (model_run) { - if (census_future.valid()) { - TimingMark census_mark(logger, "leftover-lattice census wait (model)"); - census_future.wait(); - } - p.indexer_pool.reset(); - // The model's GPU budget is 80% of the card less what the merge still holds, so from here the merge - // holds the card: a speculative probe not started yet waits for it to end (GPUMerge), and one let - // beside it is waited for. That is a probe pass of seconds: this long means it is stuck, and the - // model goes on as it would have without it. - if (p.gpu_merge) - p.gpu_merge->Holds(get_gpu_total_memory(0)); - TimingMark beside_mark(logger, "GPU work beside wait (model)"); - wait_for_gpu_work_beside(std::chrono::minutes(10)); - } if (model_run && !model_predicted && !config_.force_mr) { phase("Validating against model"); Step("Model validation"); TimingMark mark(logger, "model validation"); - const auto data_sg = experiment_.GetGemmiSpaceGroup(); - // With a reference MTZ the merohedral indexing was already resolved against it (rotation - // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. - // The validation in the model's setting further down depends on this one only through the - // setting and the indexing it settles, and both are forecast well before this one has finished - // (ModelFrameForecast). So it is started on that forecast, beside this one - on the next card - // where there is one - and kept only where the decision is the one forecast; its log is held, - // and its files wait, until then. On any other decision it is dropped, unwritten. The relabelling - // it runs on is the one AdoptModelFrame and relabel_output below make, made here on a copy. - // A model asserting the other enantiomorph is not forecast: whether the label is taken is - // only known at the end, from the anomalous map. - struct Ahead { - gemmi::Op change_of_basis_op, indexing_op; - // The reflections it validates against, relabelled as the files would be: where it is kept, - // they are the written ones, and the relabelling below is not made a second time. - std::vector merged; - Logger log = Logger::Buffered(); - std::future run; - std::promise write; // destroyed before `run`, so a run left waiting on it is released - } ahead; - // The merge is done with its engine - nothing below reads it but the finalist ledger's P1 merge, - // which runs on it beside the validation - so the device memory it holds goes back before the - // validation decides, from the card's total, what goes on its GPU. Where the ledger's merge keeps - // it, what it holds is counted against the card by every validation here, the same for each, so - // that one started on the forecast decides as the one it stands in for would. - if (!ledger_ahead.valid()) - rsm.reset(); - const size_t merge_on_card = rsm && rsm->OnGPU() ? rsm->DeviceBytes() : 0; - ModelValidationSchedule first_schedule; - first_schedule.gpu_resident_bytes = merge_on_card; - first_schedule.on_forecast = [&](const ModelFrameForecast &f) { - if (ahead.run.valid() || !config_.reference_data.empty() || !data_sg.has_value() - || !f.may_run_beside || f.enantiomorph_candidate - || f.change_of_basis_op == gemmi::Op::identity()) - return; - ahead.change_of_basis_op = f.change_of_basis_op; - ahead.indexing_op = f.indexing_op; - ahead.run = std::async(std::launch::async, - [&, merged = sm.merged, sg = experiment_.GetSpaceGroupOrP1(), cell = *result.consensus_cell, - friedel = experiment_.GetScalingSettings().GetMergeFriedel(), - rfree_fraction = experiment_.GetScalingSettings().GetRfreeFraction(), - wavelength = experiment_.GetWavelength_A(), - gate = ahead.write.get_future().share(), beside = f.gpu_bytes]() mutable { - if (const int32_t cards = get_gpu_count(); cards > 0) - pin_gpu(1 % cards); - if (!(ahead.indexing_op == gemmi::Op::identity())) { - merged = ReindexMergedIntoAsu(merged, ahead.indexing_op, sg, friedel); - cell = CellInBasis(cell, HklOperator(ahead.indexing_op)); - } - gemmi::Op to_model = ahead.change_of_basis_op; - to_model.tran = {0, 0, 0}; - const gemmi::Op cob = to_model.inverse(); - const gemmi::SpaceGroup *in = &sg; - if (const gemmi::SpaceGroup *moved = SpaceGroupInBasis(sg, cob)) { - merged = ReindexMergedIntoAsu(merged, HklOperator(cob), *moved, friedel); - cell = CellInBasis(cell, cob); - AssignRfreeFlags(merged, *moved, rfree_fraction, 500, cell, config_.nthreads); - in = moved; - } - ModelValidationSchedule schedule; - schedule.maps_gate = gate; - schedule.fit_test = false; // the verdict is the first validation's (KeepModelVerdict) - schedule.gpu_resident_bytes = merge_on_card; - schedule.gpu_beside_bytes = beside; - ahead.merged = std::move(merged); - return ValidateAgainstModel(ahead.merged, cell, *config_.model, ahead.log, - in, /*probe_indexing_ambiguity=*/false, - static_cast(config_.nthreads), wavelength, - report_shell_d_min, schedule); - }); - }; - // Whether the validation started on the forecast is the one kept: the decision is the one forecast. - // Asked by the first validation before it makes its maps, which are then left to that one, and - // again below with the verdict in hand - the same answer, as nothing it reads changes between. - const auto forecast_kept = [&](const ModelValidationResult &v) { - return ahead.run.valid() && config_.reference_data.empty() && v.model_fits - && !(v.change_of_basis_op == gemmi::Op::identity()) && !v.adopted_model_enantiomorph - && v.change_of_basis_op == ahead.change_of_basis_op && v.indexing_op == ahead.indexing_op; - }; - first_schedule.maps_superseded = forecast_kept; // The ledger's P1 merge (started above) is taken before anything below acts on what the // 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, - logger, data_sg ? &*data_sg : nullptr, - /*probe_indexing_ambiguity=*/config_.reference_data.empty(), - static_cast(config_.nthreads), experiment_.GetWavelength_A(), - report_shell_d_min, first_schedule); + ModelValidationResult validation = model_validation_run.get(); + model_validation_log.ReplayInto(logger); take_maps(result.outputs->maps, validation); take_ledger_p1(); // A model that was asked for and could not be used has to say so where anyone will see diff --git a/tests/AnisotropyAnalysisTest.cpp b/tests/AnisotropyAnalysisTest.cpp index 400ea81ee..74d5f9ba1 100644 --- a/tests/AnisotropyAnalysisTest.cpp +++ b/tests/AnisotropyAnalysisTest.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include "../image_analysis/scale_merge/AnisotropyAnalysis.h" @@ -146,6 +147,30 @@ TEST_CASE("Anisotropy tensor is recovered from a synthetic merge", "[anisotropy] CHECK_FALSE(result.refusal.empty()); } +// The French-Wilson prior from the merge alone is the analysis' own, to the bit: in a class with a free +// direction, in a cubic one, and where too few reflections leave no tensor fitted. +TEST_CASE("Anisotropy prior from the merge alone is the analysis' tensor", "[anisotropy]") { + const auto same = [](const std::vector &merged, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *sg) { + const gemmi::Mat33 a = AnisotropyTensorHKL(AnalyzeAnisotropy(merged, {}, cell, sg), cell); + const gemmi::Mat33 b = AnisotropyPriorHKL(merged, cell, sg); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + CHECK(std::memcmp(&a[i][j], &b[i][j], sizeof(double)) == 0); + }; + const gemmi::UnitCell tetragonal = Cell(79, 79, 38, 90, 90, 90); + const gemmi::SpaceGroup *p43212 = gemmi::find_spacegroup_by_name("P 43 21 2"); + const auto merged = SyntheticMerge(tetragonal, p43212, 2.0, 20.0, {-5.0, -5.0, 10.0, 0.0, 0.0, 0.0}); + same(merged, tetragonal, p43212); + same({merged.begin(), merged.begin() + 100}, tetragonal, p43212); + const gemmi::UnitCell cubic = Cell(78, 78, 78, 90, 90, 90); + const gemmi::SpaceGroup *i23 = gemmi::find_spacegroup_by_name("I 2 3"); + same(SyntheticMerge(cubic, i23, 2.5, 20.0, {8.0, -4.0, -4.0, 0.0, 0.0, 0.0}), cubic, i23); + const gemmi::UnitCell monoclinic = Cell(51, 62, 73, 90, 105, 90); + const gemmi::SpaceGroup *c2 = gemmi::find_spacegroup_by_name("C 1 2 1"); + same(SyntheticMerge(monoclinic, c2, 2.0, 20.0, {-6.0, 2.0, 4.0, 0.0, 1.5, 0.0}), monoclinic, c2); +} + // An isotropic merge must not produce a tensor, whatever the Laue class allows. TEST_CASE("Anisotropy of an isotropic synthetic merge is small", "[anisotropy]") { const gemmi::UnitCell cell = Cell(79, 79, 38, 90, 90, 90);