leonarski_fandClaude Opus 5
72cd91ccc3
A candidate can have several confirmed subgroups of the same order - 422 has both 4 and 222 - and on a twinned crystal the rival is not a harmless alternative: a P4 crystal twinned by 2[100] has 222 confirmed too, and 222 CONTAINS the twin laws, so its own merge b is already ballooned. The H test already answers to every tied parent; the systematic-b veto and rescue took whichever one the enumeration happened to list first (222 before 4, by space- group number), which disabled the veto on exactly the case it exists for. Take the smallest parent b, which is the conservative direction for both tests. The refusal message also quoted the raw parent b rather than the floored value the veto actually compared against. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
909 lines
48 KiB
C++
909 lines
48 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "SearchSpaceGroup.h"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <cstdint>
|
|
#include <iomanip>
|
|
#include <limits>
|
|
#include <map>
|
|
#include <numeric>
|
|
#include <sstream>
|
|
#include <string>
|
|
#include <tuple>
|
|
#include <unordered_map>
|
|
#include <vector>
|
|
|
|
namespace {
|
|
// A merged reflection, folded onto the +/- Friedel-equivalent it represents, used as a
|
|
// hash key to match symmetry-related reflections.
|
|
struct HKLKey {
|
|
int h = 0, k = 0, l = 0;
|
|
bool operator==(const HKLKey& o) const noexcept { return h == o.h && k == o.k && l == o.l; }
|
|
};
|
|
|
|
struct HKLKeyHash {
|
|
size_t operator()(const HKLKey& key) const noexcept {
|
|
auto mix = [](uint64_t x) {
|
|
x ^= x >> 33; x *= 0xff51afd7ed558ccdULL;
|
|
x ^= x >> 33; x *= 0xc4ceb9fe1a85ec53ULL;
|
|
x ^= x >> 33; return x;
|
|
};
|
|
return static_cast<size_t>(mix(static_cast<uint64_t>(key.h)) ^
|
|
(mix(static_cast<uint64_t>(key.k)) << 1) ^
|
|
(mix(static_cast<uint64_t>(key.l)) << 2));
|
|
}
|
|
};
|
|
|
|
HKLKey Canonicalize(int h, int k, int l, bool merge_friedel) {
|
|
if (merge_friedel && std::make_tuple(-h, -k, -l) < std::make_tuple(h, k, l))
|
|
return {-h, -k, -l};
|
|
return {h, k, l};
|
|
}
|
|
|
|
double PearsonCC(const std::vector<double>& x, const std::vector<double>& y) {
|
|
if (x.size() < 2)
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
|
|
double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
|
|
for (size_t i = 0; i < x.size(); ++i) {
|
|
sx += x[i]; sy += y[i];
|
|
sxx += x[i] * x[i]; syy += y[i] * y[i]; sxy += x[i] * y[i];
|
|
}
|
|
const double n = static_cast<double>(x.size());
|
|
const double vx = sxx - sx * sx / n;
|
|
const double vy = syy - sy * sy / n;
|
|
if (vx <= 0 || vy <= 0)
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
return (sxy - sx * sy / n) / std::sqrt(vx * vy);
|
|
}
|
|
|
|
// A reflection is extinct from lattice centering alone (independent of any screw/glide) when a
|
|
// centering translation makes its structure factor cancel. Mirrors the centering half of
|
|
// gemmi::GroupOps::is_systematically_absent, so screw absences can be judged separately.
|
|
bool CenteringAbsent(const gemmi::GroupOps& gops, const gemmi::Op::Miller& hkl) {
|
|
for (size_t i = 1; i < gops.cen_ops.size(); ++i) {
|
|
const auto& t = gops.cen_ops[i];
|
|
if ((t[0] * hkl[0] + t[1] * hkl[1] + t[2] * hkl[2]) % gemmi::Op::DEN != 0)
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// The reciprocal-space ROW a reflection lies on: its direction, reduced by the gcd and
|
|
// sign-canonicalised, so 0,0,l and 0,0,-l are one row and h,0,0 is a different one. Screw
|
|
// absences are judged against the other reflections of their own row (see below).
|
|
using AxialRow = std::array<int, 3>;
|
|
|
|
AxialRow RowOf(int h, int k, int l) {
|
|
const int g = std::gcd(std::gcd(std::abs(h), std::abs(k)), std::abs(l));
|
|
if (g > 0) { h /= g; k /= g; l /= g; }
|
|
if (std::make_tuple(-h, -k, -l) < std::make_tuple(h, k, l))
|
|
return {-h, -k, -l};
|
|
return {h, k, l};
|
|
}
|
|
|
|
// Median of an unordered set (reordered in place); 0 for an empty set.
|
|
double MedianOf(std::vector<double>& v) {
|
|
if (v.empty())
|
|
return 0.0;
|
|
const size_t mid = v.size() / 2;
|
|
std::nth_element(v.begin(), v.begin() + mid, v.end());
|
|
return v[mid];
|
|
}
|
|
|
|
std::string FormatDouble(double v, int decimals) {
|
|
std::ostringstream o;
|
|
o << std::fixed << std::setprecision(decimals) << v;
|
|
return o.str();
|
|
}
|
|
|
|
std::array<int, 9> RotKey(const gemmi::Op& op) {
|
|
std::array<int, 9> out{};
|
|
for (int i = 0; i < 3; ++i)
|
|
for (int j = 0; j < 3; ++j)
|
|
out[i * 3 + j] = op.rot[i][j];
|
|
return out;
|
|
}
|
|
|
|
// The rotation part of a space group in the reference setting (identity included), as a
|
|
// sorted list of matrices - the key that groups space groups into a candidate point group.
|
|
// It must be the rotation SET, not gemmi's PointGroup enum: P321 and P312 are both "32" yet
|
|
// have their 2-folds along different directions, and only the matrices tell them apart.
|
|
using RotationSet = std::vector<std::array<int, 9>>;
|
|
|
|
RotationSet RotationSetOf(const gemmi::SpaceGroup& sg) {
|
|
RotationSet out;
|
|
for (const auto& op : sg.operations().derive_symmorphic().sym_ops)
|
|
out.push_back(RotKey(op));
|
|
std::sort(out.begin(), out.end());
|
|
return out;
|
|
}
|
|
|
|
// Proper rotations of a crystal system's holohedry (the highest lattice symmetry it can host),
|
|
// in the reference setting. Any candidate point group must be a subgroup of this.
|
|
RotationSet HolohedryRotationSet(gemmi::CrystalSystem system) {
|
|
int number = 0;
|
|
switch (system) {
|
|
case gemmi::CrystalSystem::Triclinic: number = 1; break; // P1
|
|
case gemmi::CrystalSystem::Monoclinic: number = 3; break; // P2 (unique axis b)
|
|
case gemmi::CrystalSystem::Orthorhombic: number = 16; break; // P222
|
|
case gemmi::CrystalSystem::Tetragonal: number = 89; break; // P422
|
|
case gemmi::CrystalSystem::Trigonal: number = 155; break; // R32
|
|
case gemmi::CrystalSystem::Hexagonal: number = 177; break; // P622
|
|
case gemmi::CrystalSystem::Cubic: number = 207; break; // P432
|
|
}
|
|
const auto* sg = gemmi::find_spacegroup_by_number(number);
|
|
return sg ? RotationSetOf(*sg) : RotationSet{};
|
|
}
|
|
|
|
// A candidate point group: its proper rotations (reference setting) and a representative
|
|
// symmorphic space group (used when only the point group is wanted, or for display).
|
|
struct PointGroupInfo {
|
|
RotationSet rotation_set;
|
|
std::vector<gemmi::Op> rotations; // non-identity proper rotations
|
|
const gemmi::SpaceGroup* representative = nullptr;
|
|
};
|
|
|
|
// Enumerate candidate point groups. When a holohedry is given (from the lattice metric), keep
|
|
// only its subgroups - this both skips operators the lattice forbids and avoids accepting a
|
|
// coincidental higher symmetry; all subgroups down to P1 are still candidates.
|
|
std::vector<PointGroupInfo> EnumeratePointGroups(const std::optional<RotationSet>& holohedry) {
|
|
std::vector<PointGroupInfo> out;
|
|
std::map<RotationSet, size_t> index;
|
|
|
|
for (const auto& sg : gemmi::spacegroup_tables::main) {
|
|
if (!sg.is_sohncke() || !sg.is_reference_setting())
|
|
continue;
|
|
|
|
RotationSet rs = RotationSetOf(sg);
|
|
if (holohedry.has_value() &&
|
|
!std::includes(holohedry->begin(), holohedry->end(), rs.begin(), rs.end()))
|
|
continue;
|
|
|
|
auto it = index.find(rs);
|
|
size_t pos;
|
|
if (it == index.end()) {
|
|
PointGroupInfo info;
|
|
for (const auto& op : sg.operations().derive_symmorphic().sym_ops) {
|
|
if (op.rot == gemmi::Op::identity().rot)
|
|
continue;
|
|
info.rotations.push_back(gemmi::Op{op.rot, {0, 0, 0}, op.notation});
|
|
}
|
|
info.rotation_set = rs;
|
|
pos = out.size();
|
|
index[rs] = pos;
|
|
out.push_back(std::move(info));
|
|
} else {
|
|
pos = it->second;
|
|
}
|
|
// Prefer a symmorphic representative (the plain point-group setting).
|
|
auto& info = out[pos];
|
|
if (info.representative == nullptr ||
|
|
(!info.representative->is_symmorphic() && sg.is_symmorphic()))
|
|
info.representative = &sg;
|
|
}
|
|
return out;
|
|
}
|
|
}
|
|
|
|
SearchSpaceGroupResult SearchSpaceGroup(
|
|
const std::vector<MergedReflection>& merged,
|
|
const SearchSpaceGroupOptions& opt) {
|
|
|
|
SearchSpaceGroupResult result;
|
|
if (merged.empty())
|
|
return result;
|
|
|
|
const size_t n = merged.size();
|
|
|
|
// Flatten the reflections and mark which ones each stage may use. The correlation stage drops
|
|
// weak reflections; the absence stage must keep them - that is where the screw-axis signal is.
|
|
std::vector<int> H(n), K(n), L(n);
|
|
std::vector<double> I(n), Sigma(n), IoverSigma(n);
|
|
std::vector<HKLKey> key(n);
|
|
std::vector<char> pass_absence(n, 0), pass_cc(n, 0);
|
|
|
|
for (size_t i = 0; i < n; ++i) {
|
|
const auto& r = merged[i];
|
|
H[i] = r.h; K[i] = r.k; L[i] = r.l;
|
|
I[i] = r.I;
|
|
Sigma[i] = std::isfinite(r.sigma) && r.sigma > 0 ? r.sigma : 0.0;
|
|
key[i] = Canonicalize(r.h, r.k, r.l, opt.merge_friedel);
|
|
|
|
const bool finite = std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0 &&
|
|
std::isfinite(r.d) && r.d > 0;
|
|
const bool in_range = finite && (opt.d_min_limit_A <= 0 || r.d >= opt.d_min_limit_A);
|
|
IoverSigma[i] = finite ? r.I / r.sigma : 0.0;
|
|
pass_absence[i] = in_range;
|
|
// The correlation stage uses only genuinely-present reflections. Near-zero (systematically
|
|
// absent) reflections would otherwise form a second cluster at the origin and fake a high
|
|
// correlation for false operators - fatal on centered lattices, where half the reflections
|
|
// are extinct.
|
|
pass_cc[i] = in_range && IoverSigma[i] >= opt.present_i_over_sigma &&
|
|
(opt.min_i_over_sigma <= 0 || IoverSigma[i] >= opt.min_i_over_sigma);
|
|
}
|
|
|
|
// Resolution-normalised intensity E^2 = I / <I>(shell), from equal-count resolution shells over
|
|
// the reflections the absence test uses. Lets the absence test judge "present" by intensity
|
|
// magnitude, not by a possibly under-estimated sigma (see present_e_squared).
|
|
std::vector<double> Esq(n, 0.0);
|
|
{
|
|
std::vector<size_t> order;
|
|
order.reserve(n);
|
|
for (size_t i = 0; i < n; ++i)
|
|
if (pass_absence[i])
|
|
order.push_back(i);
|
|
std::sort(order.begin(), order.end(),
|
|
[&](size_t a, size_t b) { return merged[a].d > merged[b].d; }); // low res -> high res
|
|
const int bins = std::clamp(static_cast<int>(order.size() / 100), 1, 25);
|
|
const size_t per = (order.size() + bins - 1) / std::max(1, bins);
|
|
for (size_t b = 0; b * per < order.size(); ++b) {
|
|
const size_t lo = b * per, hi = std::min(order.size(), lo + per);
|
|
double sum = 0.0;
|
|
for (size_t j = lo; j < hi; ++j)
|
|
sum += I[order[j]];
|
|
const double mean = (hi > lo) ? sum / static_cast<double>(hi - lo) : 0.0;
|
|
for (size_t j = lo; j < hi; ++j)
|
|
Esq[order[j]] = mean > 0.0 ? I[order[j]] / mean : 0.0;
|
|
}
|
|
}
|
|
|
|
// Overlap guard (Stage A / correlation only): drop the extreme resolution-normalised-E tail, which
|
|
// on a two-lattice crystal is the one-sided overlap contamination that poisons the operator CC.
|
|
// See SearchSpaceGroupOptions::max_e_squared_for_cc. Absences (pass_absence) keep the full range.
|
|
if (opt.max_e_squared_for_cc > 0.0)
|
|
for (size_t i = 0; i < n; ++i)
|
|
if (pass_cc[i] && Esq[i] > opt.max_e_squared_for_cc)
|
|
pass_cc[i] = false;
|
|
|
|
std::unordered_map<HKLKey, int, HKLKeyHash> key_to_index;
|
|
key_to_index.reserve(n * 2);
|
|
for (size_t i = 0; i < n; ++i)
|
|
if (pass_absence[i])
|
|
key_to_index.emplace(key[i], static_cast<int>(i));
|
|
|
|
// --- Stage A: score each distinct rotation operator once ---
|
|
std::vector<uint32_t> visited(n, 0);
|
|
uint32_t epoch = 0;
|
|
|
|
auto score_operator = [&](const gemmi::Op& op) -> SpaceGroupOperatorScore {
|
|
++epoch;
|
|
std::vector<double> x, y;
|
|
for (size_t i = 0; i < n; ++i) {
|
|
if (!pass_cc[i] || visited[i] == epoch)
|
|
continue;
|
|
const auto m2 = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}});
|
|
const HKLKey k2 = Canonicalize(m2[0], m2[1], m2[2], opt.merge_friedel);
|
|
if (k2 == key[i])
|
|
continue; // reflection lies on this rotation axis
|
|
const auto it = key_to_index.find(k2);
|
|
if (it == key_to_index.end())
|
|
continue;
|
|
const int j = it->second;
|
|
if (!pass_cc[j])
|
|
continue;
|
|
x.push_back(I[i]);
|
|
y.push_back(I[j]);
|
|
visited[i] = epoch;
|
|
visited[j] = epoch;
|
|
}
|
|
SpaceGroupOperatorScore s;
|
|
s.op_triplet_hkl = op.as_hkl().triplet('h');
|
|
s.n_pairs = static_cast<int>(x.size());
|
|
s.cc = PearsonCC(x, y);
|
|
// Sigma-free disagreement over the same pairs (see SpaceGroupOptions::max_operator_h_ratio).
|
|
std::vector<double> hv;
|
|
hv.reserve(x.size());
|
|
for (size_t p = 0; p < x.size(); ++p) {
|
|
const double denom = x[p] + y[p];
|
|
if (denom > 0.0)
|
|
hv.push_back(std::fabs(x[p] - y[p]) / denom);
|
|
}
|
|
if (!hv.empty()) {
|
|
const size_t mid = hv.size() / 2;
|
|
std::nth_element(hv.begin(), hv.begin() + mid, hv.end());
|
|
s.h_stat = hv[mid];
|
|
}
|
|
s.present = s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc) &&
|
|
s.cc >= opt.min_operator_cc;
|
|
return s;
|
|
};
|
|
|
|
std::map<std::array<int, 9>, SpaceGroupOperatorScore> op_cache;
|
|
auto operator_score = [&](const gemmi::Op& op) -> const SpaceGroupOperatorScore& {
|
|
const auto rk = RotKey(op);
|
|
auto it = op_cache.find(rk);
|
|
if (it != op_cache.end())
|
|
return it->second;
|
|
return op_cache.emplace(rk, score_operator(op)).first->second;
|
|
};
|
|
|
|
// Conjugate rotations (symmetry-equivalent within the point group) relate symmetry-equivalent
|
|
// reflection sets, so on real data their CCs cluster; a noisy crystal can push one class member
|
|
// below min_operator_cc while the class is unmistakably present (e.g. one cubic 3-fold at 0.48
|
|
// among siblings at 0.53-0.66). Judge each conjugacy class by its mean CC, not its weakest
|
|
// member, so a genuine high-symmetry point group is not lost to one marginal operator. chi2_under
|
|
// (below) remains the safety net against a truly false promotion. Returns {all classes present,
|
|
// worst class-mean CC}.
|
|
auto point_group_present = [&](const std::vector<gemmi::Op>& rots) -> std::pair<bool, double> {
|
|
const size_t m = rots.size();
|
|
std::vector<int> cls(m, -1);
|
|
int n_cls = 0;
|
|
for (size_t i = 0; i < m; ++i) {
|
|
if (cls[i] >= 0)
|
|
continue;
|
|
cls[i] = n_cls;
|
|
for (size_t j = i + 1; j < m; ++j)
|
|
if (cls[j] < 0)
|
|
for (const auto& p : rots)
|
|
if ((p * rots[i] * p.inverse()).rot == rots[j].rot) {
|
|
cls[j] = n_cls;
|
|
break;
|
|
}
|
|
++n_cls;
|
|
}
|
|
bool ok = true;
|
|
double worst_mean = 1.0;
|
|
for (int c = 0; c < n_cls; ++c) {
|
|
// Average over the class members that actually have enough pairs to score; a single
|
|
// low-multiplicity / degenerate (NaN) operator in an otherwise strong class is skipped,
|
|
// not allowed to veto the class. The class must still have at least one scored member.
|
|
double sum_cc = 0.0;
|
|
int n_valid = 0;
|
|
for (size_t i = 0; i < m; ++i)
|
|
if (cls[i] == c) {
|
|
const auto& s = operator_score(rots[i]);
|
|
if (s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc)) {
|
|
sum_cc += s.cc;
|
|
++n_valid;
|
|
}
|
|
}
|
|
const double mean_cc = n_valid > 0 ? sum_cc / n_valid : 0.0;
|
|
worst_mean = std::min(worst_mean, mean_cc);
|
|
if (n_valid == 0 || mean_cc < opt.min_operator_cc)
|
|
ok = false;
|
|
}
|
|
return {ok, worst_mean};
|
|
};
|
|
|
|
std::optional<RotationSet> holohedry;
|
|
if (opt.lattice_system.has_value())
|
|
holohedry = HolohedryRotationSet(opt.lattice_system.value());
|
|
const auto point_groups = EnumeratePointGroups(holohedry);
|
|
|
|
// Reduced chi^2 of the intensities merged under a point group's rotations - how well its symmetry
|
|
// equivalents agree RELATIVE TO THEIR ERRORS. A real point group gives ~1; a false operator forces
|
|
// non-equivalent reflections together, so they disagree by many sigma and chi^2 blows up. This is
|
|
// more sensitive than R-meas to a strong pseudo-symmetry (where the intensities still correlate well
|
|
// - high operator CC - but not within their errors). Inverse-variance weighted mean per orbit, over
|
|
// the present (pass_cc) reflections.
|
|
auto chi2_under = [&](const std::vector<gemmi::Op>& rotations) -> double {
|
|
struct Acc { double sw = 0.0, swI = 0.0; int n = 0; };
|
|
std::unordered_map<HKLKey, Acc, HKLKeyHash> grp;
|
|
std::vector<HKLKey> rep(n);
|
|
for (size_t i = 0; i < n; ++i) {
|
|
if (!pass_cc[i] || !(Sigma[i] > 0.0))
|
|
continue;
|
|
HKLKey best = key[i];
|
|
for (const auto& op : rotations) {
|
|
const auto m = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}});
|
|
const HKLKey k2 = Canonicalize(m[0], m[1], m[2], opt.merge_friedel);
|
|
if (std::make_tuple(k2.h, k2.k, k2.l) < std::make_tuple(best.h, best.k, best.l))
|
|
best = k2;
|
|
}
|
|
rep[i] = best;
|
|
auto& g = grp[best];
|
|
const double w = 1.0 / (Sigma[i] * Sigma[i]);
|
|
g.sw += w; g.swI += w * I[i]; g.n += 1;
|
|
}
|
|
double chi2 = 0.0;
|
|
long dof = 0;
|
|
for (size_t i = 0; i < n; ++i) {
|
|
if (!pass_cc[i] || !(Sigma[i] > 0.0))
|
|
continue;
|
|
const auto& g = grp[rep[i]];
|
|
if (g.n < 2)
|
|
continue;
|
|
const double mean = g.swI / g.sw, dev = I[i] - mean;
|
|
chi2 += dev * dev / (Sigma[i] * Sigma[i]);
|
|
}
|
|
for (const auto& [k, g] : grp)
|
|
if (g.n >= 2)
|
|
dof += g.n - 1;
|
|
return dof > 0 ? chi2 / static_cast<double>(dof) : std::numeric_limits<double>::quiet_NaN();
|
|
};
|
|
|
|
// Extra intensity-proportional systematic error a point group's merge has to invoke to reconcile
|
|
// its symmetry equivalents: the smallest b for which sigma^2 + (b I)^2 brings the merged reduced
|
|
// chi^2 down to 1. A genuine symmetry needs almost none - its equivalents already agree within
|
|
// their errors, so the extra scatter is random and multiplicity absorbs it. A twin or pseudo-
|
|
// symmetry forces non-equivalent reflections together, and that disagreement scales with I, so b
|
|
// has to grow to swallow it (mirroring the merge error model's b / ISa collapse). This isolates
|
|
// the systematic part of the scatter, which the fixed-sigma chi^2 ratio cannot: a genuine but
|
|
// imperfectly-scaled high-symmetry merge and a twin can share a chi^2 ratio (~2) yet differ
|
|
// sharply here (a genuine cubic step b x1.04 vs a merohedral twin b x1.6).
|
|
auto merge_systematic_b = [&](const std::vector<gemmi::Op>& rotations) -> double {
|
|
struct Acc { double sw = 0.0, swI = 0.0; int n = 0; };
|
|
std::unordered_map<HKLKey, Acc, HKLKeyHash> grp;
|
|
std::vector<HKLKey> rep(n);
|
|
for (size_t i = 0; i < n; ++i) {
|
|
if (!pass_cc[i] || !(Sigma[i] > 0.0))
|
|
continue;
|
|
HKLKey best = key[i];
|
|
for (const auto& op : rotations) {
|
|
const auto m = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}});
|
|
const HKLKey k2 = Canonicalize(m[0], m[1], m[2], opt.merge_friedel);
|
|
if (std::make_tuple(k2.h, k2.k, k2.l) < std::make_tuple(best.h, best.k, best.l))
|
|
best = k2;
|
|
}
|
|
rep[i] = best;
|
|
auto& g = grp[best];
|
|
const double w = 1.0 / (Sigma[i] * Sigma[i]);
|
|
g.sw += w; g.swI += w * I[i]; g.n += 1;
|
|
}
|
|
std::vector<std::array<double, 3>> obs; // I, sigma, deviation-from-orbit-mean
|
|
for (size_t i = 0; i < n; ++i) {
|
|
if (!pass_cc[i] || !(Sigma[i] > 0.0))
|
|
continue;
|
|
const auto& g = grp[rep[i]];
|
|
if (g.n < 2)
|
|
continue;
|
|
obs.push_back({I[i], Sigma[i], I[i] - g.swI / g.sw});
|
|
}
|
|
if (obs.size() < 20)
|
|
return 0.0;
|
|
auto reduced_chi2 = [&](double b) {
|
|
double s = 0.0;
|
|
for (const auto& o : obs)
|
|
s += o[2] * o[2] / (o[1] * o[1] + (b * o[0]) * (b * o[0]));
|
|
return s / static_cast<double>(obs.size());
|
|
};
|
|
if (reduced_chi2(0.0) <= 1.0)
|
|
return 0.0;
|
|
double lo = 0.0, hi = 2.0; // b is a fraction of I; 2.0 = 200% is far past any real error model
|
|
for (int it = 0; it < 40; ++it) {
|
|
const double mid = 0.5 * (lo + hi);
|
|
(reduced_chi2(mid) > 1.0 ? lo : hi) = mid;
|
|
}
|
|
return 0.5 * (lo + hi);
|
|
};
|
|
|
|
// Operator-CC-confirmed candidates, each with its merge chi^2 and systematic-error b; chi2_ref =
|
|
// the most consistent.
|
|
struct PGCand { const PointGroupInfo* pg; int order; double min_class_cc; double chi2; double b_extra; };
|
|
int refused_order = 0;
|
|
std::string refused_pg_hm, refused_why;
|
|
std::vector<PGCand> pg_cands;
|
|
double chi2_ref = std::numeric_limits<double>::infinity();
|
|
for (const auto& pg : point_groups) {
|
|
const auto [present, min_class_cc] = point_group_present(pg.rotations);
|
|
if (!present)
|
|
continue;
|
|
const double ch = pg.rotations.empty() ? std::numeric_limits<double>::quiet_NaN()
|
|
: chi2_under(pg.rotations);
|
|
const double be = pg.rotations.empty() ? 0.0 : merge_systematic_b(pg.rotations);
|
|
pg_cands.push_back({&pg, static_cast<int>(pg.rotations.size()) + 1, min_class_cc, ch, be});
|
|
if (!pg.rotations.empty() && std::isfinite(ch))
|
|
chi2_ref = std::min(chi2_ref, ch);
|
|
}
|
|
|
|
// Choose the largest point group that is both operator-confirmed AND self-consistent (its merge
|
|
// chi^2 is not inflated past the miscalibration-widened bound below; ties -> higher min class CC).
|
|
// Identity (no operators) is always consistent, so it stays the P1 fallback.
|
|
const PointGroupInfo* best_pg = nullptr;
|
|
int best_pg_order = 0;
|
|
double best_pg_min_cc = -2.0;
|
|
for (const auto& c : pg_cands) {
|
|
// A genuine symmetry operator merges equivalent reflections, so it barely changes the reduced
|
|
// chi^2 relative to the best subgroup - across the whole rotation-test battery every correct
|
|
// point group stays within ~1.7x, even on weak or badly-integrated data (a cubic F432 chi2_ref
|
|
// 8.3 -> 1.15; a tetragonal P41212 -> 1.71). A twin law or pseudo-symmetry forces non-equivalent
|
|
// reflections together, so its ratio is markedly higher (a merohedral twin 2-fold: R3 3.02 ->
|
|
// R32 6.07, ratio 2.01). max_merge_chi2_ratio sits between the two. (An earlier log10(chi2_ref) widening
|
|
// compensated for an under-calibrated error model that inflated real-symmetry ratios with data
|
|
// weakness; the variance-floor fix removed that inflation, and the widening now only let the
|
|
// twin through, so it is gone.)
|
|
bool consistent = c.pg->rotations.empty() || !std::isfinite(c.chi2) ||
|
|
!std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * opt.max_merge_chi2_ratio;
|
|
|
|
// Systematic-error test vs the largest confirmed subgroup (by rotation-set inclusion): merging
|
|
// under a genuine operator gains multiplicity without intensity-proportional disagreement, so the
|
|
// merge error model's b barely moves; a merohedral twin forces non-equivalent reflections together
|
|
// and b balloons. It both RESCUES a genuine step whose chi^2 drifts just past the ratio bound
|
|
// (imperfectly scaled data) and VETOES a twin whose chi^2 now looks self-consistent but whose b
|
|
// balloons - the chi^2 ratio alone no longer separates them.
|
|
double parent_b = -1.0;
|
|
// Every confirmed subgroup of the largest order below this candidate. There can be more than one
|
|
// - 422 has both 4 and 222 - and on a twinned crystal the rival is not a harmless alternative: a
|
|
// P4 crystal twinned by 2[100] has its two twin 2-folds confirmed, so 222 is CC-confirmed too and
|
|
// CONTAINS the twin laws. Normalising the H test against it hides the twin among the promotion's
|
|
// own real operators (measured on the synthetic grid: ratio 8.19 against the true parent 4, 0.78
|
|
// against the rival 222). Which one is the true parent is exactly what is unknown here, so the
|
|
// promotion must answer to all of them.
|
|
std::vector<const PointGroupInfo *> parents;
|
|
if (!c.pg->rotations.empty()) {
|
|
int parent_order = 0;
|
|
for (const auto& s : pg_cands)
|
|
if (s.order < c.order && s.order >= parent_order
|
|
&& std::includes(c.pg->rotation_set.begin(), c.pg->rotation_set.end(),
|
|
s.pg->rotation_set.begin(), s.pg->rotation_set.end())) {
|
|
if (s.order > parent_order) {
|
|
parent_order = s.order;
|
|
parent_b = s.b_extra;
|
|
parents.clear();
|
|
} else {
|
|
// Tied parent. The b tests answer to the most damning of them, exactly as the H
|
|
// test below does: the SMALLEST parent b is the one that makes the veto easiest to
|
|
// trip and the rescue hardest to pass, and a rival subgroup that already contains
|
|
// the twin laws has its own b ballooned - taking it would hide the twin.
|
|
parent_b = std::min(parent_b, s.b_extra);
|
|
}
|
|
parents.push_back(s.pg);
|
|
}
|
|
}
|
|
|
|
// Sigma-free twin test: compare the disagreement H of the operators this promotion ADDS with the
|
|
// disagreement of the parent's own operators, measured on the same reflections. A real operator
|
|
// relates equal intensities and matches the parent; a twin law relates different ones and reads
|
|
// systematically higher. Skipped when either side has too few pairs to mean anything, and when
|
|
// there is no parent group to normalise against (the first step out of P1). Where several parents
|
|
// tie (see above), the promotion is judged on the most damning of them.
|
|
double h_ratio = std::numeric_limits<double>::quiet_NaN();
|
|
for (const auto *parent : parents) {
|
|
double h_new = 0.0, h_par = 0.0;
|
|
int n_new = 0, n_par = 0, pairs_new = 0, pairs_par = 0;
|
|
for (const auto &rot : c.pg->rotations) {
|
|
if (rot.rot == gemmi::Op::identity().rot)
|
|
continue;
|
|
const auto &os = operator_score(rot);
|
|
if (os.n_pairs < opt.min_pairs_per_operator)
|
|
continue;
|
|
const bool in_parent = std::binary_search(parent->rotation_set.begin(),
|
|
parent->rotation_set.end(), RotKey(rot));
|
|
if (in_parent) { h_par += os.h_stat; ++n_par; pairs_par += os.n_pairs; }
|
|
else { h_new += os.h_stat; ++n_new; pairs_new += os.n_pairs; }
|
|
}
|
|
if (n_new > 0 && n_par > 0 && pairs_new >= opt.min_pairs_for_h
|
|
&& pairs_par >= opt.min_pairs_for_h && h_par > 0.0) {
|
|
const double r = (h_new / n_new) / (h_par / n_par);
|
|
if (!std::isfinite(h_ratio) || r > h_ratio)
|
|
h_ratio = r;
|
|
}
|
|
}
|
|
// The chi^2 ratio is only trustworthy when the error model is calibrated. When even the best
|
|
// subgroup's reduced chi^2 (chi2_ref) is far above 1 - weak, low-resolution data whose merged
|
|
// sigmas are badly under-estimated - the ratio grows with point-group order for genuine high
|
|
// symmetry too and wrongly rejects it (a true weak F432 reaches ratio ~14). The systematic-b test
|
|
// re-fits its own error, so it stays valid under a broken sigma model: a genuine step's b barely
|
|
// moves (b-ratio ~1) while a twin's balloons. So once chi2_ref shows the error model is unreliable,
|
|
// a promotion is rescued on the b-test alone (subject to the balloon veto below); otherwise the
|
|
// rescue is confined to the narrow chi^2 band just past the ratio bound.
|
|
const bool miscalibrated = std::isfinite(chi2_ref) && chi2_ref > opt.chi2_ref_reliable;
|
|
if (!consistent && parent_b > 1e-4 && c.b_extra <= parent_b * opt.max_systematic_b_ratio
|
|
&& (miscalibrated || (std::isfinite(c.chi2) && std::isfinite(chi2_ref)
|
|
&& c.chi2 <= chi2_ref * opt.max_merge_chi2_rescue)))
|
|
consistent = true;
|
|
// Veto a chi^2-passing promotion whose b clearly ballooned (above the largest genuine step, below a
|
|
// twin); a genuine but imperfectly-scaled high-symmetry merge stays under the bound and is untouched.
|
|
// The parent b is floored (min_systematic_b_for_veto) so a near-zero parent on excellent data cannot
|
|
// fabricate a huge ratio out of a still-tiny absolute b (a genuine 422 at b=0.05 over a 222 parent at
|
|
// b=0.008 is not a twin - a real twin drives b to ~0.19 regardless).
|
|
bool b_refused = false;
|
|
if (consistent && parent_b > 1e-4
|
|
&& c.b_extra > std::max(parent_b, opt.min_systematic_b_for_veto) * opt.max_systematic_b_veto) {
|
|
consistent = false;
|
|
b_refused = true;
|
|
}
|
|
|
|
// The H test is a necessary condition for promotion where it can be computed: it is the only
|
|
// statistic measured to separate genuine symmetry from a merohedral twin across data amounts.
|
|
const bool h_refused = std::isfinite(h_ratio) && h_ratio > opt.max_operator_h_ratio;
|
|
if (h_refused)
|
|
consistent = false;
|
|
|
|
if (!consistent) {
|
|
// Record the highest-order refusal so the caller can say WHY it is processing lower.
|
|
if (c.order > refused_order && c.pg->representative) {
|
|
refused_order = c.order;
|
|
refused_pg_hm = c.pg->representative->point_group_hm();
|
|
if (h_refused)
|
|
refused_why = "operator disagreement H is " + FormatDouble(h_ratio, 2)
|
|
+ "x the parent's (bound " + FormatDouble(opt.max_operator_h_ratio, 2)
|
|
+ ") - the added operator relates unequal intensities, as a twin law does";
|
|
else if (b_refused)
|
|
// Against the FLOORED parent b, which is what the test compared - quoting the raw
|
|
// parent b reads as a far bigger balloon than the one that actually tripped the veto.
|
|
refused_why = "its merge's systematic error b is " + FormatDouble(c.b_extra, 3)
|
|
+ " against the subgroup's "
|
|
+ FormatDouble(std::max(parent_b, opt.min_systematic_b_for_veto), 3)
|
|
+ " (bound " + FormatDouble(opt.max_systematic_b_veto, 2)
|
|
+ "x) - the added operator forces unequal intensities together";
|
|
else if (std::isfinite(c.chi2) && std::isfinite(chi2_ref))
|
|
refused_why = "merge chi^2 is " + FormatDouble(c.chi2 / chi2_ref, 2)
|
|
+ "x the subgroup's (bound " + FormatDouble(opt.max_merge_chi2_ratio, 2) + ")";
|
|
else
|
|
refused_why = "the merge under it is not self-consistent";
|
|
}
|
|
continue;
|
|
}
|
|
if (c.order > best_pg_order || (c.order == best_pg_order && c.min_class_cc > best_pg_min_cc)) {
|
|
best_pg = c.pg;
|
|
best_pg_order = c.order;
|
|
best_pg_min_cc = c.min_class_cc;
|
|
}
|
|
}
|
|
|
|
for (const auto& [rk, s] : op_cache)
|
|
result.operator_scores.push_back(s);
|
|
std::sort(result.operator_scores.begin(), result.operator_scores.end(),
|
|
[](const auto& a, const auto& b) { return a.cc > b.cc; });
|
|
|
|
if (best_pg == nullptr) // should not happen (C1 always qualifies)
|
|
return result;
|
|
|
|
if (best_pg->representative)
|
|
result.point_group_hm = best_pg->representative->point_group_hm();
|
|
result.point_group_order = best_pg_order;
|
|
|
|
// Only report a refusal that is actually ABOVE what was adopted.
|
|
if (refused_order > best_pg_order) {
|
|
result.refused_point_group_hm = refused_pg_hm;
|
|
result.refused_reason = refused_why;
|
|
}
|
|
|
|
// --- Stage B: pick the space group within the point group ---
|
|
// Without screw/centering determination, return the symmorphic representative.
|
|
if (!opt.determine_space_group || best_pg->rotations.empty()) {
|
|
if (best_pg->representative)
|
|
result.best_space_group = *best_pg->representative;
|
|
return result;
|
|
}
|
|
|
|
for (const auto& sg : gemmi::spacegroup_tables::main) {
|
|
if (!sg.is_sohncke() || !sg.is_reference_setting() || RotationSetOf(sg) != best_pg->rotation_set)
|
|
continue;
|
|
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
SpaceGroupCandidateScore s{.space_group = sg};
|
|
double absent_sum = 0, present_sum = 0;
|
|
int present_n = 0;
|
|
|
|
// Judge centering and screw/glide absences on separate reflection sets. Lumping them lets
|
|
// a large, correct centering-absent set hide a few strong screw violations and over-claim
|
|
// screw axes (e.g. I4_132 on I432 data).
|
|
int centering_absent = 0, centering_violations = 0;
|
|
double centering_absent_sum = 0;
|
|
int present_strong = 0;
|
|
// A screw axis extinguishes only the reflections that lie ON it, so its absent class and the
|
|
// rest of that same axial row are collected apart from the general reflections and judged
|
|
// against each other, ROW BY ROW (see screw_e_squared below).
|
|
struct ScrewAbsent { AxialRow row; double e_squared; double i_over_sigma; };
|
|
std::vector<ScrewAbsent> screw_absent_refl;
|
|
std::map<AxialRow, std::vector<double>> row_present_esq;
|
|
|
|
for (size_t i = 0; i < n; ++i) {
|
|
if (!pass_absence[i])
|
|
continue;
|
|
const gemmi::Op::Miller hkl{{H[i], K[i], L[i]}};
|
|
// Present := statistically significant AND intensity-significant. The E^2 gate keeps a
|
|
// weak axial reflection with an under-estimated sigma (fake high I/sigma) from faking a
|
|
// screw-axis violation; it only relaxes "present", so it cannot over-call a screw whose
|
|
// predicted-absent class carries real intensity.
|
|
const bool present = IoverSigma[i] > opt.present_i_over_sigma &&
|
|
(opt.present_e_squared <= 0.0 || Esq[i] > opt.present_e_squared);
|
|
|
|
if (CenteringAbsent(gops, hkl)) {
|
|
s.absent_observed += 1;
|
|
absent_sum += IoverSigma[i];
|
|
centering_absent += 1;
|
|
centering_absent_sum += IoverSigma[i];
|
|
if (present) { s.absent_violations += 1; centering_violations += 1; }
|
|
} else if (gops.is_systematically_absent(hkl)) {
|
|
s.absent_observed += 1;
|
|
absent_sum += IoverSigma[i];
|
|
screw_absent_refl.push_back({RowOf(H[i], K[i], L[i]), Esq[i], IoverSigma[i]});
|
|
} else {
|
|
present_n += 1;
|
|
present_sum += IoverSigma[i];
|
|
if (present) present_strong += 1;
|
|
// A non-identity rotation of the group maps this reflection to itself, i.e. it lies on
|
|
// a rotation axis - the control class for that axis's screw absences.
|
|
if (gops.epsilon_factor_without_centering(hkl) > 1)
|
|
row_present_esq[RowOf(H[i], K[i], L[i])].push_back(Esq[i]);
|
|
}
|
|
}
|
|
|
|
// "Too strong to be systematically absent" is judged, for a SCREW, against the axial row the
|
|
// screw constrains rather than against the shell mean over all reflections. An axial row can be
|
|
// far stronger than an average reflection, and <I>(shell) falls with resolution while a
|
|
// systematically-absent reflection keeps a small non-decaying residual (background / profile
|
|
// leakage) - so at high resolution the plain E^2 cut turns those residuals into violations even
|
|
// though the reflections next to them in the same row are tens of times stronger. That cost a
|
|
// real P4_1 2_1 2 crystal its 4_1: 18 of its 47 absent 00l crossed the cut, all beyond 3.7 A,
|
|
// at 1-2% of the l=4n reflections beside them. Scaling by the row's own median E^2 removes the
|
|
// resolution dependence; the scale is floored at 1 so a row weaker than average keeps the plain
|
|
// cut, which makes this a rescue only - a screw can be recovered by it, never lost.
|
|
//
|
|
// Row by row, not pooled: a 4_1 along c and a 2_1 along a are separate conditions with separate
|
|
// control rows, and on the same crystal one row can be 20x an average reflection while another
|
|
// is half of one. Pooling them lets the weak rows set the threshold for the strong one and the
|
|
// rescue never fires (that very crystal pooled to a row median of 0.7 and stayed at P42_12).
|
|
std::map<AxialRow, double> row_median;
|
|
for (auto& [row, esq] : row_present_esq)
|
|
row_median[row] = MedianOf(esq);
|
|
|
|
const int screw_absent = static_cast<int>(screw_absent_refl.size());
|
|
int screw_violations = 0;
|
|
std::vector<double> screw_absent_esq, screw_control_esq;
|
|
for (const auto& a : screw_absent_refl) {
|
|
const auto it = row_median.find(a.row);
|
|
const double row_scale = it == row_median.end() ? 1.0 : std::max(1.0, it->second);
|
|
if (a.i_over_sigma > opt.present_i_over_sigma &&
|
|
(opt.present_e_squared <= 0.0 || a.e_squared > opt.present_e_squared * row_scale))
|
|
++screw_violations;
|
|
screw_absent_esq.push_back(a.e_squared);
|
|
if (it != row_median.end())
|
|
screw_control_esq.push_back(it->second);
|
|
}
|
|
s.absent_violations += screw_violations;
|
|
// Reported evidence: the absent class and the rows it is judged against (see the header).
|
|
s.screw_absent_median_e_squared = MedianOf(screw_absent_esq);
|
|
s.screw_row_median_e_squared = MedianOf(screw_control_esq);
|
|
|
|
if (s.absent_observed > 0)
|
|
s.absent_mean_i_over_sigma = absent_sum / s.absent_observed;
|
|
if (present_n > 0)
|
|
s.present_mean_i_over_sigma = present_sum / present_n;
|
|
|
|
// Centering is judged by class STRENGTH, not a per-reflection violation count. A real centering
|
|
// cancels structure factors, so its absent class is systematically weak - its mean signed
|
|
// I/sigma sits well below the present class - regardless of noise or obverse/reverse twinning;
|
|
// a false centering leaves the "absent" class as strong as the present one (mean ratio ~1). The
|
|
// count-of-strong-violations gate is brittle on noisy/twinned data, where enough genuinely-absent
|
|
// reflections randomly clear I/sigma>3 to trip the 10% bound though the class is 3-4x weaker (a
|
|
// true R3 at 13.5% violations, absent 1.7 vs present 6.0). The mean is well-determined here
|
|
// because a centering-absent class holds a third-to-half of all reflections. Screws keep the
|
|
// count gate: their predicted-absent class is a handful of axial reflections, too few to average
|
|
// - what they get instead is a row-relative threshold for counting a violation at all.
|
|
const double present_mean = present_n > 0 ? present_sum / present_n : 0.0;
|
|
const double centering_absent_mean =
|
|
centering_absent > 0 ? centering_absent_sum / centering_absent : 0.0;
|
|
// The centering-absent class proves itself weak in either of two floor-independent ways; a
|
|
// FALSE centering (absent as strong as present) fails both:
|
|
// (1) mean signed I/sigma well below the present class, OR
|
|
// (2) its strong-reflection RATE well below the present class's own strong rate.
|
|
// (2) is needed because weak / low-energy data carry a positive intensity floor (background /
|
|
// profile leakage) that lifts <I/s>abs to ~1.5-2.3 even for genuinely extinct reflections; when
|
|
// the present class is itself weak (small present_mean) that additive floor inflates the mean
|
|
// ratio past the bound and hides a real centering - e.g. an I-centred cubic crystal at low
|
|
// energy, whose true I-centering sat at ratio ~0.57. Normalising the violation count by the
|
|
// present class's own strong rate cancels the shared floor and stays reliable on weak data
|
|
// (both rates shrink together).
|
|
const double present_strong_rate =
|
|
present_n > 0 ? static_cast<double>(present_strong) / present_n : 0.0;
|
|
const double centering_violation_rate =
|
|
centering_absent > 0 ? static_cast<double>(centering_violations) / centering_absent : 0.0;
|
|
const bool centering_ok = centering_absent == 0
|
|
|| (present_n > 0 && centering_absent_mean <= opt.max_absent_present_ratio * present_mean)
|
|
|| (present_strong_rate > 0.0
|
|
&& centering_violation_rate <= opt.max_absent_present_ratio * present_strong_rate);
|
|
const bool screw_ok = screw_absent == 0 ||
|
|
screw_violations <= opt.max_absent_violation_fraction * screw_absent;
|
|
s.consistent = centering_ok && screw_ok;
|
|
result.candidates.push_back(std::move(s));
|
|
}
|
|
|
|
// A candidate is eligible when its absences are confirmed and there are enough of them to
|
|
// trust (the symmorphic group, with no absences, is always eligible as the fallback). Rank
|
|
// eligible candidates by how many absences they GENUINELY explain - absent_observed minus the
|
|
// violations, not the gross count. A false super-centering over-claims: F222 on a C222 crystal
|
|
// predicts every C absence (all genuinely weak) PLUS a block of C-present reflections it wrongly
|
|
// calls absent, so its gross count is larger yet its net count only equals C222's. Its diluted
|
|
// absent class (many true zeros + a strong block) also slips under the strength/rate gate, so the
|
|
// gate cannot veto it alone; netting the violations puts the two level, and the fewer-violations
|
|
// and lower-number tie-breaks then keep the honest, less-centred C222. The ranking is symmetric:
|
|
// on a genuine F222 crystal F explains strictly more weak absences and still wins.
|
|
auto net_absent = [](const SpaceGroupCandidateScore& s) {
|
|
return s.absent_observed - s.absent_violations;
|
|
};
|
|
auto eligible = [&](const SpaceGroupCandidateScore& s) {
|
|
return s.consistent && (s.absent_observed == 0 || s.absent_observed >= opt.min_absent_observed);
|
|
};
|
|
std::sort(result.candidates.begin(), result.candidates.end(),
|
|
[&](const SpaceGroupCandidateScore& a, const SpaceGroupCandidateScore& b) {
|
|
if (eligible(a) != eligible(b))
|
|
return eligible(a);
|
|
if (net_absent(a) != net_absent(b))
|
|
return net_absent(a) > net_absent(b);
|
|
if (a.absent_violations != b.absent_violations)
|
|
return a.absent_violations < b.absent_violations; // prefer the honest, less over-claiming group
|
|
// Genuinely indistinguishable (e.g. I23 vs I2_13, or an enantiomorphic pair): lower
|
|
// space-group number is the representative.
|
|
return a.space_group.number < b.space_group.number;
|
|
});
|
|
|
|
if (!result.candidates.empty() && eligible(result.candidates.front())) {
|
|
// Alternatives are only the candidates with the SAME absence signature - identical absent AND
|
|
// violation counts - as the winner: the enantiomorphic / origin-ambiguous partners the data
|
|
// truly cannot separate. A super-centering that nets the same count but over-claims differs in
|
|
// its violation count and is therefore not reported as an equal alternative.
|
|
const int sel_absent = result.candidates.front().absent_observed;
|
|
const int sel_violations = result.candidates.front().absent_violations;
|
|
for (auto& s : result.candidates) {
|
|
if (!eligible(s) || s.absent_observed != sel_absent || s.absent_violations != sel_violations)
|
|
continue;
|
|
s.selected = true;
|
|
if (!result.best_space_group.has_value())
|
|
result.best_space_group = s.space_group; // representative (lowest number)
|
|
else
|
|
result.alternatives.push_back(s.space_group);
|
|
}
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result,
|
|
size_t max_candidates_to_print) {
|
|
std::ostringstream os;
|
|
|
|
if (!result.refused_point_group_hm.empty())
|
|
os << "Higher symmetry " << result.refused_point_group_hm << " was confirmed by the operator "
|
|
"correlations but REFUSED: " << result.refused_reason << ".\n"
|
|
" Processing in the lower symmetry, which is the recoverable direction - if this is a "
|
|
"twin, merging in the higher group would average non-equivalent reflections together and "
|
|
"hide the twin law.\n";
|
|
os << "Point group: " << (result.point_group_hm.empty() ? "?" : result.point_group_hm)
|
|
<< " (from intensity correlations)\n";
|
|
|
|
os << " " << std::setw(14) << std::left << "operator" << std::right
|
|
<< std::setw(9) << "CC" << std::setw(10) << "pairs" << std::setw(9) << "symm" << "\n";
|
|
for (const auto& s : result.operator_scores) {
|
|
os << " " << std::setw(14) << std::left << s.op_triplet_hkl << std::right
|
|
<< std::setw(9) << std::fixed << std::setprecision(3) << s.cc
|
|
<< std::setw(10) << s.n_pairs
|
|
<< std::setw(9) << (s.present ? "yes" : "no") << "\n";
|
|
}
|
|
|
|
os << "\nSpace-group candidates\n";
|
|
os << " " << std::setw(10) << std::left << "SG" << std::right
|
|
<< std::setw(9) << "absent" << std::setw(7) << "viol"
|
|
<< std::setw(11) << "<I/s>abs" << std::setw(11) << "<I/s>pres"
|
|
<< std::setw(10) << "E2 screw" << std::setw(9) << "E2 row"
|
|
<< std::setw(6) << "OK" << "\n";
|
|
|
|
const size_t count = std::min(max_candidates_to_print, result.candidates.size());
|
|
for (size_t i = 0; i < count; ++i) {
|
|
const auto& c = result.candidates[i];
|
|
os << (c.selected ? "* " : " ")
|
|
<< std::setw(10) << std::left << c.space_group.short_name() << std::right
|
|
<< std::setw(9) << c.absent_observed << std::setw(7) << c.absent_violations
|
|
<< std::setw(11) << std::fixed << std::setprecision(2) << c.absent_mean_i_over_sigma
|
|
<< std::setw(11) << std::fixed << std::setprecision(2) << c.present_mean_i_over_sigma
|
|
<< std::setw(10) << std::fixed << std::setprecision(3) << c.screw_absent_median_e_squared
|
|
<< std::setw(9) << std::fixed << std::setprecision(3) << c.screw_row_median_e_squared
|
|
<< std::setw(6) << (c.consistent ? "yes" : "no") << "\n";
|
|
}
|
|
os << " absent/viol = reflections the group predicts absent, and how many are nonetheless present.\n"
|
|
" E2 screw / E2 row = median I/<I>(shell) of the SCREW-absent reflections and of the rest of\n"
|
|
" their axial rows - a real screw leaves the first far below the second. The <I/s> columns are\n"
|
|
" the centering evidence; they say little about screws (the merged sigma floor pins I/sigma).\n";
|
|
|
|
if (result.best_space_group.has_value()) {
|
|
os << "Best space group: " << result.best_space_group->short_name();
|
|
for (const auto& alt : result.alternatives)
|
|
os << " or " << alt.short_name();
|
|
if (!result.alternatives.empty())
|
|
os << " (indistinguishable from these data)";
|
|
os << "\n";
|
|
} else {
|
|
os << "Best space group: none determined\n";
|
|
}
|
|
|
|
return os.str();
|
|
}
|