// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include "MolecularReplacement.h"

#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <future>
#include <numeric>
#include <random>
#include <unordered_map>

#include <gemmi/mmread_gz.hpp>   // read_structure_gz
#include <gemmi/modify.hpp>      // remove_hydrogens
#include <gemmi/dencalc.hpp>     // DensityCalculator
#include <gemmi/it92.hpp>        // IT92 x-ray form factors
#include <gemmi/scaling.hpp>     // Scaling (bulk solvent + anisotropic B)

#include "../../common/JFJochMath.h"   // PI
#include "../../common/Logger.h"
#include "../../common/ParallelFor.h"
#include "ModelFFT.h"
#include "ModelGrid.h"
#include "ModelScaling.h"
#include "RigidBodyRefine.h"
#ifdef JFJOCH_USE_CUDA
#include "RigidBodyGPU.h"
#endif

namespace {

using Table = gemmi::IT92<float>;

// The search's resolution range: the low-resolution end is where the bulk solvent, which the search
// does not model, stops dominating; the high-resolution end is where a predicted model still agrees
// with the crystal's structure in its details.
constexpr double SEARCH_D_MIN = 4.0;
constexpr double SEARCH_D_MAX = 12.0;
constexpr int SEARCH_BINS = 10;
constexpr double ROTATION_STEP_DEG = 6.0;
// Orientations the translation function is run for, per merge, and placements per space group taken
// through the rigid body.
constexpr int TRANSLATION_ROTATIONS = 50;
constexpr int RESCORED_PER_GROUP = 3;
// The resolution the placements are refined and scored to: whether a model sits in the right place
// is a low-resolution question, the rigid body's ladder ends at 3.5 A anyway, and every R-factor of
// the comparison is made to the same limit.
constexpr double SCORE_D_MIN = 3.0;

// The decision. A translation-function peak this many standard deviations above the function's mean
// is a solution, where the model packs and one copy fits the cell; the McCoy et al. (2007) J. Appl.
// Cryst. 40, 658-674 convention for a single-copy search, taken over to this Z-score's own null on
// the Crowther & Blow function (measured: the highest peak of a wrong orientation reads 4-6).
constexpr double TF_Z_PLACED = 8.0;
constexpr double CLASH_MAX = 0.05;
// Matthews (1968) J. Mol. Biol. 33, 491-497: protein crystals have Vm of about 1.7 to 3.5 A^3/Da.
// Below VM_MIN one copy does not fit the cell; at or above VM_TWO_COPIES two copies would.
constexpr double VM_MIN = 1.5;
constexpr double VM_TWO_COPIES = 4.0;
constexpr double DALTON_PER_RESIDUE = 110.0;
// How far below every other candidate of its merge the chosen one's R-free must be.
constexpr double R_FREE_MARGIN = 0.03;
// The null (as in model validation): the same model refitted from random orientations about its
// centroid, at a coarse limit, and how far above that the real placement's R-work has to sit.
constexpr int NULL_REPLICATES = 9;
constexpr unsigned NULL_SEED = 20261011;
constexpr double NULL_D_MIN = 3.5;
constexpr double NULL_SIGMA_MIN = 3.0;

double seconds_since(std::chrono::steady_clock::time_point t0) {
    return std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
}

// A rotation drawn uniformly from SO(3), through a uniform random unit quaternion.
// Following Shoemake (1992) Graphics Gems III, 124-132
gemmi::Mat33 random_rotation(std::mt19937 &rng) {
    std::uniform_real_distribution<double> u(0.0, 1.0);
    const double u1 = u(rng), t2 = 2 * PI * u(rng), t3 = 2 * PI * u(rng);
    const double r1 = std::sqrt(1 - u1), r2 = std::sqrt(u1);
    const double x = r1 * std::sin(t2), y = r1 * std::cos(t2), z = r2 * std::sin(t3), w = r2 * std::cos(t3);
    return {1 - 2 * (y * y + z * z), 2 * (x * y - z * w),     2 * (x * z + y * w),
            2 * (x * y + z * w),     1 - 2 * (x * x + z * z), 2 * (y * z - x * w),
            2 * (x * z - y * w),     2 * (y * z + x * w),     1 - 2 * (x * x + y * y)};
}

double trace(const gemmi::Mat33 &m) {
    return m.a[0][0] + m.a[1][1] + m.a[2][2];
}

gemmi::Vec3 miller_vec(const std::array<int, 3> &h) {
    return {static_cast<double>(h[0]), static_cast<double>(h[1]), static_cast<double>(h[2])};
}

// h R for a symmetry operator's rotation: the index the operator takes h to.
gemmi::Miller rotated_index(const gemmi::Op &op, const gemmi::Miller &h) {
    const auto k = op.apply_to_hkl_without_division(h);
    return {{k[0] / gemmi::Op::DEN, k[1] / gemmi::Op::DEN, k[2] / gemmi::Op::DEN}};
}

long hkl_key(const gemmi::Miller &h) {
    return (static_cast<long>(h[0]) + 512) * 1048576L + (static_cast<long>(h[1]) + 512) * 1024L + (h[2] + 512);
}

int good_even_size(int n) {
    while (n % 2 != 0 || !gemmi::has_small_factorization(n))
        ++n;
    return n;
}

// --- placements scored against the data: R-work / R-free after the bulk solvent and the scale ---

struct PlacementScore {
    double r_work = 1, r_free = 1;
};

gemmi::AsuData<gemmi::ValueSigma<float>> working_amplitudes(const std::vector<MergedReflection> &merged,
                                                            const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
                                                            double d_min) {
    const gemmi::GroupOps gops = sg.operations();
    const gemmi::ReciprocalAsu asu(&sg);
    gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
    fobs.unit_cell_ = cell;
    fobs.spacegroup_ = &sg;
    for (const MergedReflection &r : merged) {
        if (std::isnan(r.F) || r.rfree_flag || r.d < d_min)
            continue;
        gemmi::Miller h{{r.h, r.k, r.l}};
        if (!asu.is_in(h))
            h = asu.to_asu(h, gops).first;
        fobs.v.push_back({h, {r.F, 1.0f}});
    }
    fobs.ensure_asu();
    fobs.ensure_sorted();
    return fobs;
}

// The model's Fcalc and bulk-solvent mask to d_min, scaled to the working amplitudes as model
// validation scales them (ModelValidation.cpp, fit_model), and the R-factors over the reflections
// to d_min.
PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
                               const std::vector<MergedReflection> &merged, double d_min, size_t nthreads) {
    gemmi::DensityCalculator<Table, float> dc;
    dc.d_min = d_min;
    dc.rate = 1.5;
    dc.set_grid_cell_and_spacegroup(st);
    dc.set_refmac_compatible_blur(st.models[0]);
    gemmi::Grid<float> mask;
    mask.unit_cell = dc.grid.unit_cell;
    mask.spacegroup = dc.grid.spacegroup;
    mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
    const std::vector<size_t> leaders = OrbitLeaders(mask, nthreads);
    PutMaskOnGrid(mask, st.models[0], leaders, nthreads);
    gemmi::AsuData<std::complex<float>> fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0);
    PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
    gemmi::AsuData<std::complex<float>> fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);

    gemmi::AsuData<gemmi::ValueSigma<float>> fobs = working_amplitudes(merged, st.cell, sg, d_min);
    gemmi::Scaling<float> scaling(st.cell, &sg);
    scaling.use_solvent = true;
    scaling.prepare_points(fmodel, fobs, &fmask);
    FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads);
    scaling.scale_data(fmodel, &fmask);

    const gemmi::GroupOps gops = sg.operations();
    const gemmi::ReciprocalAsu asu(&sg);
    std::unordered_map<long, std::pair<double, bool>> obs;
    for (const MergedReflection &r : merged) {
        if (std::isnan(r.F) || r.d < d_min)
            continue;
        gemmi::Miller h{{r.h, r.k, r.l}};
        if (!asu.is_in(h))
            h = asu.to_asu(h, gops).first;
        obs[hkl_key(h)] = {r.F, r.rfree_flag};
    }
    double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
    for (const auto &hv : fmodel.v) {
        const auto it = obs.find(hkl_key(hv.hkl));
        if (it == obs.end())
            continue;
        const double fo = it->second.first, fc = std::abs(hv.value);
        if (it->second.second) {
            num_f += std::fabs(fo - fc);
            den_f += fo;
        } else {
            num_w += std::fabs(fo - fc);
            den_w += fo;
        }
    }
    PlacementScore out;
    out.r_work = den_w > 0 ? num_w / den_w : 1;
    out.r_free = den_f > 0 ? num_f / den_f : 1;
    return out;
}

