feat: Orientation from Euler angles, dot, dotOuter and distanceMatrix

Sources/DiffractoOrientation/Misorientation.swift

@@ -147,6 +147,75 @@ public struct Orientation: Sendable {
         self.matrix = matrix; self.symmetry = symmetry
     }
 
+    // MARK: orix from_euler, dot, dot_outer, get_distance_matrix
+
+    /// Which way Bunge Euler angles turn (orix `from_euler(direction=)`): `labToCrystal` is the Bunge convention and orix's default,
+    /// `crystalToLab` the inverse (MTEX's).
+    public enum EulerDirection: Sendable, Equatable { case labToCrystal, crystalToLab }
+
+    /// orix `Orientation.from_euler(euler, symmetry, direction)`: Bunge (phi1, Phi, phi2) in radians, through orix's `eu2qu`
+    /// (Rowenhorst 2015, Eqs. A.5 and A.6), the inverse of `inEulerFundamentalRegion()`'s conversion. Throws `invalidRotation` for a
+    /// non-finite angle.
+    public init(euler: SIMD3<Double>, symmetry: CrystalOrientationSymmetry, direction: EulerDirection = .labToCrystal) throws {
+        guard euler.x.isFinite, euler.y.isFinite, euler.z.isFinite else { throw OrientationSymmetryError.invalidRotation }
+        let sigma = 0.5 * (euler.x + euler.z), delta = 0.5 * (euler.x - euler.z)
+        let c = cos(euler.y / 2), s = sin(euler.y / 2)
+        var q = simd_quatd(real: c * cos(sigma), imag: SIMD3(-s * cos(delta), -s * sin(delta), -c * sin(sigma)))
+        if direction == .crystalToLab { q = q.conjugate }
+        self.init(unchecked: simd_double3x3(q), symmetry: symmetry)
+    }
+
+    /// orix `Orientation.dot(other)`: the symmetry-reduced dot product, the largest `|<other * self^-1, s>|` over the proper
+    /// operators s, in [0, 1] (the cosine of half the disorientation angle). Throws `differentSymmetry` unless both share a group
+    /// (as `angle(with:)`; orix takes the product of the two groups, which this kit does not hold).
+    public func dot(_ other: Orientation) throws -> Double {
+        guard symmetry.sameGroup(as: other.symmetry) else { throw OrientationSymmetryError.differentSymmetry }
+        return Orientation.reducedDot(simd_quatd(matrix), simd_quatd(other.matrix), Orientation.operatorQuaternions(symmetry))
+    }
+
+    /// orix `Orientation.dot_outer(other)`: `result[i][j] = a[i].dot(b[j])`. Throws `differentSymmetry` unless every orientation of
+    /// both lists shares one group.
+    public static func dotOuter(_ a: [Orientation], _ b: [Orientation]) throws -> [[Double]] {
+        guard let ops = try sharedOperators(a + b) else { return a.map { _ in [] } }
+        let qa = a.map { simd_quatd($0.matrix) }, qb = b.map { simd_quatd($0.matrix) }
+        return qa.map { x in qb.map { y in reducedDot(x, y, ops) } }
+    }
+
+    /// orix `Orientation.get_distance_matrix()`: the disorientation angle (radians) between every pair, computed as orix does,
+    /// `arccos(2 d^2 - 1)` with the dot product d rounded to 12 decimals and clamped to 1 (so the diagonal is exactly 0). Throws
+    /// `differentSymmetry` unless every orientation shares one group.
+    public static func distanceMatrix(_ orientations: [Orientation]) throws -> [[Double]] {
+        guard let ops = try sharedOperators(orientations) else { return [] }
+        let q = orientations.map { simd_quatd($0.matrix) }
+        return q.map { x in
+            q.map { y in
+                let d = min((reducedDot(x, y, ops) * 1e12).rounded() / 1e12, 1)
+                return acos(2 * d * d - 1)
+            }
+        }
+    }
+
+    /// The proper operators as quaternions when every orientation shares a group; nil for an empty list.
+    private static func sharedOperators(_ orientations: [Orientation]) throws -> [simd_quatd]? {
+        guard let first = orientations.first else { return nil }
+        guard orientations.allSatisfy({ $0.symmetry.sameGroup(as: first.symmetry) }) else {
+            throw OrientationSymmetryError.differentSymmetry
+        }
+        return operatorQuaternions(first.symmetry)
+    }
+
+    private static func operatorQuaternions(_ symmetry: CrystalOrientationSymmetry) -> [simd_quatd] {
+        symmetry.operators.map { simd_quatd($0) }
+    }
+
+    /// orix's kernel: M = b * conj(a), the largest |M . s| over the operators.
+    private static func reducedDot(_ a: simd_quatd, _ b: simd_quatd, _ ops: [simd_quatd]) -> Double {
+        let m = (b * a.conjugate).vector
+        var best = 0.0
+        for s in ops { best = max(best, abs(simd_dot(m, s.vector))) }
+        return best
+    }
+
     // MARK: Euler fundamental region
 
     /// The proper subgroup's name in orix's convention for the kit's standard setting ("121" for 2/m with b unique, "312" for -31m, "32" for -3m1).

