From ac93fee58fb64410a7cdf21a21292aaaa646337f Mon Sep 17 00:00:00 2001 From: Eric Wieser Date: Wed, 8 Apr 2020 15:23:07 +0100 Subject: [PATCH] Add `LinearMatrix.from_mapping`, and replace `BladeMap` with it For the example application in spacetime algebra, `bm(x)` is at least 10x faster. --- clifford/_blademap.py | 46 ++++++++++---------- clifford/transformations.py | 85 +++++++++++++++++++++++++++++++++---- 2 files changed, 100 insertions(+), 31 deletions(-) diff --git a/clifford/_blademap.py b/clifford/_blademap.py index ba122f69..a4a67fca 100644 --- a/clifford/_blademap.py +++ b/clifford/_blademap.py @@ -1,6 +1,11 @@ -class BladeMap(object): +from .transformations import LinearMatrix + +class BladeMap: ''' - A Map Relating Blades in two different algebras + A Map Relating Blades in two different algebras. + + This is now just a thin wrapper around :func:`LinearMatrix.from_mapping`, + which is more powerful. Examples ----------- @@ -12,13 +17,15 @@ class BladeMap(object): >>> # Pauli Algebra `P` >>> P, P_blades = Cl(3, names='p') >>> locals().update(P_blades) - >>> sta_split = BladeMap([(d01, p1), - ... (d02, p2), - ... (d03, p3), - ... (d12, p12), - ... (d23, p23), - ... (d13, p13)]) - + >>> sta_split = BladeMap([ + ... (d01, p1), + ... (d02, p2), + ... (d03, p3), + ... (d12, p12), + ... (d23, p23), + ... (d13, p13), + ... (d0123, p123) + ... ]) ''' def __init__(self, blades_map, map_scalars=True): self.blades_map = blades_map @@ -29,6 +36,8 @@ def __init__(self, blades_map, map_scalars=True): s2 = self.b2[0]._newMV(dtype=int)+1 self.blades_map = [(s1, s2)] + self.blades_map + self._transformation = LinearMatrix.from_mapping(self.blades_map) + @property def b1(self): return [k[0] for k in self.blades_map] @@ -39,28 +48,19 @@ def b2(self): @property def layout1(self): - return self.b1[0].layout + return self._transformation.layout_src @property def layout2(self): - return self.b2[0].layout + return self._transformation.layout_dst def __call__(self, A): '''map an MV `A` according to blade_map''' # determine direction of map - if A.layout == self.layout1: - from_b = self.b1 - to_b = self.b2 - + if A.layout == self._transformation.layout_src: + return _transformation(A) elif A.layout == self.layout2: - from_b = self.b2 - to_b = self.b1 + return _transformation.adjoint(A) else: raise ValueError('A doesnt belong to either Algebra in this Map') - - # create empty MV, and map values - B = to_b[0]._newMV(dtype=int) - for from_obj, to_obj in zip(from_b, to_b): - B += (sum(A.value*from_obj.value)*to_obj) - return B diff --git a/clifford/transformations.py b/clifford/transformations.py index 831ad199..ba24733f 100644 --- a/clifford/transformations.py +++ b/clifford/transformations.py @@ -30,7 +30,7 @@ This module may in future become the home for optimized rotor transformations, or non-linear transformations """ -from typing import Dict, Any, Callable +from typing import Dict, Any, Callable, List, Tuple from abc import abstractmethod import numpy as np @@ -104,6 +104,22 @@ def _make_outermorphism( return t +class _Mismatch(ValueError): + def __init__(self, expected, actual): + self.expected = expected + self.actual = actual + + +def _all_same(it, first=None): + """ Helper to extract the first item from an iterable """ + for val in it: + if first is None: + first = val + elif first != val: + raise _Mismatch(first, val) + return first + + Transformation = Callable[[MultiVector], MultiVector] """ A callable mapping one MultiVector to another. """ @@ -245,17 +261,70 @@ def from_function(cls, func: Callable[[MultiVector], MultiVector], layout_src: L ] # check the layouts of the resulting blades match - for d in blades_dst: - if layout_dst is None: - layout_dst = d.layout - elif d.layout != layout_dst: - raise ValueError( - "result of func() is from the wrong layout, expected: {}, " - "got {}.".format(layout_dst, d.layout)) + try: + layout_dst = _all_same( + (d.layout for d in blades_dst), first=layout_dst) + except _Mismatch as m: + raise ValueError( + "result of func() is from the wrong layout, expected: {}, " + "got {}.".format(m.expected, m.actual)) matrix = np.array([b_dst.value for b_dst in blades_dst]) return cls(matrix, layout_src, layout_dst) + @classmethod + def from_mapping( + cls, + mapping: List[Tuple[MultiVector, MultiVector]], + layout_src: Layout = None, + layout_dst: Layout = None + ): + """ Compute a transformation that maps between the pairs in `mapping`. + + The mapping is computed via a least-squares fit. + + Parameters + ---------- + mapping : + List of ``(src_mv, dest_mv)`` pairs to map between. + layout_src : + The source layout. Inferred if not provided. + layout_dst : + The destination layout. Inferred if not provided. + """ + try: + layout_src = _all_same( + (s.layout for s, d in mapping), first=layout_src) + except _Mismatch as m: + raise ValueError( + "mapping source is from the wrong layout, expected: {}, " + "got {}.".format(m.expected, m.actual)) + try: + layout_dst = _all_same( + (d.layout for s, d in mapping), first=layout_dst) + except _Mismatch as m: + raise ValueError( + "mapping source is from the wrong layout, expected: {}, " + "got {}.".format(m.expected, m.actual)) + + N = len(mapping) + + mat_src = np.stack([s.value for s, d in mapping], axis=-1) + mat_dst = np.stack([d.value for s, d in mapping], axis=-1) + + # rewrite in the form needed by lstsq, + # `mat @ mat_src = mat_dst` => `mat_src.T @ mat.T = mat_dst.T` + mat_T, resid, rank, s = np.linalg.lstsq(mat_src.T, mat_dst.T, rcond=None) + mat = mat_T.T + + # if the input mapping was integral, and the matrix contains only + # integers, then for convenience keep the result as integers + rounded = mat.astype(np.result_type(mat_src.dtype, mat_dst.dtype), copy=False) + if rounded.dtype == mat.dtype or np.all(rounded == mat): + mat = rounded + + return cls(mat, layout_src, layout_dst) + @classmethod def from_rotor(cls, rotor: MultiVector) -> 'LinearMatrix': """ Build a linear transformation from the result of applying a rotor sandwich.