diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 3484e35e9..37d0b3e0d 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -8,6 +8,8 @@ ADD_LIBRARY(Rugnux STATIC Rugnux.h RugnuxCommandLine.cpp RugnuxCommandLine.h + ModelFFT.cpp + ModelFFT.h ModelScaling.cpp ModelScaling.h ModelValidation.cpp @@ -27,7 +29,7 @@ ADD_LIBRARY(Rugnux STATIC WriteModel.h ) -TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochWriter gemmi) +TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochWriter gemmi fftw3f) # rugnux is the single offline analysis CLI; --mode picks what it does: mx (the full pipeline, the # default), azint, scale (re-scale/merge stored reflections) or calibration (powder-ring geometry). diff --git a/rugnux/ModelFFT.cpp b/rugnux/ModelFFT.cpp new file mode 100644 index 000000000..73daab1fc --- /dev/null +++ b/rugnux/ModelFFT.cpp @@ -0,0 +1,73 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ModelFFT.h" + +#include +#include +#include +#include +#include + +#include + +#include "gemmi/fail.hpp" + +namespace { + +// The r2c plan for an (nu, nv, nw) map stored u fastest, halving w - the layout gemmi's own transform +// reads and writes. FFTW's planner is not thread-safe, so plans are made under a lock; executing one +// on new arrays (fftwf_execute_dft_r2c) is, which is how every caller uses it. +fftwf_plan PlanFor(int nu, int nv, int nw) { + static std::mutex m; + static std::map, fftwf_plan> plans; + std::lock_guard lock(m); + const std::array key{nu, nv, nw}; + if (const auto it = plans.find(key); it != plans.end()) + return it->second; + float *in = fftwf_alloc_real(static_cast(nu) * nv * nw); + fftwf_complex *out = fftwf_alloc_complex(static_cast(nu) * nv * (nw / 2 + 1)); + // The last dimension is the one r2c halves; the strides put u fastest on both sides. + fftwf_iodim dims[3] = {{nu, 1, 1}, {nv, nu, nu}, {nw, nu * nv, nu * nv}}; + fftwf_plan plan = (in != nullptr && out != nullptr) + ? fftwf_plan_guru_dft_r2c(3, dims, 0, nullptr, in, out, FFTW_ESTIMATE) + : nullptr; + fftwf_free(in); + fftwf_free(out); + if (plan == nullptr) + gemmi::fail("MapToFPhi(): no FFTW plan for the grid"); + plans.emplace(key, plan); + return plan; +} + +} // namespace + +gemmi::FPhiGrid MapToFPhi(const gemmi::Grid &map) { + if (map.axis_order == gemmi::AxisOrder::ZYX) + gemmi::fail("MapToFPhi(): ZYX order is not supported"); + gemmi::FPhiGrid hkl; + hkl.unit_cell = map.unit_cell; + hkl.spacegroup = map.spacegroup; + hkl.axis_order = map.axis_order; + hkl.half_l = true; + hkl.set_size_without_checking(map.nu, map.nv, map.nw / 2 + 1); + const float norm = static_cast(map.unit_cell.volume / map.point_count()); + + const fftwf_plan plan = PlanFor(map.nu, map.nv, map.nw); + // Buffers from fftwf_malloc: a plan may only be executed on arrays aligned as the ones it was made + // on were, which the vectors' own allocations do not promise. + float *in = fftwf_alloc_real(map.data.size()); + fftwf_complex *out = fftwf_alloc_complex(hkl.data.size()); + if (in == nullptr || out == nullptr) { + fftwf_free(in); + fftwf_free(out); + gemmi::fail("MapToFPhi(): cannot allocate the FFT buffers"); + } + std::copy(map.data.begin(), map.data.end(), in); + fftwf_execute_dft_r2c(plan, in, out); + for (size_t i = 0; i < hkl.data.size(); ++i) + hkl.data[i] = std::complex(out[i][0], out[i][1]) * norm; + fftwf_free(in); + fftwf_free(out); + return hkl; +} diff --git a/rugnux/ModelFFT.h b/rugnux/ModelFFT.h new file mode 100644 index 000000000..b98778bd2 --- /dev/null +++ b/rugnux/ModelFFT.h @@ -0,0 +1,18 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "gemmi/grid.hpp" +#include "gemmi/recgrid.hpp" + +// The structure factors of a real-space map on the model path (Fcalc from the model density, Fmask +// from the bulk-solvent mask): what gemmi::transform_map_to_f_phi(map, true) returns - the forward +// transform scaled by V/N, kept for l >= 0 on an (nu, nv, nw/2+1) grid in the map's own XYZ order, +// ready for prepare_asu_data() - computed with FFTW rather than gemmi's single-threaded pocketfft. +// The same transform, not the same bits: the two libraries order their arithmetic differently. +// +// Plans are FFTW_ESTIMATE only - MEASURE would time candidate algorithms and could pick a different +// one from run to run - and made once per grid size and kept, so the rigid body's hundred-odd +// evaluations on one grid plan once. Safe to call from several threads at once. +gemmi::FPhiGrid MapToFPhi(const gemmi::Grid &map); diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 4deecd583..d85b0430b 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -19,7 +19,7 @@ #include // MaybeGzipped #include // IT92 x-ray form factors #include // DensityCalculator -#include // transform_map_to_f_phi, get_f_phi_on_grid, transform_f_phi_grid_to_map +#include // get_f_phi_on_grid, transform_f_phi_grid_to_map (the maps) #include // SolventMasker #include // Scaling (bulk solvent + anisotropic B) #include // Ccp4 map I/O @@ -30,6 +30,7 @@ #include "../common/Logger.h" #include "../common/ParallelFor.h" // ParallelFor #include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections +#include "ModelFFT.h" #include "RigidBodyRefine.h" #include "SigmaA.h" @@ -283,7 +284,7 @@ double frame_probe_r(const gemmi::Structure &st, const gemmi::SpaceGroup *sg, dc.set_refmac_compatible_blur(st.models[0]); dc.put_model_density_on_grid(st.models[0]); gemmi::AsuData> fcalc = - gemmi::transform_map_to_f_phi(dc.grid, true) + MapToFPhi(dc.grid) .prepare_asu_data(probe_d_min, dc.blur, false, false, false); gemmi::GroupOps gops = sg->operations(); @@ -476,7 +477,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & dc.set_grid_cell_and_spacegroup(ms.st); dc.set_refmac_compatible_blur(ms.st.models[0]); dc.put_model_density_on_grid(ms.st.models[0]); - ms.fcalc = gemmi::transform_map_to_f_phi(dc.grid, true) + ms.fcalc = MapToFPhi(dc.grid) .prepare_asu_data(dc.d_min, dc.blur, false, false, false); // Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost. @@ -486,7 +487,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & mask_grid.spacegroup = dc.grid.spacegroup; mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); masker.put_mask_on_grid(mask_grid, ms.st.models[0]); - ms.fmask = gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); + ms.fmask = MapToFPhi(mask_grid).prepare_asu_data(dc.d_min, 0); }; compute_model_factors(mdl); @@ -1246,7 +1247,7 @@ std::vector ModelReferenceIntensities(const std::string &model dc.set_refmac_compatible_blur(st.models[0]); dc.put_model_density_on_grid(st.models[0]); gemmi::AsuData> fcalc = - gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false); + MapToFPhi(dc.grid).prepare_asu_data(dc.d_min, dc.blur, false, false, false); // Flat bulk solvent at the standard constants. Nothing here is fitted - there are no observations // yet - but without it the few lowest-resolution reflections are the largest and the most wrong, @@ -1260,7 +1261,7 @@ std::vector ModelReferenceIntensities(const std::string &model mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); masker.put_mask_on_grid(mask_grid, st.models[0]); gemmi::AsuData> fmask = - gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); + MapToFPhi(mask_grid).prepare_asu_data(dc.d_min, 0); std::unordered_map> mask_by_hkl; mask_by_hkl.reserve(fmask.v.size()); diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index e9a34551f..91d25843c 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -13,11 +13,11 @@ #include #include "gemmi/dencalc.hpp" // DensityCalculator -#include "gemmi/fourier.hpp" // transform_map_to_f_phi #include "gemmi/it92.hpp" // IT92 x-ray form factors #include "gemmi/scaling.hpp" // Scaling (bulk solvent + anisotropic B) #include "gemmi/solmask.hpp" // SolventMasker +#include "ModelFFT.h" // MapToFPhi #include "ModelScaling.h" // FitModelScale #include "../common/JFJochMath.h" // PI #include "../common/Logger.h" @@ -104,7 +104,7 @@ public: dc.set_refmac_compatible_blur(model_); dc.put_model_density_on_grid(model_); gemmi::AsuData> fcalc = - gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false); + MapToFPhi(dc.grid).prepare_asu_data(dc.d_min, dc.blur, false, false, false); gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac); gemmi::Grid mask_grid; @@ -113,7 +113,7 @@ public: mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); masker.put_mask_on_grid(mask_grid, model_); gemmi::AsuData> fmask = - gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); + MapToFPhi(mask_grid).prepare_asu_data(dc.d_min, 0); if (fmask.size() != fcalc.size()) return false; diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index b7e4ec822..943a559cc 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -11,8 +11,10 @@ #include #include +#include #include "../common/Logger.h" +#include "../rugnux/ModelFFT.h" #include "../rugnux/ModelValidation.h" #include "../rugnux/RigidBodyRefine.h" #include "../rugnux/SigmaA.h" @@ -587,3 +589,36 @@ TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") std::filesystem::remove(prefix + suffix); std::filesystem::remove(path); } + +// The model path's structure factors come from FFTW; they must be gemmi's own transform to float +// precision, in the same layout, so prepare_asu_data() reads the same reflections from either. +TEST_CASE("ModelValidation_MapToFPhiMatchesGemmi", "[ModelValidation]") { + gemmi::Grid map; + map.unit_cell.set(40.0, 50.0, 60.0, 90.0, 95.0, 90.0); + map.spacegroup = gemmi::find_spacegroup_by_name("P 1"); + map.set_size(20, 24, 30); + std::mt19937 rng(7); + std::uniform_real_distribution dist(-1.0f, 1.0f); + for (auto &x : map.data) + x = dist(rng); + + gemmi::FPhiGrid ref = gemmi::transform_map_to_f_phi(map, true); + gemmi::FPhiGrid ours = MapToFPhi(map); + REQUIRE(ours.nu == ref.nu); + REQUIRE(ours.nv == ref.nv); + REQUIRE(ours.nw == ref.nw); + REQUIRE(ours.half_l == ref.half_l); + REQUIRE(ours.data.size() == ref.data.size()); + double largest = 0, worst = 0; + for (size_t i = 0; i < ref.data.size(); i++) { + largest = std::max(largest, static_cast(std::abs(ref.data[i]))); + worst = std::max(worst, static_cast(std::abs(ours.data[i] - ref.data[i]))); + } + CHECK(worst <= 1e-5 * largest); + + const auto a = ref.prepare_asu_data(4.0, 0, false, false, false); + const auto b = ours.prepare_asu_data(4.0, 0, false, false, false); + REQUIRE(a.v.size() == b.v.size()); + for (size_t i = 0; i < a.v.size(); i++) + CHECK(a.v[i].hkl == b.v[i].hkl); +}