diff --git a/clifford/__init__.py b/clifford/__init__.py index 8e360349..9f6142d1 100644 --- a/clifford/__init__.py +++ b/clifford/__init__.py @@ -82,6 +82,7 @@ # Major library imports. import numpy as np import numba +import numba.numpy_support import sparse @@ -320,7 +321,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 diff --git a/clifford/_bit_helpers.py b/clifford/_bit_helpers.py index 57138135..b5af9433 100644 --- a/clifford/_bit_helpers.py +++ b/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) diff --git a/clifford/_blademap.py b/clifford/_blademap.py index 72c1e4b1..082c10bd 100644 --- a/clifford/_blademap.py +++ b/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: diff --git a/clifford/test/test_bit_helpers.py b/clifford/test/test_bit_helpers.py index 6b5e057a..7e9ec307 100644 --- a/clifford/test/test_bit_helpers.py +++ b/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]), + ) diff --git a/docs/tutorials/apollonius-cga-augmented.ipynb b/docs/tutorials/apollonius-cga-augmented.ipynb index d14093e9..3bffcbbe 100644 --- a/docs/tutorials/apollonius-cga-augmented.ipynb +++ b/docs/tutorials/apollonius-cga-augmented.ipynb @@ -20,25 +20,38 @@ "metadata": {}, "outputs": [], "source": [ - "from clifford import ConformalLayout, BasisVectorIds, MultiVector\n", + "from clifford import ConformalLayout, BasisVectorIds, MultiVector, BasisVectorMap\n", "\n", "class OurCustomLayout(ConformalLayout):\n", " def __init__(self, ndims):\n", " self.ndims = ndims\n", + " \n", + " euclidean_vectors = [str(i + 1) for i in range(ndims)]\n", + " conformal_vectors = ['m2', 'm1']\n", "\n", " # Construct our custom algebra. Note that ConformalLayout requires the e- and e+ basis vectors to be last.\n", " ConformalLayout.__init__(\n", " self,\n", " [1]*ndims + [-1] + [1, -1],\n", - " ids=BasisVectorIds([str(i + 1) for i in range(ndims)] + ['np1', 'm2', 'm1'])\n", + " ids=BasisVectorIds(euclidean_vectors + ['np1'] + conformal_vectors)\n", " )\n", " self.enp1 = self.basis_vectors_lst[ndims]\n", "\n", " # Construct a base algebra without the extra `enp1`, which would not be understood by pyganja.\n", " self.conformal_base = ConformalLayout(\n", " [1]*ndims + [1, -1],\n", - " ids=BasisVectorIds([str(i + 1) for i in range(ndims)] + ['m2', 'm1'])\n", - " )" + " ids=BasisVectorIds(euclidean_vectors + conformal_vectors)\n", + " )\n", + " \n", + " # this lets us convert between the two layouts\n", + " self.to_conformal = BasisVectorMap(self, self.conformal_base)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The code above also defines a stardard conformal $\\mathbb{R}^{N+1,1}$ layout without this new basis vector. This is primarily to support rendering with `pyganja`, which doesn't support the presence of this extra vector. `BasisVectorMap` defaults to preserving vectors by name between one algebra and another, while throwing away blades containing vectors missing from the destination algebra." ] }, { @@ -89,53 +102,19 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Visualization of custom algebras\n", - "\n", - "In order to render with `pyganja`, we'll need a helper to convert from our custom $\\mathbb{R}^{N+1,2}$ layout into a standard conformal $\\mathbb{R}^{N+1,1}$ layout. `clifford` maps indices in `.value` to basis blades via `layout._basis_blade_order.index_to_bitmap`, which we can use to convert the indices in one layout to the indices in another.\n", - "\n", - "
\n", "\n", - "Warning\n", "\n", - "This function uses the private attribute `Layout._basis_blade_order`, which may change or be removed in future.\n", - "Future releases of `clifford` may contain tools to make switching between algebras easier.\n", - "\n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def to_conformal(self, mv: MultiVector) -> MultiVector:\n", - " \"\"\" Convert a mv with `mv.layout == self` to one with `mv_new.layout == self.conformal_base`\n", - " \n", - " This is done by discarding coefficients of blades containing enp1\n", - " \"\"\"\n", - " bits_base = self.conformal_base._basis_blade_order.index_to_bitmap\n", - " # insert the np1 bit into each basis blade bitmask, by...\n", - " bits = (\n", - " # ... leaving e1...en in place ...\n", - " bits_base & ~(~0 << self.ndims) |\n", - " # and moving em1 and em2 one bit to the left, leaving a 0 in the np1 bit\n", - " (bits_base & (0b11 << self.ndims)) << 1\n", - " )\n", - " # convert the new bitmaps into indices into the augmented value\n", - " inds = self._basis_blade_order.bitmap_to_index[bits]\n", - "\n", - " # and pick only those indicess\n", - " return self.conformal_base.MultiVector(mv.value[inds])\n", - "\n", - "OurCustomLayout.to_conformal = to_conformal; del to_conformal" + "In order to render with `pyganja`, we'll need a helper to convert from our custom $\\mathbb{R}^{N+1,2}$ layout into a standard conformal $\\mathbb{R}^{N+1,1}$ layout. `clifford` maps indices in `.value` to basis blades via `layout._basis_blade_order.index_to_bitmap`, which we can use to convert the indices in one layout to the indices in another." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Finally, we'll define a plotting function, which plots the problem and solution circles in suitable colors via `pyganja`:" + "## Visualization\n", + "\n", + "Finally, we'll define a plotting function, which plots the problem and solution circles in suitable colors via `pyganja`.\n", + "Note that this all works because of our definition of the `to_conformal` `BasisVectorMap`." ] }, {