Add a helper for mapping between different algebras which share basis…

… blades

clifford/__init__.py

@@ -320,7 +320,7 @@ def val_get_right_gmt_matrix(mt: sparse.COO, x):
 from ._layout_helpers import BasisVectorIds, BasisBladeOrder  # noqa: F401
 from ._mvarray import MVArray, array  # noqa: F401
 from ._frame import Frame  # noqa: F401
-from ._blademap import BladeMap  # noqa: F401
+from ._blademap import BladeMap, BasisVectorMap  # noqa: F401
 
 
 # copied from the itertools docs

clifford/_bit_helpers.py

@@ -46,3 +46,55 @@ def count_set_bits(bitmap: int) -> int:
     def count_set_bits(x: int) -> int:
         """ Counts the number of bits set to 1 in bitmap """
         return __builtin_popcnt(x)
+
+
+def left_shift(x, shift):
+    if shift >= 0:
+        return x << shift
+    else:
+        return x >> -shift
+
+
+class BitPermuter:
+    """
+    Produce a callable that sets bit `i` of out to bit `src_bits[i]` of in.
+
+    Attributes
+    ----------
+    domain : int
+        A mask of the bits that this callable accepts
+    range : int
+        A mask of the bits that this callable produces
+    """
+    _inverse = None
+
+    def __init__(self, src_to_dst):
+        # find bits which need shifting by the same amount
+        self._mask_for_shift = {}
+        self.domain = 0
+        self.range = 0
+        for s, d in src_to_dst.items():
+            self._mask_for_shift[d - s] = self._mask_for_shift.setdefault(d - s, 0) | 1 << s
+            self.domain |= 1 << s
+            self.range |= 1 << d
+
+        # this makes the degenerate case work for numpy arrays
+        if not src_to_dst:
+            self._mask_for_shift[0] = 0
+
+        print([(shift, bin(mask)) for shift, mask in self._mask_for_shift.items()])
+
+    def __call__(self, bitmap):
+        ret = 0
+        for shift, mask in self._mask_for_shift.items():
+            ret |= left_shift(bitmap & mask, shift)
+        return ret
+
+    def inverse(self, bitmap):
+        ret = 0
+        for shift, mask in self._mask_for_shift.items():
+            ret |= left_shift(bitmap, -shift) & mask
+        return ret
+
+    def __repr__(self):
+        return "BitPermuter({})".format(self)

clifford/_blademap.py

