Features from meuns/galgebra: GA4CS examples, PGA, code-generation, and more

galgebra/.gitignore

@@ -0,0 +1 @@
+*.pyc

galgebra/ga.py

@@ -419,7 +419,7 @@ def __init__(self, bases, **kwargs):
         metric.Metric.__init__(self, bases, **kwargs)
 
         self.par_coords = None
-        self._build_bases()
+        self._build_bases(kwargs.get('sign_and_indexes', None))
         self.dot_mode = '|'
         self._build_basis_product_tables()
 
@@ -513,6 +513,9 @@ def E(self):  # Unnoromalized pseudo-scalar
     def I(self):  # Noromalized pseudo-scalar
         return self.i
 
+    def I_inv(self):
+        return self.i_inv
+
     @property
     def mv_I(self):
         # This exists for backwards compatibility. Note this is not `I()`!
@@ -700,7 +703,37 @@ def _build_basis_blade_symbol(self, base_index):
                 raise ValueError('!!!!No unique root symbol for wedge_print = False!!!!')
         return Symbol(symbol_str, commutative=False)
 
-    def _build_bases(self):
+    def build_cobases(self, coindexes=None):
+        """
+        Cobases for building Poincare duality, this is useful for defining wedge and vee without using I nor any metric.
+        """
+        # TODO: check this can be used with another GA than 3D PGA...
+        if coindexes is None:
+            raise NotImplementedError('!!!!We should provide a default implementation!!!!')
+        else:
+            self.coindexes = coindexes
+            self.coindexes_lst = [index for cograde_index in coindexes for index in cograde_index]
+
+        n = self.n
+
+        self.coblades_lst = []
+        for cograde_index in self.coindexes:
+            k = len(cograde_index[0]) if len(cograde_index) > 0 else 0
+            for cobase_index in cograde_index:
+                coblade_sign = -1 if k == n - 1 and k % 2 == 1 else 1
+                coblade = coblade_sign * self._build_basis_blade_symbol(cobase_index)
+                self.coblades_lst.append(coblade)
+
+        self.coblades_inv_lst = []
+        for grade_index in self.indexes:
+            k = len(grade_index[0]) if len(grade_index) > 0 else 0
+            for base_index in grade_index:
+                coblade_inv_sign = -1 if k == 1 and k % 2 == 1 else 1
+                coblade_inv = coblade_inv_sign * self._build_basis_blade_symbol(base_index)
+                self.coblades_inv_lst.append(coblade_inv)
+        self.coblades_inv_lst = list(reversed(self.coblades_inv_lst))
+
+    def _build_bases(self, sign_and_indexes=None):
         r"""
         The bases for the multivector (geometric) algebra are formed from
         all combinations of the bases of the vector space and the scalars.
@@ -736,14 +769,18 @@ def _build_bases(self):
         """
 
         # index list for multivector bases and blades by grade
-        basis_indexes = tuple(self.n_range)
-        self.indexes = []
-        self._all_indexes_lst = []
-        for i in range(len(basis_indexes) + 1):
-            base_tuple = tuple(combinations(basis_indexes, i))
-            self.indexes.append(base_tuple)
-            self._all_indexes_lst += list(base_tuple)
-        self.indexes = tuple(self.indexes)
+        if sign_and_indexes is None:
+            basis_indexes = tuple(self.n_range)
+            self.indexes = []
+            self._all_indexes_lst = []
+            for i in range(len(basis_indexes) + 1):
+                base_tuple = tuple(combinations(basis_indexes, i))
+                self.indexes.append(base_tuple)
+                self._all_indexes_lst += list(base_tuple)
+            self.indexes = tuple(self.indexes)
+        else:
+            self.indexes = sign_and_indexes[1]
+            self._all_indexes_lst = [index for grade_index in self.indexes for index in grade_index]
 
         # list of non-commutative symbols for multivector bases and blades
         # by grade and as a flattened list
@@ -762,9 +799,21 @@ def _build_bases(self):
 
         self.blades_to_indexes = []
         self.indexes_to_blades = []
