more optimized motor examples and test coverage

numga/multivector/extension/decompose.py

@@ -71,10 +71,30 @@ def motor_rotor(m) -> "Translator":
 	""""""
 	return m.nondegenerate()
 
-mv.motor_split = SubspaceDispatch("""Split motor wrt a given origin""")
-@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.vector())
-def motor_split(m: Motor, o: Vector) -> Tuple[Motor, Motor]:
-	o = o.dual()
+mv.motor_split = SubspaceDispatch("""
+	Split motor wrt a given origin in two components;
+	one that leaves the origin invariant and one
+	complimentary to that
+	""")
+@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.antivector() and o.dual().is_degenerate)
+def motor_split_euclidian(m: Motor, o: AntiVector) -> Tuple[Motor, Motor]:
+	"""If the dual of o is degenerate, it represents the euclidean origin
+	In this special case there is a more optimized solution
+	m = t * r
+	"""
+	return motor_translator(m), motor_rotor(m)
+@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.antivector() and o.is_blade)
+def motor_split_blade(m: Motor, o: AntiVector) -> Tuple[Motor, Motor]:
+	"""If the dual of o is a blade, we can also specialize
+	"""
+	# drop translational axes
+	R = m.subspace.reject(o.subspace.dual())
+	T = m.context.subspace.vector().product(o.subspace.dual())
+	t, r = motor_split(m, o)
+	return t.select_subspace(T), r.select_subspace(R)
+@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.antivector())
+def motor_split(m: Motor, o: AntiVector) -> Tuple[Motor, Motor]:
 	# construct shortest trajectory from o to m >> o
 	t = ((m >> o) / o).motor_square_root()
-	return t, ~t * m
+	# the residual is a rotation that leaves o invariant
+	return t, ~t * m

numga/multivector/extension/norms.py

@@ -68,4 +68,5 @@ def reverse_study_normalized(x):
 	"""Binds to 4d and 5d even grade, as well as many other objects"""
 	s = x.symmetric_reverse_product()
 	s = s.square_root().inverse()
+	# NOTE: order of s and x here is important!
 	return s * x

numga/multivector/extension/optimized.py