// The rigid-body engines for one model in one group, where there is a GPU.
struct RigidBodyEngines {
#ifdef JFJOCH_USE_CUDA
    std::unique_ptr<RigidBodyGPUPool> pool;
#endif
    RigidBodyGPUPool *get() {
#ifdef JFJOCH_USE_CUDA
        return pool.get();
#else
        return nullptr;
#endif
    }
};

RigidBodyEngines make_engines(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
                              const std::vector<MergedReflection> &merged, double d_min, size_t engines,
                              Logger &logger) {
    RigidBodyEngines out;
#ifdef JFJOCH_USE_CUDA
    const double finest = RigidBodyLadder(d_min).back();
    size_t observations = 0;
    for (const MergedReflection &r : merged)
        if (!std::isnan(r.F) && r.d >= finest)
            ++observations;
    out.pool = RigidBodyGPUPool::Create(st.models[0], st.cell, sg, d_min, observations, engines, logger);
#else
    (void) st; (void) sg; (void) merged; (void) d_min; (void) engines; (void) logger;
#endif
    return out;
}

// The placement through the rigid body, and its R-factors where the rigid body left it. A search
// placement is a few degrees and a grid step from its answer, well inside the rigid body's reach, so
// the step is taken without the R-free check model validation makes before committing one.
PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &sg,
                                const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
                                size_t nthreads) {
    const std::vector<gemmi::Position> start = ModelPositions(st.models[0]);
    Logger quiet = Logger::Buffered();
    const RigidBodyRefineResult rb = RefineRigidBody(st.models[0], st.cell, sg,
                                                     working_amplitudes(merged, st.cell, sg, d_min), d_min, quiet,
                                                     nthreads, pool);
    if (!rb.converged)
        SetModelPositions(st.models[0], start);
    return score_placement(st, sg, merged, d_min, nthreads);
}

// The placed model's R-work at NULL_D_MIN against the same model refitted from random orientations
// about its centroid, in standard deviations of the latter: the null model validation builds for a
// model that claims something (ModelValidation.cpp). An orientation the crystal cannot tell from the
// placement's - within the rigid body's reach of one of its symmetry classes - is drawn again.
double null_sigma(const gemmi::Structure &placed, const gemmi::Mat33 &orientation, const gemmi::SpaceGroup &sg,
                  const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
                  size_t nthreads) {
    const double null_d_min = std::max(d_min, NULL_D_MIN);
    std::vector<MergedReflection> coarse;
    for (const MergedReflection &r : merged)
        if (r.d >= null_d_min)
            coarse.push_back(r);
    const double real = score_placement(placed, sg, coarse, null_d_min, nthreads).r_work;

    const std::vector<gemmi::Mat33> pg = CartesianPointGroup(sg, placed.cell);
    const double reach = RigidBodyReachDeg(placed.models[0], null_d_min);
    std::vector<gemmi::Mat33> draws;
    std::mt19937 rng(NULL_SEED);
    while (static_cast<int>(draws.size()) < NULL_REPLICATES) {
        const gemmi::Mat33 r = random_rotation(rng);
        if (RotationDistanceDeg(r.multiply(orientation), orientation, pg) >= reach)
            draws.push_back(r);
    }
    std::vector<double> r_work(NULL_REPLICATES);
    std::atomic<int> next{0};
    std::vector<std::future<void>> workers;
    for (size_t t = 0; t < std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES); t++)
        workers.push_back(std::async(std::launch::async, [&] {
            for (int i = next++; i < NULL_REPLICATES; i = next++) {
                gemmi::Structure rep = placed;
                std::vector<gemmi::Position> pos = ModelPositions(rep.models[0]);
                gemmi::Vec3 c;
                for (const gemmi::Position &p : pos)
                    c += p;
                c *= 1.0 / static_cast<double>(pos.size());
                for (gemmi::Position &p : pos)
                    p = gemmi::Position(draws[i].multiply(gemmi::Vec3(p) - c) + c);
                SetModelPositions(rep.models[0], pos);
                r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads).r_work;
            }
        }));
    for (auto &w : workers)
        w.get();
    const double mean = std::accumulate(r_work.begin(), r_work.end(), 0.0) / NULL_REPLICATES;
    double s2 = 0;
    for (double x : r_work)
        s2 += (x - mean) * (x - mean);
    const double sd = std::sqrt(s2 / (NULL_REPLICATES - 1));
    return sd > 0 ? (mean - real) / sd : 0.0;
}