-        for (index, blade) in zip(self._all_indexes_lst, self._all_blades_lst):
-            self.blades_to_indexes.append((blade, index))
-            self.indexes_to_blades.append((index, blade))
+
+        if sign_and_indexes is None:
+            for (index, blade) in zip(self._all_indexes_lst, self._all_blades_lst):
+                self.blades_to_indexes.append((blade, (1, index)))
+                self.indexes_to_blades.append((index, blade))
+        else:
+            basis_indexes = tuple(self.n_range)
+            default_indexes_lst = []
+            for i in range(len(basis_indexes) + 1):
+                base_tuple = tuple(combinations(basis_indexes, i))
+                default_indexes_lst += list(base_tuple)
+            signs_lst = [sign for grade_sign in sign_and_indexes[0] for sign in grade_sign]
+            for (default_index, sign, blade) in zip(default_indexes_lst, signs_lst, self._all_blades_lst):
+                self.blades_to_indexes.append((blade, (sign, default_index)))
+                self.indexes_to_blades.append((default_index, sign * blade))
         self.blades_to_indexes_dict = OrderedDict(self.blades_to_indexes)
         self.indexes_to_blades_dict = OrderedDict(self.indexes_to_blades)
 
@@ -908,11 +957,11 @@ def geometric_product_basis_blades(self, blade12):
         # geometric (*) product for orthogonal basis
         if self.is_ortho:
             (blade1, blade2) = blade12
-            index1 = self.blades_to_indexes_dict[blade1]
-            index2 = self.blades_to_indexes_dict[blade2]
+            sign1, index1 = self.blades_to_indexes_dict[blade1]
+            sign2, index2 = self.blades_to_indexes_dict[blade2]
             blade_index = list(index1 + index2)
             repeats = []
-            sgn = 1
+            sgn = sign1 * sign2
             for i in range(1, len(blade_index)):
                 save = blade_index[i]
                 j = i
@@ -1028,7 +1077,7 @@ def reduce_basis_loop(g, blst):
                 else:
                     blst1_flg = False  # more revision needed
                 return a1, blst1, blst1_flg, blst
- 
+
         return True  # revision complete, blst in normal order
 
     #******************* Outer/wedge (^) product **********************#
@@ -1052,15 +1101,15 @@ def wedge_product_basis_blades(self, blade12):  # blade12 = blade1*blade2
         # outer (^) product of basis blades
         # this method works for both orthogonal and non-orthogonal basis
         (blade1, blade2) = blade12
-        index1 = self.blades_to_indexes_dict[blade1]
-        index2 = self.blades_to_indexes_dict[blade2]
+        sign1, index1 = self.blades_to_indexes_dict[blade1]
+        sign2, index2 = self.blades_to_indexes_dict[blade2]
         index12 = list(index1 + index2)
 
         if len(index12) > self.n:
             return 0
         (sgn, wedge12) = Ga.blade_reduce(index12)
         if sgn != 0:
-            return(sgn * self.indexes_to_blades_dict[tuple(wedge12)])
+            return(sgn * sign1 * sign2 * self.indexes_to_blades_dict[tuple(wedge12)])
         else:
             return 0
 
@@ -1070,8 +1119,8 @@ def dot_product_basis_blades(self, blade12, mode):
         # dot (|), left (<), and right (>) products
         # dot product for orthogonal basis
         (blade1, blade2) = blade12
-        index1 = self.blades_to_indexes_dict[blade1]
-        index2 = self.blades_to_indexes_dict[blade2]
+        sign1, index1 = self.blades_to_indexes_dict[blade1]
+        sign2, index2 = self.blades_to_indexes_dict[blade2]
         index = list(index1 + index2)
         grade1 = len(index1)
         grade2 = len(index2)
@@ -1087,7 +1136,7 @@ def dot_product_basis_blades(self, blade12, mode):
             if grade < 0:
                 return 0
         n = len(index)
-        sgn = 1
+        sgn = sign1 * sign2
         result = 1
         ordered = False
         while n > grade:
@@ -1207,7 +1256,7 @@ def _build_base_blade_conversions(self):
 
         # expand blade basis in terms of base basis
         for blade in self._all_blades_lst:
-            index = self.blades_to_indexes_dict[blade]
+            sign, index = self.blades_to_indexes_dict[blade]
             grade = len(index)
             if grade <= 1:
                 blade_expansion.append(blade)

galgebra/generator.py

