generalized symmetry operator to arbitrary number of axes

numga/algebra/test/test_algebra.py

@@ -179,7 +179,7 @@ def test_norm_squared():
 	kv = ga.subspace.k_vector(grade)
 	# kv = ga.subspace.from_grades([0, 4, 8])
 	# construct symbolically simplified operator for (x * ~x)
-	S = ga.operator.reverse_product(kv, kv).symmetry((0, 1), +1)
+	S = ga.operator.reverse_product(kv, kv).symmetry((0, 1))
 	print(np.unique(ga.grade(S.output.blades)))
 	assert S.output == ga.subspace.from_grades(output_grades)
 

numga/operator/factory.py

@@ -95,19 +95,19 @@ def restrict(self, i: SubSpace, o: SubSpace) -> Operator:
 	# 	sign = parity_to_sign([g in grades for g in subspace.grades()])
 	# 	return self.diagonal(subspace, sign)
 	@cache
-	def scalar_negation(self, subspace) -> "Operator":
-		sign = parity_to_sign(subspace.grades() > 0)
-		return self.diagonal(subspace, sign)
+	def scalar_negation(self, v) -> "Operator":
+		sign = parity_to_sign(v.subspace.grades() > 0)
+		return self.diagonal(v, sign)
 	@cache
-	def pseudoscalar_negation(self, subspace) -> "Operator":
+	def pseudoscalar_negation(self, v) -> "Operator":
 		# grades = (1, subspace.algebra.n_dimensions - 1)
 		if self.algebra.n_dimensions % 2 == 1:
-			grades = (1, subspace.algebra.n_dimensions - 1)
+			grades = (1, v.subspace.algebra.n_dimensions - 1)
 		else:
 			grades = (1, )
-		grades = (0, subspace.algebra.n_dimensions)
-		sign = -parity_to_sign(np.array([int(g) in grades for g in subspace.grades()]))
-		return self.diagonal(subspace, sign)
+		grades = (0, v.subspace.algebra.n_dimensions)
+		sign = -parity_to_sign(np.array([int(g) in grades for g in v.subspace.grades()]))
+		return self.diagonal(v, sign)
 
 
 	@cache
@@ -425,30 +425,30 @@ def anti_reverse_product(self, l: SubSpace, r: SubSpace) -> Operator:
 	@cache
 	def symmetric_reverse_product(self, x: SubSpace) -> Operator:
 		"""x * ~x"""
-		return self.reverse_product(x, x).symmetry((0, 1), +1)
+		return self.reverse_product(x, x).symmetry((0, 1))
 	@cache
 	def symmetric_involute_product(self, x: SubSpace) -> Operator:
 		"""x * x.involute()"""
-		return self.involute_product(x, x).symmetry((0, 1), +1)
+		return self.involute_product(x, x).symmetry((0, 1))
 	@cache
 	def symmetric_conjugate_product(self, x: SubSpace) -> Operator:
 		"""x * x.conjugate()"""
-		return self.conjugate_product(x, x).symmetry((0, 1), +1)
+		return self.conjugate_product(x, x).symmetry((0, 1))
 	@cache
 	def symmetric_study_conjugate_product(self, x: SubSpace) -> Operator:
 		"""x * x.conjugate()"""
-		return self.study_conjugate_product(x, x).symmetry((0, 1), +1)
+		return self.study_conjugate_product(x, x).symmetry((0, 1))
 	@cache
 	def symmetric_scalar_negation_product(self, x: SubSpace) -> Operator:
-		return self.scalar_negation_product(x, x).symmetry((0, 1), +1)
+		return self.scalar_negation_product(x, x).symmetry((0, 1))
 	@cache
 	def symmetric_pseudoscalar_negation_product(self, x: SubSpace) -> Operator:
-		return self.pseudoscalar_negation_product(x, x).symmetry((0, 1), +1)
+		return self.pseudoscalar_negation_product(x, x).symmetry((0, 1))
 
 	@cache
 	def squared(self, v: SubSpace) -> Operator:
 		"""x * x"""
-		return self.product(v, v).symmetry((0, 1), +1)
+		return self.product(v, v).symmetry((0, 1))
 
 	# FIXME: should we introduce dual/anti norms as well?
 	# @cache
@@ -493,7 +493,7 @@ def full_sandwich(self, R: SubSpace, v: SubSpace) -> Operator:
 		however if the motor being sandwiched is indeed satisfies the requirements of a motor,
 		the product will be grade preserving and the grade-5 part will be numericailly zero.
 		"""
-		return self.product(self.product(R, v), self.reverse(R)).symmetry((0, 2), +1)
+		return self.product(self.product(R, v), self.reverse(R)).symmetry((0, 2))
 		# return self.product(R, self.product(v, self.reverse(R))
 
 	# @cache

numga/operator/operator.py

@@ -64,10 +64,10 @@ def swapaxes(self, a, b) -> "Operator":
 			tuple(axes)
 		)
 
-	def is_diagonal(self, axes: Tuple[int, int]) -> bool:
-		"""Test if a particular axes pair represents a diagonal"""
-		sym = self.symmetry(axes, -1)
-		return np.all(sym.kernel == 0)
+	# def is_diagonal(self, axes: Tuple[int, int]) -> bool:
+	# 	"""Test if a particular axes pair represents a diagonal"""
+	# 	sym = self.symmetry(axes, -1)
+	# 	return np.all(sym.kernel == 0)
 
 	@cached_property
 	def arity(self) -> int:
@@ -105,22 +105,27 @@ def __mul__(self, s):
 	def __rmul__(self, s):
 		return self.mul(s)
 
-	def symmetry(self, axes: Tuple[int, int], sign: int) -> "Operator":
-		"""Enforce commutativity symmetry of the given sign on the given axes [a, b]
+	def symmetry(self, axes: Tuple[int]) -> "Operator":
+		"""Enforce symmetry on the given axes [a, b]
 		That is the output op will have the property
-			op(a, b) = sign * op(b, a)
+			op(a, b) = op(b, a)
 
 		We may choose to apply this to our operator,
 		to encode knowledge about inputs being numerically identical,
 		such as encountered in the sandwich product or norm calculations.
+		which usually allows simplification of the resulting expressions
 		"""
-		# FIXME: is it appropriate to view everything as real numbers at this point?
-		#  i think so; but i keep confusing myself
-		assert self.axes[axes[0]] == self.axes[axes[1]]
-		kernel2 = self.kernel + sign * np.swapaxes(self.kernel, *axes)
-		kernel = kernel2 // 2
-		if not np.all(kernel * 2 == kernel2):
-			raise Exception('inappropriate integer division')
+		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)
+
+		kernel = k // len(P)
+		if not np.all(kernel * len(k) == k):
+			kernel = k / len(P) # fallback to float division if required
 		return Operator(kernel, self.axes).squeeze()
 
 	def fuse(self, other: "Operator", axis: int) -> "Operator":