From 7a7710826c7e7f1183ab4996ddeba1fe1f7233f9 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Fri, 24 May 2013 00:05:08 -0400
Subject: [PATCH] Fix signs, Cauchy's inequality in GA.

---
 pymbolic/geometric_algebra.py | 28 ++++++----------------------
 test/test_pymbolic.py         | 21 +++++++++++++--------
 2 files changed, 19 insertions(+), 30 deletions(-)

diff --git a/pymbolic/geometric_algebra.py b/pymbolic/geometric_algebra.py
index 3fd55ed..3aa3e98 100644
--- a/pymbolic/geometric_algebra.py
+++ b/pymbolic/geometric_algebra.py
@@ -196,10 +196,7 @@ class _GAProduct(object):
 class _OuterProduct(_GAProduct):
     @staticmethod
     def generic_blade_product_weight(a_bits, b_bits, space):
-        if a_bits & b_bits:
-            return 0
-        else:
-            return canonical_reordering_sign(a_bits, b_bits)
+        return int(not a_bits & b_bits)
 
     orthogonal_blade_product_weight = generic_blade_product_weight
 
@@ -211,13 +208,12 @@ class _GeometricProduct(_GAProduct):
 
     @staticmethod
     def orthogonal_blade_product_weight(a_bits, b_bits, space):
-        result = canonical_reordering_sign(a_bits, b_bits)
         shared_bits = a_bits & b_bits
 
         if shared_bits:
-            result *= _shared_metric_coeff(shared_bits, space)
-
-        return result
+            return _shared_metric_coeff(shared_bits, space)
+        else:
+            return 1
 
 class _InnerProduct(_GAProduct):
     @staticmethod
@@ -230,9 +226,6 @@ class _InnerProduct(_GAProduct):
         shared_bits = a_bits & b_bits
 
         if shared_bits == a_bits or shared_bits == b_bits:
-            # No reordering (and hence no sign computation) needs to be
-            # performed, since all the inner product does is zap basis
-            # elements out of a blade.
             return _shared_metric_coeff(shared_bits, space)
         else:
             return 0
@@ -248,9 +241,6 @@ class _LeftContractionProduct(_GAProduct):
         shared_bits = a_bits & b_bits
 
         if shared_bits == b_bits:
-            # No reordering (and hence no sign computation) needs to be
-            # performed, since all the contraction product does is zap basis
-            # elements out of a blade.
             return _shared_metric_coeff(shared_bits, space)
         else:
             return 0
@@ -266,9 +256,6 @@ class _RightContractionProduct(_GAProduct):
         shared_bits = a_bits & b_bits
 
         if shared_bits == a_bits:
-            # No reordering (and hence no sign computation) needs to be
-            # performed, since all the contraction product does is zap basis
-            # elements out of a blade.
             return _shared_metric_coeff(shared_bits, space)
         else:
             return 0
@@ -282,9 +269,6 @@ class _ScalarProduct(_GAProduct):
     @staticmethod
     def orthogonal_blade_product_weight(a_bits, b_bits, space):
         if a_bits == b_bits:
-            # No reordering (and hence no sign computation) needs to be
-            # performed, since all the contraction product does is zap basis
-            # elements out of a blade.
             return _shared_metric_coeff(a_bits, space)
         else:
             return 0
@@ -486,12 +470,12 @@ class MultiVector(object):
         for sbits, scoeff in self.data.iteritems():
             for obits, ocoeff in other.data.iteritems():
                 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 * scoeff * ocoeff
+                    coeff = weight * canonical_reordering_sign(sbits, obits) * scoeff * ocoeff
                     new_coeff = new_data.setdefault(new_bits, 0) + coeff
                     if is_zero(new_coeff):
                         del new_data[new_bits]
diff --git a/test/test_pymbolic.py b/test/test_pymbolic.py
index cdce639..c06e9fb 100644
--- a/test/test_pymbolic.py
+++ b/test/test_pymbolic.py
@@ -228,21 +228,20 @@ def test_geometric_algebra(dims):
             (vec1*vec2, vec3, vec4),
             (vec1, vec2*vec3, vec4),
             (vec1, vec2, vec3*vec4),
-            (vec1, vec2, vec3*vec4),
             (vec1, vec2, vec3*vec4*vec5),
             (vec1, vec2*vec1, vec3*vec4*vec5),
             ]:
 
         # 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() )
 
-        # FIXME: still broken
-        #assert obj2.norm_squared() >= 0
-        #assert obj3.norm_squared() >= 0
+        assert obj2.norm_squared() >= 0
+        assert obj3.norm_squared() >= 0
 
         # Cauchy's inequality
-        # FIXME: still broken
-        #assert obj2.scalar_product(obj3) <= abs(obj2)*abs(obj3) + 1e-13
+        assert obj2.scalar_product(obj3) <= abs(obj2)*abs(obj3) + 1e-13
 
         # reverse/dual properties (Sec 2.9.5 DFW)
         assert obj3.rev().rev() == obj3
@@ -267,14 +266,20 @@ def playground():
     import numpy as np
     from pymbolic.geometric_algebra import MultiVector as MV
 
-    dims = 3
+    dims = 2
     vec1 = MV(np.random.randn(dims))
     vec2 = MV(np.random.randn(dims))
     vec3 = MV(np.random.randn(dims))
     vec4 = MV(np.random.randn(dims))
     vec5 = MV(np.random.randn(dims))
 
-    print (vec3*vec4).norm_squared()
+    a = vec3^vec4
+    print (a.rev()*a).project(0)
+    print a.scalar_product(a.rev())
+
+    #print a.norm_squared()
+    #print ((a.rev()*a).project(0) ).close_to( a.norm_squared() )
+
 
 
 
-- 
GitLab