higher order symmetry refinements

numga/backend/test/test_numpy.py

@@ -122,38 +122,66 @@ def check_inverse(x, i, atol=1e-9):
 
 @pytest.mark.parametrize(
 	'descr', [
-		# (1, 0, 0), (0, 1, 0),
-		# (2, 0, 0), (1, 1, 0), (0, 2, 0),
-		# (3, 0, 0), (2, 1, 0), (1, 2, 0),
-		# (4, 0, 0), #(3, 1, 0), (2, 2, 0),
-		(5, 0, 0),# (4, 1, 0), (3, 2, 0)
+		(1, 0, 0), (0, 1, 0),
+		(2, 0, 0), (1, 1, 0), (0, 2, 0),
+		(3, 0, 0), (2, 1, 0), (1, 2, 0),
+		(4, 0, 0), (3, 1, 0), (2, 2, 0),
+		(5, 0, 0), (4, 1, 0), (3, 2, 0)
 	],
+	# 'descr', [(4, 0, 0)],
 )
 def test_inverse_exhaustive(descr):
 	"""Test general inversion cases for all multivector grade combos in dimension < 6"""
 	np.random.seed(0)   # fix seed to prevent chasing ever changing outliers
 	ga = NumpyContext(Algebra.from_pqr(*descr))
 
-	N = 100
+	N = 10
 	print()
 	print(descr)
 
 	all_grades = np.arange(ga.algebra.n_grades)
 	import itertools
 	for r in all_grades:
 		for grades in itertools.combinations(all_grades, r+1):
-			try:
-				V = ga.subspace.from_grades(list(grades))
-				# print()
-				x = random_subspace(ga, V, (N,))
-				i = x.inverse()
-				check_inverse(x, i, atol=1e-10)
+			V = ga.subspace.from_grades(list(grades))
+			print()
+			x = random_subspace(ga, V, (N,))
+			i = x.inverse()
+			check_inverse(x, i, atol=1e-11)
 
