diff --git a/__init__.py b/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/galgebra/.gitignore b/galgebra/.gitignore new file mode 100644 index 00000000..0d20b648 --- /dev/null +++ b/galgebra/.gitignore @@ -0,0 +1 @@ +*.pyc diff --git a/galgebra/ga.py b/galgebra/ga.py index 9946d2c2..2306f920 100644 --- a/galgebra/ga.py +++ b/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) diff --git a/galgebra/generator.py b/galgebra/generator.py new file mode 100644 index 00000000..dd03d17d --- /dev/null +++ b/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)) diff --git a/galgebra/metric.py b/galgebra/metric.py index 4c9b38cc..1ea2ad76 100644 --- a/galgebra/metric.py +++ b/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): diff --git a/galgebra/mv.py b/galgebra/mv.py index 081af205..8e967f2c 100644 --- a/galgebra/mv.py +++ b/galgebra/mv.py @@ -170,6 +170,14 @@ def characterise_Mv(self): # positionally. When python 3.8 is the lowest supported version, we can # switch to using the / syntax from PEP570 + @staticmethod + def _make_blade(ga, __name, __grade, **kwargs): + if isinstance(__name, str) and isinstance(__grade, int): + return reduce(Mv.__xor__, [ga.mv('%s%d' % (__name, i), 'vector') for i in range(__grade)], + ga.mv(1, 'scalar')).obj + else: + raise TypeError("name must be a string and grade must be an int") + @staticmethod def _make_grade(ga, __name_or_coeffs, __grade, **kwargs): """ Make a pure grade multivector. """ @@ -484,7 +492,7 @@ def __mul__(self, A): if isinstance(A, Dop): return A.Mul(self, A, op='*') - if self.is_scalar(): + if self.is_scalar() or A.is_scalar(): return Mv(self.obj * A, ga=self.Ga) if self.is_blade_rep and A.is_blade_rep: @@ -548,6 +556,7 @@ def Mv_str(self): if self.i_grade == 0: return str(self.obj) self.obj = expand(self.obj) + self.char_Mv = False self.characterise_Mv() self.obj = metric.Simp.apply(self.obj) if self.is_blade_rep or self.Ga.is_ortho: @@ -728,7 +737,7 @@ def __xor__(self, A): # wedge (^) product if isinstance(A, Dop): return A.Mul(self, A, op='^') - if self.is_scalar(): + if self.is_scalar() or isinstance(A, numbers.Number) or A.is_scalar(): return self * A self = self.blade_rep() @@ -956,6 +965,20 @@ def get_coefs(self, grade): (coefs, bases) = list(zip(*cb)) return coefs + def J(self): + # TODO: If we pick our basis smartly, we can probably use J in place of Jinv... + obj = 0 + for coef, coblade in zip(self.blade_coefs(self.Ga._all_mv_blades_lst), self.Ga.coblades_lst): + obj += coef * coblade + return Mv(obj, ga=self.Ga) + + def Jinv(self): + # TODO: If we pick our basis smartly, we can probably use J in place of Jinv... + obj = 0 + for coef, coblade_inv in zip(self.blade_coefs(self.Ga._all_mv_blades_lst), self.Ga.coblades_inv_lst): + obj += coef * coblade_inv + return Mv(obj, ga=self.Ga) + def blade_coefs(self, blade_lst=None): """ For a multivector, A, and a list of basis blades, blade_lst return @@ -965,12 +988,10 @@ def blade_coefs(self, blade_lst=None): if blade_lst is None: blade_lst = self.Ga._all_mv_blades_lst - - #print 'Enter blade_coefs blade_lst =', blade_lst, type(blade_lst), [i.is_blade() for i in blade_lst] - - for blade in blade_lst: - if not blade.is_base() or not blade.is_blade(): - raise ValueError("%s expression isn't a basis blade" % blade) + else: + for blade in blade_lst: + if not blade.is_base() or not blade.is_blade(): + raise ValueError("%s expression isn't a basis blade" % blade) blade_lst = [x.obj for x in blade_lst] (coefs, bases) = metric.linear_expand(self.obj) coef_lst = [] @@ -2505,6 +2526,20 @@ def correlation(u, v, dec=3): # Compute the correlation coefficient of vectors return ulocal.dot(vlocal) / (ulocal.norm() * vlocal.norm()). evalf(dec) +def J(A): + if isinstance(A, Mv): + return A.J() + else: + raise ValueError('A not a multivector in J(A)') + + +def Jinv(A): + if isinstance(A, Mv): + return A.Jinv() + else: + raise ValueError('A not a multivector in J(A)') + + def cross(v1, v2): if v1.is_vector() and v2.is_vector() and v1.Ga.name == v2.Ga.name and v1.Ga.n == 3: return -v1.Ga.I() * (v1 ^ v2) diff --git a/test/test_chapter11.py b/test/test_chapter11.py new file mode 100644 index 00000000..5914e457 --- /dev/null +++ b/test/test_chapter11.py @@ -0,0 +1,156 @@ +from .test_utils import TestCase + +from functools import reduce +from sympy import sqrt, Rational, Symbol +from galgebra.ga import Ga +from galgebra.mv import Mv + + +class TestChapter11(TestCase): + + def test11_4(self): + """ + All planes are 3-blades. + """ + GA, e_0, e_1, e_2, e_3 = Ga.build("e*0|1|2|3", g='-1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + + p0 = Symbol('p0', real=True) + q0 = Symbol('q0', real=True) + r0 = Symbol('r0', real=True) + + p = GA.mv((p0, Symbol('p1', real=True), Symbol('p2', real=True), Symbol('p3', real=True)), 'vector') + q = GA.mv((q0, Symbol('q1', real=True), Symbol('q2', real=True), Symbol('q3', real=True)), 'vector') + r = GA.mv((r0, Symbol('r1', real=True), Symbol('r2', real=True), Symbol('r3', real=True)), 'vector') + + p_inf = p.subs({p0: 0}) + q_inf = q.subs({q0: 0}) + r_inf = r.subs({r0: 0}) + + p = p.subs({p0: 1}) + q = q.subs({q0: 1}) + r = r.subs({r0: 1}) + + self.assertEqual(p ^ q ^ r, p ^ (q - p) ^ (r - p)) + self.assertEqual(p ^ q ^ r, p ^ (q_inf - p_inf) ^ (r_inf - p_inf)) + self.assertEqual(p ^ q ^ r, ((p + q + r) / 3) ^ ((p ^ q) + (q ^ r) + (r ^ p))) + + def test11_6(self): + """ + Dual representation. + """ + Ga.dual_mode('Iinv+') + + GA_list = [ + Ga("e*0|1|2", g='-1 0 0, 0 1 0, 0 0 1'), + Ga("e*0|1|2|3", g='-1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1'), + Ga("e*0|1|2|3|5", g='-1 0 0 0 0, 0 1 0 0 0, 0 0 1 0 0, 0 0 0 1 0, 0 0 0 0 1'), + ] + + for GA in GA_list: + e_0 = GA.mv_basis[0] + e_0_inv = e_0.inv() + + Ip = GA.I() + Ip_inv = GA.I_inv() + Ir = e_0_inv < Ip + Ir_inv = Ir.inv() + self.assertEqual(Ip, e_0 ^ Ir) + self.assertEqual(Ip, e_0 * Ir) + + p = GA.mv([1] + [Symbol('p%d' % i, real=True) for i in range(1, GA.n)], 'vector') + + v = [ + GA.mv([0] + [Symbol('q%d' % i, real=True) for i in range(1, GA.n)], 'vector'), + GA.mv([0] + [Symbol('r%d' % i, real=True) for i in range(1, GA.n)], 'vector'), + GA.mv([0] + [Symbol('s%d' % i, real=True) for i in range(1, GA.n)], 'vector'), + GA.mv([0] + [Symbol('t%d' % i, real=True) for i in range(1, GA.n)], 'vector'), + ] + + # We test available finite k-flats + for k in range(1, GA.n): + A = reduce(Mv.__xor__, v[:k]) + X = (p ^ A) + self.assertNotEqual(X, 0) + M = e_0_inv < (e_0 ^ X) + + # Very slow + # d = (e_0_inv < (e_0 ^ X)) / (e_0_inv < X) + # d_inv = d.inv() + + def hat(A): + return ((-1) ** A.pure_grade()) * A + + self.assertEqual(hat(A < Ir_inv), ((-1) ** (GA.n - 1)) * (hat(A) < Ir_inv)) + + Xd = (p ^ A).dual() + self.assertEqual(Xd, (p ^ A) < Ip_inv) + self.assertEqual(Xd, p < (A < Ip_inv)) + self.assertEqual(Xd, p < ((A < Ir_inv) * e_0_inv)) + self.assertEqual(Xd, hat(A < Ir_inv) - e_0_inv * (p < hat(A < Ir_inv))) + # Very slow + # self.assertEquals(Xd, hat(A < Ir_inv) + e_0_inv * hat(M < Ir_inv)) + # self.assertEquals(Xd, (e_0_inv - d_inv) * hat(M < Ir_inv)) + + Ga.dual_mode() + + def test11_12_1(self): + """ + Compute the 2-blades corresponding to the lines gives by the data below. Which of + the lines are the same, considered as weighted oriented elements of geometry, which + are the same as offset subspaces ? + """ + GA, e_0, e_1, e_2, e_3 = Ga.build("e*0|1|2|3", g='-1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + + p = e_0 + e_1 + q = e_0 + e_2 + d = e_2 - e_1 + self.assertEqual(p ^ q, p ^ d) + self.assertEqual(p ^ q, q ^ d) + + r = e_0 + 2 * (e_2 - e_1) + e = 2 * (e_2 - e_1) + s = e_0 + 3 * (e_2 - e_1) + t = 2 * (e_0 + e_2) + self.assertEqual(2 * (p ^ q), p ^ e) + self.assertEqual(2 * (p ^ q), p ^ t) + + def test11_12_2_1(self): + """ + Let an orthonormal coordinate system be given in 3-dimensional Euclidean space. + Compute the support vector of the line with direction u = e1 + 2 e2 - e3, through the point p = e1 - 3 e2. + What is the distance of the line to the origin ? + """ + GA, e_0, e_1, e_2, e_3 = Ga.build("e*0|1|2|3", g='-1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + + L = (e_1 + 2 * e_2 - e_3) ^ (e_0 + e_1 - 3 * e_2) + d = (e_0.inv() < (e_0 ^ L)) / (e_0.inv() < L) + self.assertEqual(d.norm(), sqrt(Rational(35, 6))) + + def test11_12_2_2(self): + """ + Convert the line of the previous exercise into a parametric equation x = p + t * u; express t as a function of x + for a point x on the line. + """ + GA, e_0, e_1, e_2, e_3 = Ga.build("e*0|1|2|3", g='-1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + + u = e_1 + 2 * e_2 - e_3 + p = e_0 + e_1 - 3 * e_2 + L = u ^ p + + t = Symbol('t', real=True) + x_1 = Symbol('x_1', real=True) + x_2 = Symbol('x_2', real=True) + x_3 = Symbol('x_3', real=True) + x = e_0 + x_1 * e_1 + x_2 * e_2 + x_3 * e_3 + + # x(t) + x_t = (e_0 + e_1 - 3 * e_2) + t * (e_1 + 2 * e_2 - e_3) + + # t(x) + t_x = (x - (e_0 + e_1 - 3 * e_2)) / (e_1 + 2 * e_2 - e_3) + + for i in range(11): + t_value = Rational(i, 10) + x_value = x_t.subs({t: t_value}) + self.assertEqual(x_value ^ L, 0) + self.assertEqual(t_x.subs(list(zip([x_1, x_2, x_3], x_value.blade_coefs([e_1, e_2, e_3])))), t_value) diff --git a/test/test_chapter13.py b/test/test_chapter13.py new file mode 100644 index 00000000..ea4b800e --- /dev/null +++ b/test/test_chapter13.py @@ -0,0 +1,164 @@ +from .test_utils import TestCase + +from sympy import Symbol, S, sqrt +from galgebra.ga import Ga + +USE_DIAGONAL_G = True + + +class TestChapter13(TestCase): + + def setUp(self): + """ + Initialize CGA. + """ + if USE_DIAGONAL_G: + # This is way faster with galgebra but o and inf can't be printed... + g = [ + [ 1, 0, 0, 0, 0], + [ 0, 1, 0, 0, 0], + [ 0, 0, 1, 0, 0], + [ 0, 0, 0, 1, 0], + [ 0, 0, 0, 0, -1], + ] + + ga, e, e_1, e_2, e_3, eb = Ga.build('e e1 e2 e3 eb', g=g) + o = S.Half * (e + eb) + inf = eb - e + + else: + g = [ + [ 0, 0, 0, 0, -1], + [ 0, 1, 0, 0, 0], + [ 0, 0, 1, 0, 0], + [ 0, 0, 0, 1, 0], + [-1, 0, 0, 0, 0], + ] + + ga, o, e_1, e_2, e_3, inf = Ga.build('o e1 e2 e3 inf', g=g) + + self.ga = ga + self.o = o + self.e_1 = e_1 + self.e_2 = e_2 + self.e_3 = e_3 + self.inf = inf + + def vector(self, x, y, z): + """ + Make a vector. + """ + return x * self.e_1 + y * self.e_2 + z * self.e_3 + + def point(self, alpha, v): + """ + Make a point. + """ + return alpha * (self.o + v + S.Half * v * v * self.inf) + + def dual_plane(self, n, delta): + """ + Make a dual plane. + """ + return n + delta * self.inf + + def dual_sphere(self, alpha, c, r): + """ + Make a dual sphere. + """ + return alpha * (c - S.Half * r * r * self.inf) + + def dual_im_sphere(self, alpha, c, r): + """ + Make a dual imaginary sphere. + """ + return alpha * (c + S.Half * r * r * self.inf) + + def test_13_1_2(self): + """ + Points as null vectors. + """ + p = self.point(1, self.vector(Symbol('px', real=True), Symbol('py', real=True), Symbol('pz', real=True))) + q = self.point(1, self.vector(Symbol('qx', real=True), Symbol('qy', real=True), Symbol('qz', real=True))) + + self.assertEqual(p | q, -S.Half * q * q + p | q - S.Half * p * p) + self.assertEqual(p | q, -S.Half * (q - p) * (q - p)) + + def test_13_1_3(self): + """ + General vectors represent dual planes and spheres. + """ + alpha = Symbol('alpha', real=True) + p = self.point(alpha, self.vector(Symbol('px', real=True), Symbol('py', real=True), Symbol('pz', real=True))) + self.assertEqual(p | p, S.Zero) + self.assertEqual(self.inf | p, -alpha) + + # Dual plane + nx = Symbol('nx', real=True) + ny = Symbol('ny', real=True) + nz = Symbol('nz', real=True) + p = self.dual_plane(self.vector(nx, ny, nz), Symbol('delta', real=True)) + self.assertEqual(p | p, nx * nx + ny * ny + nz * nz) + self.assertEqual(self.inf | p, S.Zero) + + # Dual sphere + cx = Symbol('cx', real=True) + cy = Symbol('cy', real=True) + cz = Symbol('cz', real=True) + r = Symbol('r', real=True) + + c = self.point(1, self.vector(cx, cy, cz)) + self.assertEqual(c * c, S.Zero) + self.assertEqual(-self.inf | c, S.One) + + s = self.dual_sphere(alpha, c, r) + self.assertEqual(s | s, alpha * alpha * r * r) + self.assertEqual(-self.inf | s, alpha) + + im_s = self.dual_im_sphere(alpha, c, r) + self.assertEqual(im_s | im_s, -alpha * alpha * r * r) + self.assertEqual(-self.inf | im_s, alpha) + + def symbol_vector(self, name='n', normalize=False): + nx = Symbol(name + 'x', real=True) + ny = Symbol(name + 'y', real=True) + nz = Symbol(name + 'z', real=True) + return self.vector(nx, ny, nz) / (sqrt(nx * nx + ny * ny + nz * nz) if normalize else S.One) + + def test_13_2_2(self): + """ + Proper Euclidean motions as even versors : Translations. + """ + delta_1 = Symbol('delta_1', real=True) + delta_2 = Symbol('delta_2', real=True) + n = self.symbol_vector('n', normalize=True) + + Tt = self.dual_plane(n, delta_2) * self.dual_plane(n, delta_1) + + # Even versor + self.assertEqual(Tt * Tt.rev(), 1) + + # Translation + r = Tt * self.o * Tt.inv() + t = self.point(1, 2 * (delta_2 - delta_1) * n) + self.assertEqual(r, t) + + # TODO: This exponential isn't available in galgebra + #Te = (-t * self.inf * S.Half).exp() + #self.assertEqual(Te * Te.rev(), 1) + #self.assertEqual(Te, Tt) + + def test_13_2_3(self): + """ + Proper Euclidean motions as even versors : Rotations in the origin. + """ + n0 = self.symbol_vector('n0') + n1 = self.symbol_vector('n1') + R = n1 * n0 # 2 reflections + Rinv = R.inv() + + p = self.symbol_vector('p') + + p0 = R * self.point(1, p) * Rinv + p1 = self.point(1, R * p * Rinv) + self.assertEqual(p0, p1) diff --git a/test/test_chapter2.py b/test/test_chapter2.py new file mode 100755 index 00000000..efd3ac74 --- /dev/null +++ b/test/test_chapter2.py @@ -0,0 +1,347 @@ +from .test_utils import TestCase + +from functools import reduce +from itertools import product +from sympy import S, Symbol, Matrix, solve, solve_poly_system, cos, sin +from galgebra.ga import Ga +from galgebra.mv import Mv + + +class TestChapter2(TestCase): + + def test2_3_2(self): + """ + Introducing the outer product. + """ + GA = Ga('e*1|2|3') + + a = GA.mv('a', 'vector') + b = GA.mv('b', 'vector') + c = GA.mv('c', 'vector') + alpha = Symbol('alpha', real=True) + + self.assertEqual(a ^ b, -b ^ a) + self.assertEqual(a ^ (alpha * b), alpha * (a ^ b)) + self.assertEqual(GA.mv(alpha, 'scalar') ^ b, alpha * b) + self.assertEqual(a ^ (b + c), (a ^ b) + (a ^ c)) + + def test_2_4_2(self): + """ + Associativity of the outer product and calculating determinants. + """ + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3') + + a_1 = Symbol('a_1', real=True) + a_2 = Symbol('a_2', real=True) + a_3 = Symbol('a_3', real=True) + b_1 = Symbol('b_1', real=True) + b_2 = Symbol('b_2', real=True) + b_3 = Symbol('b_3', real=True) + c_1 = Symbol('c_1', real=True) + c_2 = Symbol('c_2', real=True) + c_3 = Symbol('c_3', real=True) + + a = GA.mv((a_1, a_2, a_3), 'vector') + b = GA.mv((b_1, b_2, b_3), 'vector') + c = GA.mv((c_1, c_2, c_3), 'vector') + + self.assertEqual(a ^ b ^ c, a ^ (b ^ c)) + self.assertEqual(a ^ b ^ c, (a ^ b) ^ c) + + m = Matrix([ + [a_1, b_1, c_1], + [a_2, b_2, c_2], + [a_3, b_3, c_3] + ]) + + self.assertEqual(a ^ b ^ c, m.det() * (e_1 ^ e_2 ^ e_3)) + + def test2_9_1(self): + """ + Blades and grades. Be careful ga.grade_decomposition and Mv.pure_grade can't tell if a multivector is a blade, + but if we know a multivector is a blade we can retrieve its grade. + """ + GA = Ga('e*1|2|3|4|5') + + # Check for k in [0, R.n] + for k in range(GA.n + 1): + Ak = GA.mv('A', 'blade', k) + self.assertEqual(Ak.pure_grade(), k) + grades = GA.grade_decomposition(Ak) + self.assertEqual(len(grades), 1) + self.assertEqual(list(grades.keys())[0], k) + + # Check for k and l in [0, R.n] + for k, l in product(range(GA.n + 1), range(GA.n + 1)): + Ak = GA.mv('A', 'blade', k) + Bl = GA.mv('B', 'blade', l) + C = Ak ^ Bl + self.assertEqual(C.pure_grade(), 0 if C == 0 else k + l) + grades = GA.grade_decomposition(C) + self.assertEqual(len(grades), 1) + self.assertEqual(list(grades.keys())[0], 0 if C == 0 else k + l) + + def test2_9_5(self): + """ + Reversion and grade involution. + """ + GA = Ga('e*1|2|3|4') + + for k in range(1, GA.n + 1): + a = [GA.mv('a%d' % i, 'vector') for i in range(k)] + A = reduce(Mv.__xor__, a) + A_rev = reduce(Mv.__xor__, reversed(a)) + self.assertEqual(A_rev, A.rev()) + self.assertEqual(A_rev, ((-1) ** ((k * (k - 1)) / 2)) * A) + + def test2_12_1_1(self): + """ + Compute the outer products of the following 3-space expressions, + giving the result relative to the basis {1, e_1, e_2, e_3, e_1^e_2, e_1^e_3, e_2^e_3, e_1^e_2^e_3}. + """ + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3') + self.assertEqual((e_1 + e_2) ^ (e_1 + e_3), (-e_1 ^ e_2) + (e_1 ^ e_3) + (e_2 ^ e_3)) + self.assertEqual((e_1 + e_2 + e_3) ^ (2 * e_1), -2 * (e_1 ^ e_2) - 2 * (e_1 ^ e_3)) + self.assertEqual((e_1 - e_2) ^ (e_1 - e_3), (e_1 ^ e_2) - (e_1 ^ e_3) + (e_2 ^ e_3)) + self.assertEqual( + (e_1 + e_2) ^ (0.5 * e_1 + 2 * e_2 + 3 * e_3), 1.5 * (e_1 ^ e_2) + 3 * (e_1 ^ e_3) + 3 * (e_2 ^ e_3)) + self.assertEqual((e_1 ^ e_2) ^ (e_1 + e_3), (e_1 ^ e_2 ^ e_3)) + self.assertEqual((e_1 + e_2) ^ ((e_1 ^ e_2) + (e_2 ^ e_3)), (e_1 ^ e_2 ^ e_3)) + + def test2_12_1_2(self): + """ + Given the 2-blade B = e_1 ^ (e_2 - e_3) that represents a plane, + determine if each of the following vectors lies in that plane. + """ + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3') + B = e_1 ^ (e_2 - e_3) + self.assertEqual(e_1 ^ B, 0) + self.assertNotEqual((e_1 + e_2) ^ B, 0) + self.assertNotEqual((e_1 + e_2 + e_3) ^ B, 0) + self.assertEqual((2 * e_1 - e_2 + e_3) ^ B, 0) + + def test2_12_1_3(self): + """ + What is the area of the parallelogram spanned by the vectors a = e_1 + 2*e_2 + and b = -e_1 - e_2 (relative to the area of e_1 ^ e_2) ? + """ + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3') + a = e_1 + 2 * e_2 + b = -e_1 - e_2 + B = a ^ b + self.assertEqual(B, 1 * (e_1 ^ e_2)) + + def test2_12_1_4(self): + """ + Compute the intersection of the non-homogeneous line L with position vector e_1 + and direction vector e_2, and the line M with position vector e_2 and direction vector + (e_1 + e_2), using 2-blades. + """ + GA, e_1, e_2 = Ga.build('e*1|2') + + # x is defined in the basis {e_1, e_2} + a = Symbol('a', real=True) + b = Symbol('b', real=True) + x = a * e_1 + b * e_2 + + # x lies on L and M (eq. L == 0 and M == 0) + L = (x ^ e_2) - (e_1 ^ e_2) + M = (x ^ (e_1 + e_2)) - (e_2 ^ (e_1 + e_2)) + + # Solve the linear system + R = solve([L.obj, M.obj], a, b, dict=True) # TODO: fix this... + self.assertEqual(len(R), 1) + + # Replace symbols + x = x.subs(R[0]) + + self.assertEqual(x, e_1 + 2 * e_2) + + def test2_12_1_5(self): + """ + Compute (2 + 3 * e_3) ^ (e_1 + (e_2 ^ e_3) using the grade based defining equations of section 2.9.4. + """ + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3') + + A = 2 + 3 * e_3 + B = e_1 + (e_2 ^ e_3) + + C = GA.mv(0, 'scalar') + for k, l in product(range(GA.n + 1), range(GA.n + 1)): + C += A.get_grade(k) ^ B.get_grade(l) + + self.assertEqual(C, (2 * e_1 + 3 * (e_3 ^ e_1) + 2 * (e_2 ^ e_3))) + + def test2_12_2_1(self): + """ + In R2 with Euclidean metric, choose an orthonormal basis {e_1, e_2} in the plane of a and b such that e1 is parallel to a. + Write x = a * e_1 and y = b * (cos(t) * e_1 + sin(t) * e_2), whete t is the angle from a to b. + Evaluate the outer product. What is the geometrical interpretation ? + """ + GA, e_1, e_2 = Ga.build('e*1|2', g='1 0, 0 1') + + # TODO: use alpha, beta and theta instead of a, b and t (it crashes sympy) + a = Symbol('a', real=True) + b = Symbol('b', real=True) + t = Symbol('t', real=True) + x = a * e_1 + y = b * (cos(t) * e_1 + sin(t) * e_2) + B = x ^ y + self.assertEqual(B, (a * b * sin(t) * (e_1 ^ e_2))) + + # Retrieve the parallelogram area from the 2-vector + area = B.norm() + self.assertEqual(area, (a * b * sin(t))) + + # Compute the parallelogram area using the determinant + x = [a, 0] + y = [b * cos(t), b * sin(t)] + area = Matrix([x, y]).det() + self.assertEqual(area, (a * b * sin(t))) + + def test2_12_2_2(self): + """ + """ + GA, e_1, e_2 = Ga.build('e*1|2', g='1 0, 0 1') + + a = Symbol('a', real=True) + b = Symbol('b', real=True) + t = Symbol('t', real=True) + x = a * e_1 + y = b * (cos(t) * e_1 + sin(t) * e_2) + + x_1 = x.blade_coefs([e_1])[0] + x_2 = x.blade_coefs([e_2])[0] + y_1 = y.blade_coefs([e_1])[0] + y_2 = y.blade_coefs([e_2])[0] + + cross = x_1 * y_2 - y_1 * x_2 + + self.assertEqual(cross, (a * b * sin(t))) + + def test2_12_2_4(self): + """ + Solve the linear system of equation using the outer product. + """ + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3') + + a_1 = Symbol('a_1', real=True) + a_2 = Symbol('a_2', real=True) + a_3 = Symbol('a_3', real=True) + b_1 = Symbol('b_1', real=True) + b_2 = Symbol('b_2', real=True) + b_3 = Symbol('b_3', real=True) + c_1 = Symbol('c_1', real=True) + c_2 = Symbol('c_2', real=True) + c_3 = Symbol('c_3', real=True) + a = a_1 * e_1 + a_2 * e_2 + a_3 * e_3 + b = b_1 * e_1 + b_2 * e_2 + b_3 * e_3 + c = c_1 * e_1 + c_2 * e_2 + c_3 * e_3 + + x_1 = Symbol('x_1', real=True) + x_2 = Symbol('x_2', real=True) + x_3 = Symbol('x_3', real=True) + x = x_1 * a + x_2 * b + x_3 * c + + # Solve x_1 + self.assertEqual((x ^ a), (x_2 * (b ^ a) + x_3 * (c ^ a))) + self.assertEqual((x ^ a ^ b), x_3 * (c ^ a ^ b)) + self.assertEqual((x ^ a ^ b) * (c ^ a ^ b).inv(), GA.mv(x_3, 'scalar')) + + # Solve x_2 + self.assertEqual((x ^ b), (x_1 * (a ^ b) + x_3 * (c ^ b))) + self.assertEqual((x ^ b ^ c), x_1 * (a ^ b ^ c)) + self.assertEqual((x ^ b ^ c) * (a ^ b ^ c).inv(), GA.mv(x_1, 'scalar')) + + # Solve x_3 + self.assertEqual((x ^ c), (x_1 * (a ^ c) + x_2 * (b ^ c))) + self.assertEqual((x ^ c ^ a), x_2 * (b ^ c ^ a)) + self.assertEqual((x ^ c ^ a) * (b ^ c ^ a).inv(), GA.mv(x_2, 'scalar')) + + def test2_12_2_5(self): + """ + Consider R4 with basis {e_1, e_2, e_3, e_4}. Show that the 2-vector B = (e_1 ^ e_2) + (e_3 ^ e_4) + is not a 2-blade (i.e., it cannot be written as the outer product of two vectors). + """ + GA, e_1, e_2, e_3, e_4 = Ga.build('e*1|2|3|4') + + # B + B = (e_1 ^ e_2) + (e_3 ^ e_4) + + # C is the product of a and b vectors + a_1 = Symbol('a_1', real=True) + a_2 = Symbol('a_2', real=True) + a_3 = Symbol('a_3', real=True) + a_4 = Symbol('a_4', real=True) + a = a_1 * e_1 + a_2 * e_2 + a_3 * e_3 + a_4 * e_4 + + b_1 = Symbol('b_1', real=True) + b_2 = Symbol('b_2', real=True) + b_3 = Symbol('b_3', real=True) + b_4 = Symbol('b_4', real=True) + b = b_1 * e_1 + b_2 * e_2 + b_3 * e_3 + b_4 * e_4 + + C = a ^ b + + # other coefficients are null + blades = [ + e_1 ^ e_2, e_1 ^ e_3, e_1 ^ e_4, e_2 ^ e_3, e_2 ^ e_4, e_3 ^ e_4, + ] + + C_coefs = C.blade_coefs(blades) + B_coefs = B.blade_coefs(blades) + + # try to solve the system and show there is no solution + system = [ + (C_coef) - (B_coef) for C_coef, B_coef in zip(C_coefs, B_coefs) + ] + + unknowns = [ + a_1, a_2, a_3, a_4, b_1, b_2, b_3, b_4 + ] + + # TODO: use solve if sympy fix it + result = solve_poly_system(system, unknowns) + self.assertTrue(result is None) + + def test2_12_2_6(self): + """ + Show that B = e1 ^ e2 + e3 ^ e4 of the previous exercise doesn't contain any other vector than 0. + """ + GA, e_1, e_2, e_3, e_4 = Ga.build('e*1|2|3|4') + + # B + B = (e_1 ^ e_2) + (e_3 ^ e_4) + + # x + x_1 = Symbol('x_1', real=True) + x_2 = Symbol('x_2', real=True) + x_3 = Symbol('x_3', real=True) + x_4 = Symbol('x_4', real=True) + x = x_1 * e_1 + x_2 * e_2 + x_3 * e_3 + x_4 * e_4 + + # Solve x ^ B = 0 + system = (x ^ B).blade_coefs() + + unknowns = [ + x_1, x_2, x_3, x_4 + ] + + # TODO: use solve if sympy fix it + result = solve_poly_system(system, unknowns) + result = result[0] + + self.assertEqual(result[0], S.Zero) + self.assertEqual(result[1], S.Zero) + self.assertEqual(result[2], S.Zero) + self.assertEqual(result[3], S.Zero) + + def test2_12_2_9(self): + """ + Prove Ak ^ Bl = (-1**kl) Bl ^ Ak. + """ + for GA in [Ga('e*1|2'), Ga('e*1|2|3'), Ga('e*1|2|3|4')]: # , Ga('e*1|2|3|4|5')]: + for k, l in product(range(GA.n + 1), range(GA.n + 1)): + Ak = GA.mv('A', 'blade', k) + Bl = GA.mv('B', 'blade', l) + self.assertEqual(Ak ^ Bl, (-1) ** (k * l) * (Bl ^ Ak)) diff --git a/test/test_chapter3.py b/test/test_chapter3.py new file mode 100644 index 00000000..4bd4bd01 --- /dev/null +++ b/test/test_chapter3.py @@ -0,0 +1,468 @@ +from .test_utils import TestCase + +from itertools import product +from sympy import S, Symbol, cos, sin, sqrt +from galgebra.ga import Ga +from galgebra.mv import cross + + +class TestChapter3(TestCase): + + def test_3_1_2(self): + """ + Definition of the scalar product. + """ + GA = Ga('e*1|2|3|4') + A_blades = [(k, GA.mv('A', 'blade', k)) for k in range(GA.n + 1)] + B_blades = [(l, GA.mv('B', 'blade', l)) for l in range(GA.n + 1)] + + # GAlgebra doesn't define any scalar product but rely on the geometric product instead + for (k, A), (l, B) in product(A_blades, B_blades): + if k == l: + self.assertNotEqual((A * B).scalar(), S.Zero) + else: + self.assertEqual((A * B).scalar(), S.Zero) + + def test_3_2_2(self): + """ + Computing the contraction explicitly. + """ + GA = Ga('e*1|2|3') + A_blades = [(k, GA.mv('A', 'blade', k)) for k in range(GA.n + 1)] + B_blades = [(l, GA.mv('B', 'blade', l)) for l in range(GA.n + 1)] + C_blades = [(m, GA.mv('C', 'blade', m)) for m in range(GA.n + 1)] + + # scalar and blades of various grades + k, A = A_blades[0] + for l, B in B_blades: + self.assertEqual(A < B, A * B) + + k, A = A_blades[0] + for l, B in B_blades: + self.assertEqual(B < A, 0 if l > 0 else (A * B).scalar()) + + # vectors + k, A = A_blades[1] + l, B = B_blades[1] + self.assertEqual(A < B, (A * B).scalar()) + + # vector and the outer product of 2 blades of various grades (scalars, vectors, 2-vectors...) + k, A = A_blades[1] + for (l, B), (m, C) in product(B_blades, C_blades): + self.assertEqual(A < (B ^ C), ((A < B) ^ C) + (-1) ** l * (B ^ (A < C))) + + # vector and the outer product of 2 blades of various grades (scalars, vectors, 2-vectors...) + for (k, A), (l, B), (m, C) in product(A_blades, B_blades, C_blades): + self.assertEqual((A ^ B) < C, A < (B < C)) + + # distributive properties + for (k, A), (l, B), (m, C) in product(A_blades, B_blades, C_blades): + self.assertEqual((A + B) < C, (A < C) + (B < C)) + + for (k, A), (l, B), (m, C) in product(A_blades, B_blades, C_blades): + self.assertEqual(A < (B + C), (A < B) + (A < C)) + + alpha = Symbol('alpha', real=True) + for (k, A), (l, B) in product(A_blades, B_blades): + self.assertEqual((alpha * A) < B, alpha * (A < B)) + self.assertEqual((alpha * A) < B, A < (alpha * B)) + + a = GA.mv('a', 'blade', 1) + for (k, A_minus1), (l, B) in product(A_blades[:-1], B_blades): + A = A_minus1 ^ a + self.assertEqual(A < B, (A_minus1 ^ a) < B) + self.assertEqual(A < B, A_minus1 < (a < B)) + + def test_3_4(self): + """ + The other contraction. + """ + GA = Ga('e*1|2|3') + A_grade_and_blades = [(k, GA.mv('A', 'blade', k)) for k in range(GA.n + 1)] + B_grade_and_blades = [(l, GA.mv('B', 'blade', l)) for l in range(GA.n + 1)] + + for (k, A), (l, B) in product(A_grade_and_blades, B_grade_and_blades): + self.assertEqual(B > A, ((-1) ** (k * (l - 1))) * (A < B)) + + for (k, A), (l, B) in product(A_grade_and_blades, B_grade_and_blades): + C = B > A + C_grades = GA.grade_decomposition(C).keys() + self.assertEqual(len(C_grades), 1) + self.assertTrue(C == 0 or C_grades[0] == l - k) + + def test_3_5_2(self): + """ + The inverse of a blade. + """ + Ga.dual_mode("Iinv+") + + for GA in [Ga('e*1|2'), Ga('e*1|2|3'), Ga('e*1|2|3|4')]: # , Ga('e*1|2|3|4|5')]: + A_grade_and_blades = [(k, GA.mv('A', 'blade', k)) for k in range(GA.n + 1)] + for k, A in A_grade_and_blades: + inv_A = A.inv() + rev_A = A.rev() + rev_sign = ((-1) ** (k * (k - 1) / 2)) + norm2 = A.norm2() + self.assertEqual(rev_A, rev_sign * A) + self.assertEqual(inv_A, rev_A / norm2) + # We compute the scalar product using the geometric product + self.assertEqual(inv_A, rev_A / (A * rev_A).scalar()) + self.assertEqual(inv_A, rev_sign * (A / norm2)) + + for k, A in A_grade_and_blades: + self.assertEqual(A < A.inv(), 1) + + Ga.dual_mode() + + def test_3_5_3(self): + """ + Orthogonal complement and duality. + """ + Ga.dual_mode("Iinv+") + + # some blades by grades for each space + spaces = [([GA.mv('A', 'blade', k) for k in range(GA.n + 1)], GA) for GA in + [Ga('e*1|2'), Ga('e*1|2|3'), Ga('e*1|2|3|4')]] + + # dualization + for blades, GA in spaces: + for A in blades: + self.assertEqual(A.dual(), A < GA.I_inv()) + + # dualization sign + for blades, R in spaces: + for A in blades: + self.assertEqual(A.dual().dual(), ((-1) ** (GA.n * (GA.n - 1) / 2)) * A) + + # undualization + for blades, R in spaces: + for A in blades: + self.assertEqual(A, A.dual() < GA.I()) + + Ga.dual_mode() + + def test_3_5_4(self): + """ + The duality relationships. + """ + Ga.dual_mode("Iinv+") + + GA = Ga('e*1|2|3') + A_blades = [GA.mv('A', 'blade', k) for k in range(GA.n + 1)] + B_blades = [GA.mv('B', 'blade', l) for l in range(GA.n + 1)] + + for A, B in product(A_blades, B_blades): + self.assertEqual((A ^ B).dual(), A < B.dual()) + + for A, B in product(A_blades, B_blades): + self.assertEqual((A < B).dual(), A ^ B.dual()) + + Ga.dual_mode() + + def test_3_6(self): + """ + Orthogonal projection of subspaces. + """ + GA = Ga('e*1|2|3') + X_blades = [GA.mv('X', 'blade', k) for k in range(GA.n + 1)] + B_blades = [GA.mv('B', 'blade', l) for l in range(GA.n + 1)] + + # projection of X on B + def P(X, B): + return (X < B.inv()) < B + + # a projection should be idempotent + for X, B in product(X_blades, B_blades): + self.assertEqual(P(X, B), P(P(X, B), B)) + + # with the contraction + for X, B in product(X_blades, B_blades): + self.assertEqual(X < B, P(X, B) < B) + + def test3_7_2(self): + """ + The cross product incorporated. + """ + Ga.dual_mode("Iinv+") + + GA = Ga('e*1|2|3') + a = GA.mv('a', 'vector') + b = GA.mv('b', 'vector') + + A = a < GA.I() + B = b < GA.I() + + # Cross product definition + self.assertEqual(GA.I_inv(), -GA.I()) + self.assertEqual(cross(a, b), (a ^ b).dual()) + + # About velocities + self.assertEqual(cross(a, b), (b ^ a) < GA.I()) + self.assertEqual(cross(a, b), (a ^ b) < GA.I_inv()) + self.assertEqual(cross(a, b), -(b ^ a).dual()) + self.assertEqual(cross(a, b), -b < a.dual()) + self.assertEqual(cross(a, b), b < A) + + # Intersecting planes + self.assertEqual(cross(a, b), ((A < GA.I_inv()) ^ (B < GA.I_inv())) < GA.I_inv()) + self.assertEqual(cross(a, b), (B < GA.I_inv()) < ((A < GA.I_inv()) < GA.I())) + self.assertEqual(cross(a, b), (B < GA.I_inv()) < A) + self.assertEqual(cross(a, b), B.dual() < A) + + Ga.dual_mode() + + def test_3_10_1_1(self): + """ + Let a = e_1 + e_2 and b = e_2 + e_3 in a 3D Euclidean space R(3,0) with an orthonormal basis {e_1, e_2, e_3}. + Compute the following expressions, giving the results relative to the basis + {1, e_1, e_2, e_3, e_1 ^ e_2, e_2 ^ e_3, e_3 ^ e_1, e_1 ^ e_2 ^ e_3}. + """ + Ga.dual_mode("Iinv+") + + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3', g='1 0 0, 0 1 0, 0 0 1') + a = e_1 + e_2 + b = e_2 + e_3 + + # a + self.assertEqual(e_1 < a, (e_1 < e_1) + (e_1 < e_2)) + self.assertEqual(e_1 < e_1, e_1 | e_1) + self.assertEqual(e_1 < e_1, 1) + self.assertEqual(e_1 < e_2, e_1 | e_2) + self.assertEqual(e_1 < e_2, 0) + self.assertEqual(e_1 < a, 1) + + # b + self.assertEqual(e_1 < (a ^ b), ((e_1 < a) ^ b) - (a ^ (e_1 < b))) + self.assertEqual((e_1 < a) ^ b, b) # reusing drill a + self.assertEqual(e_1 < b, (e_1 < e_2) + (e_1 < e_3)) + self.assertEqual(e_1 < b, 0) + self.assertEqual(e_1 < (a ^ b), b) + self.assertEqual(e_1 < (a ^ b), e_2 + e_3) + + # c + self.assertEqual((a ^ b) < e_1, a < (b < e_1)) + self.assertEqual((a ^ b) < e_1, a < ((e_2 < e_1) + (e_3 < e_1))) + self.assertEqual((a ^ b) < e_1, 0) + + # d + self.assertEqual(((2 * a) + b) < (a + b), ((2 * a) < a) + ((2 * a) < b) + (b < a) + (b < b)) + self.assertEqual(((2 * a) + b) < (a + b), ((2 * (a < a)) + (2 * (a < b)) + (b < a) + (b < b))) + self.assertEqual(((2 * a) + b) < (a + b), 4 + 2 + 1 + 2) + + # e + self.assertEqual(a < (e_1 ^ e_2 ^ e_3), a < ((e_1 ^ e_2) ^ e_3)) + self.assertEqual(a < (e_1 ^ e_2 ^ e_3), ((a < (e_1 ^ e_2)) ^ e_3) + ((e_1 ^ e_2) ^ (a < e_3))) + self.assertEqual(a < (e_1 ^ e_2 ^ e_3), (((e_1 < (e_1 ^ e_2)) + (e_2 < (e_1 ^ e_2))) ^ e_3) + ( + (e_1 ^ e_2) ^ ((e_1 < e_3) + (e_2 < e_3)))) + self.assertEqual(((e_1 < (e_1 ^ e_2)) + (e_2 < (e_1 ^ e_2))) ^ e_3, + (((e_1 < e_1) ^ e_2) - (e_1 ^ (e_1 < e_2)) + ((e_2 < e_1) ^ e_2) - (e_1 ^ (e_2 < e_2))) ^ e_3) + self.assertEqual(((e_1 < (e_1 ^ e_2)) + (e_2 < (e_1 ^ e_2))) ^ e_3, (e_2 - e_1) ^ e_3) + self.assertEqual(((e_1 < (e_1 ^ e_2)) + (e_2 < (e_1 ^ e_2))) ^ e_3, (e_2 ^ e_3) - (e_1 ^ e_3)) + self.assertEqual(e_1 < e_3, 0) + self.assertEqual(e_2 < e_3, 0) + self.assertEqual(((e_1 ^ e_2) ^ ((e_1 < e_3) + (e_2 < e_3))), 0) + self.assertEqual(a < (e_1 ^ e_2 ^ e_3), (e_2 ^ e_3) + (e_3 ^ e_1)) + + # f + self.assertEqual(a < GA.I_inv(), a < (-e_1 ^ e_2 ^ e_3)) # equals because we fixed the metric + self.assertEqual(a < GA.I_inv(), a < ((e_1 ^ e_2) ^ (-e_3))) # almost last drill + self.assertEqual(a < GA.I_inv(), -(e_2 ^ e_3) - (e_3 ^ e_1)) + + # g + self.assertEqual((a ^ b) < GA.I_inv(), -((b ^ a) < GA.I_inv())) + self.assertEqual((a ^ b) < GA.I_inv(), -(b < (a < GA.I_inv()))) + self.assertEqual((a ^ b) < GA.I_inv(), -((e_2 + e_3) < (-(e_2 ^ e_3) + (e_1 ^ e_3)))) + self.assertEqual((a ^ b) < GA.I_inv(), ((e_2 + e_3) < (e_2 ^ e_3)) - ((e_2 + e_3) < (e_1 ^ e_3))) + self.assertEqual((a ^ b) < GA.I_inv(), + ((e_2 < (e_2 ^ e_3)) + (e_3 < (e_2 ^ e_3)) - (e_2 < (e_1 ^ e_3)) - (e_3 < (e_1 ^ e_3)))) + self.assertEqual(e_2 < (e_2 ^ e_3), e_3) + self.assertEqual(e_3 < (e_2 ^ e_3), -e_2) + self.assertEqual(e_2 < (e_1 ^ e_3), 0) + self.assertEqual(e_3 < (e_1 ^ e_3), -e_1) + self.assertEqual((a ^ b) < GA.I_inv(), e_3 - e_2 + e_1) + + # h + self.assertEqual(a < (b < GA.I_inv()), (a ^ b) < GA.I_inv()) + self.assertEqual(a < (b < GA.I_inv()), e_3 - e_2 + e_1) + + Ga.dual_mode() + + def test_3_10_1_2(self): + """ + Compute the cosine of the angle between the following subspaces given on an orthonormal basis of a Euclidean space. + """ + Ga.dual_mode("Iinv+") + + GA, e_1, e_2, e_3, e_4 = Ga.build('e*1|2|3|4', g='1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + alpha = Symbol('alpha', real=True) + theta = Symbol('theta', real=True) + + # e_1 and alpha * e_1 + A = e_1 + B = alpha * e_1 + cosine = (A | B.rev()) / (A.norm() * B.norm()) + self.assertEqual(A | B.rev(), e_1 | (alpha * e_1)) + self.assertEqual(A.norm() * B.norm(), alpha) + self.assertEqual(cosine, (e_1 | (alpha * e_1)) / alpha) + self.assertEqual(cosine, (alpha * (e_1 | e_1)) / alpha) + self.assertEqual(cosine, 1) + + # (e_1 + e_2) ^ e_3 and e_1 ^ e_3 + A = (e_1 + e_2) ^ e_3 + B = e_1 ^ e_3 + num = A | B.rev() + den = A.norm() * B.norm() + cosine = num / den + self.assertEqual(num, ((e_1 ^ e_3) + (e_2 ^ e_3)) | (e_3 ^ e_1)) + self.assertEqual(num, ((e_1 ^ e_3) | (e_3 ^ e_1)) + ((e_2 ^ e_3) | (e_3 ^ e_1))) + self.assertEqual(((e_1 ^ e_3) | (e_3 ^ e_1)), 1) + self.assertEqual(((e_2 ^ e_3) | (e_3 ^ e_1)), 0) + self.assertEqual(num, 1) + self.assertEqual(den, ((e_1 + e_2) ^ e_3).norm() * (e_1 ^ e_3).norm()) + self.assertEqual(den, ((e_1 ^ e_3) + (e_2 ^ e_3)).norm()) + den2 = ((e_1 ^ e_3) + (e_2 ^ e_3)).norm2() + self.assertEqual(den2, ((e_1 ^ e_3) + (e_2 ^ e_3)) | ((e_1 ^ e_3) + (e_2 ^ e_3)).rev()) + self.assertEqual(den2, ((e_1 ^ e_3) + (e_2 ^ e_3)) | ((e_3 ^ e_1) + (e_3 ^ e_2))) + self.assertEqual(den2, + ((e_1 ^ e_3) | (e_3 ^ e_1)) + ((e_1 ^ e_3) | (e_3 ^ e_2)) + ((e_2 ^ e_3) | (e_3 ^ e_1)) + ( + (e_2 ^ e_3) | (e_3 ^ e_2))) + self.assertEqual(den2, 1 + 0 + 0 + 1) + self.assertEqual(cosine, 1 / sqrt(2)) + + # (cos(theta) * e_1 + sin(theta) * e_2) ^ e_3 and e_2 ^ e_3 + A = (cos(theta) * e_1 + sin(theta) * e_2) ^ e_3 + B = e_2 ^ e_3 + num = A | B.rev() + den = A.norm() * B.norm() + cosine = num / den + self.assertEqual(num, (cos(theta) * (e_1 ^ e_3) + sin(theta) * (e_2 ^ e_3)) | (e_3 ^ e_2)) + self.assertEqual(num, ((cos(theta) * (e_1 ^ e_3)) | (e_3 ^ e_2)) + ((sin(theta) * (e_2 ^ e_3)) | (e_3 ^ e_2))) + self.assertEqual((cos(theta) * (e_1 ^ e_3)) | (e_3 ^ e_2), 0) + self.assertEqual((sin(theta) * (e_2 ^ e_3)) | (e_3 ^ e_2), sin(theta)) + self.assertEqual(num, sin(theta)) + self.assertEqual(den, ((cos(theta) * e_1 + sin(theta) * e_2) ^ e_3).norm() * (e_2 ^ e_3).norm()) + self.assertEqual(den, ((cos(theta) * e_1 + sin(theta) * e_2) ^ e_3).norm()) + den2 = ((cos(theta) * e_1 + sin(theta) * e_2) ^ e_3).norm2() + self.assertEqual(den2, ((cos(theta) * e_1 + sin(theta) * e_2) ^ e_3) | ( + (cos(theta) * e_1 + sin(theta) * e_2) ^ e_3).rev()) + self.assertEqual(den2, (cos(theta) * (e_1 ^ e_3) + sin(theta) * (e_2 ^ e_3)) | ( + cos(theta) * (e_1 ^ e_3) + sin(theta) * (e_2 ^ e_3)).rev()) + self.assertEqual(den2, (cos(theta) * (e_1 ^ e_3) + sin(theta) * (e_2 ^ e_3)) | ( + cos(theta) * (e_3 ^ e_1) + sin(theta) * (e_3 ^ e_2))) + self.assertEqual(den2, ((cos(theta) * (e_1 ^ e_3)) | (cos(theta) * (e_3 ^ e_1))) + ( + (cos(theta) * (e_1 ^ e_3)) | (sin(theta) * (e_3 ^ e_2))) + ( + (sin(theta) * (e_2 ^ e_3)) | (cos(theta) * (e_3 ^ e_1))) + ( + (sin(theta) * (e_2 ^ e_3)) | (sin(theta) * (e_3 ^ e_2)))) + self.assertEqual(den2, cos(theta) * cos(theta) + sin(theta) * sin(theta)) + self.assertEqual(den2, 1) + self.assertEqual(cosine, sin(theta)) + + # e_1 ^ e_2 and e_3 ^ e_4 + A = e_1 ^ e_2 + B = e_3 ^ e_4 + num = A | B.rev() + den = A.norm() * B.norm() + cosine = num / den + self.assertEqual(num, (e_1 ^ e_2) | (e_4 ^ e_3)) + self.assertEqual(num, 0) + self.assertEqual(cosine, 0) + + Ga.dual_mode() + + def test3_10_2_1(self): + """ + In 2-D Euclidean space R(2,0) with orthogonal basis {e1, e2}, let us determine the value of the contraction + e1 < (e1 ^ e2) by means of its implicit definition with A = e1 and B = e1 ^ e2. Let X range over the basis of + the blades {1, e1, e2, e1 ^ e2}. This produces four equations, each of which gives you information on the + coefficient of the corresponding basis element in the final result. + Show that e1 < (e1 ^ e2) = (0)1 + 0(e1) + 1(e2) + 0(e1 ^ e2). + """ + GA, e_1, e_2 = Ga.build('e*1|2', g='1 0, 0 1') + A = e_1 + B = e_1 ^ e_2 + + X = GA.mv(1) + self.assertEqual(((X ^ A) * B).scalar(), (X * (A < B)).scalar()) + self.assertEqual(((X ^ A) * B).scalar(), 0) + + X = e_1 + self.assertEqual(((X ^ A) * B).scalar(), (X * (A < B)).scalar()) + self.assertEqual(((X ^ A) * B).scalar(), 0) + + X = e_2 + self.assertEqual(((X ^ A) * B).scalar(), (X * (A < B)).scalar()) + self.assertEqual(((X ^ A) * B).scalar(), 1) + + X = e_1 ^ e_2 + self.assertEqual(((X ^ A) * B).scalar(), (X * (A < B)).scalar()) + self.assertEqual(((X ^ A) * B).scalar(), 0) + + def test3_10_2_2(self): + """ + Change the metric such that e2 . e2 == 0. Show that you can't determine the coefficient of e2 in the value + of e1 < (e1 ^ e2) like the previous exercise. The use the explicit definition of the contraction to show the + contraction is still well defined, and equal to e1 < (e1 ^ e2) == e2. + """ + GA, e_1, e_2 = Ga.build('e*1|2', g='1 0, 0 0') + + # We can't use the scalar product anymore (because of the metric) + A = e_1 + B = e_1 ^ e_2 + X = e_2 + self.assertEqual(((X ^ A) * B).scalar(), (X * (A < B)).scalar()) + self.assertEqual(((X ^ A) * B).scalar(), 0) # Which is false + + e_1_grades = list(GA.grade_decomposition(e_1).keys()) + self.assertEqual(len(e_1_grades), 1) + self.assertEqual(e_1_grades[0], 1) + + # We use the explicit definition of the contraction instead + self.assertEqual(e_1 < (e_1 ^ e_2), ((e_1 < e_1) ^ e_2) + ((-1) ** e_1_grades[0]) * (e_1 ^ (e_1 < e_2))) + # We can't use the definition a < b = a . b using GAlgebra so we solved it ourselves... + self.assertEqual(e_1 < (e_1 ^ e_2), (GA.mv(1) ^ e_2) + ((-1) ** e_1_grades[0]) * (e_1 ^ GA.mv(0))) + self.assertEqual(e_1 < (e_1 ^ e_2), e_2) # Which is true + + def test3_10_2_3(self): + """ + Derive the following dualities for right contraction : + C > (B ^ A) = (C > B) > A and C > (B > A) = (C > B) ^ A when A included in C. + """ + GA = Ga('e*1|2|3') + + A_blades = [GA.mv('X', 'blade', k) for k in range(GA.n + 1)] + B_blades = [GA.mv('B', 'blade', l) for l in range(GA.n + 1)] + C_blades = [GA.mv('B', 'blade', m) for m in range(GA.n + 1)] + + for A, B, C in product(A_blades, B_blades, C_blades): + M = C > (B ^ A) + self.assertEqual(M, ((A.rev() ^ B.rev()) < C.rev()).rev()) + self.assertEqual(M, (A.rev() < (B.rev() < C.rev())).rev()) + self.assertEqual(M, (B.rev() < C.rev()).rev() > A) + self.assertEqual(M, (C > B) > A) + + M = C > (B > A) + # TODO + + def test3_10_2_11(self): + """ + Derive the notorious bac-cab formula for the cross product (a x (b x c) = b (a . c) - c (a . b)), + directly from its definition (3.28). What is the corresponding formula using ^ and < ? + """ + Ga.dual_mode("Iinv+") + + GA = Ga('e*1|2|3') + a = GA.mv('a', 'vector') + b = GA.mv('b', 'vector') + c = GA.mv('c', 'vector') + + xx = cross(a, cross(b, c)) + self.assertEqual(xx, (a ^ (b ^ c).dual()).dual()) + self.assertEqual(xx, a < (b ^ c).dual().dual()) + self.assertEqual(xx, -a < (b ^ c)) + self.assertEqual(xx, ((-a < b) ^ c) - (b ^ (-a < c))) + self.assertEqual(xx, (b ^ (a < c)) - (c ^ (a < b))) + self.assertTrue((a < c).is_scalar()) + self.assertTrue((a < b).is_scalar()) + self.assertEqual(xx, b * (a < c) - c * (a < b)) + + Ga.dual_mode() diff --git a/test/test_chapter6.py b/test/test_chapter6.py new file mode 100644 index 00000000..8ecd238e --- /dev/null +++ b/test/test_chapter6.py @@ -0,0 +1,297 @@ +from .test_utils import TestCase + +from itertools import product +from sympy import Symbol, solve_poly_system +from galgebra.ga import Ga + + +class TestChapter6(TestCase): + + def test6_1_4_1(self): + """ + The geometric product for vectors on a basis. + """ + GA, e_1, e_2, e_3 = Ga.build('e*1|2|3', g='1 0 0, 0 1 0, 0 0 1') + + self.assertEqual(e_1 * e_1, 1) + self.assertEqual(e_1 * e_2, e_1 ^ e_2) + self.assertEqual(e_1 * e_3, e_1 ^ e_3) + self.assertEqual(e_2 * e_1, -e_1 ^ e_2) + self.assertEqual(e_2 * e_2, 1) + self.assertEqual(e_2 * e_3, e_2 ^ e_3) + self.assertEqual(e_3 * e_1, -e_1 ^ e_3) + self.assertEqual(e_3 * e_2, -e_2 ^ e_3) + self.assertEqual(e_3 * e_3, 1) + + e_12 = e_1 * e_2 + self.assertEqual(e_12 * e_12, -1) + + e_13 = e_1 * e_3 + self.assertEqual(e_13 * e_13, -1) + + e_23 = e_2 * e_3 + self.assertEqual(e_23 * e_23, -1) + + def test6_1_4_2(self): + """ + The geometric product for vectors on a basis. + """ + GA, e_1, e_2 = Ga.build('e*1|2', g='1 0, 0 1') + ONE = GA.mv(1, 'scalar') + e_12 = e_1 ^ e_2 + + self.assertEqual(ONE * ONE, 1) + self.assertEqual(ONE * e_1, e_1) + self.assertEqual(ONE * e_2, e_2) + self.assertEqual(ONE * e_12, e_12) + + self.assertEqual(e_1 * ONE, e_1) + self.assertEqual(e_1 * e_1, 1) + self.assertEqual(e_1 * e_2, e_12) + self.assertEqual(e_1 * e_12, e_2) + + self.assertEqual(e_2 * ONE, e_2) + self.assertEqual(e_2 * e_1, -e_12) + self.assertEqual(e_2 * e_2, 1) + self.assertEqual(e_2 * e_12, -e_1) + + self.assertEqual(e_12 * ONE, e_12) + self.assertEqual(e_12 * e_1, -e_2) + self.assertEqual(e_12 * e_2, e_1) + self.assertEqual(e_12 * e_12, -1) + + a_1 = Symbol('a_1', real=True) + a_2 = Symbol('a_2', real=True) + a = GA.mv((a_1, a_2), 'vector') + + b_1 = Symbol('b_1', real=True) + b_2 = Symbol('b_2', real=True) + b = GA.mv((b_1, b_2), 'vector') + + self.assertEqual(a * b, a_1 * b_1 + a_2 * b_2 + (a_1 * b_2 - a_2 * b_1) * (e_1 ^ e_2)) + + def test6_3_1(self): + """ + The subspace products from symmetry. + """ + GA = Ga('e*1|2|3') + a = GA.mv('a', 'vector') + B_grade_and_blades = [(k, GA.mv('B%d' % k, 'blade', k)) for k in range(GA.n + 1)] + + def hat(k, M): + return ((-1) ** k) * M + + for k, B in B_grade_and_blades: + self.assertEqual(a ^ B, 0.5 * (a * B + hat(k, B) * a)) + self.assertEqual(B ^ a, 0.5 * (B * a + a * hat(k, B))) + self.assertEqual(a < B, 0.5 * (a * B - hat(k, B) * a)) + self.assertEqual(B > a, 0.5 * (B * a - a * hat(k, B))) + self.assertEqual(a * B, (a < B) + (a ^ B)) + + for k, B in B_grade_and_blades: + self.assertEqual(hat(k, B) * a, (hat(k, B) > a) + (hat(k, B) ^ a)) + self.assertEqual(hat(k, B) * a, -(a < B) + (a ^ B)) + self.assertEqual(hat(k, B) * a, a * B - 2 * (a < B)) + self.assertEqual(hat(k, B) * a, -(a * B) + 2 * (a ^ B)) + + for k, B in B_grade_and_blades: + self.assertEqual(a * B, hat(k, B) * a + 2 * (a < B)) + self.assertEqual(a * B, -hat(k, B) * a + 2 * (a ^ B)) + + M = GA.mv('m', 'mv') + self.assertEqual(a * M, (a < M) + (a ^ M)) + + def test6_3_2(self): + """ + The subspace products as selected grades. + """ + GA = Ga('e*1|2|3') + + # This test should work for blades and other graded multivectors + A_grade_and_blades = [(k, GA.mv('A%d' % k, 'grade', k)) for k in range(GA.n + 1)] + B_grade_and_blades = [(l, GA.mv('B%d' % l, 'grade', l)) for l in range(GA.n + 1)] + + for (k, A), (l, B) in product(A_grade_and_blades, B_grade_and_blades): + self.assertEqual(A ^ B, (A * B).get_grade(k + l)) + self.assertEqual(A < B, 0 if k > l else (A * B).get_grade(l - k)) + self.assertEqual(A > B, 0 if l > k else (A * B).get_grade(k - l)) + # TODO: scalar product + + def test6_6_2_3(self): + """ + The outer product can be defined as the completely antisymmetric summed average of all permutations + of the geometric product of its factors, with a sign for each term depending on oddness or evenness of + the permutation. Derive x ^ y ^ z = 1/3! * (xyz - yxz + yzx - zyx + zxy - xzy) + """ + GA = Ga('e*1|2|3') + x = GA.mv('x', 'vector') + y = GA.mv('y', 'vector') + z = GA.mv('z', 'vector') + + def hat(M): + M_grades = list(GA.grade_decomposition(M).keys()) + self.assertEqual(len(M_grades), 1) + return ((-1) ** M_grades[0]) * M + + self.assertEqual(x * y * z, ((x < y) + (x ^ y)) * z) + self.assertEqual(x * y * z, ((x < y) * z) + ((x ^ y) * z)) + self.assertEqual(x * y * z, ((x < y) * z) - (z < hat(x ^ y)) + (z ^ hat(x ^ y))) + self.assertEqual(hat(x ^ y), x ^ y) + self.assertEqual(x * y * z, ((x < y) * z) - (z < (x ^ y)) + (z ^ x ^ y)) + self.assertEqual(x * y * z, ((x < y) * z) - (z < (x ^ y)) + (x ^ y ^ z)) + + self.assertEqual(y * x * z, ((y < x) + (y ^ x)) * z) + self.assertEqual(y * x * z, ((y < x) * z) + ((y ^ x) * z)) + self.assertEqual(y * x * z, ((y < x) * z) - (z < hat(y ^ x)) + (z ^ hat(y ^ x))) + self.assertEqual(hat(y ^ x), y ^ x) + self.assertEqual(y * x * z, ((y < x) * z) - (z < (y ^ x)) + (z ^ y ^ x)) + self.assertEqual(y * x * z, ((x < y) * z) + (z < (x ^ y)) - (x ^ y ^ z)) + + self.assertEqual(y * z * x, ((y < z) + (y ^ z)) * x) + self.assertEqual(y * z * x, ((y < z) * x) - (x < hat(y ^ z)) + (x ^ hat(y ^ z))) + self.assertEqual(y * z * x, ((y < z) * x) - (x < (y ^ z)) + (x ^ y ^ z)) + + self.assertEqual(z * y * x, ((z < y) + (z ^ y)) * x) + self.assertEqual(z * y * x, ((z < y) * x) - (x < hat(z ^ y)) + (x ^ hat(z ^ y))) + self.assertEqual(z * y * x, ((z < y) * x) - (x < (z ^ y)) - (x ^ y ^ z)) + self.assertEqual(z * y * x, ((y < z) * x) + (x < (y ^ z)) - (x ^ y ^ z)) + + self.assertEqual(z * x * y, ((z < x) + (z ^ x)) * y) + self.assertEqual(z * x * y, ((z < x) * y) + ((z ^ x) * y)) + self.assertEqual(z * x * y, ((z < x) * y) - (y < hat(z ^ x)) + (y ^ hat(z ^ x))) + self.assertEqual(z * x * y, ((z < x) * y) - (y < (z ^ x)) + (y ^ z ^ x)) + self.assertEqual(z * x * y, ((z < x) * y) - (y < (z ^ x)) + (x ^ y ^ z)) + + self.assertEqual(x * z * y, ((x < z) + (x ^ z)) * y) + self.assertEqual(x * z * y, ((x < z) * y) + ((x ^ z) * y)) + self.assertEqual(x * z * y, ((x < z) * y) - (y < hat(x ^ z)) + (y ^ hat(x ^ z))) + self.assertEqual(x * z * y, ((x < z) * y) - (y < (x ^ z)) + (y ^ x ^ z)) + self.assertEqual(x * z * y, ((z < x) * y) + (y < (z ^ x)) - (x ^ y ^ z)) + + self.assertEqual(x * y * z - y * x * z, 2 * ((x ^ y ^ z) - (z < (x ^ y)))) + self.assertEqual(y * z * x - z * y * x, 2 * ((x ^ y ^ z) - (x < (y ^ z)))) + self.assertEqual(z * x * y - x * z * y, 2 * ((x ^ y ^ z) - (y < (z ^ x)))) + + self.assertEqual(z < (x ^ y), ((z < x) ^ y) - (x ^ (z < y))) + self.assertEqual(x < (y ^ z), ((x < y) ^ z) - (y ^ (x < z))) + self.assertEqual(y < (z ^ x), ((y < z) ^ x) - (z ^ (y < x))) + + self.assertEqual(((z < x) ^ y) - (x ^ (z < y)) + ((x < y) ^ z) - (y ^ (x < z)) + ((y < z) ^ x) - (z ^ (y < x)), + 0) + + self.assertEqual(x * y * z - y * x * z + y * z * x - z * y * x + z * x * y - x * z * y, 6 * (x ^ y ^ z)) + + def test6_6_2_4(self): + """ + The parts of a certain grade of a geometric product of blades are not necessarily blades. + Show that in a 4D space with orthogonal basis, a counterexample is the grade 2 of the product + e1 * (e1 + e2) * (e2 + e3) * (e1 + e4). + """ + GA, e_1, e_2, e_3, e_4 = Ga.build('e*1|2|3|4', g='1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + R = (e_1 * (e_1 + e_2) * (e_2 + e_3) * (e_1 + e_4)).get_grade(2) + + # A 2-blade + a_1 = Symbol('a_1', real=True) + a_2 = Symbol('a_2', real=True) + a_3 = Symbol('a_3', real=True) + a_4 = Symbol('a_4', real=True) + a = GA.mv((a_1, a_2, a_3, a_4), 'vector') + + b_1 = Symbol('b_1', real=True) + b_2 = Symbol('b_2', real=True) + b_3 = Symbol('b_3', real=True) + b_4 = Symbol('b_4', real=True) + b = GA.mv((b_1, b_2, b_3, b_4), 'vector') + S = a ^ b + + # Try to solve the system and show there is no solution + system = [ + S_coef - R_coef for S_coef, R_coef in zip(S.blade_coefs(), R.blade_coefs()) + ] + + unknowns = [ + a_1, a_2, a_3, a_4, b_1, b_2, b_3, b_4 + ] + + # TODO: use solve if sympy fix it + result = solve_poly_system(system, unknowns) + self.assertTrue(result is None) + + def test6_6_2_6(self): + """ + Prove (X ^ A) * B = X * (A < B) using the grade based definition of ^, * and <. + """ + GA = Ga('e*1|2|3') + X_grade_blades = [(j, GA.mv('X', 'blade', j)) for j in range(GA.n + 1)] + A_grade_blades = [(k, GA.mv('A', 'blade', k)) for k in range(GA.n + 1)] + B_grade_blades = [(l, GA.mv('B', 'blade', l)) for l in range(GA.n + 1)] + + for (j, X), (k, A), (l, B) in product(X_grade_blades, A_grade_blades, B_grade_blades): + self.assertEqual((X * A).get_grade(j + k), X ^ A) + self.assertEqual((A * B).get_grade(l - k), 0 if k > l else A < B) + self.assertTrue(((X * A).get_grade(j + k) * B).get_grade(0) != 0 if j + k == l else True) + self.assertTrue((X * (A * B).get_grade(l - k)).get_grade(0) != 0 if j == l - k else True) + self.assertEqual(((X * A).get_grade(j + k) * B).get_grade(0), (X * (A * B).get_grade(l - k)).get_grade(0)) + + def test6_6_2_7(self): + """ + In the formula (x < (1/A)) * A, show we can replace the geometric product by a contraction, so that + it is in fact the projection (x < (1/A)) < A. + """ + + GA_list = [ + Ga('e*1|2', g='1 0, 0 1'), + Ga('e*1|2|3', g='1 0 0, 0 1 0, 0 0 1'), + Ga('e*1|2|3|4', g='1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1'), + # Ga('e*1|2|3|4|5', g='1 0 0 0 0, 0 1 0 0 0, 0 0 1 0 0, 0 0 0 1 0, 0 0 0 0 1'), + ] + + for GA in GA_list: + x = GA.mv('x', 'vector') + A_grade_and_blades = [(k, GA.mv('A', 'blade', k)) for k in range(GA.n + 1)] + for k, A in A_grade_and_blades: + self.assertEqual((x < A.inv()) * A, (x < A.inv()) < A) + + def test6_6_2_9(self): + """ + In a 4D space with orthonormal basis {e1, e2, e3, e4}, project the 2-blade X = (e1 + e2) ^ (e3 + e4) onto + the 2-blade A = (e1 ^ e3). Then determine the rejection as the difference of X and its projection. Show + that is not a blade. + """ + GA, e_1, e_2, e_3, e_4 = Ga.build('e*1|2|3|4', g='1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + X = (e_1 + e_2) ^ (e_3 + e_4) + A = (e_1 ^ e_3) + + # projection of X onto B + def proj(X, B): + return (X < B.inv()) < B + + P = proj(X, A) + R = X - P + + # A 2-blade + a_1 = Symbol('a_1', real=True) + a_2 = Symbol('a_2', real=True) + a_3 = Symbol('a_3', real=True) + a_4 = Symbol('a_4', real=True) + a = GA.mv((a_1, a_2, a_3, a_4), 'vector') + + b_1 = Symbol('b_1', real=True) + b_2 = Symbol('b_2', real=True) + b_3 = Symbol('b_3', real=True) + b_4 = Symbol('b_4', real=True) + b = GA.mv((b_1, b_2, b_3, b_4), 'vector') + S = a ^ b + + # Try to solve the system and show there is no solution + system = [ + S_coef - R_coef for S_coef, R_coef in zip(S.blade_coefs(), R.blade_coefs()) + ] + + unknowns = [ + a_1, a_2, a_3, a_4, b_1, b_2, b_3, b_4 + ] + + # TODO: use solve if sympy fix it + result = solve_poly_system(system, unknowns) + self.assertTrue(result is None) diff --git a/test/test_chapter7.py b/test/test_chapter7.py new file mode 100644 index 00000000..937ee7d2 --- /dev/null +++ b/test/test_chapter7.py @@ -0,0 +1,83 @@ +from .test_utils import TestCase + +from sympy import Symbol, pi, cos, sin, solve, sqrt, Rational, Mod +from galgebra.ga import Ga + + +class TestChapter7(TestCase): + + def test7_9_1(self): + """ + Drills. + """ + Ga.dual_mode("Iinv+") + + GA, e_1, e_2, e_3 = Ga.build("e*1|2|3", g='1 0 0, 0 1 0, 0 0 1') + + # .1 Compute R1 + R1 = ((e_1 ^ e_2) * (-pi / 4)).exp() + self.assertEqual(R1 * e_1 * R1.rev(), e_2) + + # .2 Compute R2 + R2 = ((e_3 ^ e_1) * (pi / 4)).exp() + + # .3 Compute R2 R1 + R2R1 = R2 * R1 + self.assertEqual(R2R1 * (e_1 ^ e_2) * R2R1.rev(), -e_2 ^ e_3) + + # .4 Compute the axis and angle of R2 R1 + a_1 = Symbol('a_1', real=True) + a_2 = Symbol('a_2', real=True) + a_3 = Symbol('a_3', real=True) + a = GA.mv((a_1, a_2, a_3), 'vector') + + theta = Symbol('theta', real=True) + A = a.dual() + R = cos(theta / 2) + A * sin(theta / 2) + + system = (R2R1 - R).blade_coefs() + unknowns = [a_1, a_2, a_3, theta] + + results = solve(system, unknowns, dict=True) + self.assertEqual(a.subs(results[0]), (-e_1 - e_2 + e_3) / sqrt(3)) + self.assertEqual(Mod(theta.subs(results[0]), 2 * pi), Mod(Rational(2, 3) * pi, 2 * pi)) + self.assertEqual(a.subs(results[1]), -(-e_1 - e_2 + e_3) / sqrt(3)) + self.assertEqual(Mod(theta.subs(results[1]), 2 * pi), Mod(-Rational(2, 3) * pi, 2 * pi)) + + # Check .4 against .3 + a = a.subs(results[0]) + theta = theta.subs(results[0]) + A = a.dual() + R = cos(theta / 2) + A * sin(theta / 2) + self.assertEqual(R * (e_1 ^ e_2) * R.rev(), -e_2 ^ e_3) + + # .5 Compute the product of rotors + GA, e_1, e_2, e_3, e_4 = Ga.build("e*1|2|3|4", g='1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1') + + R1 = ((e_1 ^ e_4) * (-pi / 4)).exp() + R2 = ((e_2 ^ e_3) * (-pi / 4)).exp() + R = R2 * R1 + + self.assertEqual(R * (e_1 ^ e_2) * R.rev(), -e_3 ^ e_4) + + # .6 Reflect in a plane + A = (e_1 + e_2) ^ e_3 + B = e_1 ^ e_4 + b = B.dual() + + def hat(k, M): + return ((-1) ** k) * M + + self.assertEqual(b * hat(2, A) * b.inv(), (-e_1 + e_2) ^ e_3) + + # .7 Reflect the dual plane reflector e1 in the plane e1 ^ e3 + GA, e_1, e_2, e_3 = Ga.build("e*1|2|3", g='1 0 0, 0 1 0, 0 0 1') + + a = e_1 + B = e_1 ^ e_3 + b = B.dual() + + self.assertEqual(-b * a * b.inv(), e_1) + + # Reset + Ga.dual_mode() diff --git a/test/test_chapter8.py b/test/test_chapter8.py new file mode 100644 index 00000000..6c078c25 --- /dev/null +++ b/test/test_chapter8.py @@ -0,0 +1,39 @@ +from .test_utils import TestCase, com + +from sympy import S +from galgebra.ga import Ga + + +class TestChapter8(TestCase): + + def test8_2_1(self): + """ + The commutator product. + """ + GA = Ga("e*1|2|3", g=[1, 1, 1]) + + A = GA.mv('A', 'mv') + B = GA.mv('B', 'mv') + C = GA.mv('C', 'mv') + + self.assertEqual(com(com(A, B), C) - com(A, com(B, C)), com(B, com(C, A))) + self.assertEqual(com(com(A, B), C) + com(com(C, A), B) + com(com(B, C), A), S.Zero) + + B = GA.mv('B', 'blade', 2) + X = GA.mv('A', 'mv') + self.assertEqual(B.pure_grade(), 2) + self.assertEqual(com(X, B).pure_grade(), -2) # Not pure + + E = com(X, B).rev() + self.assertEqual(E, (B.rev() * X.rev() - X.rev() * B.rev()) / 2) + self.assertEqual(E, (X.rev() * B - B * X.rev()) / 2) + self.assertEqual(E, com(X.rev(), B)) + + alpha = GA.mv('alpha', 'scalar') + self.assertEqual(alpha * X, alpha ^ X) + + a = GA.mv('a', 'vector') + self.assertEqual(a * X, (a < X) + (a ^ X)) + + A = GA.mv('A', 'grade', 2) + self.assertEqual(A * X, (A < X) + com(A, X) + (A ^ X)) diff --git a/test/test_generator.py b/test/test_generator.py new file mode 100644 index 00000000..23cb7408 --- /dev/null +++ b/test/test_generator.py @@ -0,0 +1,93 @@ +from .test_utils import TestCase + +import importlib + +from galgebra.ga import Ga +from galgebra.mv import J, Jinv +from galgebra.generator import format_geometric_algebra, flatten, expand + + +class TestGenerator(TestCase): + + def setUp(self): + Ga.dual_mode("Iinv+") + GA = Ga('e*1|2|3|4', g=[0, 1, 1, 1]) + GA.build_cobases() + with open('flat_ga.py', 'wt') as flat_GA_file: + flat_GA_file.write(format_geometric_algebra(GA)) + self.GA = GA + + flat_GA = importlib.import_module('flat_ga') + self.flat_GA = flat_GA + + def test_add(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(X + Y, expand(GA, flatten(flat_GA, X) + flatten(flat_GA, Y))) + + def test_sub(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(X - Y, expand(GA, flatten(flat_GA, X) - flatten(flat_GA, Y))) + + def test_mul(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(X * Y, expand(GA, flatten(flat_GA, X) * flatten(flat_GA, Y))) + + def test_and(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(Jinv(J(X) ^ J(Y)), expand(GA, flatten(flat_GA, X) & flatten(flat_GA, Y))) + + def test_xor(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(X ^ Y, expand(GA, flatten(flat_GA, X) ^ flatten(flat_GA, Y))) + + def test_lshift(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(X << Y, expand(GA, flatten(flat_GA, X) << flatten(flat_GA, Y))) + + def test_rshift(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(X >> Y, expand(GA, flatten(flat_GA, X) >> flatten(flat_GA, Y))) + + def test_meet_and_join(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + Y = GA.mv('y', 'mv') + + self.assertEqual(Jinv(J(X) ^ J(Y)), expand(GA, flatten(flat_GA, X).meet(flatten(flat_GA, Y)))) + self.assertEqual(X ^ Y, expand(GA, flatten(flat_GA, X).join(flatten(flat_GA, Y)))) + + def test_rev(self): + GA = self.GA + flat_GA = self.flat_GA + X = GA.mv('x', 'mv') + + self.assertEqual(X.rev(), expand(GA, flatten(flat_GA, X).rev())) diff --git a/test/test_pga.py b/test/test_pga.py new file mode 100644 index 00000000..277979de --- /dev/null +++ b/test/test_pga.py @@ -0,0 +1,796 @@ +from .test_utils import TestCase, com + +from sympy import acos, asin, cos, sin, pi, Rational, S, sqrt, Symbol, symbols +from galgebra.ga import Ga +from galgebra.mv import J, Jinv + + +class TestPGA(TestCase): + + def setUp(self): + """ + Setup 3D Projective Geometric Algebra aka PGA. + """ + indexes = ( + ((),), + ((0,), (1,), (2,), (3,)), + ((0, 1), (0, 2), (0, 3), (1, 2), (3, 1), (2, 3)), + ((0, 2, 1), (0, 1, 3), (0, 3, 2), (1, 2, 3)), + ((0, 1, 2, 3),) + ) + + # Internal products use lexical ordering for computing expressions + signs = ( + (1,), + (1, 1, 1, 1), + (1, 1, 1, 1, -1, 1), + (-1, 1, -1, 1), + (1,) + ) + + # TODO: build this from indexes... + coindexes = ( + ((0, 1, 2, 3),), + ((1, 2, 3), (0, 3, 2), (0, 1, 3), (0, 2, 1)), + ((2, 3), (3, 1), (1, 2), (0, 3), (0, 2), (0, 1)), + ((3,), (2,), (1,), (0,)), + ((),), + ) + + PGA, e_0, e_1, e_2, e_3 = Ga.build('e*0|1|2|3', g=[0, 1, 1, 1], sign_and_indexes=(signs, indexes)) + + # TODO: move this somewhere useful... + PGA.build_cobases(coindexes=coindexes) + + self.PGA = PGA + self.e_0 = e_0 + self.e_1 = e_1 + self.e_2 = e_2 + self.e_3 = e_3 + self.e_01 = e_0 ^ e_1 + self.e_02 = e_0 ^ e_2 + self.e_03 = e_0 ^ e_3 + self.e_12 = e_1 ^ e_2 + self.e_31 = e_3 ^ e_1 + self.e_23 = e_2 ^ e_3 + self.e_032 = e_0 ^ e_3 ^ e_2 + self.e_013 = e_0 ^ e_1 ^ e_3 + self.e_021 = e_0 ^ e_2 ^ e_1 + self.e_123 = e_1 ^ e_2 ^ e_3 + self.e_0123 = e_0 ^ e_1 ^ e_2 ^ e_3 + + self.e_13 = e_1 ^ e_3 + + def point(self, x, y, z): + """ + Make a point. + """ + return x * self.e_032 + y * self.e_013 + z * self.e_021 + self.e_123 + + def direction(self, x, y, z): + """ + Make an ideal point (direction). + """ + return x * self.e_032 + y * self.e_013 + z * self.e_021 + + def plane(self, a, b, c, d): + """ + Make a plane. + """ + return a * self.e_1 + b * self.e_2 + c * self.e_3 + d * self.e_0 + + def norm(self, X): + assert len(self.PGA.grade_decomposition(X)) == 1 + squared_norm = X | X + if not squared_norm.is_scalar(): + raise ValueError("X | X isn't a scalar") + squared_norm = squared_norm.scalar() + if squared_norm == S.Zero: + raise ValueError("X k-vector is null") + return sqrt(abs(squared_norm)) + + def norm2(self, X): + assert len(self.PGA.grade_decomposition(X)) == 1 + squared_norm = X | X + if not squared_norm.is_scalar(): + raise ValueError("X | X isn't a scalar") + squared_norm = squared_norm.scalar() + return squared_norm + + def ideal_norm(self, X): + assert len(self.PGA.grade_decomposition(X)) == 1 + squared_norm = 0 + for c in X.blade_coefs(): + squared_norm += c * c + if squared_norm == S.Zero: + raise ValueError("X k-vector is null") + return sqrt(squared_norm) + + def test_multiplication_table(self): + """ + Reproduce multiplication table. + """ + # 1 + self.assertEqual(S.One * self.e_0, self.e_0) + self.assertEqual(S.One * self.e_1, self.e_1) + self.assertEqual(S.One * self.e_2, self.e_2) + self.assertEqual(S.One * self.e_3, self.e_3) + self.assertEqual(S.One * self.e_01, self.e_01) + self.assertEqual(S.One * self.e_02, self.e_02) + self.assertEqual(S.One * self.e_03, self.e_03) + self.assertEqual(S.One * self.e_12, self.e_12) + self.assertEqual(S.One * self.e_31, self.e_31) + self.assertEqual(S.One * self.e_23, self.e_23) + self.assertEqual(S.One * self.e_021, self.e_021) + self.assertEqual(S.One * self.e_013, self.e_013) + self.assertEqual(S.One * self.e_032, self.e_032) + self.assertEqual(S.One * self.e_123, self.e_123) + self.assertEqual(S.One * self.e_0123, self.e_0123) + + # e0 + self.assertEqual(self.e_0 * self.e_0, S.Zero) + self.assertEqual(self.e_0 * self.e_1, self.e_01) + self.assertEqual(self.e_0 * self.e_2, self.e_02) + self.assertEqual(self.e_0 * self.e_3, self.e_03) + self.assertEqual(self.e_0 * self.e_01, S.Zero) + self.assertEqual(self.e_0 * self.e_02, S.Zero) + self.assertEqual(self.e_0 * self.e_03, S.Zero) + self.assertEqual(self.e_0 * self.e_12, -self.e_021) + self.assertEqual(self.e_0 * self.e_31, -self.e_013) + self.assertEqual(self.e_0 * self.e_23, -self.e_032) + self.assertEqual(self.e_0 * self.e_021, S.Zero) + self.assertEqual(self.e_0 * self.e_013, S.Zero) + self.assertEqual(self.e_0 * self.e_032, S.Zero) + self.assertEqual(self.e_0 * self.e_123, self.e_0123) + self.assertEqual(self.e_0 * self.e_0123, S.Zero) + + # e1 + self.assertEqual(self.e_1 * self.e_0, -self.e_01) + self.assertEqual(self.e_1 * self.e_1, S.One) + self.assertEqual(self.e_1 * self.e_2, self.e_12) + self.assertEqual(self.e_1 * self.e_3, -self.e_31) + self.assertEqual(self.e_1 * self.e_01, -self.e_0) + self.assertEqual(self.e_1 * self.e_02, self.e_021) + self.assertEqual(self.e_1 * self.e_03, -self.e_013) + self.assertEqual(self.e_1 * self.e_12, self.e_2) + self.assertEqual(self.e_1 * self.e_31, -self.e_3) + self.assertEqual(self.e_1 * self.e_23, self.e_123) + self.assertEqual(self.e_1 * self.e_021, self.e_02) + self.assertEqual(self.e_1 * self.e_013, -self.e_03) + self.assertEqual(self.e_1 * self.e_032, self.e_0123) + self.assertEqual(self.e_1 * self.e_123, self.e_23) + self.assertEqual(self.e_1 * self.e_0123, self.e_032) + + # e2 + self.assertEqual(self.e_2 * self.e_0, -self.e_02) + self.assertEqual(self.e_2 * self.e_1, -self.e_12) + self.assertEqual(self.e_2 * self.e_2, S.One) + self.assertEqual(self.e_2 * self.e_3, self.e_23) + self.assertEqual(self.e_2 * self.e_01, -self.e_021) + self.assertEqual(self.e_2 * self.e_02, -self.e_0) + self.assertEqual(self.e_2 * self.e_03, self.e_032) + self.assertEqual(self.e_2 * self.e_12, -self.e_1) + self.assertEqual(self.e_2 * self.e_31, self.e_123) + self.assertEqual(self.e_2 * self.e_23, self.e_3) + self.assertEqual(self.e_2 * self.e_021, -self.e_01) + self.assertEqual(self.e_2 * self.e_013, self.e_0123) + self.assertEqual(self.e_2 * self.e_032, self.e_03) + self.assertEqual(self.e_2 * self.e_123, self.e_31) + self.assertEqual(self.e_2 * self.e_0123, self.e_013) + + # e3 + self.assertEqual(self.e_3 * self.e_0, -self.e_03) + self.assertEqual(self.e_3 * self.e_1, self.e_31) + self.assertEqual(self.e_3 * self.e_2, -self.e_23) + self.assertEqual(self.e_3 * self.e_3, S.One) + self.assertEqual(self.e_3 * self.e_01, self.e_013) + self.assertEqual(self.e_3 * self.e_02, -self.e_032) + self.assertEqual(self.e_3 * self.e_03, -self.e_0) + self.assertEqual(self.e_3 * self.e_12, self.e_123) + self.assertEqual(self.e_3 * self.e_31, self.e_1) + self.assertEqual(self.e_3 * self.e_23, -self.e_2) + self.assertEqual(self.e_3 * self.e_021, self.e_0123) + self.assertEqual(self.e_3 * self.e_013, self.e_01) + self.assertEqual(self.e_3 * self.e_032, -self.e_02) + self.assertEqual(self.e_3 * self.e_123, self.e_12) + self.assertEqual(self.e_3 * self.e_0123, self.e_021) + + # e01 + self.assertEqual(self.e_01 * self.e_0, S.Zero) + self.assertEqual(self.e_01 * self.e_1, self.e_0) + self.assertEqual(self.e_01 * self.e_2, -self.e_021) + self.assertEqual(self.e_01 * self.e_3, self.e_013) + self.assertEqual(self.e_01 * self.e_01, S.Zero) + self.assertEqual(self.e_01 * self.e_02, S.Zero) + self.assertEqual(self.e_01 * self.e_03, S.Zero) + self.assertEqual(self.e_01 * self.e_12, self.e_02) + self.assertEqual(self.e_01 * self.e_31, -self.e_03) + self.assertEqual(self.e_01 * self.e_23, self.e_0123) + self.assertEqual(self.e_01 * self.e_021, S.Zero) + self.assertEqual(self.e_01 * self.e_013, S.Zero) + self.assertEqual(self.e_01 * self.e_032, S.Zero) + self.assertEqual(self.e_01 * self.e_123, -self.e_032) + self.assertEqual(self.e_01 * self.e_0123, S.Zero) + + # e02 + self.assertEqual(self.e_02 * self.e_0, S.Zero) + self.assertEqual(self.e_02 * self.e_1, self.e_021) + self.assertEqual(self.e_02 * self.e_2, self.e_0) + self.assertEqual(self.e_02 * self.e_3, -self.e_032) + self.assertEqual(self.e_02 * self.e_01, S.Zero) + self.assertEqual(self.e_02 * self.e_02, S.Zero) + self.assertEqual(self.e_02 * self.e_03, S.Zero) + self.assertEqual(self.e_02 * self.e_12, -self.e_01) + self.assertEqual(self.e_02 * self.e_31, self.e_0123) + self.assertEqual(self.e_02 * self.e_23, self.e_03) + self.assertEqual(self.e_02 * self.e_021, S.Zero) + self.assertEqual(self.e_02 * self.e_013, S.Zero) + self.assertEqual(self.e_02 * self.e_032, S.Zero) + self.assertEqual(self.e_02 * self.e_123, -self.e_013) + self.assertEqual(self.e_02 * self.e_0123, S.Zero) + + # e03 + self.assertEqual(self.e_03 * self.e_0, S.Zero) + self.assertEqual(self.e_03 * self.e_1, -self.e_013) + self.assertEqual(self.e_03 * self.e_2, self.e_032) + self.assertEqual(self.e_03 * self.e_3, self.e_0) + self.assertEqual(self.e_03 * self.e_01, S.Zero) + self.assertEqual(self.e_03 * self.e_02, S.Zero) + self.assertEqual(self.e_03 * self.e_03, S.Zero) + self.assertEqual(self.e_03 * self.e_12, self.e_0123) + self.assertEqual(self.e_03 * self.e_31, self.e_01) + self.assertEqual(self.e_03 * self.e_23, -self.e_02) + self.assertEqual(self.e_03 * self.e_021, S.Zero) + self.assertEqual(self.e_03 * self.e_013, S.Zero) + self.assertEqual(self.e_03 * self.e_032, S.Zero) + self.assertEqual(self.e_03 * self.e_123, -self.e_021) + self.assertEqual(self.e_03 * self.e_0123, S.Zero) + + # e12 + self.assertEqual(self.e_12 * self.e_0, -self.e_021) + self.assertEqual(self.e_12 * self.e_1, -self.e_2) + self.assertEqual(self.e_12 * self.e_2, self.e_1) + self.assertEqual(self.e_12 * self.e_3, self.e_123) + self.assertEqual(self.e_12 * self.e_01, -self.e_02) + self.assertEqual(self.e_12 * self.e_02, self.e_01) + self.assertEqual(self.e_12 * self.e_03, self.e_0123) + self.assertEqual(self.e_12 * self.e_12, -S.One) + self.assertEqual(self.e_12 * self.e_31, self.e_23) + self.assertEqual(self.e_12 * self.e_23, -self.e_31) + self.assertEqual(self.e_12 * self.e_021, self.e_0) + self.assertEqual(self.e_12 * self.e_013, self.e_032) + self.assertEqual(self.e_12 * self.e_032, -self.e_013) + self.assertEqual(self.e_12 * self.e_123, -self.e_3) + self.assertEqual(self.e_12 * self.e_0123, -self.e_03) + + # e31 + self.assertEqual(self.e_31 * self.e_0, -self.e_013) + self.assertEqual(self.e_31 * self.e_1, self.e_3) + self.assertEqual(self.e_31 * self.e_2, self.e_123) + self.assertEqual(self.e_31 * self.e_3, -self.e_1) + self.assertEqual(self.e_31 * self.e_01, self.e_03) + self.assertEqual(self.e_31 * self.e_02, self.e_0123) + self.assertEqual(self.e_31 * self.e_03, -self.e_01) + self.assertEqual(self.e_31 * self.e_12, -self.e_23) + self.assertEqual(self.e_31 * self.e_31, -S.One) + self.assertEqual(self.e_31 * self.e_23, self.e_12) + self.assertEqual(self.e_31 * self.e_021, -self.e_032) + self.assertEqual(self.e_31 * self.e_013, self.e_0) + self.assertEqual(self.e_31 * self.e_032, self.e_021) + self.assertEqual(self.e_31 * self.e_123, -self.e_2) + self.assertEqual(self.e_31 * self.e_0123, -self.e_02) + + # e23 + self.assertEqual(self.e_23 * self.e_0, -self.e_032) + self.assertEqual(self.e_23 * self.e_1, self.e_123) + self.assertEqual(self.e_23 * self.e_2, -self.e_3) + self.assertEqual(self.e_23 * self.e_3, self.e_2) + self.assertEqual(self.e_23 * self.e_01, self.e_0123) + self.assertEqual(self.e_23 * self.e_02, -self.e_03) + self.assertEqual(self.e_23 * self.e_03, self.e_02) + self.assertEqual(self.e_23 * self.e_12, self.e_31) + self.assertEqual(self.e_23 * self.e_31, -self.e_12) + self.assertEqual(self.e_23 * self.e_23, -S.One) + self.assertEqual(self.e_23 * self.e_021, self.e_013) + self.assertEqual(self.e_23 * self.e_013, -self.e_021) + self.assertEqual(self.e_23 * self.e_032, self.e_0) + self.assertEqual(self.e_23 * self.e_123, -self.e_1) + self.assertEqual(self.e_23 * self.e_0123, -self.e_01) + + # e021 + self.assertEqual(self.e_021 * self.e_0, S.Zero) + self.assertEqual(self.e_021 * self.e_1, self.e_02) + self.assertEqual(self.e_021 * self.e_2, -self.e_01) + self.assertEqual(self.e_021 * self.e_3, -self.e_0123) + self.assertEqual(self.e_021 * self.e_01, S.Zero) + self.assertEqual(self.e_021 * self.e_02, S.Zero) + self.assertEqual(self.e_021 * self.e_03, S.Zero) + self.assertEqual(self.e_021 * self.e_12, self.e_0) + self.assertEqual(self.e_021 * self.e_31, self.e_032) + self.assertEqual(self.e_021 * self.e_23, -self.e_013) + self.assertEqual(self.e_021 * self.e_021, S.Zero) + self.assertEqual(self.e_021 * self.e_013, S.Zero) + self.assertEqual(self.e_021 * self.e_032, S.Zero) + self.assertEqual(self.e_021 * self.e_123, self.e_03) + self.assertEqual(self.e_021 * self.e_0123, S.Zero) + + # e013 + self.assertEqual(self.e_013 * self.e_0, S.Zero) + self.assertEqual(self.e_013 * self.e_1, -self.e_03) + self.assertEqual(self.e_013 * self.e_2, -self.e_0123) + self.assertEqual(self.e_013 * self.e_3, self.e_01) + self.assertEqual(self.e_013 * self.e_01, S.Zero) + self.assertEqual(self.e_013 * self.e_02, S.Zero) + self.assertEqual(self.e_013 * self.e_03, S.Zero) + self.assertEqual(self.e_013 * self.e_12, -self.e_032) + self.assertEqual(self.e_013 * self.e_31, self.e_0) + self.assertEqual(self.e_013 * self.e_23, self.e_021) + self.assertEqual(self.e_013 * self.e_021, S.Zero) + self.assertEqual(self.e_013 * self.e_013, S.Zero) + self.assertEqual(self.e_013 * self.e_032, S.Zero) + self.assertEqual(self.e_013 * self.e_123, self.e_02) + self.assertEqual(self.e_013 * self.e_0123, S.Zero) + + # e032 + self.assertEqual(self.e_032 * self.e_0, S.Zero) + self.assertEqual(self.e_032 * self.e_1, -self.e_0123) + self.assertEqual(self.e_032 * self.e_2, self.e_03) + self.assertEqual(self.e_032 * self.e_3, -self.e_02) + self.assertEqual(self.e_032 * self.e_01, S.Zero) + self.assertEqual(self.e_032 * self.e_02, S.Zero) + self.assertEqual(self.e_032 * self.e_03, S.Zero) + self.assertEqual(self.e_032 * self.e_12, self.e_013) + self.assertEqual(self.e_032 * self.e_31, -self.e_021) + self.assertEqual(self.e_032 * self.e_23, self.e_0) + self.assertEqual(self.e_032 * self.e_021, S.Zero) + self.assertEqual(self.e_032 * self.e_013, S.Zero) + self.assertEqual(self.e_032 * self.e_032, S.Zero) + self.assertEqual(self.e_032 * self.e_123, self.e_01) + self.assertEqual(self.e_032 * self.e_0123, S.Zero) + + # e123 + self.assertEqual(self.e_123 * self.e_0, -self.e_0123) + self.assertEqual(self.e_123 * self.e_1, self.e_23) + self.assertEqual(self.e_123 * self.e_2, self.e_31) + self.assertEqual(self.e_123 * self.e_3, self.e_12) + self.assertEqual(self.e_123 * self.e_01, self.e_032) + self.assertEqual(self.e_123 * self.e_02, self.e_013) + self.assertEqual(self.e_123 * self.e_03, self.e_021) + self.assertEqual(self.e_123 * self.e_12, -self.e_3) + self.assertEqual(self.e_123 * self.e_31, -self.e_2) + self.assertEqual(self.e_123 * self.e_23, -self.e_1) + self.assertEqual(self.e_123 * self.e_021, -self.e_03) + self.assertEqual(self.e_123 * self.e_013, -self.e_02) + self.assertEqual(self.e_123 * self.e_032, -self.e_01) + self.assertEqual(self.e_123 * self.e_123, -S.One) + self.assertEqual(self.e_123 * self.e_0123, self.e_0) + + # e0123 + self.assertEqual(self.e_0123 * self.e_0, S.Zero) + self.assertEqual(self.e_0123 * self.e_1, -self.e_032) + self.assertEqual(self.e_0123 * self.e_2, -self.e_013) + self.assertEqual(self.e_0123 * self.e_3, -self.e_021) + self.assertEqual(self.e_0123 * self.e_01, S.Zero) + self.assertEqual(self.e_0123 * self.e_02, S.Zero) + self.assertEqual(self.e_0123 * self.e_03, S.Zero) + self.assertEqual(self.e_0123 * self.e_12, -self.e_03) + self.assertEqual(self.e_0123 * self.e_31, -self.e_02) + self.assertEqual(self.e_0123 * self.e_23, -self.e_01) + self.assertEqual(self.e_0123 * self.e_021, S.Zero) + self.assertEqual(self.e_0123 * self.e_013, S.Zero) + self.assertEqual(self.e_0123 * self.e_032, S.Zero) + self.assertEqual(self.e_0123 * self.e_123, -self.e_0) + self.assertEqual(self.e_0123 * self.e_0123, S.Zero) + + def test_J(self): + """ + Check we can join and meet using J and Jinv. + """ + PGA = self.PGA + for k in range(PGA.n + 1): + X = PGA.mv('x', 'grade', k) + self.assertEqual(X, Jinv(J(X))) + X = PGA.mv('x', 'mv') + self.assertEqual(X, Jinv(J(X))) + + def test_geometry_incidence_planes_meet_into_points(self): + """ + Planes meet into points. + """ + x = Symbol('x') + y = Symbol('y') + z = Symbol('z') + p1 = self.plane(1, 0, 0, -x) + p2 = self.plane(0, 1, 0, -y) + p3 = self.plane(0, 0, 1, -z) + R = self.point(x, y, z) + self.assertProjEqual(R, p1 ^ p2 ^ p3) + self.assertProjEqual(R, p1 ^ p3 ^ p2) + self.assertProjEqual(R, p2 ^ p1 ^ p3) + self.assertProjEqual(R, p2 ^ p3 ^ p1) + self.assertProjEqual(R, p3 ^ p1 ^ p2) + self.assertProjEqual(R, p3 ^ p2 ^ p1) + + def test_geometry_incidence_planes_meet_into_points_2(self): + """ + Planes meet into points. + """ + P1 = self.point(*symbols('x1 y1 z1', real=True)) + P2 = self.point(*symbols('x2 y2 z2', real=True)) + P3 = self.point(*symbols('x3 y3 z3', real=True)) + P4 = self.point(*symbols('x4 y4 z4', real=True)) + p123 = Jinv(J(P1) ^ J(P2) ^ J(P3)) + p124 = Jinv(J(P1) ^ J(P2) ^ J(P4)) + p234 = Jinv(J(P2) ^ J(P3) ^ J(P4)) + p314 = Jinv(J(P3) ^ J(P1) ^ J(P4)) + + # TODO: find a way to meet planes faster... + # self.assertProjEquals(p123 ^ p124 ^ p314, P1) + # self.assertProjEquals(p123 ^ p124 ^ p234, P2) + # self.assertProjEquals(p123 ^ p234 ^ p314, P3) + # self.assertProjEquals(p124 ^ p234 ^ p314, P4) + + def test_geometry_incidence_points_join_into_planes(self): + """ + Points join into planes. + """ + x1, y1, z1 = symbols('x1 y1 z1', real=True) + P1 = self.point(x1, y1, z1) + + x2, y2, z2 = symbols('x2 y2 z2', real=True) + P2 = self.point(x2, y2, z2) + + x3, y3, z3 = symbols('x3 y3 z3', real=True) + P3 = self.point(x3, y3, z3) + + pp = Jinv(J(P1) ^ J(P2) ^ J(P3)) + pm = Jinv(J(P1) ^ J(P3) ^ J(P2)) + self.assertEqual(pp, -pm) + + a = Symbol('a') + b = Symbol('b') + c = Symbol('c') + + coefs = pp.blade_coefs([self.e_0, self.e_1, self.e_2, self.e_3]) + p = coefs[0] + a * coefs[1] + b * coefs[2] + c * coefs[3] + self.assertEqual(p.subs({a: x1, b: y1, c: z1}), S.Zero) + self.assertEqual(p.subs({a: x2, b: y2, c: z2}), S.Zero) + self.assertEqual(p.subs({a: x3, b: y3, c: z3}), S.Zero) + + coefs = pm.blade_coefs([self.e_0, self.e_1, self.e_2, self.e_3]) + p = coefs[0] + a * coefs[1] + b * coefs[2] + c * coefs[3] + self.assertEqual(p.subs({a: x1, b: y1, c: z1}), S.Zero) + self.assertEqual(p.subs({a: x2, b: y2, c: z2}), S.Zero) + self.assertEqual(p.subs({a: x3, b: y3, c: z3}), S.Zero) + + def test_geometry_incidence_join_and_plane_side(self): + """ + Join and plane side. + """ + P1 = self.point(1, 0, 0) + P2 = self.point(0, 1, 0) + P3 = self.point(0, 0, 1) + pp = Jinv(J(P1) ^ J(P2) ^ J(P3)) + + px = Symbol('px') + py = Symbol('py') + pz = Symbol('pz') + coefs = pp.blade_coefs([self.e_0, self.e_1, self.e_2, self.e_3]) + p = coefs[0] + px * coefs[1] + py * coefs[2] + pz * coefs[3] + + self.assertTrue(p.subs({px: Rational(1, 3) - 0.01, py: Rational(1, 3) - 0.01, pz: Rational(1, 3) - 0.01}) < 0.0) + self.assertTrue(p.subs({px: Rational(1, 3) + 0.01, py: Rational(1, 3) + 0.01, pz: Rational(1, 3) + 0.01}) > 0.0) + + def test_geometry_incidence_points_join_into_lines(self): + """ + Points join into lines. + """ + x1, y1, z1 = symbols('x1 y1 z1') + P1 = self.point(x1, y1, z1) + + x2, y2, z2 = symbols('x2 y2 z2') + P2 = self.point(x2, y2, z2) + + lp = Jinv(J(P1) ^ J(P2)) + lm = Jinv(J(P2) ^ J(P1)) + self.assertEqual(lp, -lm) + + # TODO: add more tests... + + def test_metric_distance_of_points(self): + """ + We can measure distance between normalized points using the joining line norm. + """ + x0, y0, z0 = symbols('x0 y0 z0') + P0 = self.point(x0, y0, z0) + + x1, y1, z1 = symbols('x1 y1 z1') + P1 = self.point(x1, y1, z1) + + d = sqrt(abs((x1 - x0) ** 2 + (y1 - y0) ** 2 + (z1 - z0) ** 2)) + + lp = Jinv(J(P0) ^ J(P1)) + self.assertEqual(self.norm(lp), d) + + lm = Jinv(J(P1) ^ J(P0)) + self.assertEqual(self.norm(lm), d) + + def test_metric_angle_of_intersecting_planes(self): + """ + We can measure the angle between normalized planes using the meeting line norm. + """ + x = Symbol('x') + y = Symbol('y') + p1 = self.plane(1, 0, 0, -x) + p2 = self.plane(0, 1, 0, -y) + + self.assertEqual(asin(self.norm(p1 ^ p2)), pi / 2) + + def test_metric_distance_between_parallel_planes(self): + """ + We can measure the distance between two parallel and normalized planes using the meeting line norm. + """ + nx, ny, nz, x, y = symbols('nx ny nz x y', real=True) + + n_norm = sqrt(nx * nx + ny * ny + nz * nz) + p1 = self.plane(nx / n_norm, ny / n_norm, nz / n_norm, -x) + p2 = self.plane(nx / n_norm, ny / n_norm, nz / n_norm, -y) + + self.assertEqual(self.ideal_norm(p1 ^ p2), sqrt((x - y) ** 2)) + + def test_metric_oriented_distance_between_point_and_plane(self): + """ + We can measure the distance between a normalized point and a normalized plane. + """ + x0, y0, z0 = symbols('x0 y0 z0') + P0 = self.point(x0, y0, z0) + + d0 = Symbol('d0') + p0 = self.plane(1, 0, 0, -d0) + + self.assertEqual(Jinv(J(P0) ^ J(p0)), d0 - x0) + self.assertEqual(Jinv(J(p0) ^ J(P0)), x0 - d0) + + # TODO: We could assume d0 - x0 is not null for simplifying further + self.assertEqual(self.ideal_norm(P0 ^ p0), sqrt((d0 - x0) ** 2)) + + def test_metric_oriented_distance_between_point_and_line(self): + """ + We can measure the distance between a normalized point and a normalized line. + """ + x0, y0, z0 = symbols('x0 y0 z0') + P0 = self.point(x0, y0, z0) + + x1, y1, z1 = symbols('x1 y1 z1') + P1 = self.point(x1, y1, z1) + + x2, y2, z2 = symbols('x2 y2 z2') + P2 = self.point(x2, y2, z2) + + l0 = Jinv(J(P1) ^ J(P2)) + + d = self.norm(Jinv(J(P0) ^ J(l0))) + + self.assertEqual(d.subs({x0: 2, y0: 0, z0: 0, x1: 0, y1: 2, z1: 0, x2: 1, y2: 2, z2: 0}), 2) + self.assertEqual(d.subs({x0: 3, y0: 0, z0: 0, x1: 0, y1: 2, z1: 0, x2: 1, y2: 2, z2: 0}), 2) + self.assertEqual(d.subs({x0: 2, y0: 0, z0: 0, x1: 0, y1: 1, z1: 0, x2: 1, y2: 1, z2: 0}), 1) + self.assertEqual(d.subs({x0: 2, y0: 0, z0: 0, x1: 0, y1: 2, z1: 0, x2: 2, y2: 0, z2: 0}), 0) + + def test_metric_common_normal_line(self): + """ + We can find the common normal line of two normalized lines. + """ + x0, y0, z0 = symbols('x0 y0 z0', real=True) + P0 = self.point(x0, y0, z0) + + x1, y1, z1 = symbols('x1 y1 z1', real=True) + P1 = self.point(x1, y1, z1) + + x2, y2, z2 = symbols('x2 y2 z2', real=True) + P2 = self.point(x2, y2, z2) + + x3, y3, z3 = symbols('x3 y3 z3', real=True) + P3 = self.point(x3, y3, z3) + + l0 = Jinv(J(P0) ^ J(P1)) + l1 = Jinv(J(P2) ^ J(P3)) + ln = com(l0, l1) + + ln /= self.norm(ln) + l0 /= self.norm(l0) + l1 /= self.norm(l1) + + dot0 = l0 | ln + dot1 = l1 | ln + self.assertTrue(dot0.is_scalar()) + self.assertTrue(dot1.is_scalar()) + + self.assertEqual(acos(dot0.scalar()), S.Half * pi) + self.assertEqual(acos(dot1.scalar()), S.Half * pi) + + def test_metric_angle_between_lines(self): + """ + We can measure the angle between to normalized lines. + """ + x1, y1 = symbols('x1 y1', real=True) + P0 = self.point(0, 0, 0) + P1 = self.point(x1, y1, 0) + P2 = self.point(-y1, x1, 0) + + l0 = Jinv(J(P0) ^ J(P1)) # TODO: this feels weird... but ganja does the same + l1 = Jinv(J(P2) ^ J(P0)) # + + l0 /= self.norm(l0) + l1 /= self.norm(l1) + + dot = l0 | l1 + self.assertTrue(dot.is_scalar()) + + self.assertEqual(acos(dot.scalar()), pi / 2) + + @staticmethod + def rotor_cs(alpha, l): + return cos(-alpha / 2) + sin(-alpha / 2) * l # TODO: this feels weird... + + @staticmethod + def rotor_exp(alpha, l): + return (-alpha / 2 * l).exp() # TODO: this feels weird... + + def test_motors_rotator(self): + """ + Rotate anything around a normalized line. + """ + Or = self.point(+0, +0, +0) + Xp = self.point(+1, +0, +0) + Yp = self.point(+0, +1, +0) + Zp = self.point(+0, +0, +1) + Xm = self.point(-1, +0, +0) + Ym = self.point(+0, -1, +0) + Zm = self.point(+0, +0, -1) + + lXp = Jinv(J(Or) ^ J(Xp)) + lXm = Jinv(J(Or) ^ J(Xm)) + lYp = Jinv(J(Or) ^ J(Yp)) + lYm = Jinv(J(Or) ^ J(Ym)) + lZp = Jinv(J(Or) ^ J(Zp)) + lZm = Jinv(J(Or) ^ J(Zm)) + + d = Symbol('d') + pXp = self.plane(+1, 0, 0, d) + pXm = self.plane(-1, 0, 0, d) + pYp = self.plane(0, +1, 0, d) + pYm = self.plane(0, -1, 0, d) + pZp = self.plane(0, 0, +1, d) + pZm = self.plane(0, 0, -1, d) + + # Around X+ + RXp = self.rotor_cs(pi / 2, lXp) + self.assertEqual(RXp, self.rotor_exp(pi / 2, lXp)) + # Points + self.assertProjEqual(RXp * Yp * RXp.rev(), Zp) + self.assertProjEqual(RXp * Zp * RXp.rev(), Ym) + self.assertProjEqual(RXp * Ym * RXp.rev(), Zm) + self.assertProjEqual(RXp * Zm * RXp.rev(), Yp) + # Lines + self.assertProjEqual(RXp * lYp * RXp.rev(), lZp) + self.assertProjEqual(RXp * lZp * RXp.rev(), lYm) + self.assertProjEqual(RXp * lYm * RXp.rev(), lZm) + self.assertProjEqual(RXp * lZm * RXp.rev(), lYp) + # Planes + self.assertProjEqual(RXp * pYp * RXp.rev(), pZp) + self.assertProjEqual(RXp * pZp * RXp.rev(), pYm) + self.assertProjEqual(RXp * pYm * RXp.rev(), pZm) + self.assertProjEqual(RXp * pZm * RXp.rev(), pYp) + + # Around X- + RXm = self.rotor_cs(pi / 2, lXm) + self.assertEqual(RXm, self.rotor_exp(pi / 2, lXm)) + # Points + self.assertProjEqual(RXm * Yp * RXm.rev(), Zm) + self.assertProjEqual(RXm * Zm * RXm.rev(), Ym) + self.assertProjEqual(RXm * Ym * RXm.rev(), Zp) + self.assertProjEqual(RXm * Zp * RXm.rev(), Yp) + # Lines + self.assertProjEqual(RXm * lYp * RXm.rev(), lZm) + self.assertProjEqual(RXm * lZm * RXm.rev(), lYm) + self.assertProjEqual(RXm * lYm * RXm.rev(), lZp) + self.assertProjEqual(RXm * lZp * RXm.rev(), lYp) + # Planes + self.assertProjEqual(RXm * pYp * RXm.rev(), pZm) + self.assertProjEqual(RXm * pZm * RXm.rev(), pYm) + self.assertProjEqual(RXm * pYm * RXm.rev(), pZp) + self.assertProjEqual(RXm * pZp * RXm.rev(), pYp) + + # Around Y+ + RYp = self.rotor_cs(pi / 2, lYp) + self.assertEqual(RYp, self.rotor_exp(pi / 2, lYp)) + # Points + self.assertProjEqual(RYp * Xp * RYp.rev(), Zm) + self.assertProjEqual(RYp * Zm * RYp.rev(), Xm) + self.assertProjEqual(RYp * Xm * RYp.rev(), Zp) + self.assertProjEqual(RYp * Zp * RYp.rev(), Xp) + # Lines + self.assertProjEqual(RYp * lXp * RYp.rev(), lZm) + self.assertProjEqual(RYp * lZm * RYp.rev(), lXm) + self.assertProjEqual(RYp * lXm * RYp.rev(), lZp) + self.assertProjEqual(RYp * lZp * RYp.rev(), lXp) + # Planes + self.assertProjEqual(RYp * pXp * RYp.rev(), pZm) + self.assertProjEqual(RYp * pZm * RYp.rev(), pXm) + self.assertProjEqual(RYp * pXm * RYp.rev(), pZp) + self.assertProjEqual(RYp * pZp * RYp.rev(), pXp) + + # Around Y- + RYm = self.rotor_cs(pi / 2, lYm) + self.assertEqual(RYm, self.rotor_exp(pi / 2, lYm)) + # Points + self.assertProjEqual(RYm * Xp * RYm.rev(), Zp) + self.assertProjEqual(RYm * Zp * RYm.rev(), Xm) + self.assertProjEqual(RYm * Xm * RYm.rev(), Zm) + self.assertProjEqual(RYm * Zm * RYm.rev(), Xp) + # Lines + self.assertProjEqual(RYm * lXp * RYm.rev(), lZp) + self.assertProjEqual(RYm * lZp * RYm.rev(), lXm) + self.assertProjEqual(RYm * lXm * RYm.rev(), lZm) + self.assertProjEqual(RYm * lZm * RYm.rev(), lXp) + # Planes + self.assertProjEqual(RYm * pXp * RYm.rev(), pZp) + self.assertProjEqual(RYm * pZp * RYm.rev(), pXm) + self.assertProjEqual(RYm * pXm * RYm.rev(), pZm) + self.assertProjEqual(RYm * pZm * RYm.rev(), pXp) + + # Around Z+ + RZp = self.rotor_cs(pi / 2, lZp) + self.assertEqual(RZp, self.rotor_exp(pi / 2, lZp)) + # Points + self.assertProjEqual(RZp * Xp * RZp.rev(), Yp) + self.assertProjEqual(RZp * Yp * RZp.rev(), Xm) + self.assertProjEqual(RZp * Xm * RZp.rev(), Ym) + self.assertProjEqual(RZp * Ym * RZp.rev(), Xp) + # Lines + self.assertProjEqual(RZp * lXp * RZp.rev(), lYp) + self.assertProjEqual(RZp * lYp * RZp.rev(), lXm) + self.assertProjEqual(RZp * lXm * RZp.rev(), lYm) + self.assertProjEqual(RZp * lYm * RZp.rev(), lXp) + # Planes + self.assertProjEqual(RZp * pXp * RZp.rev(), pYp) + self.assertProjEqual(RZp * pYp * RZp.rev(), pXm) + self.assertProjEqual(RZp * pXm * RZp.rev(), pYm) + self.assertProjEqual(RZp * pYm * RZp.rev(), pXp) + + # Around Z- + RZm = self.rotor_cs(pi / 2, lZm) + self.assertEqual(RZm, self.rotor_exp(pi / 2, lZm)) + # Points + self.assertProjEqual(RZm * Xp * RZm.rev(), Ym) + self.assertProjEqual(RZm * Ym * RZm.rev(), Xm) + self.assertProjEqual(RZm * Xm * RZm.rev(), Yp) + self.assertProjEqual(RZm * Yp * RZm.rev(), Xp) + # Lines + self.assertProjEqual(RZm * lXp * RZm.rev(), lYm) + self.assertProjEqual(RZm * lYm * RZm.rev(), lXm) + self.assertProjEqual(RZm * lXm * RZm.rev(), lYp) + self.assertProjEqual(RZm * lYp * RZm.rev(), lXp) + # Planes + self.assertProjEqual(RZm * pXp * RZm.rev(), pYm) + self.assertProjEqual(RZm * pYm * RZm.rev(), pXm) + self.assertProjEqual(RZm * pXm * RZm.rev(), pYp) + self.assertProjEqual(RZm * pYp * RZm.rev(), pXp) + + def translator(self, d, l): + return 1 + (d / 2) * (l * self.e_0123) + + def test_motors_translator(self): + """ + Translate anything by a normalized line. + """ + x0, y0, z0 = symbols('x0 y0 z0', real=True) + Px = self.point(x0, y0, z0) + + P0 = self.point(0, 0, 0) + P1 = self.point(1, 0, 0) + l = Jinv(J(P0) ^ J(P1)) + + d = Symbol('d') + T = self.translator(d, l) + self.assertProjEqual(T * Px * T.rev(), self.point(x0 + d, y0, z0)) # TODO : like ganja.js but weird... diff --git a/test/test_utils.py b/test/test_utils.py new file mode 100644 index 00000000..0ab99ecc --- /dev/null +++ b/test/test_utils.py @@ -0,0 +1,76 @@ +import unittest + +from sympy import expand, S, simplify +from galgebra.ga import Ga +from galgebra.mv import Mv +from galgebra.metric import Simp + + +def com(A, B): + """ + I like free functions... + """ + return Ga.com(A, B) + + +class TestCase(unittest.TestCase): + + def assertEqual(self, first, second): + """ + Compare two expressions are equals. + """ + if isinstance(first, Mv): + first = first.obj + + if isinstance(second, Mv): + second = second.obj + + # We need to help sympy a little... + first = Simp.apply(expand(first)) + second = Simp.apply(expand(second)) + + # Check + diff = simplify(first - second) + + self.assertTrue(diff == 0, "\n%s\n==\n%s\n%s" % (first, second, diff)) + + def assertProjEqual(self, first, second): + """ + Compare two points, two planes or two lines up to a scalar. + """ + assert isinstance(first, Mv) + assert isinstance(second, Mv) + + # TODO: this should use Mv methods and not the derived test case methods... + first /= self.norm(first) + second /= self.norm(second) + + # We need to help sympy a little... + X = Simp.apply(expand(first.obj)) + Y = Simp.apply(expand(second.obj)) + + # We can't easily retrieve the sign, so we test both + diff = simplify(X.obj - Y.obj) + if diff != S.Zero: + diff = simplify(X.obj + Y.obj) + + self.assertTrue(diff == S.Zero, "\n%s\n==\n%s" % (X, Y)) + + def assertNotEqual(self, first, second): + """ + Compare two expressions are not equals. + """ + if isinstance(first, Mv): + first = first.obj + + if isinstance(second, Mv): + second = second.obj + + # We need to help sympy a little... + first = Simp.apply(expand(first)) + second = Simp.apply(expand(second)) + + # Check + diff = simplify(first - second) + + self.assertTrue(diff != 0, "\n%s\n!=\n%s\n%s" % (first, second, diff))