Add LinearMatrix.from_mapping, and replace BladeMap with it

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

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.