bool same_rotations_and_centring(const gemmi::GroupOps &a, const gemmi::GroupOps &b) {
    if (a.sym_ops.size() != b.sym_ops.size() || a.cen_ops.size() != b.cen_ops.size())
        return false;
    for (const gemmi::Op &op : a.sym_ops)
        if (b.find_by_rotation(op.rot) == nullptr)
            return false;
    for (const gemmi::Op::Tran &c : a.cen_ops)
        if (std::find(b.cen_ops.begin(), b.cen_ops.end(), c) == b.cen_ops.end())
            return false;
    return true;
}

bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
    if (!a.is_enantiomorphic())
        return false;
    gemmi::GroupOps e = a.operations();
    e.change_basis_forward(a.change_of_hand_op());
    const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e);
    return other != nullptr && other->number == b.number;
}

}  // namespace

// --- the model ---

SearchModel PrepareSearchModel(const std::string &path) {
    SearchModel out;
    try {
        out.st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
    } catch (const std::exception &e) {
        out.failure_reason = std::string("cannot read the model: ") + e.what();
        return out;
    }
    if (out.st.models.empty()) {
        out.failure_reason = "the model has no atoms";
        return out;
    }
    out.st.models.resize(1);
    gemmi::remove_hydrogens(out.st);
    gemmi::Model &model = out.st.models[0];
    for (gemmi::Chain &ch : model.chains) {
        for (const gemmi::Residue &r : ch.residues)
            for (const gemmi::Atom &a : r.atoms)
                if (a.b_iso > 100.0f) {
                    out.failure_reason = "the B-factor column holds values above 100, so it is not a pLDDT: "
                                         "not a predicted model";
                    return out;
                }
        out.residues_read += static_cast<int>(ch.residues.size());
        // The residues a prediction is confident of: pLDDT 70 and above, as Oeffner et al. (2022)
        // Acta Cryst. D78, 1303-1314 trim an AlphaFold model for molecular replacement.
        std::erase_if(ch.residues, [](const gemmi::Residue &r) { return r.atoms.empty() || r.atoms[0].b_iso < 70.0f; });
        out.residues_kept += static_cast<int>(ch.residues.size());
    }
    std::erase_if(model.chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); });
    if (out.residues_kept == 0) {
        out.failure_reason = "no residue of the model has a pLDDT of 70 or more";
        return out;
    }
    // B from pLDDT through the coordinate error it predicts: rmsd = 1.5 exp(4 (0.7 - pLDDT/100)) A
    // and B = 8 pi^2 rmsd^2 / 3. Following Oeffner et al. (2022) Acta Cryst. D78, 1303-1314
    std::vector<gemmi::Position> pos;
    for (gemmi::Chain &ch : model.chains)
        for (gemmi::Residue &r : ch.residues)
            for (gemmi::Atom &a : r.atoms) {
                const double rmsd = 1.5 * std::exp(4.0 * (0.7 - a.b_iso / 100.0));
                a.b_iso = static_cast<float>(8.0 * PI * PI * rmsd * rmsd / 3.0);
                pos.push_back(a.pos);
            }
    gemmi::Vec3 c;
    for (const gemmi::Position &p : pos)
        c += p;
    c *= 1.0 / static_cast<double>(pos.size());
    for (gemmi::Position &p : pos) {
        p = gemmi::Position(gemmi::Vec3(p) - c);
        out.radius_A = std::max(out.radius_A, p.length());
    }
    SetModelPositions(model, pos);
    out.mass_da = DALTON_PER_RESIDUE * out.residues_kept;
    return out;
}

// --- rotations ---

std::vector<gemmi::Mat33> UniformRotations(double step_deg) {
    const double step = step_deg * PI / 180.0;
    const double golden = PI * (3.0 - std::sqrt(5.0));
    const int shells = std::max(1, static_cast<int>(std::lround(PI / step)));
    std::vector<gemmi::Mat33> out{gemmi::Mat33()};
    for (int i = 1; i <= shells; i++) {
        const double kappa = PI * i / shells;
        const double reach = step / (2.0 * std::sin(kappa / 2.0));   // the axis spacing this shell needs
        int n = std::max(1, static_cast<int>(std::lround(4.0 * PI / (reach * reach))));
        // At 180 deg an axis and its opposite are the same rotation: one hemisphere.
        const bool half = i == shells;
        if (half)
            n = (n + 1) / 2;
        for (int k = 0; k < n; k++) {
            const double z = half ? 1.0 - (k + 0.5) / n : 1.0 - (2.0 * k + 1.0) / n;
            const double r = std::sqrt(std::max(0.0, 1.0 - z * z));
            const double x = r * std::cos(golden * k), y = r * std::sin(golden * k);
            const double s = std::sin(kappa), c1 = 1.0 - std::cos(kappa);
            out.push_back({1 + c1 * (x * x - 1), -s * z + c1 * x * y, s * y + c1 * x * z,
                           s * z + c1 * x * y, 1 + c1 * (y * y - 1), -s * x + c1 * y * z,
                           -s * y + c1 * x * z, s * x + c1 * y * z, 1 + c1 * (z * z - 1)});
        }
    }
    return out;
}