@@ -17,6 +17,8 @@
 [2] May the Forque Be with You; Dynamics in PGA
 
 """
+import numpy as np
+
 # FIXME: move this to numpy backend module?
 from numga.multivector.multivector import AbstractMultiVector as mv
 @mv.normalized.register(
@@ -54,3 +56,75 @@ def normalize_3dpga_motor_numpy(m: "Motor") -> "Motor":
 
 # reset cache
 mv.normalized.cache = {}
+@mv.exp.register(
+	lambda s:
+	s.named_str == 'xw,yw,zw' and
+	s.algebra.description.description_str == 'x+y+z+w0',
+	position=0      # try and match this before all else
+)
+def exponentiate_degenerate_bivector_numpy(b: "BiVector") -> "Motor":
+	"""Optimized 3d pga exponentiation, for the given memory layout
+	"""
+	M = np.empty(b.shape + (4,))
+	M[..., 0] = 1
+	M[..., 1:] = b.values
+	return b.context.multivector(values=M, subspace=b.context.subspace.translator())
+
+@mv.exp.register(
+	lambda s:
+	s.named_str == 'xy,xz,yz' and
+	s.algebra.description.description_str == 'x+y+z+w0',
+	position=0      # try and match this before all else
+)
+def exponentiate_nondegenerate_bivector_numpy(b: "BiVector") -> "Motor":
+	"""Optimized 3d pga exponentiation, for the given memory layout
+	"""
+	xy, xz, yz = np.moveaxis(b.values, -1, 0)
+	l = xy*xy + xz*xz + yz*yz
+	a = np.sqrt(l)
+	c = np.cos(a)
+	s = np.where(a > 1e-20, np.sin(a) / a, 1)
+
+	M = np.empty(b.shape + (4,))
+	M[..., 0] = c
+	M[..., 1] = s*xy
+	M[..., 2] = s*xz
+	M[..., 3] = s*yz
+	return b.context.multivector(values=M, subspace=b.context.subspace.rotor())
+
+@mv.exp.register(
+	lambda s:
+	s.named_str == 'xy,xz,yz,xw,yw,zw' and
+	s.algebra.description.description_str == 'x+y+z+w0',
+	position=0      # try and match this before all else
+)
+def exponentiate_bivector_numpy(b: "BiVector") -> "Motor":
+	"""Optimized 3d pga exponentiation, for the given memory layout
+
+	Notes
+	-----
+	Terms involving xz are negated compared to the reference [1],
+	since our normalized blade order is to sort the axes alphabetically,
+	and our subspace object does not currently store signs of basis blades
+	"""
+	xy, xz, yz, xw, yw, zw = np.moveaxis(b.values, -1, 0)
+	l = xy*xy + xz*xz + yz*yz
+	a = np.sqrt(l)
+	m = xy*zw - xz*yw + yz*xw
+	c = np.cos(a)
+	s = np.where(a > 1e-20, np.sin(a) / a, 1)
+	t = m * (c - s) / l
+
+	M = np.empty(b.shape + (8,))
+	M[..., 0] = c
+	M[..., 1] = s*xy
+	M[..., 2] = s*xz
+	M[..., 3] = s*yz
+	M[..., 4] = s*xw + t*yz
+	M[..., 5] = s*yw - t*xz
+	M[..., 6] = s*zw + t*xy
+	M[..., 7] = m*s
+	return b.context.multivector(values=M, subspace=b.context.subspace.motor())
+
+# reset cache
+mv.exp.cache = {}

numga/multivector/extension/test/test_decompose.py

@@ -14,57 +14,74 @@
 	(5, 0, 0), (4, 0, 1), (4, 1, 0), (3, 1, 1), #(2, 2, 1),
 ])
 def test_invariant_decomposition(descr):
+	np.random.seed(1)
 	print()
 	print(descr)
 	algebra = Algebra.from_pqr(*descr)
 	ga = NumpyContext(algebra)
 	for i in range(10):
 		b = random_motor(ga).select[2]
 
-		print('decompose')
 		l, r = b.decompose_invariant()
-		# if ga.algebra.pseudo_scalar_squared == 0:
-		# 	assert r.subspace.is_degenerate
-		q = l.inner(r)
-		npt.assert_allclose(q.values, 0, atol=1e-9)
-		q = l.squared() # test they square to a scalar
-		npt.assert_allclose(q.values[1:], 0, atol=1e-9)
-		q = r.squared()
-		npt.assert_allclose(q.values[1:], 0, atol=1e-9)
+		assert_close(l.inner(r), 0)
+		assert_close(l.squared().select.nonscalar(), 0)
+		assert_close(r.squared().select.nonscalar(), 0)
 
-		npt.assert_allclose(b.values, (l+r).values, atol=1e-9)
+		assert_close(b, l + r)
 
 		L = l.exp()
 		R = r.exp()
 		# test that component bivectors commute under exponentiation
-		npt.assert_allclose((L * R).values, (R * L).values)
+		assert_close(L * R, R * L)
 
 
 @pytest.mark.parametrize('descr', [
 	(1, 0, 1), (2, 0, 1), (3, 0, 1), (4, 0, 1),
 ])
 def test_motor_decompose_euclidean(descr):
-
 	ga = NumpyContext(descr)
 	m = random_motor(ga, (10,))
 	t = m.motor_translator()
-	print(t)
+	# print(t)
 	r = m.motor_rotor()
 
-	assert np.allclose((t * r - m).values, 0, atol=1e-9)
+	assert_close(t * r, m)
 	# alternate method of constructing translator
 	tt = (m * ~r).select.translator()
-	assert np.allclose((tt - t).values, 0, atol=1e-9)
-
-
-def test_motor_decompose_elliptical():
-	ga = NumpyContext((4,0,0))
-	m = random_motor(ga, (10,))
-
-	t, r = m.motor_split(ga.multivector.a)
+	assert_close(tt, t)
 
-	assert np.allclose((t * r - m).values, 0, atol=1e-9)
 
-	assert np.allclose((t.symmetric_reverse_product() - 1).values, 0, atol=1e-9)
-	assert np.allclose((r.symmetric_reverse_product() - 1).values, 0, atol=1e-9)
+@pytest.mark.parametrize('descr', [
+	(2, 0, 0), (3, 0, 0), (4, 0, 0), (5, 0, 0),
+	(1, 0, 1), (2, 0, 1), (3, 0, 1), (4, 0, 1),
+	(1, 1, 0), (2, 1, 0), (3, 1, 0), (4, 1, 0),
+	# (3, 0, 0)
+])
+def test_motor_split(descr):
+	"""test splitting of motors relative to a specific point"""
+	np.random.seed(10)
+	ga = NumpyContext(descr)
+	N = 10
+	m = random_motor(ga, (N,))
+
+	# test for both canonical origin, or some arbitrary point
+	o = ga.multivector.basis()[-1].dual()
+	os = [o, random_motor(ga, (N,)) >> o]
+
+	for o in os:
+		t, r = m.motor_split(o)
+		# print(m)
+		# print(t.subspace)
+		# print(t.subspace.is_subalgebra)
+		# print(t.subspace)
+		# print(r.subspace.is_subalgebra)
+		# check that it is a valid decomposition of the original motor
+		assert_close(t * r, m)
+		# check that r leaves o invariant
+		assert_close(r >> o, o, atol=1e-6)
+		# in d<6, both of these should be simple motors
+		assert_close(t * ~t, 1, atol=1e-6)
+		assert_close(r * ~r, 1, atol=1e-6)
+		# and their respective bivectors are orthogonal
+		assert_close(t.motor_log().inner(r.motor_log()), 0)
 

numga/multivector/extension/test/test_logexp.py

@@ -1,4 +1,4 @@
-
+import numpy as np
 import pytest
 
 from numga.backend.numpy.context import NumpyContext
@@ -10,12 +10,13 @@
 
 @pytest.mark.parametrize('descr', [
 	(2, 0, 0), (1, 0, 1), (1, 1, 0), #(0, 1, 1),
-	(3, 0, 0), (2, 0, 1), (2, 1, 0), (1, 1, 1), (1, 0, 2), (1, 2, 0),
+	(3, 0, 0), (2, 0, 1), (2, 1, 0), (1, 1, 1), (1, 0, 2), #(1, 2, 0),
 	(4, 0, 0), (3, 0, 1), (3, 1, 0), (2, 1, 1), (2, 0, 2), #(2, 2, 0),
 	(5, 0, 0), (4, 0, 1), (4, 1, 0), (3, 1, 1),
 ])
 def test_bisect(descr):
 	"""Test roundtrip qualities of bisection based log and exp"""
+	np.random.seed(0)
 	print()
 	print(descr)
 	algebra = Algebra.from_pqr(*descr)
@@ -25,12 +26,12 @@ def test_bisect(descr):
 
 	# check exact inverse props
 	r = m.motor_log().exp()
-	npt.assert_allclose(m.values, r.values, atol=1e-9)
+	assert_close(m, r)
 
 	# check for low iteration count
 	from numga.multivector.extension.logexp import exp_bisect, motor_log_bisect
 	r = exp_bisect(motor_log_bisect(m, 2), 2)
-	npt.assert_allclose(m.values, r.values, atol=1e-9)
+	assert_close(m, r)
 
 
 def test_interpolate():
@@ -89,6 +90,4 @@ def test_exp_high_dim():
 		r = h.squared() * i
 	# else:
 	# 	r = h.squared().solve(q)
-	res = r.symmetric_reverse_product() - 1
-	print(res)
-	npt.assert_allclose(res.values, 0, atol=1e-9)
+	assert_close(r.symmetric_reverse_product(), 1)

numga/multivector/extension/test/test_norms.py

@@ -27,27 +27,25 @@ def test_study(descr):
 	s = m.symmetric_reverse_product()
 	assert s.subspace.is_study
 
+	mn = m.normalized()
+	assert_close(mn * ~mn, 1)
+
 	print('inverse')
 	Is = s.inverse()
-	r = Is * s - 1
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(Is * s, 1)
 
 	print('square root')
 	rs = s.square_root()
-	r = rs * rs - s
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(rs * rs, s)
 	zero = s * 0
 	rs = zero.square_root()
-	npt.assert_allclose(rs.values, 0, atol=1e-9)
+	assert_close(rs, 0)
 
 	print('inverse square root')
 	Irs = s.inverse().square_root()
-
 	Irs2 = s.square_root().inverse()
-	npt.assert_allclose(Irs2.values, Irs.values, atol=1e-6)
-
-	r = s * Irs * Irs - 1
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(Irs2, Irs, atol=1e-6)
+	assert_close(s * Irs * Irs, 1)
 
 
 def test_motor_normalize_6d(descr=[6,0,0]):

numga/multivector/extension/test/test_optimized.py

@@ -3,11 +3,12 @@
 import numpy as np
 import numpy.testing as npt
 from numga.backend.numpy.context import NumpyContext as Context
+from numga.multivector.test.util import *
 
 
-def test_normalize_extended():
+def test_normalize():
 	ga = Context('x+y+z+w0')
-	m = ga.multivector.motor(np.random.normal(size=(20, 8)))
+	m = random_non_motor(ga, (10,))
 	print()
 	print(m)
 	print(m.norm())
@@ -23,4 +24,29 @@ def test_normalize_extended():
 	print(n2)
 	print(n2.norm())
 
-	npt.assert_allclose(n1.values, n2.values, atol=1e-12)
+	assert_close(n1, n2)
+
+
+def test_exp():
+	np.random.seed(0)
+	from time import time
+	from numga.multivector.extension import optimized
+	from numga.multivector.extension.logexp import exp_bisect
+	from numga.backend.numpy.operator import NumpyEinsumOperator as Operator
+
+	ga = Context('x+y+z+w0', otype=Operator)
+	B = ga.subspace.bivector()
+
+	subspaces = [B, B.degenerate(), B.nondegenerate()]
+	for s in subspaces:
+		print(s)
+		b = random_subspace(ga, s, (10, 10, ))
+		b = random_subspace(ga, s, ( ))
+		t = time()
+		m1 = exp_bisect(b)
+		print(time() - t)
+		t = time()
+		m2 = b.exp()
+		print(time() - t)
+
+		assert_close(m1, m2, atol=1e-8)