rugnux: headline the twin law whose pairs agree best, at its own fraction

TWIN_LAW= named the first operator the refused group added and quoted beside it the
mean H over all the operators the refusal added, or failing that the L-test's fraction -
neither the headlined law's own number. A 3-fold crystal twinned at 0.4 by one of its
three laws was headlined at 0.20 while TWIN_LAW_1, the same law, read H 0.10, i.e. 0.40
(the deposition's 0.39, xtriage's 0.38). The headline is now the law of the TWIN_LAW_n
list whose pairs agree best, at 0.5 - its H, so the two lines cannot disagree.

The TWIN_LAW_n operators are printed as Miller-index operators (as TWIN_ZONE and TWIN_LAW
print theirs): find_twin_laws returns real-space operations, and the two differ for every
law that is not symmetric - a pseudo-merohedral law of a centred lattice read literally as
hkl relates unrelated reflections.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK

rugnux/ResultReport.cpp

@@ -2608,8 +2608,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
             if (!result.twinning.twin_law_operator.empty()) {
                 Add(s, KeyText("TWIN_LAW", fmt::format("{} fraction {:.2f}", result.twinning.twin_law_operator,
                                                        result.twinning.twin_law_fraction)));
-                Add(s, Prose(fmt::format(" Twin law identified: operator {} at a fraction of about {:.2f} (a lower bound from\n"
-                                         " the operator's own disagreement). The mmCIF carries it as _pdbx_reflns_twin; refine\n"
+                Add(s, Prose(fmt::format(" Twin law identified: operator {} at a fraction of about {:.2f} (0.5 - H of the\n"
+                                         " pairs it relates, as TWIN_LAW_n reads it). The mmCIF carries it as _pdbx_reflns_twin; refine\n"
                                          " against these data in this space group with that law.",
                                          result.twinning.twin_law_operator, result.twinning.twin_law_fraction)));
             }

rugnux/RugnuxScaleMerge.cpp

@@ -4505,7 +4505,11 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
         // would not do: the call replaces the whole struct.)
         result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group;
         // A refusal that stood on the twin-immune zone names the law: the first operator the
-        // refused group adds, at the fraction its own disagreement implies (a lower bound).
+        // refused group adds, at the fraction its own disagreement implies (a lower bound). Both are
+        // replaced below by the law whose own pairs agree best, once every law of the lattice has been
+        // read: the fraction here is the mean H over ALL the operators the refused group added, or the
+        // L-test's, and neither is the headlined law's own - a 3-fold crystal twinned at 0.4 by one of
+        // its three laws was headlined at 0.20 beside its law's own H reading 0.40.
         if (twin_sg && result.space_group_search
             && result.space_group_search->refused_point_group_representative
             && result.twin_zone_refusal_stood) {
@@ -4628,11 +4632,14 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
             };
             std::optional<OperatorScorer> law_correlation;
             std::string laws;
+            std::string best_law;          // the law whose pairs agree best, and its H
+            double best_law_h = 1.0;
             for (const auto &op : twin_ops) {
-                std::string t = op.triplet();
-                std::replace(t.begin(), t.end(), 'x', 'h');
-                std::replace(t.begin(), t.end(), 'y', 'k');
-                std::replace(t.begin(), t.end(), 'z', 'l');
+                // find_twin_laws returns real-space operations; the Miller-index operator is the
+                // transpose (an hkl triplet, as TWIN_ZONE and TWIN_LAW print theirs), and the two differ
+                // for every law that is not symmetric - a pseudo-merohedral one of a centred lattice,
+                // say, whose real-space triplet relates unrelated reflections when read as hkl.
+                const std::string t = op.as_hkl().triplet();
                 laws += (laws.empty() ? "" : " ; ") + t;
                 // What each law reads on the merge: the disagreement H of the reflections it
                 // relates (Yeates (1997) Methods Enzymol. 276, 344-358: a twin of fraction a reads
@@ -4649,16 +4656,27 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                 law_opts.nthreads = static_cast<size_t>(std::max(1, config_.nthreads));
                 if (!law_correlation)
                     law_correlation.emplace(expanded_merge(), law_opts);
-                if (const auto s = (*law_correlation)(m); s.has_value())
+                if (const auto s = (*law_correlation)(m); s.has_value()) {
                     result.twin_law_evidence.push_back(fmt::format(
                             "{}: {} pairs, H {:.3f} (twin fraction {:.2f} if only twinning makes the pairs "
                             "agree; untwinned crystals read up to about 0.25), R {:.3f}, CC {:.3f}",
                             t, s->n_pairs, s->h_stat, std::max(0.0, 0.5 - s->h_stat), s->r_stat, s->cc));
-                else
+                    if (s->h_stat < best_law_h) {
+                        best_law_h = s->h_stat;
+                        best_law = t;
+                    }
+                } else
                     result.twin_law_evidence.push_back(fmt::format("{}: too few pairs to read", t));
             }
             for (const auto &e : result.twin_law_evidence)
                 logger.Info("Twin law {}", e);
+            // The twin law a refusal identified (above) is the one of these laws whose pairs agree
+            // best - a twin law makes its own pairs agree, the lattice's other laws relate unrelated
+            // intensities - at the fraction its own H implies, the number TWIN_LAW_n prints for it.
+            if (!result.twinning.twin_law_operator.empty() && !best_law.empty()) {
+                result.twinning.twin_law_operator = best_law;
+                result.twinning.twin_law_fraction = std::max(0.0, 0.5 - best_law_h);
+            }
             // A model can resolve it too, but whether it does is only known once the model has
             // been fitted (validation step below), so the warning stands here exactly as in a run
             // without a model, and is withdrawn there if the model's indexing probe decided it.