std::vector<gemmi::Mat33> CartesianPointGroup(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell) {
    std::vector<gemmi::Mat33> out;
    for (const gemmi::Op &op : sg.operations().sym_ops) {
        gemmi::Mat33 f;
        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                f.a[i][j] = static_cast<double>(op.rot[i][j]) / gemmi::Op::DEN;
        out.push_back(cell.orth.mat.multiply(f).multiply(cell.frac.mat));
    }
    return out;
}

std::vector<gemmi::Mat33> RotationAsu(const std::vector<gemmi::Mat33> &all,
                                      const std::vector<gemmi::Mat33> &point_group) {
    std::vector<gemmi::Mat33> out;
    for (const gemmi::Mat33 &r : all) {
        const double t = trace(r);
        if (std::all_of(point_group.begin(), point_group.end(),
                        [&](const gemmi::Mat33 &s) { return t >= trace(s.multiply(r)) - 1e-9; }))
            out.push_back(r);
    }
    return out;
}

double RotationDistanceDeg(const gemmi::Mat33 &a, const gemmi::Mat33 &b,
                           const std::vector<gemmi::Mat33> &point_group) {
    double best = -1.0;
    for (const gemmi::Mat33 &s : point_group)
        best = std::max(best, trace(a.multiply(s.multiply(b).transpose())));
    return std::acos(std::clamp((best - 1.0) / 2.0, -1.0, 1.0)) * 180.0 / PI;
}

// --- the molecular transform ---

MolecularTransform::MolecularTransform(const SearchModel &model, double d_min, size_t nthreads) {
    // On a P1 box six times the model's radius: the transform of a body of radius r turns its phase
    // by pi r / box per grid step, and trilinear interpolation reads it to a few percent at pi/6.
    // Interpolating the molecular transform of one model in place of a Fourier summation per
    // orientation, following Navaza (1994) Acta Cryst. A50, 157-163
    box_ = std::max(6.0 * model.radius_A, 120.0);
    gemmi::Structure st = model.st;
    st.cell = gemmi::UnitCell(box_, box_, box_, 90, 90, 90);
    st.spacegroup_hm = "P 1";
    st.setup_cell_images();
    gemmi::DensityCalculator<Table, float> dc;
    dc.d_min = d_min;
    dc.rate = 1.5;
    dc.set_grid_cell_and_spacegroup(st);
    dc.set_refmac_compatible_blur(st.models[0]);
    PutModelDensityOnGrid(dc, st.models[0], {}, nthreads);
    grid_ = MapToFPhi(dc.grid);
    blur_ = dc.blur;
}

std::complex<double> MolecularTransform::At(const gemmi::Vec3 &s) const {
    const double u[3] = {s.x * box_, s.y * box_, s.z * box_};
    const int i0[3] = {static_cast<int>(std::floor(u[0])), static_cast<int>(std::floor(u[1])),
                       static_cast<int>(std::floor(u[2]))};
    const double f[3] = {u[0] - i0[0], u[1] - i0[1], u[2] - i0[2]};
    const auto get = [&](int i, int j, int l) {
        bool conj = false;
        if (l < 0) {   // the half-l grid holds l >= 0; F(-h) = conj F(h)
            i = -i; j = -j; l = -l;
            conj = true;
        }
        i %= grid_.nu; if (i < 0) i += grid_.nu;
        j %= grid_.nv; if (j < 0) j += grid_.nv;
        const std::complex<float> v = grid_.data[i + static_cast<size_t>(grid_.nu) * (j + static_cast<size_t>(grid_.nv) * l)];
        return conj ? std::conj(std::complex<double>(v)) : std::complex<double>(v);
    };
    std::complex<double> out = 0;
    for (int dx = 0; dx < 2; dx++)
        for (int dy = 0; dy < 2; dy++)
            for (int dz = 0; dz < 2; dz++) {
                const double w = (dx ? f[0] : 1 - f[0]) * (dy ? f[1] : 1 - f[1]) * (dz ? f[2] : 1 - f[2]);
                out += w * get(i0[0] + dx, i0[1] + dy, i0[2] + dz);
            }
    return out * std::exp(blur_ * s.length_sq() / 4.0);
}

// --- the observations ---

SearchReflections MakeSearchReflections(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
                                        const gemmi::SpaceGroup &sg, double d_min, double d_max) {
    const gemmi::GroupOps gops = sg.operations();
    const gemmi::ReciprocalAsu asu(&sg);
    struct Row { gemmi::Miller h; double i_eps; double s2; double mult; };
    std::vector<Row> rows;
    for (const MergedReflection &r : merged) {
        if (!std::isfinite(r.I) || r.d < d_min || r.d > d_max)
            continue;
        gemmi::Miller h{{r.h, r.k, r.l}};
        if (!asu.is_in(h))
            h = asu.to_asu(h, gops).first;
        if (gops.is_systematically_absent(h))
            continue;
        // The reflections the whole sphere holds for this one: its images under the point group and
        // their Friedel mates, each counted once.
        std::vector<gemmi::Miller> images;
        for (const gemmi::Op &op : gops.sym_ops)
            for (int sign : {1, -1}) {
                gemmi::Miller k = rotated_index(op, h);
                for (int &x : k)
                    x *= sign;
                if (std::find(images.begin(), images.end(), k) == images.end())
                    images.push_back(k);
            }
        rows.push_back({h, r.I / gops.epsilon_factor_without_centering(h), cell.calculate_1_d2(h),
                        static_cast<double>(images.size())});
    }
    std::stable_sort(rows.begin(), rows.end(), [](const Row &a, const Row &b) { return a.s2 < b.s2; });
    SearchReflections out;
    const size_t n = rows.size();
    out.bins = static_cast<int>(std::min<size_t>(SEARCH_BINS, std::max<size_t>(n / 100, 1)));
    std::vector<double> sum(out.bins, 0.0);
    std::vector<int> count(out.bins, 0);
    for (size_t i = 0; i < n; i++) {
        const int b = static_cast<int>(i * out.bins / n);
        out.hkl.push_back(rows[i].h);
        out.multiplicity.push_back(rows[i].mult);
        out.bin.push_back(b);
        sum[b] += rows[i].i_eps;
        count[b]++;
    }
    for (size_t i = 0; i < n; i++) {
        const double mean = sum[out.bin[i]] / count[out.bin[i]];
        out.e2.push_back(mean > 0 ? rows[i].i_eps / mean : 0.0);
    }
    return out;
}

