test: name why orix Misorientation.reduce misses the smallest angle

Sources/DiffractoOrientation/Misorientation.swift

@@ -13,9 +13,14 @@
 //
 //  DEVIATION (orix `Misorientation.reduce`): orix returns the equivalent inside its misorientation fundamental zone, which is documented as
 //  "the smallest rotation angle", and is not for every pair of groups: in the pins, 1 of 41 draws of (m-3m, 6/mmm) and 3 of 41 of (m-3, -3)
-//  come back with a larger angle than another equivalent (orix 2.18 against 0.92 rad on one draw). The kit returns the smallest angle
-//  (brute force over the equivalents, ties by the canonical quaternion of `CrystalOrientationSymmetry.reduce`); the pins compare the kit with
-//  orix where orix is minimal and with the minimum where it is not.
+//  come back with a larger angle than another equivalent (orix 2.18 against 0.92 rad on one draw). The cause: orix's loop takes
+//  `end * M * start` and stops at the first one inside `OrientationRegion.from_symmetry(start, end)`, whose distinguished points are
+//  start . end, while the smallest angle of `end * M * start` lies in the Voronoi cell of the identity among end . start, the inverse set.
+//  The two agree when end . start = start . end (the product set is a group, for example when one group normalises the other) and differ
+//  otherwise (432 with 622, 23 with 3): the zone then holds no equivalent on some draws and the loop falls through to its last pair, or it stops at a larger angle.
+//  The same orix loop with the zone built as `from_symmetry(end, start)` returns the smallest angle on every draw (asserted by the pin script).
+//  The kit returns the smallest angle (brute force over the equivalents, ties by the canonical quaternion of `CrystalOrientationSymmetry.reduce`);
+//  the pins compare the kit with orix's loop on the ordered zone on every draw.
 //
 
 import Foundation

Tests/DiffractoKitParityTests/Orientation/OrixMisorientationTests.swift

@@ -6,13 +6,15 @@ import DiffractoCore
 import DiffractoOrientation
 
 /// Pins: tools/pins/orix/misorientation-pins.py (orix 0.15.0). Tolerances: angles atol 1e-4 (orix's test_misorientation tolerance; measured 1e-15),
