From 42a2af88260cf28ce9986ce799c04b4b531f3ec1 Mon Sep 17 00:00:00 2001 From: Benjamin Date: Fri, 27 Feb 2026 21:21:16 +0000 Subject: [PATCH 1/6] Replace N! permutation loop in RotationAnalyser with symmetry-reduced approach Use bsym to reduce vertex permutations to symmetry-inequivalent representatives, then compose the best-fit Procrustes rotation with the point group's proper rotation matrices to recover all orientation candidates. A single PointGroupAnalyzer call provides both the reduced permutations and the 3x3 rotation matrices. Gives ~50x speedup for octahedral geometries (15 CSM evaluations instead of 720). The composition is exact for undistorted geometries and numerically indistinguishable from brute force for moderate distortions. The docstring documents an orbit-expansion alternative for cases requiring exact rotation matrices. Also fixes the OrientationDict type annotation for all_rotational_distances (np.ndarray, not float). --- polyhedral_analysis/rotation_analyser.py | 252 +++++++++++++++++------ tests/test_rotation_analyser.py | 24 +++ 2 files changed, 215 insertions(+), 61 deletions(-) diff --git a/polyhedral_analysis/rotation_analyser.py b/polyhedral_analysis/rotation_analyser.py index 20b6fa1..7c5abc3 100644 --- a/polyhedral_analysis/rotation_analyser.py +++ b/polyhedral_analysis/rotation_analyser.py @@ -1,23 +1,89 @@ +from __future__ import annotations + from polyhedral_analysis.csm import continuous_symmetry_measure -from itertools import permutations -import math import numpy as np from typing import TypedDict from polyhedral_analysis.coordination_polyhedron import CoordinationPolyhedron +def _analyse_point_group(reference_points: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Analyse the point group of a reference geometry. + + Performs a single symmetry analysis and returns both the + symmetry-reduced permutation indices and the proper rotation matrices. + + Args: + reference_points: An Nx3 array of reference vertex coordinates. + + Returns: + A tuple of (reduced_permutations, proper_rotations) where: + + - reduced_permutations is an MxN integer array of symmetry-inequivalent + permutation representatives. + - proper_rotations is a Kx3x3 array of proper (det > 0) rotation + matrices of the point group. + """ + from pymatgen.core import Molecule + from pymatgen.symmetry.analyzer import PointGroupAnalyzer + from pymatgen.util.coord import coord_list_mapping + from bsym import ConfigurationSpace # type: ignore[import] + from bsym import PointGroup as BsymPointGroup # type: ignore[import] + from bsym.symmetry_operation import SymmetryOperation # type: ignore[import] + + n = len(reference_points) + mol = Molecule(['H'] * n, reference_points - reference_points.mean(axis=0)) + pga = PointGroupAnalyzer(mol) + symm_ops = pga.get_symmetry_operations() + + # Extract proper 3x3 rotation matrices + proper_rotations = np.array([ + op.rotation_matrix for op in symm_ops + if np.linalg.det(op.rotation_matrix) > 0 + ]) + + # Convert symmetry operations to permutation representations for bsym + coords = mol.cart_coords + seen: set[tuple[int, ...]] = set() + mappings: list[list[int]] = [] + for op in symm_ops: + mapped_coords = op.operate_multi(coords) + new_mol = Molecule(mol.species, mapped_coords) + mapping = [x + 1 for x in coord_list_mapping( + new_mol.cart_coords, coords)] + key = tuple(mapping) + if key not in seen: + seen.add(key) + mappings.append(mapping) + + sym_ops = [SymmetryOperation.from_vector(m) for m in mappings] + point_group = BsymPointGroup(symmetry_operations=sym_ops) + config_space = ConfigurationSpace( + objects=list(range(n)), + symmetry_group=point_group, + ) + unique_configs = config_space.unique_configurations( + site_distribution={i: 1 for i in range(n)}, + ) + reduced_permutations = np.array( + [c.tolist() for c in unique_configs], + dtype=np.intp, + ) + + return reduced_permutations, proper_rotations + + class OrientationDict(TypedDict): orientation_index: int reference_geometry_index: int rotational_distance: float symmetry_measure: float - all_rotational_distances: float + all_rotational_distances: np.ndarray -class RotationAnalyser(object): +class RotationAnalyser: """Class for analysing rotational orientation of polyhedra. Attributes: - reference_points (np.array): A Mx3xN numpy array of points around the + reference_points (np.ndarray): A Mx3xN numpy array of points around the origin ([0.0, 0.0, 0.0]) that define sets of centre-vertex vectors for specific reference orientations, e.g. for a single orientation of a tetrahedron, reference points will be a size 1x3x4 @@ -35,16 +101,13 @@ def __init__(self, """Initialise a RotationAnalyser instance. Args: - reference_points (np.array): Either a 3xN or Mx3xN numpy array - of points around the origin ([0.0, 0.0, 0.0]) that define sets of - centre-vertex vectors, and are used to classify the rotational - orientation of polyhedra. - - If a 3xN array is passed in, this will be converted to a 1x3xN array. - - Returns: - None + reference_points: Either a 3xN or Mx3xN numpy array of points + around the origin ([0.0, 0.0, 0.0]) that define sets of + centre-vertex vectors, and are used to classify the rotational + orientation of polyhedra. + If a 3xN array is passed in, this will be converted to a + 1x3xN array. """ if len(reference_points.shape) == 2: reference_points = np.array([reference_points]) @@ -52,69 +115,136 @@ def __init__(self, raise ValueError( 'Reference points must all contain the same number of coordinates') self.reference_points = reference_points + self._reduced_permutations: np.ndarray | None = None + self._proper_rotations: np.ndarray | None = None + + def _ensure_symmetry_analysis(self) -> None: + """Run point group analysis if not already done. + + Populates both ``_reduced_permutations`` and ``_proper_rotations`` + from a single symmetry analysis. + """ + if self._reduced_permutations is None: + self._reduced_permutations, self._proper_rotations = ( + _analyse_point_group(self.reference_points[0])) + + @property + def reduced_permutations(self) -> np.ndarray: + """Symmetry-inequivalent permutation indices for the reference geometry. + + Computed lazily on first access. + """ + self._ensure_symmetry_analysis() + return self._reduced_permutations # type: ignore[return-value] + + @property + def proper_rotations(self) -> np.ndarray: + """Proper rotation matrices (Kx3x3) of the reference geometry's point group. + + Computed lazily on first access. + """ + self._ensure_symmetry_analysis() + return self._proper_rotations # type: ignore[return-value] def discrete_orientation(self, points: np.ndarray) -> OrientationDict: """Find the discrete "closest orientation" for an input polyhedron of points. - For example, a tetrahedon has 12 pure rotation symmetry operations. A distorted - tetrahedron with vertices approximately aligned along the unordered vectors:: + For example, a tetrahedron has 12 proper rotation symmetry operations. + A distorted tetrahedron with vertices approximately aligned along the + unordered vectors:: [[ 1.0, -1.0, 1.0], [-1.0, -1.0, -1.0], [ 1.0, 1.0, -1.0], [-1.0, 1.0, 1.0]] - can be in one of 12 orientations. This method compares the input points to - all permutations of self.reference_points that can be generated by proper rotations, and returns - an index for orientation that minimises the rotational distance between the reordered - reference points and the input points. - - The algorithm is: - 1. Perform continuous symmetry measure analysis for all permutations of reference points. - 2. Collect operations that have the minimum continous symmetry measure value. - These give the set of all propert and improper rotations. - 3. Find all proper rotations using det(M_rot)>0. - 4. Calculate the rotational distance. + can be in one of 12 orientations. This method compares the input + points to the reference orientations and returns the orientation that + minimises the rotational distance. + + The algorithm uses a two-phase approach to avoid the N! cost of + evaluating all vertex permutations: + + 1. Find the best vertex alignment using symmetry-reduced permutations + of the reference points (N!/|G| CSM evaluations instead of N!). + 2. Compose the best-fit Procrustes rotation R_best with all proper + rotation matrices S_i of the reference geometry's point group to + give the rotation matrix for each orientation: R_i = S_i @ R_best. + 3. Calculate rotational distances from the composed matrices. + 4. Return the orientation with the minimum rotational distance. + + The composed rotation matrices in step 2 are exact for undistorted + geometries. For distorted geometries they are approximate, since the + true Procrustes rotation for each symmetry-equivalent permutation + differs slightly from the composed matrix. In practice the results + are numerically indistinguishable from the brute-force approach for + moderate distortions. + + If exact rotation matrices are needed for highly distorted + geometries, an alternative is a two-phase orbit-expansion approach: + use reduced permutations to identify the minimum-CSM orbit (phase 1), + then expand that orbit to all |G| equivalent permutations and compute + CSM for each to obtain exact rotation matrices (phase 2). This costs + N!/|G| + |G| CSM evaluations per reference orientation rather than + N!/|G|, but avoids the composition approximation. Args: - points (np.array(float)): Nx3 numpy array of points describing the coordinates of the - input polyhedron, centered around zero. + points: Nx3 numpy array of points describing the coordinates of + the input polyhedron, centred around zero. Returns: - (dict): Dictionary describing the orientation, with keys: - - ``orientation_index`` (int): Index of this particular orientation. - - ``reference_geometry_index`` (int): Index of the reference geometry - the closest discrete orientation is equivalent to. - - ``rotational_distance`` (float): Angle of rotation from the relevant - reference orientation, in radians. - - ``all_rotational_distances`` (np.array(float)): Array of angles of rotation - from each reference orientation. - - ``symmetry_measure`` (float): The continuous symmetry measure (CSM) - for this polyhedron (J. Am. Chem. Soc., 114, 7843-7851 (1992)) - + Dictionary describing the orientation, with keys: + + - ``orientation_index`` (int): Index of this particular orientation. + - ``reference_geometry_index`` (int): Index of the reference geometry + the closest discrete orientation is equivalent to. + - ``rotational_distance`` (float): Angle of rotation from the + relevant reference orientation, in radians. + - ``all_rotational_distances`` (np.ndarray): Array of angles of + rotation from each reference orientation. + - ``symmetry_measure`` (float): The continuous symmetry measure + (CSM) for this polyhedron. """ - points -= np.mean(points, axis=0, dtype=float) - sm = [] - for rp in self.reference_points: - sm.extend([continuous_symmetry_measure(points, np.array(p)) for p in permutations(rp)]) - min_sm = min(s.symmetry_measure for s in sm) - pure_rot_sm = [s for s in sm if math.isclose( - s.symmetry_measure, min_sm)] - proper_rot_sm = [s for s in pure_rot_sm if np.linalg.det( - s.rotation_matrix) > 0] - proper_rot_matrices = np.array( - [s.rotation_matrix for s in proper_rot_sm]) - trace = np.trace(proper_rot_matrices, axis1=1, axis2=2) - rot_distance = np.arccos(np.clip((trace - 1.0) / 2.0, -1.0, 1.0)) - index = int(np.argmin(rot_distance)) - reference_geometry_index = index // int( - len(proper_rot_sm) / len(self.reference_points)) - return OrientationDict(orientation_index=index, - reference_geometry_index=reference_geometry_index, - rotational_distance=rot_distance[index], - symmetry_measure=proper_rot_sm[index].symmetry_measure, - all_rotational_distances=rot_distance) + points = points - np.mean(points, axis=0, dtype=float) + # Phase 1: find best alignment using reduced permutations + best_csm = float('inf') + best_rotation = np.eye(3) + best_ref_index = 0 + for i, rp in enumerate(self.reference_points): + for perm in self.reduced_permutations: + result = continuous_symmetry_measure(points, rp[perm]) + if result.symmetry_measure < best_csm: + best_csm = result.symmetry_measure + best_rotation = result.rotation_matrix + best_ref_index = i + # Phase 2: compose best-fit rotation with proper symmetry operations + n_proper = len(self.proper_rotations) + all_rot_distances = np.empty(len(self.reference_points) * n_proper) + for i, rp in enumerate(self.reference_points): + # Find best rotation for this specific reference orientation + if i == best_ref_index: + R_i = best_rotation + else: + R_i = min( + (continuous_symmetry_measure(points, rp[perm]) + for perm in self.reduced_permutations), + key=lambda r: r.symmetry_measure, + ).rotation_matrix + composed = self.proper_rotations @ R_i + traces = np.trace(composed, axis1=1, axis2=2) + start = i * n_proper + all_rot_distances[start:start + n_proper] = np.arccos( + np.clip((traces - 1.0) / 2.0, -1.0, 1.0)) + index = int(np.argmin(all_rot_distances)) + reference_geometry_index = index // n_proper + return OrientationDict( + orientation_index=index, + reference_geometry_index=reference_geometry_index, + rotational_distance=all_rot_distances[index], + symmetry_measure=best_csm, + all_rotational_distances=all_rot_distances, + ) def polyhedron_orientation(self, polyhedron: CoordinationPolyhedron) -> OrientationDict: diff --git a/tests/test_rotation_analyser.py b/tests/test_rotation_analyser.py index ed38c5c..a53b696 100644 --- a/tests/test_rotation_analyser.py +++ b/tests/test_rotation_analyser.py @@ -52,6 +52,30 @@ def test_polyhedron_orientation(self): self.ra.discrete_orientation.assert_called_once_with('foo') self.assertEqual(orientation, 'bar') + def test_discrete_orientation_identifies_correct_reference(self): + """A slightly distorted copy of reference orientation 0 should be identified as such.""" + rng = np.random.RandomState(42) + distorted = self.ra.reference_points[0] + rng.normal(0, 0.05, self.ra.reference_points[0].shape) + result = self.ra.discrete_orientation(distorted) + self.assertEqual(result['reference_geometry_index'], 0) + self.assertAlmostEqual(result['rotational_distance'], 0.0, places=1) + self.assertLess(result['symmetry_measure'], 1.0) + + def test_discrete_orientation_identifies_inverted_reference(self): + """A slightly distorted copy of reference orientation 1 should be identified as such.""" + rng = np.random.RandomState(123) + distorted = self.ra.reference_points[1] + rng.normal(0, 0.05, self.ra.reference_points[1].shape) + result = self.ra.discrete_orientation(distorted) + self.assertEqual(result['reference_geometry_index'], 1) + self.assertAlmostEqual(result['rotational_distance'], 0.0, places=1) + + def test_discrete_orientation_perfect_geometry_has_zero_csm(self): + """A perfect (undistorted) reference geometry should have CSM close to zero.""" + points = self.ra.reference_points[0].copy() + result = self.ra.discrete_orientation(points) + self.assertAlmostEqual(result['symmetry_measure'], 0.0, places=5) + self.assertAlmostEqual(result['rotational_distance'], 0.0, places=5) + if __name__ == '__main__': unittest.main() From ff5b606a08ba48816eb4973bde3ce67c90cb91a9 Mon Sep 17 00:00:00 2001 From: Benjamin Date: Fri, 27 Feb 2026 22:14:26 +0000 Subject: [PATCH 2/6] Update changelog with rotation analyser optimisation and cleanup --- CHANGELOG.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ef3e8fc..a2b9963 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,15 @@ - Cached `symmetry_measure` property on `CoordinationPolyhedron` to avoid redundant CSM calculations on repeated access. - `edge_graph` and `symmetry_measure` now return defensive copies to prevent callers from corrupting the internal cache. +- Replaced N! permutation loop in `RotationAnalyser` with symmetry-reduced approach using bsym, giving ~50x speedup for octahedral geometries. + +### Bug fixes + +- Fixed `OrientationDict.all_rotational_distances` type annotation (`np.ndarray`, not `float`). + +### Other changes + +- Removed redundant `assert isinstance` checks across the codebase. ## 0.4.0 From 2cbf4012941725b2fb791bea0e7f3250bc26c02b Mon Sep 17 00:00:00 2001 From: Benjamin Date: Mon, 2 Mar 2026 07:12:23 +0000 Subject: [PATCH 3/6] Trim README to match actual functionality Remove repeated feature description from intro paragraph and correct "stereographic projection plotting" to "orientation distribution plotting" (the plotting module uses a KDE contour plot, not a stereographic projection). --- README.md | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 40a5d7a..c7e9fc8 100644 --- a/README.md +++ b/README.md @@ -4,9 +4,7 @@ [![Coverage Status](https://coveralls.io/repos/github/bjmorgan/polyhedral-analysis/badge.svg?branch=main)](https://coveralls.io/github/bjmorgan/polyhedral-analysis?branch=main) [![Documentation Status](https://readthedocs.org/projects/polyhedral-analysis/badge/?version=latest)](http://polyhedral-analysis.readthedocs.io/en/latest/?badge=latest) -`polyhedral-analysis` is a Python module for analysing coordination polyhedra in crystal structures and molecular dynamics trajectories. It identifies coordination environments around central atoms and computes geometric properties such as distortion measures, bond lengths, angles, volumes, and connectivity. - -Built on top of [pymatgen](https://pymatgen.org), the package works with any structure that pymatgen can read. +`polyhedral-analysis` is a Python module for analysing coordination polyhedra in crystal structures and molecular dynamics trajectories. Built on top of [pymatgen](https://pymatgen.org), it works with any structure that pymatgen can read. ## Features @@ -16,7 +14,7 @@ Built on top of [pymatgen](https://pymatgen.org), the package works with any str - Bond lengths, angles, volumes, and edge connectivity graphs - Corner-, edge-, and face-sharing neighbour analysis - Off-centre displacement and radial distortion parameters -- Vertex vector orientation analysis with stereographic projection plotting +- Vertex vector orientation analysis with orientation distribution plotting - Trajectory analysis for tracking polyhedral distortions over molecular dynamics runs ## Installation From f3531fab8667636ef01798260214b13377bedce6 Mon Sep 17 00:00:00 2001 From: Benjamin Date: Mon, 2 Mar 2026 07:16:27 +0000 Subject: [PATCH 4/6] Address review comments on RotationAnalyser - Track per-reference CSM so returned symmetry_measure matches the chosen reference geometry, not the global best. - Fix docstring array shapes: Mx3xN -> MxNx3, 3xN -> Nx3, etc. - Document that all reference orientations must share the same point group (symmetry analysis is performed once from the first). - Simplify Phase 2 loop by storing per-reference best rotations in Phase 1. --- polyhedral_analysis/rotation_analyser.py | 52 +++++++++++------------- 1 file changed, 23 insertions(+), 29 deletions(-) diff --git a/polyhedral_analysis/rotation_analyser.py b/polyhedral_analysis/rotation_analyser.py index 7c5abc3..58786e7 100644 --- a/polyhedral_analysis/rotation_analyser.py +++ b/polyhedral_analysis/rotation_analyser.py @@ -82,12 +82,16 @@ class OrientationDict(TypedDict): class RotationAnalyser: """Class for analysing rotational orientation of polyhedra. + All reference orientations must be rotated copies of the same geometry + (i.e. share the same point group). The symmetry analysis is performed + once from the first reference orientation and reused for all others. + Attributes: - reference_points (np.ndarray): A Mx3xN numpy array of points around the - origin ([0.0, 0.0, 0.0]) that define sets of centre-vertex + reference_points (np.ndarray): An MxNx3 numpy array of points around + the origin ([0.0, 0.0, 0.0]) that define sets of centre-vertex vectors for specific reference orientations, e.g. for a single - orientation of a tetrahedron, reference points will be a size 1x3x4 - array, e.g.:: + orientation of a tetrahedron, reference points will be a size + 1x4x3 array, e.g.:: [[[ 1.0, -1.0, 1.0], [-1.0, -1.0, -1.0], @@ -101,13 +105,13 @@ def __init__(self, """Initialise a RotationAnalyser instance. Args: - reference_points: Either a 3xN or Mx3xN numpy array of points + reference_points: Either an Nx3 or MxNx3 numpy array of points around the origin ([0.0, 0.0, 0.0]) that define sets of centre-vertex vectors, and are used to classify the rotational orientation of polyhedra. - If a 3xN array is passed in, this will be converted to a - 1x3xN array. + If an Nx3 array is passed in, this will be converted to a + 1xNx3 array. """ if len(reference_points.shape) == 2: reference_points = np.array([reference_points]) @@ -207,31 +211,21 @@ def discrete_orientation(self, (CSM) for this polyhedron. """ points = points - np.mean(points, axis=0, dtype=float) - # Phase 1: find best alignment using reduced permutations - best_csm = float('inf') - best_rotation = np.eye(3) - best_ref_index = 0 + n_refs = len(self.reference_points) + n_proper = len(self.proper_rotations) + # Phase 1: find best alignment per reference orientation + best_csm_per_ref = np.full(n_refs, float('inf')) + best_rotation_per_ref: list[np.ndarray] = [np.eye(3)] * n_refs for i, rp in enumerate(self.reference_points): for perm in self.reduced_permutations: result = continuous_symmetry_measure(points, rp[perm]) - if result.symmetry_measure < best_csm: - best_csm = result.symmetry_measure - best_rotation = result.rotation_matrix - best_ref_index = i + if result.symmetry_measure < best_csm_per_ref[i]: + best_csm_per_ref[i] = result.symmetry_measure + best_rotation_per_ref[i] = result.rotation_matrix # Phase 2: compose best-fit rotation with proper symmetry operations - n_proper = len(self.proper_rotations) - all_rot_distances = np.empty(len(self.reference_points) * n_proper) - for i, rp in enumerate(self.reference_points): - # Find best rotation for this specific reference orientation - if i == best_ref_index: - R_i = best_rotation - else: - R_i = min( - (continuous_symmetry_measure(points, rp[perm]) - for perm in self.reduced_permutations), - key=lambda r: r.symmetry_measure, - ).rotation_matrix - composed = self.proper_rotations @ R_i + all_rot_distances = np.empty(n_refs * n_proper) + for i in range(n_refs): + composed = self.proper_rotations @ best_rotation_per_ref[i] traces = np.trace(composed, axis1=1, axis2=2) start = i * n_proper all_rot_distances[start:start + n_proper] = np.arccos( @@ -242,7 +236,7 @@ def discrete_orientation(self, orientation_index=index, reference_geometry_index=reference_geometry_index, rotational_distance=all_rot_distances[index], - symmetry_measure=best_csm, + symmetry_measure=float(best_csm_per_ref[reference_geometry_index]), all_rotational_distances=all_rot_distances, ) From 5b60773165ec6041113877173a34471cb51e80d2 Mon Sep 17 00:00:00 2001 From: Benjamin Date: Mon, 2 Mar 2026 07:54:03 +0000 Subject: [PATCH 5/6] Address code review: validation, cache design, and test coverage - Validate coord_list_mapping produces a bijection - Guard against empty proper_rotations array - Use det > 0.5 tolerance for proper rotation filtering - Fix list aliasing in best_rotation_per_ref initialisation - Bundle lazy cache into single _symmetry_data tuple field - Add PEP 8 blank line between OrientationDict and RotationAnalyser - Add tests for _analyse_point_group outputs - Add brute-force equivalence regression test - Add tests for output shape, index consistency, input mutation - Add known-rotation test - Replace deprecated np.random.RandomState with default_rng --- polyhedral_analysis/rotation_analyser.py | 38 ++++++----- tests/test_rotation_analyser.py | 82 +++++++++++++++++++++++- 2 files changed, 101 insertions(+), 19 deletions(-) diff --git a/polyhedral_analysis/rotation_analyser.py b/polyhedral_analysis/rotation_analyser.py index 58786e7..692f529 100644 --- a/polyhedral_analysis/rotation_analyser.py +++ b/polyhedral_analysis/rotation_analyser.py @@ -38,8 +38,11 @@ def _analyse_point_group(reference_points: np.ndarray) -> tuple[np.ndarray, np.n # Extract proper 3x3 rotation matrices proper_rotations = np.array([ op.rotation_matrix for op in symm_ops - if np.linalg.det(op.rotation_matrix) > 0 + if np.linalg.det(op.rotation_matrix) > 0.5 ]) + if len(proper_rotations) == 0: + raise RuntimeError( + 'No proper rotations found for the reference geometry.') # Convert symmetry operations to permutation representations for bsym coords = mol.cart_coords @@ -48,8 +51,12 @@ def _analyse_point_group(reference_points: np.ndarray) -> tuple[np.ndarray, np.n for op in symm_ops: mapped_coords = op.operate_multi(coords) new_mol = Molecule(mol.species, mapped_coords) - mapping = [x + 1 for x in coord_list_mapping( - new_mol.cart_coords, coords)] + mapping_zero = coord_list_mapping(new_mol.cart_coords, coords) + if len(set(mapping_zero)) != n: + raise RuntimeError( + 'coord_list_mapping did not produce a bijection; ' + 'symmetry operation maps multiple vertices to the same target.') + mapping = [x + 1 for x in mapping_zero] key = tuple(mapping) if key not in seen: seen.add(key) @@ -79,6 +86,7 @@ class OrientationDict(TypedDict): symmetry_measure: float all_rotational_distances: np.ndarray + class RotationAnalyser: """Class for analysing rotational orientation of polyhedra. @@ -119,18 +127,18 @@ def __init__(self, raise ValueError( 'Reference points must all contain the same number of coordinates') self.reference_points = reference_points - self._reduced_permutations: np.ndarray | None = None - self._proper_rotations: np.ndarray | None = None + self._symmetry_data: tuple[np.ndarray, np.ndarray] | None = None - def _ensure_symmetry_analysis(self) -> None: + def _ensure_symmetry_analysis(self) -> tuple[np.ndarray, np.ndarray]: """Run point group analysis if not already done. - Populates both ``_reduced_permutations`` and ``_proper_rotations`` - from a single symmetry analysis. + Returns: + A tuple of (reduced_permutations, proper_rotations) from + a single symmetry analysis. """ - if self._reduced_permutations is None: - self._reduced_permutations, self._proper_rotations = ( - _analyse_point_group(self.reference_points[0])) + if self._symmetry_data is None: + self._symmetry_data = _analyse_point_group(self.reference_points[0]) + return self._symmetry_data @property def reduced_permutations(self) -> np.ndarray: @@ -138,8 +146,7 @@ def reduced_permutations(self) -> np.ndarray: Computed lazily on first access. """ - self._ensure_symmetry_analysis() - return self._reduced_permutations # type: ignore[return-value] + return self._ensure_symmetry_analysis()[0] @property def proper_rotations(self) -> np.ndarray: @@ -147,8 +154,7 @@ def proper_rotations(self) -> np.ndarray: Computed lazily on first access. """ - self._ensure_symmetry_analysis() - return self._proper_rotations # type: ignore[return-value] + return self._ensure_symmetry_analysis()[1] def discrete_orientation(self, points: np.ndarray) -> OrientationDict: @@ -215,7 +221,7 @@ def discrete_orientation(self, n_proper = len(self.proper_rotations) # Phase 1: find best alignment per reference orientation best_csm_per_ref = np.full(n_refs, float('inf')) - best_rotation_per_ref: list[np.ndarray] = [np.eye(3)] * n_refs + best_rotation_per_ref: list[np.ndarray] = [np.eye(3) for _ in range(n_refs)] for i, rp in enumerate(self.reference_points): for perm in self.reduced_permutations: result = continuous_symmetry_measure(points, rp[perm]) diff --git a/tests/test_rotation_analyser.py b/tests/test_rotation_analyser.py index a53b696..c11fedc 100644 --- a/tests/test_rotation_analyser.py +++ b/tests/test_rotation_analyser.py @@ -1,6 +1,8 @@ +import itertools import unittest import numpy as np -from polyhedral_analysis.rotation_analyser import RotationAnalyser +from polyhedral_analysis.rotation_analyser import RotationAnalyser, _analyse_point_group +from polyhedral_analysis.csm import continuous_symmetry_measure from unittest.mock import Mock from polyhedral_analysis.coordination_polyhedron import CoordinationPolyhedron @@ -33,6 +35,33 @@ def test_rotation_analyser_init_raises_error_for_misformed_reference_points(self with self.assertRaises(ValueError): ra = RotationAnalyser(reference_points=reference_points) +TETRAHEDRON_VERTICES = np.array([[1.0, -1.0, 1.0], + [-1.0, -1.0, -1.0], + [1.0, 1.0, -1.0], + [-1.0, 1.0, 1.0]]) + + +class TestAnalysePointGroup(unittest.TestCase): + + def test_tetrahedron_reduced_permutation_count(self): + """Td has |G|=24, so 4!/24 = 1 reduced permutation.""" + reduced_perms, _ = _analyse_point_group(TETRAHEDRON_VERTICES) + self.assertEqual(len(reduced_perms), 1) + self.assertEqual(reduced_perms.shape[1], 4) + + def test_tetrahedron_proper_rotation_count(self): + """Td has 12 proper rotations.""" + _, proper_rots = _analyse_point_group(TETRAHEDRON_VERTICES) + self.assertEqual(len(proper_rots), 12) + + def test_proper_rotations_are_orthogonal_with_det_plus_one(self): + """All proper rotation matrices should be orthogonal with det = +1.""" + _, proper_rots = _analyse_point_group(TETRAHEDRON_VERTICES) + for R in proper_rots: + np.testing.assert_allclose(R @ R.T, np.eye(3), atol=1e-10) + self.assertAlmostEqual(np.linalg.det(R), 1.0, places=10) + + class TestRotationAnalyser(unittest.TestCase): def setUp(self): @@ -54,7 +83,7 @@ def test_polyhedron_orientation(self): def test_discrete_orientation_identifies_correct_reference(self): """A slightly distorted copy of reference orientation 0 should be identified as such.""" - rng = np.random.RandomState(42) + rng = np.random.default_rng(42) distorted = self.ra.reference_points[0] + rng.normal(0, 0.05, self.ra.reference_points[0].shape) result = self.ra.discrete_orientation(distorted) self.assertEqual(result['reference_geometry_index'], 0) @@ -63,7 +92,7 @@ def test_discrete_orientation_identifies_correct_reference(self): def test_discrete_orientation_identifies_inverted_reference(self): """A slightly distorted copy of reference orientation 1 should be identified as such.""" - rng = np.random.RandomState(123) + rng = np.random.default_rng(123) distorted = self.ra.reference_points[1] + rng.normal(0, 0.05, self.ra.reference_points[1].shape) result = self.ra.discrete_orientation(distorted) self.assertEqual(result['reference_geometry_index'], 1) @@ -76,6 +105,53 @@ def test_discrete_orientation_perfect_geometry_has_zero_csm(self): self.assertAlmostEqual(result['symmetry_measure'], 0.0, places=5) self.assertAlmostEqual(result['rotational_distance'], 0.0, places=5) + def test_discrete_orientation_matches_brute_force(self): + """Reduced approach should give the same minimum CSM as brute-force N! permutations.""" + rng = np.random.default_rng(42) + distorted = self.ra.reference_points[0] + rng.normal(0, 0.05, self.ra.reference_points[0].shape) + result = self.ra.discrete_orientation(distorted) + # Brute force over all 4! = 24 permutations for reference 0 + points = distorted - np.mean(distorted, axis=0) + best_csm = float('inf') + for perm in itertools.permutations(range(4)): + csm_result = continuous_symmetry_measure( + points, self.ra.reference_points[0][list(perm)]) + if csm_result.symmetry_measure < best_csm: + best_csm = csm_result.symmetry_measure + self.assertAlmostEqual(result['symmetry_measure'], best_csm, places=10) + + def test_all_rotational_distances_shape(self): + """all_rotational_distances should have length n_refs * n_proper.""" + points = self.ra.reference_points[0].copy() + result = self.ra.discrete_orientation(points) + n_refs = len(self.ra.reference_points) + n_proper = len(self.ra.proper_rotations) + self.assertEqual(len(result['all_rotational_distances']), n_refs * n_proper) + + def test_orientation_index_indexes_into_all_rotational_distances(self): + """orientation_index should correspond to the minimum in all_rotational_distances.""" + points = self.ra.reference_points[0].copy() + result = self.ra.discrete_orientation(points) + idx = result['orientation_index'] + self.assertEqual( + result['rotational_distance'], + result['all_rotational_distances'][idx]) + + def test_discrete_orientation_does_not_mutate_input(self): + """The input points array should not be modified.""" + points = self.ra.reference_points[0] + 0.5 # offset from origin + original = points.copy() + self.ra.discrete_orientation(points) + np.testing.assert_array_equal(points, original) + + def test_proper_rotation_of_reference_gives_zero_distance(self): + """Applying a proper rotation to the reference should yield zero rotational distance.""" + R = self.ra.proper_rotations[1] # non-identity proper rotation + rotated = (R @ self.ra.reference_points[0].T).T + result = self.ra.discrete_orientation(rotated) + self.assertAlmostEqual(result['symmetry_measure'], 0.0, places=5) + self.assertAlmostEqual(result['rotational_distance'], 0.0, places=5) + if __name__ == '__main__': unittest.main() From cd64321b26c619b1671b6950a55ac541a62f207f Mon Sep 17 00:00:00 2001 From: Benjamin Date: Mon, 2 Mar 2026 07:56:43 +0000 Subject: [PATCH 6/6] Add docstring, input validation, and test for discrete_orientation - Add Google-style docstring to polyhedron_orientation - Validate points shape in discrete_orientation - Add test for mismatched vertex count --- polyhedral_analysis/rotation_analyser.py | 18 ++++++++++++++++++ tests/test_rotation_analyser.py | 6 ++++++ 2 files changed, 24 insertions(+) diff --git a/polyhedral_analysis/rotation_analyser.py b/polyhedral_analysis/rotation_analyser.py index 692f529..4da1440 100644 --- a/polyhedral_analysis/rotation_analyser.py +++ b/polyhedral_analysis/rotation_analyser.py @@ -216,6 +216,11 @@ def discrete_orientation(self, - ``symmetry_measure`` (float): The continuous symmetry measure (CSM) for this polyhedron. """ + expected_n = self.reference_points.shape[1] + if points.shape != (expected_n, 3): + raise ValueError( + f'Expected points with shape ({expected_n}, 3), ' + f'got {points.shape}.') points = points - np.mean(points, axis=0, dtype=float) n_refs = len(self.reference_points) n_proper = len(self.proper_rotations) @@ -248,5 +253,18 @@ def discrete_orientation(self, def polyhedron_orientation(self, polyhedron: CoordinationPolyhedron) -> OrientationDict: + """Find the discrete orientation of a coordination polyhedron. + + Convenience wrapper around :meth:`discrete_orientation` that + extracts vertex vectors from the polyhedron relative to its + central atom. + + Args: + polyhedron: The coordination polyhedron to analyse. + + Returns: + Dictionary describing the orientation. + See :meth:`discrete_orientation` for details. + """ points = polyhedron.vertex_vectors(reference='central_atom') return self.discrete_orientation(points) diff --git a/tests/test_rotation_analyser.py b/tests/test_rotation_analyser.py index c11fedc..0d5ae5e 100644 --- a/tests/test_rotation_analyser.py +++ b/tests/test_rotation_analyser.py @@ -120,6 +120,12 @@ def test_discrete_orientation_matches_brute_force(self): best_csm = csm_result.symmetry_measure self.assertAlmostEqual(result['symmetry_measure'], best_csm, places=10) + def test_discrete_orientation_rejects_wrong_vertex_count(self): + """discrete_orientation should raise ValueError for mismatched vertex count.""" + wrong_shape = np.zeros((3, 3)) # 3 vertices instead of 4 + with self.assertRaises(ValueError): + self.ra.discrete_orientation(wrong_shape) + def test_all_rotational_distances_shape(self): """all_rotational_distances should have length n_refs * n_proper.""" points = self.ra.reference_points[0].copy()