-				print(V.simplicity, list(grades), list(np.unique(i.subspace.grades())), end='')
-				print()
-				print(i)
-			except:
-				pass
+			print(V.simplicity, list(grades), list(np.unique(i.subspace.grades())), end='')
+			# print()
+			# print(i)
+
+
+@pytest.mark.parametrize(
+	'descr', [
+		(1, 0, 0), (0, 1, 0),
+		(2, 0, 0), (1, 1, 0), (0, 2, 0),
+		(3, 0, 0), (2, 1, 0), (1, 2, 0),
+	],
+)
+def test_inverse_hitzer(descr):
+	"""Test non recursive hitzer inversion cases for all multivector grade combos in dimension < 4"""
+	np.random.seed(0)   # fix seed to prevent chasing ever changing outliers
+	ga = NumpyContext(Algebra.from_pqr(*descr))
+
+	N = 10
+	print()
+	print(descr)
+
+	all_grades = np.arange(ga.algebra.n_grades)
+	import itertools
+	for r in all_grades:
+		for grades in itertools.combinations(all_grades, r+1):
+			V = ga.subspace.from_grades(list(grades))
+			print()
+			x = random_subspace(ga, V, (N,))
+
+			i = x.inverse_hitzer()
+			check_inverse(x, i, atol=1e-8)
+
+			print(list(grades), list(np.unique(i.subspace.grades())), end='')
 
 
 @pytest.mark.parametrize(
@@ -253,6 +281,10 @@ def test_inverse_simplicifation_failure():
 	assert z.subspace == ga.subspace.from_grades([0, 5])
 	assert np.allclose(z.select[5].values, 0)
 
+	# second-order optimized hitzer term does reduce to scalar
+	op = ga.operator.inverse_factor_completed2(V)
+	assert op.output.equals.scalar()
+
 	V = ga.subspace.from_grades([2])
 	assert V.simplicity == 2    # need two steps; but can do without the extra zeros
 	x = random_subspace(ga, V, (1,))
@@ -278,3 +310,69 @@ def test_inverse_6d():
 
 	# m = random_motor(ga, (1,))
 	# check_inverse(m, m.inverse())
+
+
+# def test_inverse_hitzer2():
+# 	"""see if higher order hitzer terms add value"""
+# 	ga = NumpyContext(Algebra.from_pqr(5,0,0))
+# 	V = ga.subspace.from_grades([1,2,5])
+# 	x = random_subspace(ga, V, (1,))
+# 	q = x.inverse_factor()
+# 	print(q)
+# def test_inverse_hitzer():
+# 	"""see if higher order hitzer terms add value"""
+# 	np.random.seed(0)
+# 	ga = NumpyContext(Algebra.from_pqr(3,0,0))
+# 	V = ga.subspace.from_grades([0,1])
+# 	x = random_subspace(ga, V, (1,))
+# 	y = x#.symmetric_reverse_product()
+# 	q = y.inverse_factor()
+# 	z = x.reverse() * x.scalar_negation() * x
+# 	print(z)
+# 	print(q)
+
+
+@pytest.mark.parametrize(
+	'descr', [
+		# (1, 0, 0), (0, 1, 0),
+		# (2, 0, 0), (1, 1, 0), (0, 2, 0),
+		# (3, 0, 0), (2, 1, 0), (1, 2, 0),
+		(5, 0, 0), #(3, 1, 0), (2, 2, 0),
+		# (5, 0, 0), #(4, 1, 0), (3, 2, 0)
+	],
+)
+def test_inverse_factor_exhaustive(descr):
+	"""Test general inversion cases for all multivector grade combos in dimension < 6"""
+	np.random.seed(0)   # fix seed to prevent chasing ever changing outliers
+	ga = NumpyContext(Algebra.from_pqr(*descr))
+
+	N = 1
+	print()
+	print(descr)
+
+	all_grades = np.arange(ga.algebra.n_grades)
+	import itertools
+	for r in all_grades:
+		for grades in itertools.combinations(all_grades, r+1):
+			try:
+				V = ga.subspace.from_grades(list(grades))
+				print()
+				print('grades: ', grades)
+				x = random_subspace(ga, V, (N,))
+				op = ga.operator.inverse_factor_completed2(V)
+				print(V)
+				print('inv-fac-2', op.output)
+				# print(x.symmetric_reverse_product().inverse_factor().subspace)
+				op=ga.operator(op)
+				y = op(x,x,x,x)
+				print(y)
+				i = x.inverse()
+				assert i.values.size > 0
+				check_inverse(x, i, atol=1e-10)
+				print(i)
+				# check_inverse(x, x.inverse_la(), atol=1e-5)
+
+				# print(V.simplicity, list(grades), list(np.unique(i.subspace.grades())), end='')
+				# print(i)
+			except:
+				pass

numga/operator/factory.py

@@ -48,6 +48,8 @@ class OperatorFactory:
 
 	Caching of al created operators happens on the basis of the subspace arguments,
 	and the subspace class implement a flyweight pattern, to make this efficient.
+	NOTE: we also may bother the cache with operator intermediate arguments
+	these are not flyweighted probably a bad idea?
 	"""
 
 	def __init__(self, algebra: Algebra):
@@ -543,18 +545,39 @@ def inverse_transform(self, R: SubSpace, v: SubSpace) -> Operator:
 	def inverse_factor(self, x: SubSpace) -> Operator:
 		"""Compute inverse_factor(x) = conj(x) * involute(x) * reverse(x)
 
+		For many multivectors, this gives:
 		inverse(x) = (x * inverse_factor(x))<0>
 		"""
-		p = self.product(self.conjugate(x), self.involute(x))
-		q = self.product(p, x.reverse())
-		return q.symmetry((0,1,2))
+		p = self.product(x, self.reverse(x)).symmetry((0, 1))
+		q = self.product(self.reverse(x), self.conjugate(p))
+		return q.symmetry((0, 1, 2))
+
+	@cache
+	def inverse_factor_completed(self, x: SubSpace) -> Operator:
+		"""Compute x * inverse_factor(x)
+		"""
+		p = self.symmetric_reverse_product(x)
+		q = self.product(p, self.conjugate(p)).symmetry((0, 1, 2, 3))
+		return q
+
+	@cache
+	def inverse_factor_completed_alt(self, x: SubSpace) -> Operator:
+		"""Fusing with scalar negation gives different reduction possibilities
+		"""
+		p = self.symmetric_reverse_product(x)
+		q = self.product(p, self.pseudoscalar_negation(p)).symmetry((0, 1, 2, 3))
+		return q
+
+	def compose_symmetry_ops(self, x: SubSpace, op1, op2):
+		"""constructed fused higher order hitzer ops"""
+
 
 	@cache
 	def inertia(self, l: SubSpace, r: SubSpace) -> Operator:
 		"""Compute inertia operator; l.regressive(l x r)"""
 		# FIXME: it is tempting to enforce symmetry here, but it does not improve sparsity,
 		#  and does force us to use a float kernel, so nevermind
-		return self.regressive(l, self.commutator(l, r))#.symmetry((0, 1), +1)
+		return self.regressive(l, self.commutator(l, r))#.symmetry((0, 1))
 
 
 	# # projections according to erik lengyel

numga/operator/operator.py

@@ -106,7 +106,7 @@ def __rmul__(self, s):
 		return self.mul(s)
 
 	def symmetry(self, axes: Tuple[int]) -> "Operator":
-		"""Enforce symmetry on the given axes [a, b]
+		"""Enforce symmetry on the given input axes [a, b]
 		That is the output op will have the property
 			op(a, b) = op(b, a)
 
@@ -118,14 +118,13 @@ def symmetry(self, axes: Tuple[int]) -> "Operator":
 		match([self.axes[a] for a in axes])
 
 		import itertools
-		k = 0
 		P = list(itertools.permutations(axes, len(axes)))
-		for a in P:
-			k = k + np.moveaxis(self.kernel, axes, a)
-
+		k = sum(np.moveaxis(self.kernel, axes, a) for a in P)
 		kernel = k // len(P)
-		if not np.all(kernel * len(k) == k):
-			kernel = k / len(P) # fallback to float division if required
+		if not np.all(kernel * len(P) == k):
+			# inherit nonzero pattern without explicit symmetric structure
+			mask = k.any(axis=tuple(range(self.arity)), keepdims=True)
+			kernel = self.kernel * mask
 		return Operator(kernel, self.axes).squeeze()
 
 	def fuse(self, other: "Operator", axis: int) -> "Operator":
@@ -207,13 +206,17 @@ def restrict_grade(self, k: int) -> "Operator":
 
 	def squeeze(self) -> "Operator":
 		"""Drop known zero terms from output subspace"""
-		indices = np.flatnonzero(self.kernel.any(axis=tuple(range(self.arity))))
+		k = self.kernel
+		# if np.issubdtype(k, np.floating):
+		# 	eps = 1e-6   # need to consider an epsilon to make chained float kernel simplifications work
+		# 	k = np.abs(k) > eps
+		# get nonzero output indices
+		indices = np.flatnonzero(k.any(axis=tuple(range(self.arity))))
 		output: SubSpace = self.output.slice_subspace(indices)
 
 		axes = list(self.axes)
 		axes[-1] = output
-		kernel = self.kernel.take(indices, axis=-1)
-		return Operator(kernel, tuple(axes))
+		return Operator(self.kernel.take(indices, axis=-1), tuple(axes))
 
 	def grade_selection(self, formula) -> "Operator":
 		"""Apply a grade selection formula"""

numga/subspace/subspace.py

@@ -13,7 +13,7 @@
 # FIXME: some texts seem to prefer element over blade; not sure which is more ideomatic
 class SubSpace(FlyweightMixin):
 	"""Describes a subspace of a full geometric algebra,
-	by means of an array of bit-blades, denoting the blades present in the subalgebra
+	by means of an array of bit-blades, denoting the blades present in the subspace
 	"""
 	def __init__(self, algebra: Algebra, blades: np.array):
 		# parent algebra
@@ -124,6 +124,8 @@ def pretty_str(self):
 	# FIXME: these compete with getattr operator construction syntax, perhaps? not working atm anyway tho...
 	#  perhaps forego subspace construction in favor of operator construction; pretty easy to write V.wedge(B).subspace, after all.
 	#  alternatively, write V.operator.wedge(B), if operator is desired?
+	#  note we are hung up on requiring context previously
+	#  just optionally pass in context? only invoked when constructing operators
 	@cache
 	def complement(self) -> "SubSpace":
 		"""Complementary set of blades, such that self * self.complement() == I"""