// --- the rotation function ---

std::vector<double> RotationFunction(const MolecularTransform &ft, const SearchReflections &refl,
                                     const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
                                     const std::vector<gemmi::Mat33> &rotations, size_t nthreads) {
    const gemmi::GroupOps gops = sg.operations();
    const size_t nops = gops.sym_ops.size(), n = refl.hkl.size();
    const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
    // Copy j of the model contributes |F_model(s)|^2 at s = (h R_j) as a Cartesian reciprocal vector.
    std::vector<gemmi::Vec3> s(n * nops);
    for (size_t i = 0; i < n; i++)
        for (size_t j = 0; j < nops; j++)
            s[i * nops + j] = frac_t.multiply(miller_vec(rotated_index(gops.sym_ops[j], refl.hkl[i])));
    double w_mean = 0;
    for (double e2 : refl.e2)
        w_mean += e2;
    w_mean /= static_cast<double>(n);
    double w_ss = 0;
    for (double e2 : refl.e2)
        w_ss += (e2 - w_mean) * (e2 - w_mean);

    std::vector<double> out(rotations.size(), 0.0);
    ParallelFor(static_cast<int>(rotations.size()), nthreads, [&](int r) {
        // The model turned by R has F_R(s) = F_model(R^T s).
        const gemmi::Mat33 rt = rotations[r].transpose();
        std::vector<double> ic(n, 0.0), bin_sum(refl.bins, 0.0);
        std::vector<int> bin_n(refl.bins, 0);
        for (size_t i = 0; i < n; i++) {
            for (size_t j = 0; j < nops; j++)
                ic[i] += std::norm(ft.At(rt.multiply(s[i * nops + j])));
            bin_sum[refl.bin[i]] += ic[i];
            bin_n[refl.bin[i]]++;
        }
        double x_mean = 0;
        for (size_t i = 0; i < n; i++) {
            ic[i] /= bin_sum[refl.bin[i]] / bin_n[refl.bin[i]];
            x_mean += ic[i];
        }
        x_mean /= static_cast<double>(n);
        double sxy = 0, sxx = 0;
        for (size_t i = 0; i < n; i++) {
            sxy += (ic[i] - x_mean) * (refl.e2[i] - w_mean);
            sxx += (ic[i] - x_mean) * (ic[i] - x_mean);
        }
        out[r] = sxx > 0 && w_ss > 0 ? sxy / std::sqrt(sxx * w_ss) : 0.0;
    });
    return out;
}

// --- the translation function ---

