Simplify the implementation of ReciprocalFrame

This is fewer lines and more comments

Author: eric-wieser

galgebra/ga.py

@@ -5,10 +5,10 @@
 import operator
 import copy
 from collections import OrderedDict
-from itertools import combinations
+import itertools
 import functools
 from functools import reduce
-from typing import Tuple, TypeVar, Callable, Mapping
+from typing import Tuple, TypeVar, Callable, Mapping, Sequence
 
 from sympy import (
     diff, Rational, Symbol, S, Mul, Add,
@@ -27,6 +27,10 @@
 zero = S(0)
 
 
+# needed to avoid clashes below
+_mv = mv
+
+
 def all_same(items):
     return all(x == items[0] for x in items)
 
@@ -832,7 +836,7 @@ def _build_bases(self):
         # index list for multivector bases and blades by grade
         basis_indexes = tuple(self.n_range)
         self.indexes = GradedTuple(
-            tuple(combinations(basis_indexes, i))
+            tuple(itertools.combinations(basis_indexes, i))
             for i in range(len(basis_indexes) + 1)
         )
 
@@ -1978,7 +1982,8 @@ def connection(self, rbase, key_base, mode, left):
             self.connect[mode_key].append((key, C))
         return C
 
-    def ReciprocalFrame(self, basis, mode='norm'):
+    # need _mv as mv would refer to the method!
+    def ReciprocalFrame(self, basis: Sequence[_mv.Mv], mode: str = 'norm') -> Tuple[_mv.Mv]:
         """
         Compute the reciprocal frame of a set of vectors
 
@@ -1990,51 +1995,43 @@ def ReciprocalFrame(self, basis, mode='norm'):
           normalization coefficient should be appended to the returned tuple.
           One can divide by this coefficient to normalize the vectors.
         """
-        dim = len(basis)
-
-        indexes = tuple(range(dim))
-        index = [()]
-
-        for i in indexes[-2:]:
-            index.append(tuple(combinations(indexes, i + 1)))
-
-        MFbasis = []
-
-        for igrade in index[-2:]:
-            grade = []
-            for iblade in igrade:
-                blade = self.mv(S(1), 'scalar')
-                for ibasis in iblade:
-                    blade ^= basis[ibasis]
-                blade = blade.trigsimp()
-                grade.append(blade)
-            MFbasis.append(grade)
-        E = MFbasis[-1][0]
-        E_sq = trigsimp((E * E).scalar())
 
-        duals = copy.copy(MFbasis[-2])
+        def wedge_reduce(mvs):
+            """ wedge together a list of multivectors """
+            if not mvs:
+                return self.mv(S(1), 'scalar')
+            return functools.reduce(operator.xor, mvs).trigsimp()
+
+        E = wedge_reduce(basis)
+
+        # elements are such that `basis[i] ^ co_basis[i] == E`
+        co_basis = [
+            sign * wedge_reduce(basis_subset)
+            for sign, basis_subset in zip(
+                # alternating signs
+                itertools.cycle([S(1), S(-1)]),
+                # tuples with one basis missing
+                itertools.combinations(basis, len(basis) - 1),
+            )
+        ]
 
-        duals.reverse()
-        sgn = S(1)
-        rbasis = []
-        for dual in duals:
-            recpv = (sgn * dual * E).trigsimp()
-            rbasis.append(recpv)
-            sgn = -sgn
+        # take the dual without normalization
+        r_basis = [(co_base * E).trigsimp() for co_base in co_basis]
 
+        # normalize
+        E_sq = trigsimp((E * E).scalar())
         if mode == 'norm':
-            for i in range(dim):
-                rbasis[i] = rbasis[i] / E_sq
+            r_basis = [r_base / E_sq for r_base in r_basis]
         else:
             if mode != 'append':
                 # galgebra 0.5.0
                 warnings.warn(
                     "Mode {!r} not understood, falling back to {!r}  but this "
                     "is deprecated".format(mode, 'append'),
                     DeprecationWarning, stacklevel=2)
-            rbasis.append(E_sq)
+            r_basis.append(E_sq)
 
-        return tuple(rbasis)
+        return tuple(r_basis)
 
     def Mlt(self, *args, **kwargs):
         return lt.Mlt(args[0], self, *args[1:], **kwargs)