From 469c5e1ca040c1c82898ae8996d5620368e5f4d9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 24 May 2013 01:32:02 -0400 Subject: [PATCH] GA: More build-out, more identities to test. --- pymbolic/geometric_algebra.py | 88 +++++++++++++++++++++++++++++++++-- test/test_pymbolic.py | 73 +++++++++++++++++++++-------- 2 files changed, 138 insertions(+), 23 deletions(-) diff --git a/pymbolic/geometric_algebra.py b/pymbolic/geometric_algebra.py index 3aa3e98..7d0a63c 100644 --- a/pymbolic/geometric_algebra.py +++ b/pymbolic/geometric_algebra.py @@ -240,7 +240,7 @@ class _LeftContractionProduct(_GAProduct): def orthogonal_blade_product_weight(a_bits, b_bits, space): shared_bits = a_bits & b_bits - if shared_bits == b_bits: + if shared_bits == a_bits: return _shared_metric_coeff(shared_bits, space) else: return 0 @@ -255,7 +255,7 @@ class _RightContractionProduct(_GAProduct): def orthogonal_blade_product_weight(a_bits, b_bits, space): shared_bits = a_bits & b_bits - if shared_bits == a_bits: + if shared_bits == b_bits: return _shared_metric_coeff(shared_bits, space) else: return 0 @@ -450,7 +450,7 @@ class MultiVector(object): # }}} - # {{{ products + # {{{ multiplicative operators def _generic_product(self, other, product_class): """ @@ -472,7 +472,6 @@ class MultiVector(object): new_bits = sbits ^ obits weight = bpw(sbits, obits, self.space) - print self.space.blade_bits_to_str(sbits), self.space.blade_bits_to_str(obits) if not is_zero(weight): # These are nonzero by definition. coeff = weight * canonical_reordering_sign(sbits, obits) * scoeff * ocoeff @@ -540,8 +539,57 @@ class MultiVector(object): return self._generic_product(other, _ScalarProduct).as_scalar() + def x(self, other): + """Return the commutator product. + + See (1.55) in [HS]. + """ + return (self*other - other*self)/2 + + def __truediv__(self, other): + """Return ``self*(1/other)``. + """ + if not isinstance(other, MultiVector): + other = MultiVector(other, self.space) + + return self*other.inv() + + def __rtruediv__(self, other): + """Return ``other * (1/self)``. + """ + other = MultiVector(other, self.space) + + return other * self.inv() + # }}} + def inv(self): + """Return the *multiplicative inverse* of the blade *self*. + """ + + nsqr = self.norm_squared() + if len(self.data) == 0: + raise ZeroDivisionError + if len(self.data) > 1: + if self.get_pure_grade() in [0, 1, self.space.dimensions]: + return MultiVector(dict( + (bits, coeff/nsqr) for bits, coeff in self.data.iteritems()), + self.space) + + else: + raise NotImplementedError("division by non-blades") + + (bits, coeff), = self.data.iteritems() + + # (1.54) in [HS] + grade = bit_count(bits) + if grade*(grade-1)//2 % 2: + coeff = -coeff + + coeff = coeff/nsqr + + return MultiVector({bits: coeff}, self.space) + def rev(self): """Return the *reverse* of *self*, i.e. the multivector obtained by reversing the order of all component blades. @@ -570,6 +618,11 @@ class MultiVector(object): return MultiVector(new_data, self.space) + def dual(self): + """Return the dual of *self*, see (2.26) in [HS].""" + + return self | self.I.rev() + def norm_squared(self): return self.rev().scalar_product(self) @@ -614,6 +667,15 @@ class MultiVector(object): # {{{ grade manipulation + def gen_blades(self, grade=None): + if grade is None: + for bits, coeff in self.data.iteritems(): + yield MultiVector({bits: coeff}, self.space) + else: + for bits, coeff in self.data.iteritems(): + if bit_count(bits) == grade: + yield MultiVector({bits: coeff}, self.space) + def project(self, grade): new_data = {} for bits, coeff in self.data.iteritems(): @@ -642,6 +704,24 @@ class MultiVector(object): return result + def odd(self): + """Extract the odd-grade blades.""" + new_data = {} + for bits, coeff in self.data.iteritems(): + if bit_count(bits) % 2: + new_data[bits] = coeff + + return MultiVector(new_data, self.space) + + def even(self): + """Extract the even-grade blades.""" + new_data = {} + for bits, coeff in self.data.iteritems(): + if bit_count(bits) % 2 == 0: + new_data[bits] = coeff + + return MultiVector(new_data, self.space) + def as_scalar(self): result = 0 for bits, coeff in self.data.iteritems(): diff --git a/test/test_pymbolic.py b/test/test_pymbolic.py index c06e9fb..43cf32c 100644 --- a/test/test_pymbolic.py +++ b/test/test_pymbolic.py @@ -223,7 +223,12 @@ def test_geometric_algebra(dims): from operator import xor as outer assert reduce(outer, vecs).close_to(0) - for obj1, obj2, obj3 in [ + assert ( vec1.inv()*vec1 ).close_to( 1 ) + assert ( vec1*vec1.inv() ).close_to( 1 ) + assert ( (1/vec1)*vec1 ).close_to( 1 ) + assert ( vec1/vec1 ).close_to( 1 ) + + for a, b, c in [ (vec1, vec2, vec3), (vec1*vec2, vec3, vec4), (vec1, vec2*vec3, vec4), @@ -232,32 +237,62 @@ def test_geometric_algebra(dims): (vec1, vec2*vec1, vec3*vec4*vec5), ]: + # Associativity + assert ((a*b)*c).close_to( + a*(b*c)) + assert ((a^b)^c).close_to( + a^(b^c)) + assert ((a*b)*c).close_to( + a*(b*c)) + # scalar product - assert ( (obj3*obj2).project(0) ) .close_to( obj2.scalar_product(obj3) ) - assert ( (obj3.rev()*obj2).project(0) ) .close_to( obj2.rev().scalar_product(obj3) ) - assert ( (obj2.rev()*obj2).project(0) ) .close_to( obj2.norm_squared() ) + assert ( (c*b).project(0) ) .close_to( b.scalar_product(c) ) + assert ( (c.rev()*b).project(0) ) .close_to( b.rev().scalar_product(c) ) + assert ( (b.rev()*b).project(0) ) .close_to( b.norm_squared() ) - assert obj2.norm_squared() >= 0 - assert obj3.norm_squared() >= 0 + assert b.norm_squared() >= 0 + assert c.norm_squared() >= 0 # Cauchy's inequality - assert obj2.scalar_product(obj3) <= abs(obj2)*abs(obj3) + 1e-13 + assert b.scalar_product(c) <= abs(b)*abs(c) + 1e-13 + + # contractions + + # (3.18) in [DFM] + assert abs(b.scalar_product(a ^ c) - (b >> a).scalar_product(c)) < 1e-13 + # duality, (3.20) in [DFM] + assert ((a ^ b) << c) .close_to( a << (b << c) ) + + # inverse + for div in list(b.gen_blades()) + [vec1, vec1.I]: + assert ( div.inv()*div ).close_to( 1 ) + assert ( div*div.inv() ).close_to( 1 ) + assert ( (1/div)*div ).close_to( 1 ) + assert ( div/div ).close_to( 1 ) + assert ((c/div)*div ).close_to( c ) + assert ((c*div)/div).close_to( c ) - # reverse/dual properties (Sec 2.9.5 DFW) - assert obj3.rev().rev() == obj3 - assert (obj2^obj3).rev() .close_to( (obj3.rev() ^ obj2.rev()) ) + # reverse properties (Sec 2.9.5 [DFM]) + assert c.rev().rev() == c + assert (b^c).rev() .close_to( (c.rev() ^ b.rev()) ) + + # dual properties + # (2.26) in [HS] + assert c.dual() .close_to( c|c.I.rev() ) + assert c.dual() .close_to( c*c.I.rev() ) # involution properties (Sec 2.9.5 DFW) - assert obj3.invol().invol() == obj3 - assert (obj2^obj3).invol() .close_to( (obj2.invol() ^ obj3.invol()) ) + assert c.invol().invol() == c + assert (b^c).invol() .close_to( (b.invol() ^ c.invol()) ) + + # commutator properties + + # Jacobi identity (1.56c) in [HS] + assert ( a.x(b.x(c)) + b.x(c.x(a)) + c.x(a.x(b)) ).close_to(0) + # (1.57) in [HS] + assert a.x(b*c) .close_to( a.x(b)*c + b*a.x(c)) + - # Associativity - assert ((obj1*obj2)*obj3).close_to( - obj1*(obj2*obj3)) - assert ((obj1^obj2)^obj3).close_to( - obj1^(obj2^obj3)) - assert ((obj1*obj2)*obj3).close_to( - obj1*(obj2*obj3)) -- GitLab