@@ -0,0 +1,122 @@
+# -*- coding: utf-8 -*-
+
+from sympy import Symbol
+
+from .ga import Ga
+from .mv import J, Jinv
+
+
+def create_multivector(GA, name):
+    blades = [1] + GA.blades_lst
+    mv = GA.mv(0, 'scalar')
+    for blade_index, blade in enumerate(blades):
+        mv += Symbol('{name}[{i}]'.format(name=name, i=blade_index)) * blade
+    return mv
+
+
+_CLASS_TEMPLATE = '''# -*- coding: utf-8 -*-
+
+
+class FlatMv(object):
+    def __init__(self, coefs):
+        assert len(coefs) == {class_blade_count}
+        self.coefs = coefs
+
+    def __getitem__(self, index):
+        return self.coefs[index]
+'''
+
+_BINARY_OPERATOR_TEMPLATE = '''
+    def __{op_name}__(self, other):
+        x = self.coefs
+        y = other.coefs
+        return FlatMv([
+            {op_list}
+        ])
+'''
+
+_UNARY_METHOD_TEMPLATE = '''
+    def {method_name}(self):
+        x = self.coefs
+        return FlatMv([
+            {method_list}
+        ])
+'''
+
+_BINARY_METHOD_TEMPLATE = '''
+    def {method_name}(self, other):
+        x = self.coefs
+        y = other.coefs
+        return FlatMv([
+            {method_list}
+        ])
+'''
+
+
+def format_class(class_blade_count):
+    return _CLASS_TEMPLATE.format(class_blade_count=class_blade_count)
+
+
+def format_op_name(name):
+    return name
+
+
+def format_op_list(mv):
+    return ',\n            '.join(str(blade_coef) for blade_coef in mv.blade_coefs())
+
+
+def format_binary_operator(name, mv):
+    return _BINARY_OPERATOR_TEMPLATE.format(op_name=format_op_name(name), op_list=format_op_list(mv))
+
+
+def format_method_name(name):
+    return name
+
+
+def format_method_list(mv):
+    return ',\n            '.join(str(blade_coef) for blade_coef in mv.blade_coefs())
+
+
+def format_unary_method(name, mv):
+    return _UNARY_METHOD_TEMPLATE.format(method_name=format_method_name(name), method_list=format_method_list(mv))
+
+
+def format_binary_method(name, mv):
+    return _BINARY_METHOD_TEMPLATE.format(method_name=format_method_name(name), method_list=format_method_list(mv))
+
+
+def format_geometric_algebra(GA):
+
+    X = create_multivector(GA, 'x')
+    Y = create_multivector(GA, 'y')
+
+    flat_geometric_algebra = format_class(len(GA._all_blades_lst))
+    flat_geometric_algebra += format_binary_operator('add', X + Y)
+    flat_geometric_algebra += format_binary_operator('sub', X - Y)
+    flat_geometric_algebra += format_binary_operator('mul', X * Y)
+    flat_geometric_algebra += format_binary_operator('and', Jinv(J(X) ^ J(Y)))
+    flat_geometric_algebra += format_binary_operator('xor', X ^ Y)
+    flat_geometric_algebra += format_binary_operator('lshift', X << Y)
+    flat_geometric_algebra += format_binary_operator('rshift', X >> Y)
+    flat_geometric_algebra += format_binary_method('meet', Jinv(J(X) ^ J(Y)))
+    flat_geometric_algebra += format_binary_method('join', X ^ Y)
+    flat_geometric_algebra += format_unary_method('rev', X.rev())
+
+    return flat_geometric_algebra
+
+
+def flatten(flat_ga_module, mv):
+    return flat_ga_module.FlatMv([blade_coef for blade_coef in mv.blade_coefs()])
+
+
+def expand(GA, flat_mv):
+    assert len(flat_mv.coefs) == len(GA._all_blades_lst)
+    mv = GA.mv(0, 'scalar')
+    for blade_coef, blade in zip(flat_mv.coefs, GA._all_blades_lst):
+        mv += blade_coef * blade
+    return mv
+
+
+if __name__ == "__main__":
+    GA = Ga('e*1|2|3', g=[1, 1, 1])
+    print(format_geometric_algebra(GA))

galgebra/metric.py

@@ -273,7 +273,8 @@ class Metric(object):
                   'gsym': (None, 'String s to use "det("+s+")" function in reciprocal basis'),
                   'sig': ('e', 'Signature of metric, default is (n,0) a Euclidean metric'),
                   'Isq': ('-', "Sign of square of pseudo-scalar, default is '-'"),
-                  'wedge': (True, 'Use ^ symbol to print basis blades')}
+                  'wedge': (True, 'Use ^ symbol to print basis blades'),
+                  'sign_and_indexes': (None, 'bases indexes')}
 
     @staticmethod
     def dot_orthogonal(V1, V2, g=None):