From 7e5e4c54e42d2bb0ba2769b22bcc462c40e54cae Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 19 Oct 2026 13:43:28 +0200 Subject: [PATCH] rugnux: the model's first validation starts beside the analyses instead of after them Validation reads of the final merge only its French-Wilson amplitudes, which the analyses made last, from the anisotropy tensor. That tensor is a fit to the merged intensities alone (the unmerged observations only give its error bar), so AnisotropyPriorHKL computes it from the merge without the rest of the analysis, bit for bit what AnisotropyTensorHKL of the full analysis gives (test). The first ValidateAgainstModel now runs from the start of the analyses on a copy of the merge given those same amplitudes, on the caller's card (get_gpu), with its log held and replayed where it was printed before; nothing it decides is acted on earlier. The model's preparation (census wait, indexer pool, GPU hold, merge engine release) moves with it. Measured, --simulate-live 10 --fast-live, 4axt model, RTX 5080, analyses + validation steps: strong lysozyme 0.91 -> 0.63 s, reference lysozyme 1.02 -> 0.66 s (validation beside the analyses 0.60-0.62 s vs 0.39-0.45 s alone, CPU shared; analyses 0.53 -> 0.62 s). The analyses were 80% waiting for the anisotropy (0.43 of 0.53 s), so starting on the full analysis' tensor would have gained ~0.1 s. p.hkl, p_maps.mtz identical, p.mtz data identical, all validation numbers identical. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4 --- common/CUDAWrapper.cpp | 4 + common/CUDAWrapper.cu | 7 + common/CUDAWrapper.h | 4 + .../scale_merge/AnisotropyAnalysis.cpp | 27 ++ .../scale_merge/AnisotropyAnalysis.h | 6 + rugnux/RugnuxScaleMerge.cpp | 256 ++++++++++-------- tests/AnisotropyAnalysisTest.cpp | 25 ++ 7 files changed, 211 insertions(+), 118 deletions(-) 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);