TranslationPeak TranslationFunction(const MolecularTransform &ft, const SearchReflections &refl,
                                    const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
                                    const gemmi::Mat33 &rotation) {
    // The full-symmetry translation function as one Fourier synthesis, following Crowther & Blow
    // (1967) Acta Cryst. 23, 544-548: with the model at t, copy j adds G_j(h) exp(2 pi i (h R_j).t)
    // exp(2 pi i h.t_j) to F(h), so |F(h)|^2 summed against E^2 - 1 over the sphere is, for each pair
    // of copies j != k, a Fourier series in t with coefficient (E^2 - 1) G_j conj(G_k) exp(2 pi i
    // h.(t_j - t_k)) at the index h R_j - h R_k. Every pair's terms go on one grid, and one FFT gives
    // the function at every t. The pairs j = k are the rotation function's, constant in t.
    const gemmi::GroupOps gops = sg.operations();
    const size_t nops = gops.sym_ops.size();
    TranslationPeak out;
    if (nops < 2)
        return out;   // P1: every position is the same
    const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
    const gemmi::Mat33 rt = rotation.transpose();
    // Per copy j: k = h R_j, and g = G_j(h) exp(2 pi i h.t_j), so a pair's coefficient is g_j conj(g_k).
    struct Term { gemmi::Miller h; std::vector<gemmi::Miller> k; std::vector<std::complex<double>> g; double w; };
    std::vector<Term> terms;
    std::vector<double> bin_sum(refl.bins, 0.0);
    std::vector<int> bin_n(refl.bins, 0);
    std::vector<int> term_bin;
    for (size_t i = 0; i < refl.hkl.size(); i++) {
        if (gops.is_systematically_absent(refl.hkl[i]))
            continue;
        Term t{refl.hkl[i], {}, {}, refl.multiplicity[i] * (refl.e2[i] - 1.0)};
        double sum = 0;
        for (const gemmi::Op &op : gops.sym_ops) {
            t.k.push_back(rotated_index(op, refl.hkl[i]));
            const double ht = (t.h[0] * op.tran[0] + t.h[1] * op.tran[1] + t.h[2] * op.tran[2])
                              / static_cast<double>(gemmi::Op::DEN);
            t.g.push_back(ft.At(rt.multiply(frac_t.multiply(miller_vec(t.k.back())))) * std::polar(1.0, 2.0 * PI * ht));
            sum += std::norm(t.g.back());
        }
        bin_sum[refl.bin[i]] += sum;
        bin_n[refl.bin[i]]++;
        term_bin.push_back(refl.bin[i]);
        terms.push_back(std::move(t));
    }
    if (terms.empty())
        return out;
    int max_m[3] = {0, 0, 0};
    for (const Term &t : terms)
        for (size_t j = 0; j < nops; j++)
            for (size_t k = 0; k < nops; k++)
                for (int a = 0; a < 3; a++)
                    max_m[a] = std::max(max_m[a], std::abs(t.k[j][a] - t.k[k][a]));
    const int nu = good_even_size(2 * max_m[0] + 2), nv = good_even_size(2 * max_m[1] + 2),
              nw = good_even_size(2 * max_m[2] + 2);
    const int nwh = nw / 2 + 1;
    std::vector<std::complex<double>> c(static_cast<size_t>(nu) * nv * nwh, 0.0);
    for (size_t i = 0; i < terms.size(); i++) {
        const Term &t = terms[i];
        const double norm = bin_sum[term_bin[i]] / bin_n[term_bin[i]];
        for (size_t j = 0; j < nops; j++)
            for (size_t k = 0; k < nops; k++) {
                if (j == k)
                    continue;
                const int m[3] = {t.k[j][0] - t.k[k][0], t.k[j][1] - t.k[k][1], t.k[j][2] - t.k[k][2]};
                // Only l >= 0 is stored: the pair (k, j) carries the conjugate term at -m.
                if (m[2] < 0)
                    continue;
                const std::complex<double> v = t.w / norm * t.g[j] * std::conj(t.g[k]);
                const int u = (m[0] % nu + nu) % nu, vv = (m[1] % nv + nv) % nv;
                c[u + static_cast<size_t>(nu) * (vv + static_cast<size_t>(nv) * m[2])] += v;
            }
    }
    // MapFromFPhi sums F exp(-2 pi i m.x): with F = conj(C) that is the conjugate of the series, which
    // is real.
    gemmi::FPhiGrid<float> hkl;
    hkl.unit_cell = cell;
    hkl.spacegroup = &gemmi::get_spacegroup_p1();
    hkl.axis_order = gemmi::AxisOrder::XYZ;
    hkl.half_l = true;
    hkl.set_size_without_checking(nu, nv, nwh);
    for (size_t i = 0; i < c.size(); i++)
        hkl.data[i] = std::complex<float>(std::conj(c[i]));
    const gemmi::Grid<float> map = MapFromFPhi(hkl);
    double sum = 0, sum2 = 0;
    size_t best = 0;
    for (size_t i = 0; i < map.data.size(); i++) {
        sum += map.data[i];
        sum2 += static_cast<double>(map.data[i]) * map.data[i];
        if (map.data[i] > map.data[best])
            best = i;
    }
    const double count = static_cast<double>(map.data.size());
    const double mean = sum / count, sd = std::sqrt(std::max(0.0, sum2 / count - mean * mean));
    out.z = sd > 0 ? (map.data[best] - mean) / sd : 0.0;
    const size_t u = best % map.nu, v = (best / map.nu) % map.nv, w = best / (static_cast<size_t>(map.nu) * map.nv);
    out.t = gemmi::Fractional(static_cast<double>(u) / map.nu, static_cast<double>(v) / map.nv,
                              static_cast<double>(w) / map.nw);
    return out;
}

// --- placements ---

gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t,
                                  const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg) {
    gemmi::Structure st = model.st;
    const gemmi::Position shift = cell.orthogonalize(t);
    std::vector<gemmi::Position> pos = ModelPositions(st.models[0]);
    for (gemmi::Position &p : pos)
        p = gemmi::Position(rotation.multiply(p) + shift);
    SetModelPositions(st.models[0], pos);
    st.cell = cell;
    st.spacegroup_hm = sg.xhm();
    st.setup_cell_images();
    return st;
}

double ClashFraction(const gemmi::Structure &placed) {
    const gemmi::UnitCell &cell = placed.cell;
    const gemmi::SpaceGroup *sg = placed.find_spacegroup();
    std::vector<gemmi::Fractional> ca;
    for (const gemmi::Chain &ch : placed.models[0].chains)
        for (const gemmi::Residue &r : ch.residues)
            for (const gemmi::Atom &a : r.atoms)
                if (a.name == "CA")
                    ca.push_back(cell.fractionalize(a.pos));
    if (ca.empty() || sg == nullptr)
        return 0.0;
    gemmi::Fractional c0;
    for (const gemmi::Fractional &f : ca)
        c0 = c0 + f;
    c0 = gemmi::Fractional(c0 * (1.0 / static_cast<double>(ca.size())));
    double radius = 0;
    for (const gemmi::Fractional &f : ca)
        radius = std::max(radius, cell.orthogonalize_difference(f - c0).length());
    constexpr double CLASH_A = 3.0;
    std::vector<char> clash(ca.size(), 0);
    for (const gemmi::Op &op : sg->operations().all_ops_sorted()) {
        std::vector<gemmi::Fractional> mate;
        for (const gemmi::Fractional &f : ca) {
            const std::array<double, 3> x = op.apply_to_xyz({{f.x, f.y, f.z}});
            mate.push_back(gemmi::Fractional(x[0], x[1], x[2]));
        }
        gemmi::Fractional cm;
        for (const gemmi::Fractional &f : mate)
            cm = cm + f;
        cm = gemmi::Fractional(cm * (1.0 / static_cast<double>(mate.size())));
        const gemmi::Fractional nearest(std::round(c0.x - cm.x), std::round(c0.y - cm.y), std::round(c0.z - cm.z));
        for (int dx = -1; dx <= 1; dx++)
            for (int dy = -1; dy <= 1; dy++)
                for (int dz = -1; dz <= 1; dz++) {
                    const gemmi::Fractional shift(nearest.x + dx, nearest.y + dy, nearest.z + dz);
                    if (op == gemmi::Op::identity() && shift.x == 0 && shift.y == 0 && shift.z == 0)
                        continue;
                    if (cell.orthogonalize_difference(cm + shift - c0).length() > 2 * radius + CLASH_A)
                        continue;
                    for (size_t a = 0; a < ca.size(); a++) {
                        if (clash[a])
                            continue;
                        for (const gemmi::Fractional &m : mate)
                            if (cell.orthogonalize_difference(m + shift - ca[a]).length_sq() < CLASH_A * CLASH_A) {
                                clash[a] = 1;
                                break;
                            }
                    }
                }
    }
    return static_cast<double>(std::count(clash.begin(), clash.end(), 1)) / static_cast<double>(ca.size());
}