-/// equivalent sets and reduced matrices atol 1e-9 where the smallest angle is unique, the Euler angles atol 1e-9. Where orix's `reduce` is not
-/// the smallest angle (documented in Misorientation.swift) the kit is compared with the minimum, not with orix.
+/// equivalent sets and reduced matrices atol 1e-9 where the smallest angle is unique, the Euler angles atol 1e-9. Every draw is compared with
+/// orix's `reduce` loop run with its zone built as `from_symmetry(end, start)` (`orderedZone*`, the cause in Misorientation.swift); where orix's
+/// shipped `reduce` is not the smallest angle, the test asserts that it is larger than the minimum.
 private struct Pins: Decodable {
     struct Item: Decodable {
         let matrix, reduced: [[Double]]
         let equivalentCount: Int
-        let angle, minAngle: Double
+        let angle, minAngle, orderedZoneAngle: Double
+        let orderedZoneReduced: [[Double]]
         let orixIsMinimal, unique: Bool
         let equivalents: [[[Double]]]?
         let grainExchangeCount: Int?
@@ -64,6 +66,12 @@ struct OrixMisorientationTests {
                 let r = m.reduce()
                 worstAngle = max(worstAngle, abs(r.angleRad - item.minAngle))
                 #expect(abs(r.angleRad - item.minAngle) < 1e-9, Comment(rawValue: "\(p.source)>\(p.target) minimum angle"))
+                #expect(abs(r.angleRad - item.orderedZoneAngle) < 1e-4, Comment(rawValue: "\(p.source)>\(p.target) orix ordered-zone angle"))
+                if item.unique {
+                    let g = gap(r.matrix, mat(item.orderedZoneReduced))
+                    worstMatrix = max(worstMatrix, g)
+                    #expect(g < 1e-9, Comment(rawValue: "\(p.source)>\(p.target) orix ordered-zone reduced matrix"))
+                }
                 if item.orixIsMinimal {
                     #expect(abs(r.angleRad - item.angle) < 1e-4, Comment(rawValue: "\(p.source)>\(p.target) orix angle"))
                     if item.unique {

docs/open-items.md

@@ -32,6 +32,8 @@ A plain list, worst first. Each line: the item, the cause where it is known, wha
 - **Ring fit.** `fitRing` at R 600 px on 2048² takes tens of seconds (a coarse-then-fine annulus, or Metal). A one-ring profile
   model that tolerates a neighbouring ring is missing (the {111} bias of `RingFitNeighbourTailTests` is the model's, not a bug).
 - **orix API gaps.** `Orientation.from_euler(symmetry:)`, `dot`, `dot_outer`, `get_distance_matrix` and the orientation
-  (Rodrigues) fundamental zone. orix's `reduce()` disagrees with the smallest angle (parity P36); not reported upstream yet.
+  (Rodrigues) fundamental zone. orix's `Misorientation.reduce()` misses the smallest angle when the two groups' product set is
+  not symmetric (432 with 622, 23 with 3): it builds its zone from start . end while its equivalents need end . start (cause and
+  pins in `Misorientation.swift`); the kit returns the smallest angle. Not reported upstream yet.
 - **Not yet used by a known consumer** (0.10.0): `PatternCentre.find`, `PatternRadialProfile.about`, `LinearScale.fit` and
   `Misorientation`. Revisit at 0.12.

tools/pins/orix/misorientation-pins.py

@@ -7,6 +7,13 @@
     The kit's `Misorientation(matrix, source: Gl, target: Gr)` is the same object: the matrix maps source-crystal vectors to
     target-crystal vectors, the equivalents are g_target R g_source. Asserted here: the matrix of each orix equivalent equals the
     product of the matrices of its factors.
+  * orix's `reduce()` loop with its zone built as `OrientationRegion.from_symmetry(end, start)` (`orderedZone*` fields). orix's
+    `reduce()` takes `end * M * start` for every symmetry pair and stops at the first one inside `from_symmetry(start, end)`, whose
+    distinguished points are `start.outer(end)` = start . end. The smallest angle of `end * M * start` lies in the Voronoi cell of
+    the identity among end . start, the inverse set. The two sets are the same when end . start = start . end (the product set is
+    a group, for example when one group normalises the other); otherwise orix's zone is the inverse of the right one, holds no equivalent on some draws,
+    and the loop falls through to the last pair (or stops at a larger angle inside the wrong zone). The same loop with
+    `OrientationRegion.from_symmetry(end, start)` returns the smallest angle on every draw (asserted here); the kit is compared with it.
   * `Orientation.in_euler_fundamental_region()` for the eleven Laue classes (and -3m in both settings).
 The symmetry groups are built from the kit's own generators (CrystalOrientationSymmetry.init, copied below), so orix and the kit
 share the settings: unique axis b for 2/m, hexagonal a = x, b at 120 degrees, -31m (twofold perpendicular to a) by default.
@@ -21,6 +28,7 @@
 from pathlib import Path
 import numpy as np
 from orix.quaternion import Misorientation, Orientation, Rotation, Symmetry, symmetry as S
+from orix.quaternion.orientation_region import OrientationRegion
 from scipy.spatial.transform import Rotation as SR
 
 KIT = Path(__file__).resolve().parents[3]
@@ -76,6 +84,16 @@ def ang(R): return float(np.arccos(np.clip((np.trace(R) - 1) / 2, -1, 1)))
 
 def mats(r): return r.to_matrix().reshape(-1, 3, 3)
 
+def reduce_ordered_zone(m, fz):
+    """orix 0.15.0 `Misorientation.reduce()` line for line (same loop order, same `<` test), with `fz` built from (end, start)."""
+    start, end = m._symmetry
+    for sym_start in start:
+        for sym_end in end:
+            reduced = sym_end * m * sym_start
+            if reduced < fz:
+                return reduced
+    raise AssertionError('no equivalent inside the ordered zone')
+
 def main():
     rng = np.random.default_rng(SEED)
     out = {'orix': md.version('orix'), 'seed': SEED, 'pairs': [], 'euler': [], 'beta': None}
@@ -88,6 +106,7 @@ def main():
         gl = sa.proper_subgroup
         gr = sb.proper_subgroup
         la = mats(gl); lb = mats(gr)
+        ordered = OrientationRegion.from_symmetry(gr, gl)
         rec = {'source': a, 'target': b, 'sourceOrder': len(la), 'targetOrder': len(lb), 'items': []}
         Ms = list(SR.random(40, random_state=rng).as_matrix()) + [np.eye(3)]
         keepSets = 2          # only the first few items carry the full list of equivalents (file size)
@@ -107,8 +126,11 @@ def main():
             gap = [x for x in np.diff(angles) if x > 1e-9]
             unique = bool(len(angles) < 2 or (angles[1] - angles[0]) > 1e-6)
             orix_min = abs(float(red.angle[0]) - angles[0]) < 1e-8
+            oz = reduce_ordered_zone(m, ordered)
+            assert abs(float(oz.angle[0]) - angles[0]) < 1e-8, (a, b, index)
             item = {'matrix': M.tolist(), 'equivalentCount': len(uniq), 'angle': float(red.angle[0]), 'minAngle': angles[0], 'orixIsMinimal': bool(orix_min),
-                    'unique': unique, 'reduced': mats(red)[0].tolist()}
+                    'unique': unique, 'reduced': mats(red)[0].tolist(),
+                    'orderedZoneAngle': float(oz.angle[0]), 'orderedZoneReduced': mats(oz)[0].tolist()}
             if index < keepSets or index == len(Ms) - 1:
                 item['equivalents'] = [u.tolist() for u in uniq]
             if a == b:
@@ -148,6 +170,7 @@ def main():
     out['beta'] = {'R0': R0.tolist(), 'equivalentCount': len(uniq), 'reducedAngle': float(red.angle[0]),
                    'reducedAngleClosedForm': float(np.arctan2(2, 3))}
     (FIX / 'orix-misorientation-pins.json').write_text(json.dumps(out, sort_keys=True, separators=(',', ':'), allow_nan=False) + '\n')
+    print('orix reduce() with its zone in the order of its equivalents is the minimum angle on every draw')
     print('orix reduce() is not the minimum angle in', {f"{p['source']}>{p['target']}": sum(not i['orixIsMinimal'] for i in p['items']) for p in out['pairs']})
     print('ok', [(p['source'], p['target'], p['items'][0]['equivalentCount']) for p in out['pairs']],
           [(e['symbol'], e['orixProperName'], e['region']) for e in out['euler']], out['beta']['equivalentCount'], out['beta']['reducedAngle'])