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):