std::vector<gemmi::SpaceGroup> SameRotationGroups(const gemmi::SpaceGroup &sg) {
    const gemmi::GroupOps ops = sg.operations();
    std::vector<gemmi::SpaceGroup> out{sg};
    for (const gemmi::SpaceGroup &g : gemmi::spacegroup_tables::main) {
        if (g.number == sg.number && g.xhm() == sg.xhm())
            continue;
        const gemmi::GroupOps gops = g.operations();
        if (!same_rotations_and_centring(ops, gops) || gemmi::find_spacegroup_by_ops(gops) != &g)
            continue;
        out.push_back(g);
    }
    return out;
}

// --- the search and the decision ---

MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector<MRMerge> &merges,
                                                   Logger &logger, size_t nthreads) {
    MolecularReplacementResult result;
    result.model_path = model_path;
    auto t0 = std::chrono::steady_clock::now();
    const SearchModel model = PrepareSearchModel(model_path);
    result.residues_read = model.residues_read;
    result.residues_kept = model.residues_kept;
    result.mass_da = model.mass_da;
    if (!model.failure_reason.empty()) {
        result.failure_reason = model.failure_reason;
        return result;
    }
    if (merges.empty() || merges[0].merged == nullptr) {
        result.failure_reason = "no merge to search";
        return result;
    }
    double data_d_min = 1e9;
    for (const MergedReflection &r : *merges[0].merged)
        if (r.d > 0)
            data_d_min = std::min<double>(data_d_min, r.d);
    result.d_min_search = std::max(SEARCH_D_MIN, data_d_min);
    result.d_max_search = SEARCH_D_MAX;
    result.d_min_score = std::max(SCORE_D_MIN, data_d_min);
    result.t_model = seconds_since(t0);
    logger.Info("Model prior: {} of {} residues kept at pLDDT >= 70 ({:.1f} kDa, radius {:.1f} A)",
                model.residues_kept, model.residues_read, model.mass_da / 1000.0, model.radius_A);

    t0 = std::chrono::steady_clock::now();
    const MolecularTransform ft(model, result.d_min_search, nthreads);
    result.t_transform = seconds_since(t0);
    const std::vector<gemmi::Mat33> sphere = UniformRotations(ROTATION_STEP_DEG);

    // One merge: the rotation function over its point group, then for each of its groups the
    // translation function and the best placements through packing, the rigid body and R-free.
    const auto search_merge = [&](size_t m) {
        const MRMerge &merge = merges[m];
        const gemmi::UnitCell cell(merge.cell);
        const double d_min = result.d_min_score;
        // --- rotation ---
        t0 = std::chrono::steady_clock::now();
        const SearchReflections refl = MakeSearchReflections(*merge.merged, cell, merge.space_group,
                                                             result.d_min_search, result.d_max_search);
        const std::vector<gemmi::Mat33> pg = CartesianPointGroup(merge.space_group, cell);
        const std::vector<gemmi::Mat33> rotations = RotationAsu(sphere, pg);
        const std::vector<double> rf = RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads);
        double rf_mean = 0, rf_sd = 0;
        for (double x : rf)
            rf_mean += x;
        rf_mean /= static_cast<double>(rf.size());
        for (double x : rf)
            rf_sd += (x - rf_mean) * (x - rf_mean);
        rf_sd = std::sqrt(rf_sd / static_cast<double>(rf.size()));
        // The highest orientations, one per peak: a rotation within two steps of one already taken
        // (in any of its symmetry classes) belongs to that one's peak.
        std::vector<size_t> order(rf.size());
        std::iota(order.begin(), order.end(), 0);
        std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) { return rf[a] > rf[b]; });
        std::vector<size_t> peaks;
        for (size_t i : order) {
            if (std::none_of(peaks.begin(), peaks.end(), [&](size_t p) {
                    return RotationDistanceDeg(rotations[i], rotations[p], pg) < 2 * ROTATION_STEP_DEG;
                }))
                peaks.push_back(i);
            if (static_cast<int>(peaks.size()) == TRANSLATION_ROTATIONS)
                break;
        }
        result.t_rotation += seconds_since(t0);
        logger.Info("Model prior: {} in {}: {} reflections {:.1f}-{:.1f} A, {} orientations, best RF Z {:.1f}",
                    merge.name, merge.space_group.xhm(), refl.hkl.size(), result.d_max_search, result.d_min_search,
                    rotations.size(), rf_sd > 0 ? (rf[peaks.front()] - rf_mean) / rf_sd : 0.0);

        // --- translation, for every group and every peak orientation ---
        t0 = std::chrono::steady_clock::now();
        const std::vector<gemmi::SpaceGroup> groups = SameRotationGroups(merge.space_group);
        std::vector<TranslationPeak> tf_all(groups.size() * peaks.size());
        ParallelFor(static_cast<int>(tf_all.size()), nthreads, [&](int i) {
            tf_all[i] = TranslationFunction(ft, refl, cell, groups[i / peaks.size()], rotations[peaks[i % peaks.size()]]);
        });
        result.t_translation += seconds_since(t0);

        for (size_t g = 0; g < groups.size(); g++) {
            const gemmi::SpaceGroup &sg = groups[g];
            MRCandidate cand;
            cand.merge = static_cast<int>(m);
            cand.space_group = sg;
            const gemmi::GroupOps gops = sg.operations();
            cand.vm = cell.volume / static_cast<double>(gops.order()) / model.mass_da;
            cand.copies_likely = cand.vm >= VM_TWO_COPIES ? 2 : 1;
            const TranslationPeak *tf = &tf_all[g * peaks.size()];
            std::vector<size_t> by_tf(peaks.size());
            std::iota(by_tf.begin(), by_tf.end(), 0);
            std::stable_sort(by_tf.begin(), by_tf.end(), [&](size_t a, size_t b) { return tf[a].z > tf[b].z; });

            // --- the best few through packing, the rigid body and R-free ---
            t0 = std::chrono::steady_clock::now();
            RigidBodyEngines engines;
            bool have = false;
            gemmi::Mat33 orientation;
            // Where one copy cannot be the content of the cell - too big for it, or small enough for two -
            // the group is not decided on here, and only its best translation-function placement is
            // checked for packing.
            const bool one_copy = cand.vm >= VM_MIN && cand.copies_likely == 1;
            for (int k = 0; k < std::min<int>(one_copy ? RESCORED_PER_GROUP : 1, static_cast<int>(by_tf.size())); k++) {
                const size_t i = by_tf[k];
                auto placed = std::make_shared<gemmi::Structure>(
                        PlaceSearchModel(model, rotations[peaks[i]], tf[i].t, cell, sg));
                PlacementScore score{NAN, NAN};
                if (one_copy) {
                    if (k == 0)
                        engines = make_engines(*placed, sg, *merge.merged, d_min, NULL_REPLICATES + 1, logger);
                    score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads);
                }
                const double clash = ClashFraction(*placed);
                if (!have || score.r_free < cand.r_free) {   // NaN (not refined) only ever on k = 0
                    have = true;
                    cand.rf_z = rf_sd > 0 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0;
                    cand.tf_z = tf[i].z;
                    cand.clash_fraction = clash;
                    cand.r_work = score.r_work;
                    cand.r_free = score.r_free;
                    cand.placed = placed;
                    orientation = rotations[peaks[i]];
                }
                logger.Info("Model prior:   {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}%  R-work {:.3f} R-free {:.3f}",
                            sg.xhm(), rf_sd > 0 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash,
                            score.r_work, score.r_free);
            }
            result.t_rescore += seconds_since(t0);

            if (cand.vm < VM_MIN)
                cand.verdict = "MATTHEWS";
            else if (cand.copies_likely > 1)
                cand.verdict = "MULTI_COPY";
            else if (cand.clash_fraction > CLASH_MAX)
                cand.verdict = "CLASH";
            else if (!(cand.tf_z >= TF_Z_PLACED))
                cand.verdict = "NOT_PLACED";
            else
                cand.verdict = "PLACED";
            if (cand.verdict == "PLACED") {
                t0 = std::chrono::steady_clock::now();
                cand.null_sigma = null_sigma(*cand.placed, orientation, sg, *merge.merged, d_min, engines.get(),
                                             nthreads);
                result.t_null += seconds_since(t0);
                if (!(cand.null_sigma >= NULL_SIGMA_MIN))
                    cand.verdict = "NULL";
            }
            result.candidates.push_back(std::move(cand));
        }
    };

    // --- the decision ---
    // Within one merge the candidates are scored on the same reflections, so their R-free values can
    // be compared directly: the placed candidate with the lowest R-free is chosen where it leads every
    // other candidate of that merge that packs by R_FREE_MARGIN.
    const auto choose_in = [&](int merge) -> int {
        int best = -1;
        for (size_t i = 0; i < result.candidates.size(); i++) {
            const MRCandidate &c = result.candidates[i];
            if (c.merge == merge && c.verdict == "PLACED" && (best < 0 || c.r_free < result.candidates[best].r_free))
                best = static_cast<int>(i);
        }
        if (best < 0)
            return -1;
        for (size_t i = 0; i < result.candidates.size(); i++) {
            const MRCandidate &c = result.candidates[i];
            if (static_cast<int>(i) == best || c.merge != merge || c.clash_fraction > CLASH_MAX)
                continue;
            if (!(c.r_free - result.candidates[best].r_free >= R_FREE_MARGIN)) {
                result.decision = fmt::format("UNDECIDED: {} and {} within {:.2f} in R-free",
                                              result.candidates[best].space_group.xhm(), c.space_group.xhm(),
                                              R_FREE_MARGIN);
                return -2;
            }
        }
        return best;
    };
    // A merge whose every group the model cannot be in at all - one copy does not fit the cell, or
    // every placement clashes with its own symmetry mates.
    const auto cannot_hold = [&](int merge) {
        bool any = false;
        for (const MRCandidate &c : result.candidates)
            if (c.merge == merge) {
                any = true;
                if (c.verdict != "MATTHEWS" && c.verdict != "CLASH")
                    return false;
            }
        return any;
    };
    // The twin hypothesis' merge is searched only where it can decide something: where the answer's
    // merge cannot hold the model in any of its groups.
    search_merge(0);
    const int answer = choose_in(0);
    if (answer >= 0) {
        result.chosen = answer;
    } else if (answer == -1 && merges.size() > 1 && cannot_hold(0)) {
        search_merge(1);
        const int twin = choose_in(1);
        if (twin >= 0)
            result.chosen = twin;
        else if (twin == -1)
            result.decision = "NOT_PLACED: no group of either merge holds a placed model";
    } else if (answer == -1) {
        result.decision = "NOT_PLACED: no group of the answer's merge holds a placed model";
    }
    if (result.chosen >= 0) {
        const MRCandidate &c = result.candidates[result.chosen];
        const gemmi::SpaceGroup &data = merges[0].space_group;
        if (c.merge == 0 && c.space_group.xhm() == data.xhm())
            result.decision = "AGREES";
        else if (c.merge == 0 && enantiomorphs(data, c.space_group))
            result.decision = "ENANTIOMORPH";
        else if (c.merge == 0)
            result.decision = "SPACE_GROUP";
        else
            result.decision = "TWIN_SUBGROUP";
    }
    result.ok = true;
    logger.Info("Model prior: {} (model {:.1f} s, transform {:.2f} s, rotation {:.2f} s, translation {:.2f} s, "
                "rescoring {:.2f} s, null {:.2f} s)",
                result.chosen >= 0 ? result.decision + " - " + result.candidates[result.chosen].space_group.xhm()
                                   : result.decision,
                result.t_model, result.t_transform, result.t_rotation, result.t_translation, result.t_rescore,
                result.t_null);
    return result;
}