diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 7b6a77718..4cb751bbf 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -15,6 +15,7 @@ * The rugnux results report names the groups the data cannot separate, the enantiomorph state and any refused higher point group as `SPACE_GROUP_ALTERNATIVES`, `SPACE_GROUP_ENANTIOMORPH`, `SPACE_GROUP_REFUSED_POINT_GROUP` and `SPACE_GROUP_REFUSED_REASON`; `REPORT_VERSION` is 5. * `rugnux` handles symmetry better: the lattice, the point group, the setting and the systematic absences. * `rugnux` fits the direction of the goniometer axis and not its length, so the cell chosen by the first pass is the one its own refinement scored. +* `rugnux` settles two first-pass cells whose volumes differ by a whole-number factor on which of them accounts for more of the found spots, so a crystal with a real pseudo-translation keeps its full-length axis. * `rugnux` reports the detector geometry it measured - the direct beam, the tilt and the beam centre - and what a single sweep can and cannot determine. * `rugnux` reports the resolution the CC1/2 fit reached and the strong-direction diffraction limit. * The rugnux anisotropy warning says which reciprocal axis each of the two diffraction limits it quotes belongs to. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 3df31354d..9d5348647 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1784,6 +1784,86 @@ namespace { } return best; } + + // === Is the larger of two integer-related cells real, or a spurious axis multiple? === + // + // The validation-FRAME count cannot answer it: a spurious multiple indexes every frame its true + // sub-cell does, so both reach the same total and the count saturates. The rule that then decided + // the pair was unconditionally against the larger cell, so on a crystal with a real + // pseudo-translation the true cell could not win in any scheme order. + // + // Ask the two candidates the same question at the granularity where it does not saturate: how + // many of the validation frames' SPOTS does each account for? That comparison leans towards the + // smaller cell by construction, and needs no threshold to do so. Acceptance is a fractional-Miller + // test, so multiplying an axis by n multiplies that axis's residual by n: every spot the larger + // cell places on the sub-lattice the smaller cell places n times more accurately, and the spots in + // the tolerance margin the larger cell simply loses. The only thing that can pay for that loss is + // the class of reflections the larger cell ADDS, which is all but empty for a spurious multiple + // and holds the superstructure's satellite rows for a real one. So the larger cell wins the count + // only when the extra periodicity is really there, and the bar it has to clear is set by the + // data rather than by a constant. + // + // What that added class holds is reported alongside, because it is the physics the count is a + // consequence of. Multiplying an axis by n is a CENTRING CONDITION in the larger cell, and the + // index-n sublattices of its Miller lattice are the kernels of the linear forms f = (u,v,w) mod n + // with gcd(u,v,w,n) = 1, so enumerate them and report the emptiest: that is the condition the + // larger cell would have to satisfy to BE its sub-cell. It is reported and NOT thresholded, and + // that is deliberate. Measured on real data the number does not have two populations to separate: + // an absent class reads a few tenths of a percent where nothing but the crystal is in the spot + // list and several percent where the finder also picked up things that are not this crystal's + // reflections, a real superstructure layer reads tens of percent, and a pair whose volumes are in + // an integer ratio without either cell being the other's supercell reads the chance value (n-1)/n + // - so a floor placed anywhere in that range is calibrated on whichever crystals were to hand. + // + // Following the reduction of an axis multiple to a centring condition, and its measurement as an + // absence test, from this project's own lattice-arbiter analysis. + std::vector> IndexedMillerIndices(const CrystalLattice &latt, + const std::vector &cloud, float tol) { + const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2(); + const float tol_sq = tol * tol; + std::vector> hkl; + hkl.reserve(cloud.size()); + for (const Coord &s : cloud) { + const float fh = a * s, fk = b * s, fl = c * s; // Coord operator* = dot = Miller index + const float h = std::round(fh), k = std::round(fk), l = std::round(fl); + const float dh = fh - h, dk = fk - k, dl = fl - l; + if (dh * dh + dk * dk + dl * dl < tol_sq) + hkl.push_back({static_cast(h), static_cast(k), static_cast(l)}); + } + return hkl; + } + + struct AddedClassEvidence { + double occupancy = 1.0; // share of the larger cell's indexed spots OFF the sub-lattice + int u = 0, v = 0, w = 0; // the emptiest form, i.e. the condition that would collapse it + }; + + AddedClassEvidence AddedClassOccupancy(const std::vector> &hkl, int n) { + AddedClassEvidence out; + if (n < 2 || hkl.empty()) + return out; + for (int u = 0; u < n; u++) + for (int v = 0; v < n; v++) + for (int w = 0; w < n; w++) { + if (u == 0 && v == 0 && w == 0) + continue; + // A form sharing a factor with n has a kernel of index < n - it describes a + // different, coarser condition than the one the volume ratio asks about, and + // admitting it could only make the class look emptier than it is. + if (std::gcd(std::gcd(std::gcd(u, v), w), n) != 1) + continue; + size_t off = 0; + for (const auto &m : hkl) + if (((u * m[0] + v * m[1] + w * m[2]) % n + n) % n != 0) + ++off; + const double occ = static_cast(off) / static_cast(hkl.size()); + if (occ < out.occupancy) { + out.occupancy = occ; + out.u = u; out.v = v; out.w = w; + } + } + return out; + } } ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) { @@ -2177,6 +2257,31 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // many - this is the discriminator that lets the best scheme win. const std::vector validation = select_equally_spaced_image_ordinals( images_to_process, std::min(images_to_process, 60)); + + const float index_tol = experiment_.GetIndexingSettings().GetTolerance(); + + // The validation frames' spots in reciprocal space, under a candidate's own refined geometry + // and axis. The same construction the rotation indexer makes on its accumulated frames, but + // over the validation set and over each frame's WHOLE spot list - the arbiter counts spots and + // wants the deep sample, where the first pass wants a bounded one. A superstructure layer is + // faintest exactly where the first pass's accumulation cap cuts. + auto validation_cloud = [&](const RotationIndexerResult &r) { + std::vector cloud; + if (!r.axis) + return cloud; + for (const int ordinal : validation) { + const auto it = spot_cache.find(ordinal); + if (it == spot_cache.end()) + continue; + const float ang = rot_angle(ordinal).value_or( + r.axis->GetAngle_deg(static_cast(ordinal)) + r.axis->GetWedge_deg() / 2.0f); + const auto rot = r.axis->GetTransformationAngle(ang); + for (const auto &sp : it->second) + cloud.emplace_back(rot * sp.ReciprocalCoord(r.geom)); + } + return cloud; + }; + auto validation_settings = config_.spot_finding; validation_settings.indexing = true; validation_settings.quick_integration = false; @@ -2363,6 +2468,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // ratio separates a real axis multiplication from a centering coincidence. Volumes are // primitive (see above), so a pure setting difference is a ratio of 1 and never fires. bool integer_subcell = false; + bool harmonic_settled = false; double harmonic_ratio = 0.0; if (bp.result.has_value() && vol > 1.0 && bp.vol > 1.0) { const bool tied = static_cast(score) >= bp.score * 0.9f - 0.5f; @@ -2370,6 +2476,47 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const double nearest = std::round(ratio); const bool integer_multiple = nearest >= 2.0 && std::abs(ratio - nearest) < 0.15; integer_subcell = tied && !clearly_more && integer_multiple && vol < bp.vol; + + // Ask the data which of the pair is right instead of assuming the larger is + // spurious: which cell accounts for more of the validation frames' spots (see + // IndexedMillerIndices)? The measurement is made on their FULL spot lists, which + // reach far deeper into each frame's intensity distribution than the first pass's + // own accumulation cap - and a superstructure layer is faintest exactly where that + // cap cuts. + if (tied && !clearly_more && integer_multiple) { + const bool cand_is_larger = vol > bp.vol; + const RotationIndexerResult &larger = cand_is_larger ? result : *bp.result; + const RotationIndexerResult &smaller = cand_is_larger ? *bp.result : result; + const auto hkl = IndexedMillerIndices(larger.lattice, validation_cloud(larger), + index_tol); + const size_t n_smaller = IndexedMillerIndices(smaller.lattice, + validation_cloud(smaller), + index_tol).size(); + const auto ev = AddedClassOccupancy(hkl, static_cast(nearest)); + const bool larger_is_real = hkl.size() > n_smaller; + logger.Info("Axis-harmonic arbiter: '{}' ({:.0f} A^3) vs '{}' ({:.0f} A^3), {:.2f}x - " + "the larger cell accounts for {} validation spots against the smaller " + "cell's {}, and {:.1f}% of them lie off its ({},{},{}) sub-lattice, so " + "the class it adds is {}", + name, vol, bp.name, bp.vol, ratio, hkl.size(), n_smaller, + 100.0 * ev.occupancy, ev.u, ev.v, ev.w, + larger_is_real ? "occupied enough to pay for the longer axis - the " + "larger cell is real" + : "not occupied enough to pay for the longer axis - the " + "larger cell is a multiple"); + // Adopt the candidate exactly when its role matches the verdict. + integer_subcell = (cand_is_larger == larger_is_real); + // Only a win for the larger cell settles the pair. Declining it is a fall back + // to the default, and the warning below - which names the beam centre - is then + // still the most useful thing the run can say. It is the right diagnosis in + // particular for the failure this arbiter cannot see: a centre error along the + // spindle translates the derotated cloud rigidly by HALF a lattice spacing, and + // a half-offset lattice is indexable only on a doubled axis. Such a cell fits no + // index-n sublattice at all, so its added class reads the chance value (n-1)/n - + // the largest the occupancy can be - and it accounts for far fewer spots than + // the true cell, which is what this test sees it by. + harmonic_settled = larger_is_real; + } // sqrt(3) as well as an integer: a hexagonal supercell that lengthens one axis by // sqrt(3) is the other harmonic seen on real data and its volume ratio is not a // whole number. This is a REPORT and not a decision - it is deliberately taken @@ -2388,7 +2535,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // between the deposited cell and a halved axis, and the halved one is the answer that // indexed MORE frames (100.00 % against 99.23 %). So the frame count cannot arbitrate // this pair, and until this line nothing in the run said what the pair meant. - if (harmonic_ratio > 0.0) + // ... unless the arbiter above established the larger cell on the data, which is + // exactly the question this warning says the run cannot answer. + if (harmonic_ratio > 0.0 && !harmonic_settled) logger.Warning("Scheme '{}' (primitive cell volume {:.0f} A^3) and '{}' ({:.0f} A^3) " "differ by {:.2f}x - one is an axis harmonic of the other. Which of " "the two is right is settled by the beam centre to a fraction of a " @@ -2397,10 +2546,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b "(--estimate-beam-center)", name, vol, bp.name, bp.vol, harmonic_ratio); if (!bp.result.has_value() || clearly_more || integer_subcell) { - if (integer_subcell) + if (integer_subcell && vol < bp.vol) logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at " "equal frame count - adopting the smaller true cell", name, vol, std::round(bp.vol / vol), bp.name, bp.vol); + else if (integer_subcell) + logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x supercell of '{}' (vol {:.0f}) at " + "equal frame count and accounts for more of the spots - adopting the " + "larger true cell", name, vol, std::round(vol / bp.vol), bp.name, bp.vol); bp.score = score; bp.vol = vol; bp.name = name;