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

galgebra/.gitignore

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

galgebra/ga.py

@@ -416,7 +416,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.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
+
     def X(self):
         return self.mv(sum([coord*base for (coord, base) in zip(self.coords, self.basis)]))
 
@@ -660,7 +663,54 @@ def parametric(self, coords):
     def basis_vectors(self):
         return tuple(self.basis)
 
-    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:
+            basis_indexes = tuple(self.n_range)
+            self.coindexes = []
+            self.coindexes_lst = []
+            for i in reversed(basis_indexes):
+                base_tuple = tuple(reversed(tuple(combinations(basis_indexes, i + 1))))
+                self.coindexes.append(base_tuple)
+                self.coindexes_lst += list(base_tuple)
+            self.coindexes.append(())
+            self.coindexes = tuple(self.coindexes)
+        else:
+            self.coindexes = coindexes
+            self.coindexes_lst = [index for cograde_index in coindexes for index in cograde_index]
+
+        def format_symbol_name(symbol_index):
+            return (''.join([str(self.basis[i]) + '^' for i in symbol_index]))[:-1]
+
+        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 * Symbol(format_symbol_name(cobase_index), commutative=False)
+                self.coblades_lst.append(coblade)
+        self.coblades_lst0 = self.coblades_lst + [S(1),]
+
+        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 * Symbol(format_symbol_name(base_index), commutative=False)
+                self.coblades_inv_lst.append(coblade_inv)
+        self.coblades_inv_lst = list(reversed(self.coblades_inv_lst))
+        self.coblades_inv_lst0 = self.coblades_inv_lst + [S(1),]
+
+        self.mv_blades_lst0 = []
+        for obj in self.blades_lst0:
+            self.mv_blades_lst0.append(self.mv(obj))
+
+    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.
@@ -696,14 +746,18 @@ def build_bases(self):
         """
 
         # index list for multivector bases and blades by grade
-        basis_indexes = tuple(self.n_range)
-        self.indexes = [()]
-        self.indexes_lst = []
-        for i in basis_indexes:
-            base_tuple = tuple(combinations(basis_indexes, i + 1))
-            self.indexes.append(base_tuple)
-            self.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.indexes_lst = []
+            for i in basis_indexes:
+                base_tuple = tuple(combinations(basis_indexes, i + 1))
+                self.indexes.append(base_tuple)
+                self.indexes_lst += list(base_tuple)
+            self.indexes = tuple(self.indexes)
+        else:
+            self.indexes = sign_and_indexes[1]
+            self.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
@@ -742,9 +796,20 @@ def build_bases(self):
 
         self.blades_to_indexes = []
         self.indexes_to_blades = []
-        for (index, blade) in zip(self.indexes_lst, self.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.indexes_lst, self.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 basis_indexes:
+                base_tuple = tuple(combinations(basis_indexes, i + 1))
+                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.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)
 
@@ -885,11 +950,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
@@ -1030,15 +1095,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
 
@@ -1048,8 +1113,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)
@@ -1065,7 +1130,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:
@@ -1185,7 +1250,7 @@ def base_blade_conversions(self):
 
         # expand blade basis in terms of base basis
         for blade in self.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.blades_lst0))
+    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.blades_lst0)
+    mv = GA.mv(0, 'scalar')
+    for blade_coef, blade in zip(flat_mv.coefs, GA.blades_lst0):
+        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

@@ -288,7 +288,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):