Tests/DiffractoOrientationTests/OrixOrientationAPITests.swift

@@ -0,0 +1,103 @@
+import DiffractoKitTestSupport
+import Foundation
+import Testing
+import simd
+import DiffractoCore
+import DiffractoOrientation
+
+/// Pins: tools/pins/orix/orientation-api-pins.py (orix 0.15.0). Tolerances of orix test_orientation.py: `from_euler` atol 1e-4
+/// (`test_from_euler_symmetry`), the rest np.allclose defaults (rtol 1e-5, atol 1e-8) element by element; the measured worst gaps
+/// are printed (below 1e-15 when written).
+private struct Pins: Decodable {
+    struct EulerItem: Decodable { let euler: [Double]; let lab2crystal, crystal2lab: [[Double]] }
+    struct Set: Decodable {
+        let symbol: String
+        let A, B: [[[Double]]]
+        let dot: [Double]
+        let dotOuter, distanceMatrix: [[Double]]
+    }
+    let euler: [EulerItem]
+    let sets: [Set]
+}
+private func mat(_ r: [[Double]]) -> simd_double3x3 { .init(rows: r.map { SIMD3($0[0], $0[1], $0[2]) }) }
+private func gap(_ a: simd_double3x3, _ b: simd_double3x3) -> Double {
+    var g = 0.0
+    for c in 0 ..< 3 { for r in 0 ..< 3 { g = max(g, abs(a[c][r] - b[c][r])) } }
+    return g
+}
+/// np.allclose defaults: |x - want| <= 1e-8 + 1e-5 |want|.
+private func close(_ x: Double, _ want: Double) -> Bool { abs(x - want) <= 1e-8 + 1e-5 * abs(want) }
+private func symmetry(_ symbol: String) -> CrystalOrientationSymmetry {
+    if symbol == "-31m" { return CrystalOrientationSymmetry(laueClass: .threeBarM, setting: .threeBarOneM) }
+    return pinnedSymmetry(symbol)!
+}
+
+struct OrixOrientationAPITests {
+    private func pins() throws -> Pins {
+        let file = URL(fileURLWithPath: #filePath).deletingLastPathComponent().appendingPathComponent("Fixtures/orix-orientation-api-pins.json")
+        return try JSONDecoder().decode(Pins.self, from: Data(contentsOf: file))
+    }
+
+    @Test func fromEulerMatchesOrixInBothDirections() throws {
+        let cubic = CrystalOrientationSymmetry(laueClass: .mThreeBarM)
+        var worst = 0.0
+        for item in try pins().euler {
+            let e = SIMD3(item.euler[0], item.euler[1], item.euler[2])
+            let lab = try Orientation(euler: e, symmetry: cubic)
+            let crystal = try Orientation(euler: e, symmetry: cubic, direction: .crystalToLab)
+            #expect(lab.symmetry.laueClass == .mThreeBarM && crystal.symmetry.laueClass == .mThreeBarM)
+            let g = max(gap(lab.matrix, mat(item.lab2crystal)), gap(crystal.matrix, mat(item.crystal2lab)))
+            worst = max(worst, g)
+            #expect(g < 1e-4, Comment(rawValue: "\(item.euler)"))     // orix test_from_euler_symmetry: atol 1e-4
+        }
+        print("ORIX from_euler worst: \(worst)")
+    }
+
+    @Test func fromEulerRefusesNonFiniteAngles() {
+        let cubic = CrystalOrientationSymmetry(laueClass: .mThreeBarM)
+        #expect(throws: OrientationSymmetryError.invalidRotation) { try Orientation(euler: SIMD3(.nan, 0, 0), symmetry: cubic) }
+        #expect(throws: OrientationSymmetryError.invalidRotation) { try Orientation(euler: SIMD3(0, .infinity, 0), symmetry: cubic) }
+    }
+
+    @Test func dotDotOuterAndDistanceMatrixMatchOrixForEveryClass() throws {
+        var worstDot = 0.0, worstDistance = 0.0
+        for set in try pins().sets {
+            let s = symmetry(set.symbol)
+            let a = try set.A.map { try Orientation(mat($0), symmetry: s) }
+            let b = try set.B.map { try Orientation(mat($0), symmetry: s) }
+            for (i, want) in set.dot.enumerated() {
+                let d = try a[i].dot(b[i])
+                worstDot = max(worstDot, abs(d - want))
+                #expect(close(d, want), Comment(rawValue: "\(set.symbol) dot \(i)"))
+            }
+            let outer = try Orientation.dotOuter(a, b)
+            #expect(outer.count == a.count && outer.allSatisfy { $0.count == b.count })
+            for i in a.indices { for j in b.indices {
+                worstDot = max(worstDot, abs(outer[i][j] - set.dotOuter[i][j]))
+                #expect(close(outer[i][j], set.dotOuter[i][j]), Comment(rawValue: "\(set.symbol) dotOuter \(i) \(j)"))
+            } }
+            let distance = try Orientation.distanceMatrix(a)
+            #expect(distance.count == a.count && distance.allSatisfy { $0.count == a.count })
+            for i in a.indices { for j in a.indices {
+                worstDistance = max(worstDistance, abs(distance[i][j] - set.distanceMatrix[i][j]))
+                #expect(close(distance[i][j], set.distanceMatrix[i][j]), Comment(rawValue: "\(set.symbol) distance \(i) \(j)"))
+            } }
+            // A[1] and B[0] are symmetry equivalents of A[0]: the dot is 1 and the rounded, clamped distance exactly 0.
+            #expect(distance[0][0] == 0 && distance[0][1] == 0 && distance[1][0] == 0, Comment(rawValue: "\(set.symbol) equivalent distance"))
+            let selfDot = try a[0].dot(a[1])
+            #expect(abs(outer[0][0] - 1) < 1e-12 && abs(selfDot - 1) < 1e-12, Comment(rawValue: "\(set.symbol) equivalent dot"))
+            // The distance is the disorientation angle of `angle(with:)`.
+            #expect(close(distance[1][2], try a[1].angle(with: a[2])))
+        }
+        print("ORIX dot / dot_outer worst: \(worstDot), distance matrix worst: \(worstDistance)")
+    }
+
+    @Test func dotRefusesDifferentSymmetries() throws {
+        let cubic = try Orientation(matrix_identity_double3x3, symmetry: CrystalOrientationSymmetry(laueClass: .mThreeBarM))
+        let hexagonal = try Orientation(matrix_identity_double3x3, symmetry: CrystalOrientationSymmetry(laueClass: .sixOverMmm))
+        #expect(throws: OrientationSymmetryError.differentSymmetry) { try cubic.dot(hexagonal) }
+        #expect(throws: OrientationSymmetryError.differentSymmetry) { try Orientation.dotOuter([cubic], [hexagonal]) }
+        #expect(throws: OrientationSymmetryError.differentSymmetry) { try Orientation.distanceMatrix([cubic, hexagonal]) }
+        #expect(try Orientation.distanceMatrix([]).isEmpty)
+    }
+}

docs/open-items.md

@@ -31,9 +31,10 @@ A plain list, worst first. Each line: the item, the cause where it is known, wha
   find it.
 - **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 `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.
+- **orix API gaps.** The orientation (Rodrigues) fundamental zone, and `dot` / `dot_outer` / `get_distance_matrix` across two
+  different groups (`from_euler`, `dot`, `dot_outer` and `get_distance_matrix` for one group are held). orix's
+  `Misorientation.reduce()` misses the smallest angle when the two groups' product set is not a group (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.

docs/parity.md

@@ -296,8 +296,8 @@ Single-pattern and calibration helpers the apps asked for that no Python package
 | `Rotation.random_vonmises`, `von_mises`, `Rotation.get_random_sample` | none | not held | not planned (random orientation spreads are for plotting and tests) | `test_quaternion/test_rotation.py` (no value test) |  |
 | `Orientation.reduce`, `Orientation.angle_with`, `Misorientation.reduce` (symmetry reduction, disorientation) with the same group on both sides | `CrystalOrientationSymmetry.reduce(_:)`, `disorientation(from:to:)` | **held** | study | np.allclose defaults (rtol 1e-5, atol 1e-8); `atol=1e-3` (`test_set_symmetry`) in `test_quaternion/test_orientation.py` | P27, orix 0.15.0 executed, all 11 classes; `reduce` Float angle and disorientation within orix's allclose (atol 1e-8 + rtol 1e-5 of the angle); orix's `Quaternion.from_matrix` is imprecise near a half-turn (1.3e-5 rad), so the pins build quaternions with SciPy |
 | `Misorientation((s1, s2))` with two different symmetries (phase boundaries, orientation relationships), `equivalent`, `map_into_symmetry_reduced_zone` / `reduce` | `Misorientation` (`equivalent(grainExchange:)`, `reduce()`, 0.10.0), `OrientationVariantSet`, `OrientationRelationship.variants(of:matrix:precipitate:)`, `.rotation(_:_:precipitate:matrix:)` | **held** | study | np.allclose defaults (rtol 1e-5, atol 1e-8); `test_quaternion/test_misorientation.py` (`atol=1e-4`) | P36, orix 0.15.0 executed: 12 group pairs (m-3m with 2/m, 6/mmm and m-3m; 2/m with m-3m; 6/mmm with 2/m; 4/mmm with mmm; m-3 with -3; -31m with -3m1; -1 with m-3m; mmm; 6/m with 4/m; -3m1 with 6/mmm), 41 draws each: the equivalents as sets (6e-16), the distinct count, the smallest angle (4e-15), the reduced matrix where it is unique, the grain-exchange count. **DEVIATION**: orix's `reduce()` is the equivalent inside its misorientation fundamental zone, which is documented as the smallest angle and is not for every pair: it returns a larger angle on 1 of 41 draws of (m-3m, 6/mmm) and 3 of 41 of (m-3, -3) (2.18 against 0.92 rad on one); the kit returns the smallest angle (compared with orix where orix is minimal, with brute force where it is not). beta'' in Al: 24 distinct equivalents, reduced angle arctan(2/3) = 33.690068 degrees, 12 variants (planted truth, `OrientationVariantsTests`: four needle variants per <001> axis, variant-to-variant angles 0/90/120/180, nearest variant of 48 noisy orientations exact); no Python reference computes the variant set |
-| `Orientation` (`from_euler(symmetry=)`, `from_matrix`, `from_axes_angles`, `symmetry`, `unit`, `in_euler_fundamental_region`, `dot`, `dot_outer`) | `Orientation(_:symmetry:)` with `reduce()`, `angle(with:)`, `inEulerFundamentalRegion()` (0.10.0; `OrientationResult` still holds Euler angles and a score) | **partly** | study | np.allclose defaults (rtol 1e-5, atol 1e-8); `atol=1e-4` (`test_from_euler_symmetry`) in `test_quaternion/test_orientation.py` | P36, orix 0.15.0 executed: `in_euler_fundamental_region` for the 11 classes and -3m in both settings (kit settings, 61 draws each; orix's proper-subgroup names 1, 121, 222, 4, 422, 3, 312, 32, 6, 622, 23, 432), Euler angles within 2e-15. Not held: `from_euler(symmetry=)`, `from_matrix`, `from_axes_angles`, `dot`, `dot_outer` |
-| `Orientation.get_distance_matrix`, `Misorientation.get_distance_matrix`, `angle_with_outer` | none | not held | useful | `atol=1e-4` in `test_quaternion/test_misorientation.py`; `atol=1e-6` in `test_quaternion/test_orientation.py` |  |
+| `Orientation` (`from_euler(symmetry=, direction=)`, `from_matrix`, `from_axes_angles`, `symmetry`, `unit`, `in_euler_fundamental_region`, `dot`, `dot_outer`, `get_distance_matrix`) | `Orientation(_:symmetry:)` with `reduce()`, `angle(with:)`, `inEulerFundamentalRegion()` (0.10.0; `OrientationResult` still holds Euler angles and a score); `Orientation(euler:symmetry:direction:)`, `dot(_:)`, `Orientation.dotOuter(_:_:)`, `Orientation.distanceMatrix(_:)` | **partly** | study | np.allclose defaults (rtol 1e-5, atol 1e-8); `atol=1e-4` (`test_from_euler_symmetry`) in `test_quaternion/test_orientation.py` | P36, orix 0.15.0 executed: `in_euler_fundamental_region` for the 11 classes and -3m in both settings (kit settings, 61 draws each; orix's proper-subgroup names 1, 121, 222, 4, 422, 3, 312, 32, 6, 622, 23, 432), Euler angles within 2e-15. `from_euler` (36 Bunge triplets with the gimbal-lock and half-turn edges, both directions; matrix within 3.3e-16, held at orix's atol 1e-4), `dot`, `dot_outer` and `get_distance_matrix` (the 12 class/setting entries, six and five orientations each, two of them symmetry equivalents of the first: dot 1, distance exactly 0; within 8.9e-16; the distance matrix rounds the dot product to 12 decimals as orix's kernel does) held at np.allclose defaults (`OrixOrientationAPITests`, `tools/pins/orix/orientation-api-pins.py`). `dot`, `dotOuter` and `distanceMatrix` refuse two different groups (`differentSymmetry`, as `angle(with:)`): orix's two-group `dot` searches `sym1.outer(sym2)`, which with equivalents `S * M` is the inverse of the set it needs (the cause of the `Misorientation.reduce` DEVIATION in the row above), and for 432 with 622 or 23 with 3 it differs from the largest dot over equivalent pairs (by up to 0.08 in a 200-draw check); not held. Not held: `from_matrix` beyond the matrix initialiser, `from_axes_angles` |
+| `Misorientation.get_distance_matrix`, `angle_with_outer` (`Orientation.get_distance_matrix` is held, row above) | none | not held | useful | `atol=1e-4` in `test_quaternion/test_misorientation.py`; `atol=1e-6` in `test_quaternion/test_orientation.py` |  |
 | `Symmetry` as a set of 32 point groups (`order`, `is_proper`, `subgroups`, `proper_subgroup`, `contains_inversion`, `diads`, `from_generators`, `system`) and the 11 Laue groups (`laue`, `laue_proper_subgroup`) | `PointSymmetry`, `LaueClass` (the 11 Laue classes as integer rotations acting on hkl), `CrystalOrientationSymmetry.operators` (proper operators as matrices) | **partly** | study | equality of operator sets, no numeric tolerance (`test_quaternion/test_symmetry.py`) | operation counts and rotation sets per class vs pymatgen/spglib (P2, P3, P17), and vs orix as asserted operator sets inside the P26/P27/P32 generators (the Swift tests compare the kit with the resulting numbers, not with the operator sets); the 21 non-centrosymmetric point groups, `subgroups` and `from_generators` are not held |
 | `get_point_group(space_group_number, proper=False)` | none (`PhaseModel.orientationClass(of:)` gives the Laue class of a detected structure) | not held | study | exact equality for all 230 numbers (`test_quaternion/test_symmetry.py::test_get_point_group`) | needs the space group number (ROADMAP item 2); `PhaseModel.detectedOrientationClassIsSpglibsLaueClassOrNil` (P29) checks 15 groups against spglib's Laue class |
 | `Symmetry.fundamental_sector` (`FundamentalSector`: `vertices`, `center`, `edges`), `fundamental_zone`, `euler_fundamental_region` | internal sector tables generated from orix (`OrixIPFSectors.swift`, P26); `CrystalOrientationSymmetry.sampleFundamentalZone(count:)` for the directions | **partly** | study | `atol=1e-3`, `1e-4` for the 32 sectors in `test_quaternion/test_symmetry.py`; `atol=1e-3` in `test_vector3d/test_fundamental_sector.py` | the sectors enter only through the colour key (P26) and the sector-corner pins; they are not a public type and `euler_fundamental_region` is held as a private table behind `Orientation.inEulerFundamentalRegion()` (P36) |