@@ -1,4 +1,149 @@
-class BladeMap(object):
+from ._layout import Layout
+from ._multivector import MultiVector
+from ._bit_helpers import BitPermuter
+
+import functools
+import operator
+
+
+class _BaseTransformation:
+    """ Base class for transforming multivectors between one layout and another """
+    def __init__(self, src_layout: Layout, dst_layout: Layout = None):
+        self.src_layout = src_layout
+        if dst_layout is None:
+            dst_layout = src_layout
+        self.dst_layout = dst_layout
+
+    def __repr__(self):
+        return "<{} from {} to {}>".format(type(self).__name__, self.src_layout, self.dst_layout)
+
+    @property
+    def inverse(self):
+        """ The inverse of this mapping """
+        self_rev = type(self).__new__()
+        self_rev.__inverse_init__(self)
+        return self
+
+    def __inverse_init__(self, orig):
+        """ a hook to be implemented by subclasses """
+        self.src_layout = orig.dst_layout
+        self.dst_layout = orig.src_layout
+
+
+class _IndexTransformation(_BaseTransformation):
+    """ A transformation that moves around multivector component indices """
+    def __init__(self, src_layout: Layout, dst_layout: Layout, src_inds, dst_inds):
+        super().__init__(src_layout, dst_layout)
+        self._src_inds = src_inds
+        self._dst_inds = dst_inds
+
+    def __call__(self, mv_src: MultiVector) -> MultiVector:
+        """ Apply this transformation to a MultiVector """
+        if mv_src.layout != self.src_layout:
+            raise ValueError("Multivector must belong to {}".format(self.src_layout))
+
+        if self._dst_inds == slice(None):
+            # optimization to avoid a copy
+            return self.dst_layout.MultiVector(mv_src.value[self._src_inds])
+        else:
+            mv_dst = self.dst_layout.MultiVector(dtype=mv_src.value.dtype)
+            mv_dst.value[self._dst_inds] = mv_src.value[self._src_inds]
+            return mv_dst
+
+    def __inverse_init__(self, orig):
+        super().__inverse_init__(orig)
+        self._dst_inds = orig._src_inds
+        self._src_inds = orig._dst_inds
+
+
+class BasisVectorMap(_IndexTransformation):
+    """
+    A map that converts blades from one algebra to another via converting the
+    basis vectors.
+
+    Note that does not support any scale factors, so is most useful for adding
+    and removing basis vectors from multivectors, rather than for transforming
+    the basis.
+
+    Blades from the source algebra containing basis vectors not present in the
+    destination algebra are dropped completely.
+    """
+    def __init__(
+        self,
+        src_layout: Layout, dst_layout: Layout,
+        src_vectors: list = None, dst_vectors: list = None
+    ):
+        all_src_vectors = src_layout._basis_vector_ids.values
+        all_dst_vectors = dst_layout._basis_vector_ids.values
+
+        # handle default arguments
+        if dst_vectors is None and src_vectors is None:
+            # use the common vectors
+            src_vectors = list(set(all_src_vectors) & set(all_dst_vectors))
+            dst_vectors = src_vectors
+        elif dst_vectors is None:
+            dst_vectors = src_vectors
+
+        try:
+            src_bits = [all_src_vectors.index(f) for f in src_vectors]
+        except IndexError:
+            raise ValueError("Unknown ids {!r} in src_vectors".format(set(src_vectors) - set(all_src_vectors))) from None
+        try:
+            dst_bits = [all_dst_vectors.index(f) for f in dst_vectors]
+        except IndexError:
+            raise ValueError("Unknown ids {!r} in dst_vectors".format(set(dst_vectors) - set(all_dst_vectors))) from None
+
+        # work out where each bit in `src` goes
+        src_bit_to_dst_bit = {}
+        for d, s in zip(dst_bits, src_bits):
+            if s in src_bit_to_dst_bit:
+                raise ValueError(
+                    "Bit src[{}] maps to both dst[{}] and dst[{}]".format(s, src_bit_to_dst_bit[s], d))
+            src_bit_to_dst_bit[s] = d
+
+        permute = BitPermuter(src_bit_to_dst_bit)
+
+        if len(src_bit_to_dst_bit) == len(all_dst_vectors):
+            # every destination element has a source (some source bits may be discarded)
+            src_inds = src_layout._basis_blade_order.bitmap_to_index[
+                permute.inverse(dst_layout._basis_blade_order.index_to_bitmap)
+            ]
+            dst_inds = slice(None)
+        elif len(src_bit_to_dst_bit) == len(all_src_vectors):
+            # every source bit maps to a destination (some destination bits will be 0)
+            src_inds = slice(None)
+            dst_inds = dst_layout._basis_blade_order.bitmap_to_index[
+                permute(src_layout._basis_blade_order.index_to_bitmap)
+            ]
+        else:
+            # a combination of the above two cases
+            src_inds = src_layout._basis_blade_order.index_to_bitmap & ~permute.domain == 0
+            dst_inds = dst_layout._basis_blade_order.bitmap_to_index[
+                permute.inverse(src_layout._basis_blade_order.index_to_bitmap[src_inds])
+            ]
+
+        super().__init__(src_layout, dst_layout, src_inds, dst_inds)
+        self._src_vectors = src_vectors
+        self._dst_vectors = dst_vectors
+
+    def __inverse_init__(self, orig):
+        super().__inverse_init__(orig)
+        self._dst_vectors = orig._src_vectors
+        self._src_vectors = orig._dst_vectors
+
+    def __repr__(self):
+        return "<{} from {} to {} with {}>".format(
+            type(self).__name__,
+            self.src_layout,
+            self.dst_layout,
+            ', '.join(
+                '{}->{}'.format(s, d)
+                for s, d in zip(self._src_vectors, self._dst_vectors)
+            )
+        )
+
+
+class BladeMap(_BaseTransformation):
     '''
     A Map Relating Blades in two different algebras
 
@@ -29,6 +174,9 @@ def __init__(self, blades_map, map_scalars=True):
             s2 = self.b2[0]._newMV(dtype=int)+1
             self.blades_map = [(s1, s2)] + self.blades_map
 
+        first_src, first_dest = blades_map[0]
+        super().__init__(self, first_src.layout, first_dest.layout)
+
     @property
     def b1(self):
         return [k[0] for k in self.blades_map]
@@ -39,21 +187,21 @@ def b2(self):
 
     @property
     def layout1(self):
-        return self.b1[0].layout
+        return self.src_layout
 
     @property
     def layout2(self):
-        return self.b2[0].layout
+        return self.dst_layout
 
     def __call__(self, A):
         '''map an MV `A` according to blade_map'''
 
         # determine direction of map
-        if A.layout == self.layout1:
+        if A.layout == self.src_layout:
             from_b = self.b1
             to_b = self.b2
 
-        elif A.layout == self.layout2:
+        elif A.layout == self.dst_layout:
             from_b = self.b2
             to_b = self.b1
         else:

clifford/test/test_bit_helpers.py

@@ -1,5 +1,8 @@
 """ Tests of clifford._bit_helpers """
-from clifford._bit_helpers import count_set_bits, set_bit_indices
+import pytest
+import numpy as np
+
+from clifford._bit_helpers import count_set_bits, set_bit_indices, BitPermuter
 
 
 def test_count_bits():
@@ -14,3 +17,78 @@ def test_bit_indices():
     assert list(set_bit_indices(0b1)) == [0]
     assert list(set_bit_indices(0b101)) == [0, 2]
     assert list(set_bit_indices(0b101010)) == [1, 3, 5]
+
+
+
+class TestPermuter:
+    @pytest.fixture(scope='class')
+    def permuter(self):
+        r"""
+        Performs the transformation::
+
+            src 0  1  2  x  x  x
+                |  |  |
+                |  \-----------\
+                |     |        |
+                \--------\     |
+                      |  |     |
+                   /--/  |     |
+                   V     V     V
+            dst x  1  x  3  x  5
+        """
+        return BitPermuter({0: 3, 1: 5, 2: 1})
+
+    @pytest.fixture(scope='class')
+    def null_permuter(self):
+        return BitPermuter({})
+
+    def test_call(self, permuter):
+        assert permuter(0b000) == 0b000000
+
+        assert permuter(0b001) == 0b001000
+        assert permuter(0b010) == 0b100000
+        assert permuter(0b100) == 0b000010
+
+        # out of range bits are ignored
+        assert permuter(0b100111) == 0b101010
+
+    def test_attributes(self, permuter):
+        assert permuter.domain == 0b111
+        assert permuter.range == 0b101010
+        assert permuter(permuter.domain) == permuter.range
+
+    def test_attributes_null(self, null_permuter):
+        assert null_permuter.domain == 0
+        assert null_permuter.range == 0
+
+    def test_inverse(self, permuter):
+        # test the inverse
+        assert permuter.inverse(permuter(0b001)) == 0b001
+        assert permuter.inverse(permuter(0b010)) == 0b010
+        assert permuter.inverse(permuter(0b100)) == 0b100
+
+        assert permuter(permuter.inverse(0b001000)) == 0b001000
+        assert permuter(permuter.inverse(0b100000)) == 0b100000
+        assert permuter(permuter.inverse(0b000010)) == 0b000010
+
+        assert permuter.inverse(permuter.range) == permuter.domain
+
+    def test_arrays(self, permuter):
+        np.testing.assert_equal(
+            permuter(np.array([0b001, 0b010, 0b100])),
+            np.array([0b001000, 0b100000, 0b000010]),
+        )
+        np.testing.assert_equal(
+            permuter.inverse(np.array([0b001000, 0b100000, 0b000010])),
+            np.array([0b001, 0b010, 0b100]),
+        )
+
+    def test_arrays_null(self, null_permuter):
+        np.testing.assert_equal(
+            null_permuter(np.array([0b001, 0b010, 0b100])),
+            np.array([0, 0, 0]),
+        )
+        np.testing.assert_equal(
+            null_permuter.inverse(np.array([0b001000, 0b100000, 0b000010])),
+            np.array([0, 0, 0]),
+        )