diff --git a/.github/workflows/environments/requirements-test-3.10.txt b/.github/workflows/environments/requirements-test-3.10.txt index 0297cd1..1854c61 100644 --- a/.github/workflows/environments/requirements-test-3.10.txt +++ b/.github/workflows/environments/requirements-test-3.10.txt @@ -4,6 +4,8 @@ exceptiongroup==1.3.0 # via pytest iniconfig==2.3.0 # via pytest +more-itertools==10.8.0 + # via -r requirements-test.in numpy==2.2.6 # via # -r requirements-test.in diff --git a/.github/workflows/environments/requirements-test-3.11.txt b/.github/workflows/environments/requirements-test-3.11.txt index 580469b..80be02c 100644 --- a/.github/workflows/environments/requirements-test-3.11.txt +++ b/.github/workflows/environments/requirements-test-3.11.txt @@ -2,6 +2,8 @@ # uv pip compile --python-version=3.11 requirements-test.in --output-file=requirements-test-3.11.txt iniconfig==2.3.0 # via pytest +more-itertools==10.8.0 + # via -r requirements-test.in numpy==2.3.4 # via # -r requirements-test.in diff --git a/.github/workflows/environments/requirements-test-3.12.txt b/.github/workflows/environments/requirements-test-3.12.txt index 300da18..f70586b 100644 --- a/.github/workflows/environments/requirements-test-3.12.txt +++ b/.github/workflows/environments/requirements-test-3.12.txt @@ -2,6 +2,8 @@ # uv pip compile --python-version=3.12 requirements-test.in --output-file=requirements-test-3.12.txt iniconfig==2.3.0 # via pytest +more-itertools==10.8.0 + # via -r requirements-test.in numpy==2.3.4 # via # -r requirements-test.in diff --git a/.github/workflows/environments/requirements-test-3.13.txt b/.github/workflows/environments/requirements-test-3.13.txt index 1212de8..93a4fea 100644 --- a/.github/workflows/environments/requirements-test-3.13.txt +++ b/.github/workflows/environments/requirements-test-3.13.txt @@ -2,6 +2,8 @@ # uv pip compile --python-version=3.13 requirements-test.in --output-file=requirements-test-3.13.txt iniconfig==2.3.0 # via pytest +more-itertools==10.8.0 + # via -r requirements-test.in numpy==2.3.4 # via # -r requirements-test.in diff --git a/.github/workflows/environments/requirements-test-3.9.txt b/.github/workflows/environments/requirements-test-3.9.txt index ad76788..5418305 100644 --- a/.github/workflows/environments/requirements-test-3.9.txt +++ b/.github/workflows/environments/requirements-test-3.9.txt @@ -4,6 +4,8 @@ exceptiongroup==1.3.0 # via pytest iniconfig==2.1.0 # via pytest +more-itertools==10.8.0 + # via -r requirements-test.in numpy==2.0.2 # via # -r requirements-test.in diff --git a/.github/workflows/environments/requirements-test.in b/.github/workflows/environments/requirements-test.in index f651ca5..20ab1a5 100644 --- a/.github/workflows/environments/requirements-test.in +++ b/.github/workflows/environments/requirements-test.in @@ -1,3 +1,4 @@ +more_itertools numpy pytest scipy diff --git a/doc/package-rowan.rst b/doc/package-rowan.rst index 14fb936..e58b2a9 100644 --- a/doc/package-rowan.rst +++ b/doc/package-rowan.rst @@ -33,6 +33,7 @@ rowan rowan.isfinite rowan.isinf rowan.isnan + rowan.SymmetricallyEquivalentQuaternions .. rubric:: Details diff --git a/pyproject.toml b/pyproject.toml index e5a10ab..a05ad3f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,6 +24,7 @@ classifiers = [ ] dependencies = [ "numpy>=1.21", + "more_itertools" ] [project.optional-dependencies] diff --git a/rowan/__init__.py b/rowan/__init__.py index 69b90ff..c5e6811 100644 --- a/rowan/__init__.py +++ b/rowan/__init__.py @@ -52,6 +52,7 @@ to_matrix, vector_vector_rotation, ) +from .mapping.symmetries import SymmetricallyEquivalentQuaternions # Get the version __version__ = "1.3.2" @@ -94,4 +95,5 @@ "to_euler", "to_matrix", "vector_vector_rotation", + "SymmetricallyEquivalentQuaternions", ] diff --git a/rowan/grids/__init__.py b/rowan/grids/__init__.py new file mode 100644 index 0000000..1cb9fdd --- /dev/null +++ b/rowan/grids/__init__.py @@ -0,0 +1,59 @@ +r"""Low-dispersion grids on manifolds related to rotation. + +Finding a uniform of `n` rotations (in any representation) is a generalized Tammes +problem on the manifold associated with that representation. Grids of rotation axes +reduce to the standard Tammes problem on the 2-sphere, and grids of quaternions are a +generalization to the 3-sphere. Closed form solutions for general `n` are not possible, +but we can create nearly uniform (low-dispersion) grids using Fibonacci lattices. +""" + +import numpy as np + +"""The Golden ratio, 1.61803398875...""" +GOLDEN_RATIO = (1 + np.sqrt(5)) / 2 + + +"""The angle that divides 2π radians into the golden ratio. +This takes the longer of the two segments π(√5 - 1), where π[(√5 - 1) + (3-√5)] = 2π +""" +GOLDEN_ANGLE = np.pi * (np.sqrt(5) - 1) + + +def spherical_fibonacci_lattice(n): + """A Fibonacci lattice of `n` points on the 2-sphere. + + Args: + n (int): The number of points in the lattice. + """ # TODO: example image! + t = np.arange(n) + + z = 1 - (t / (n - 1)) * 2 + radius = np.sqrt(1 - z**2) + theta = (GOLDEN_ANGLE * t) % (2 * np.pi) + + x = np.cos(theta) * radius + y = np.sin(theta) * radius + + return np.column_stack([x, y, z]).T + + +def quaternion_fibonacci_lattice(n): + """A near-uniform grid of `n` quaternions. + + This is equivalent to a Fibonacci lattice on the 3-sphere. See + `this paper `_ for a derivation. + + Args: + n (int): The number of points in the lattice. + """ + PSI = 1.533751168755204288118041413 # Solution to ψ**4 = ψ + 4 + s = np.arange(n) + 1 / 2 + t = s / n + d = 2 * np.pi * s + r0, r1 = (np.sqrt(t), np.sqrt(1 - t)) + α, β = (d / np.sqrt(2), d / PSI) + + # Allocate as rows and then transpose, rather than stacking columns + result = np.empty((4, n)) + result[...] = r0 * np.sin(α), r0 * np.cos(α), r1 * np.sin(β), r1 * np.cos(β) + return result.T diff --git a/rowan/mapping/__init__.py b/rowan/mapping/__init__.py index 4de0f71..8fd7ec6 100644 --- a/rowan/mapping/__init__.py +++ b/rowan/mapping/__init__.py @@ -59,7 +59,12 @@ from ..functions import from_matrix, rotate, to_matrix from ..geometry import angle -__all__ = ["kabsch", "davenport", "procrustes", "icp"] +__all__ = [ + "kabsch", + "davenport", + "procrustes", + "icp", +] def kabsch(X, Y, require_rotation=True): diff --git a/rowan/mapping/symmetries.py b/rowan/mapping/symmetries.py new file mode 100644 index 0000000..a8f058e --- /dev/null +++ b/rowan/mapping/symmetries.py @@ -0,0 +1,196 @@ +""".""" + +from itertools import combinations, permutations, product +from typing import Iterable + +import numpy as np +from more_itertools import distinct_permutations +from numpy.typing import ArrayLike + +# from ..functions import from_axis_angle +from rowan.functions import from_axis_angle + + +def sign_changes(even: bool = True, d: int = 3): + """Get all even (or odd) sign changes of a vector in ℝ3.""" + all_changes = np.array([*product([1, -1], repeat=d)]) + if even == "all": + return all_changes + return all_changes[np.prod(all_changes, axis=1) == (1 if even else -1)] + + +def generate_tetrahedral_group(): + """Generates the 24 quaternions of the chiral tetrahedral group .""" + quats = [ + # 8 (4+4) quaternions from permutationsof (±1, 0, 0, 0) + (1, 0, 0, 0), + (-1, 0, 0, 0), + (0, 1, 0, 0), + (0, -1, 0, 0), + (0, 0, 1, 0), + (0, 0, -1, 0), + (0, 0, 0, 1), + (0, 0, 0, -1), + *product([-0.5, 0.5], repeat=4), # 16 quaternions from 1/2 * (±1, ±1, ±1, ±1) + ] + return np.array(quats) + + +def generate_octahedral_group(): + """Generates the 48 quaternions of the octahedral group.""" + c = 1 / np.sqrt(2) + + return np.array( + [ + # 24 elements of the tetrahedral group + *generate_tetrahedral_group(), + # 24 quaternions from distinct permutations of 1/√2 * (±1, ±1, 0, 0) + *distinct_permutations([c, c, 0, 0]), + *distinct_permutations([c, -c, 0, 0]), + *distinct_permutations([-c, -c, 0, 0]), + ] + ) + + +def generate_cyclic_group(n: int, axis: ArrayLike = (0, 0, 1)): + """Generates the 2n quaternions of the cyclic group . + + Args: + n (int): The index of the cyclic group C_n. The order of the + resulting group will be 2n. + axis (np.ndarray): The axis of rotation for the group. + + Returns: + np.ndarray: A (2n, 4) array of quaternions. + """ + # return np.array([from_axis_angle(axis, 2 * np.pi * k / n) for k in range(n)]) + return from_axis_angle(axis, np.linspace(0, 2 * np.pi, n, endpoint=False)).round(15) + + +# def generate_dicyclic_group(n: int, axis: ArrayLike = (0, 0, 1)): +# """Generates the 4n quaternions of the dicyclic group . + +# Args: +# n (int): The index of the cyclic group D_n. The order of the +# resulting group will be 4n. +# axis (np.ndarray): The axis of rotation for the group. + +# Returns: +# np.ndarray: A (4n, 4) array of quaternions. +# """ +# return np.concatenate( +# ( +# # 2n quaternions from the cyclic group +# generate_cyclic_group(n, axis), +# # hetas = np.linspace(0, np.pi, n, endpoint=False) +# # rv = np.pi * np.vstack([np.zeros(n), np.cos(thetas), np.sin(thetas)]).T +# # g2 = np.roll(rv, axis, axis=1) +# # from_axis_angle(axis, np.linspace(0, np.pi, n, endpoint=False)) # TODO +# ) +# ) + + +def even_permutations(it: Iterable): + """Return permutations containing an even number of transpositions from an iterable. + + Args: + it (typing.Iterable): The iterable to permute. + + Returns: + typing.Generator: An (N!/2, len(it)) generator of permutations. + """ + return ( + p for p in permutations(it) if not sum(a > b for a, b in combinations(p, 2)) % 2 + ) + + +def generate_icosahedral_group(): + """Generates the 48 quaternions of the octahedral group.""" + φ = (1 + np.sqrt(5)) / 2 + + return np.array( + [ + # 24 elements of the tetrahedral group + *generate_tetrahedral_group(), + # 96 quaternions from even permutations of 1/2 * (0, ±1, ±1/φ, ±φ) + *( + q + for a, b, c in product([1 / 2, -1 / 2], repeat=3) + for q in even_permutations([0, a, b / φ, c * φ]) + ), + ] + ) + + +class SymmetricallyEquivalentQuaternions(np.ndarray): + """Equivalent quaternions for a particular point group. + + This class is a read-only subclass of `np.ndarray` that holds the set of + quaternions representing the symmetry operations of a specific point group + ('T', 'O', or 'I'). + """ + + def __new__(cls, data, group): # noqa: D102 + obj = np.asarray(data).view(cls) + obj.flags.writeable = False + obj._group = group # noqa: SLF001 + return obj + + def __array_finalize__(self, obj): + """Finalize the array, ensuring attributes are passed on.""" + if obj is None: + return + self._group = getattr(obj, "_group", None) + + @classmethod + def create_group(cls, group: str): + """Create the set of symmetrically equivalent quaternions for a rotation group. + + Args: + group ({'T', 'O', 'I'}): Schönflies notation for the group. + + Returns: + (..., 4) :class:`numpy.ndarray`: + Eqivalent quaternions ``q`` for the provided point group. + + Example:: + >>> from rowan import SymmetricallyEquivalentQuaternions + >>> tetrahedral_quats = SymmetricallyEquivalentQuaternions.create_group("T") + >>> tetrahedral_quats.shape + (24, 4) + >>> assert np.array_equal(tetrahedral_quats[0], [1,0,0,0]) + >>> SymmetricallyEquivalentQuaternions.create_group("T")[:] + [[ 1. 0. 0. 0. ] + [-1. 0. 0. 0. ] + [ 0. 1. 0. 0. ] + ... + [ 0.5 0.5 -0.5 0.5] + [ 0.5 0.5 0.5 -0.5] + [ 0.5 0.5 0.5 0.5]] + """ + if group == "T": + return cls(data=generate_tetrahedral_group(), group=group) + if group == "O": + return cls(data=generate_octahedral_group(), group=group) + if group == "I": + return cls(data=generate_icosahedral_group(), group=group) + if group[0] == "C": + return cls( + data=generate_cyclic_group(int(group[1:]), axis=[0, 0, 1]), group=group + ) + msg = ( + f"Unknown group '{group}' does not match valid options " + "{{'T', 'O', 'I'}}" + # "{{'T', 'O', 'I', 'Cn'}}" + ) + raise ValueError(msg) + + def __str__(self): # noqa: D105 + arrstr = np.array2string(self, floatmode="maxprec", threshold=24, precision=9) + return ( + f"SymmetricallyEquivalentQuaternions['{self._group}'] " + f"(order {len(self)}):\n{arrstr}" + ) + + def __repr__(self): # noqa: D105 + return self.__str__() diff --git a/tests/test_mapping.py b/tests/test_mapping.py index 81687f9..f8c6723 100644 --- a/tests/test_mapping.py +++ b/tests/test_mapping.py @@ -1,16 +1,37 @@ """Test algorithms for point-cloud mapping.""" import unittest +from itertools import product import numpy as np +import scipy -from rowan import from_axis_angle, from_matrix, mapping, random, rotate +from rowan import ( + SymmetricallyEquivalentQuaternions, + from_axis_angle, + from_matrix, + mapping, + random, + rotate, +) zero = np.array([0, 0, 0, 0]) one = np.array([1, 0, 0, 0]) half = np.array([0.5, 0.5, 0.5, 0.5]) +def _cyclic_permutations(a): + def make_permutation_generator(a, forward=True): + """Return a list of functions that roll an array by 0,1,2...n.""" + n = len(a) + sign = -1 if forward else 1 + return [ + lambda arr, k=j: arr[(sign * k) :] + arr[: (sign * k)] for j in range(n) + ] + + return [op(a) for op in make_permutation_generator(a)] + + class TestMapping(unittest.TestCase): """Test mapping functions.""" @@ -191,6 +212,97 @@ def test_icp_mismatched(self): # a heuristic since we can't guarantee exact matches assert np.mean(norms) < 0.5 + def test_symmetric_quaternions_tetrahedral(self): + """Verify our symmetrically equivalent quaternions are correct.""" + tet = [ + (-0.5, -0.5, 0.5), + (-0.5, 0.5, -0.5), + (0.5, -0.5, -0.5), + (0.5, 0.5, 0.5), + ] # Verts in sorted order + + quats = SymmetricallyEquivalentQuaternions.create_group("T") + assert len(quats) == 24 + for quat in quats: + res = sorted(rotate(quat, tet).tolist()) + np.testing.assert_array_equal(res, tet, err_msg=f"q={quat}") + + # Compare our array of generated quaternions to Scipy + np.testing.assert_array_equal( + sorted(quats.tolist()), + sorted( + [ + *scipy.spatial.transform.Rotation.create_group("T") + .as_quat() + .tolist(), + *( + -scipy.spatial.transform.Rotation.create_group("T").as_quat() + ).tolist(), + ] + ), + ) + + def test_symmetric_quaternions_octahedral(self): + """Verify our symmetrically equivalent quaternions are correct.""" + cube = [*product([-0.5, 0.5], repeat=3)] # Verts in sorted order + + quats = SymmetricallyEquivalentQuaternions.create_group("O") + assert len(quats) == 48 + for quat in quats: + res = sorted(rotate(quat, cube).tolist()) + np.testing.assert_allclose(res, cube, err_msg=f"q={quat}", atol=4e-16) + + # Compare our array of generated quaternions to Scipy + np.testing.assert_allclose( + sorted(quats.tolist()), + sorted( + [ + *scipy.spatial.transform.Rotation.create_group("O") + .as_quat() + .tolist(), + *( + -scipy.spatial.transform.Rotation.create_group("O").as_quat() + ).tolist(), + ] + ), + atol=4e-16, + ) + + def test_symmetric_quaternions_icosahedral(self): + """Verify our symmetrically equivalent quaternions are correct.""" + # φ = (1 + np.sqrt(5)) / 2 + # icos = np.sort( + # [ + # *_cyclic_permutations([0, -1, -φ]), + # *_cyclic_permutations([0, -1, φ]), + # *_cyclic_permutations([0, 1, -φ]), + # *_cyclic_permutations([0, 1, φ]), + # ], + # axis=0, + # ) + + quats = SymmetricallyEquivalentQuaternions.create_group("I") + assert len(quats) == 120 + + # TODO: we get the correct vertices, but the ordering is difficult to ensure. + # Is there a better way to test for equivalence here? + # for i, quat in enumerate(quats): + # res = np.sort(rotate(quat, icos).round(14), axis=0) + # np.testing.assert_allclose(res, icos, err_msg=f"q{i}={quat}", atol=4e-16) + + # Compare our array of generated quaternions to Scipy + np.testing.assert_allclose( + np.sort(quats, axis=0), + np.sort( + [ + *scipy.spatial.transform.Rotation.create_group("I").as_quat(), + *(-scipy.spatial.transform.Rotation.create_group("I").as_quat()), + ], + axis=0, + ), + atol=4e-16, + ) + if __name__ == "__main__": unittest.main()