diff --git a/numga/backend/jax/operator.py b/numga/backend/jax/operator.py index c100e57..77a1182 100644 --- a/numga/backend/jax/operator.py +++ b/numga/backend/jax/operator.py @@ -35,10 +35,6 @@ class JaxOperator(AbstractConcreteOperator): def shape(self): """Shape of the kernel, modulo subspace axis""" return self.kernel.shape[:-(self.arity+1)] - # @property - # def kernel_shape(self): - # """Shape of the kernel, modulo subspace axis""" - # return self.kernel.shape[self.arity:] def broadcast_allocate(self, inputs: Tuple[JaxMultiVector], output=None) -> JaxMultiVector: """Allocation for a set of inputs, with multivector components as last axis, @@ -124,19 +120,20 @@ class JaxDenseOperator(JaxOperator): def __init__(self, *args, **kwargs): super(JaxDenseOperator, self).__init__(*args, **kwargs) - # put kernel on device - # self.jax_kernel = np.array(self.kernel) # precompute reshape operations self.shapes = self.broadcasting_shapes + # contraction over all kernel input axes + self.sum_axes = tuple(-(a + 2) for a in range(self.arity)) @partial(jax.jit, static_argnums=(0,)) def __call__(self, *inputs: Tuple[JaxMultiVector]) -> JaxMultiVector: - # kernel, shapes = self.precompute - shape = jnp.broadcast_shapes(*(i.shape for i in inputs)) return self.context.multivector( values=jnp.sum( - math.prod((i.values.reshape(i.shape + s) for i, s in zip(inputs, self.shapes)), start=self.kernel), - axis=range(len(shape), len(shape)+self.arity)# FIXME: negative indexing rather than len(shape)? + math.prod( + (i.values.reshape(i.shape + s) for i, s in zip(inputs, self.shapes)), + start=self.kernel + ), + axis=self.sum_axes ), subspace=self.output ) diff --git a/numga/backend/numpy/multivector.py b/numga/backend/numpy/multivector.py index 678c5e5..3206d5c 100644 --- a/numga/backend/numpy/multivector.py +++ b/numga/backend/numpy/multivector.py @@ -71,11 +71,12 @@ def flatten(self): def la_inverse(self): """Inverse of x such that x * x.inverse() == 1 == x.inverse() * x""" - op = self.operator.product(self.subspace, self.subspace) + inverse_subspace = self.operator.inverse_factor(self.subspace).output + op = self.operator.product(self.subspace, inverse_subspace) k = op.partial({0: self}).kernel idx, = np.flatnonzero(op.output.blades == 0) # grab index of scalar of output; zero or raises r = np.linalg.solve( # use least squares to solve for inverse np.einsum('...ji,...ki->...jk', k, k), # k.T * k k[..., idx], #equal to k.T * unit_scalar ) - return self.copy(r) + return self.context.multivector(values=r, subspace=inverse_subspace) diff --git a/numga/backend/numpy/operator.py b/numga/backend/numpy/operator.py index 757a6d8..25da505 100644 --- a/numga/backend/numpy/operator.py +++ b/numga/backend/numpy/operator.py @@ -97,18 +97,6 @@ def transpose(self): axes=(self.operator.axes[1], self.operator.axes[0]), ) ) - # def quadratic(self, l, r=None): - # """Compute quadratic form of unary operator; l*O*r""" - # assert self.arity == 1 - # r = l if r is None else r - # output = self.broadcast_allocate((l, r), self.algebra.subspace.scalar()) - # np.einsum( - # '...ij,...i,...j->...', - # self.kernel, l.values, r.values, - # out=output.values[..., 0], - # optimize=True - # ) - # return output class NumpyEinsumOperator(NumpyOperator): diff --git a/numga/backend/test/test_numpy.py b/numga/backend/test/test_numpy.py index cf8dee3..0930ebd 100644 --- a/numga/backend/test/test_numpy.py +++ b/numga/backend/test/test_numpy.py @@ -4,22 +4,30 @@ from numga.backend.numpy.operator import NumpySparseOperator import pytest +from numga.multivector.test.util import random_motor, random_subspace + def test_basic(): print() ga = NumpyContext('x+y+z+w0') Q, V, B = ga.subspace.even_grade(), ga.subspace.vector(), ga.subspace.bivector() - q, v = ga.multivector(Q), ga.multivector(V) + q, v = random_motor(ga), random_subspace(ga, V) print(q) print(v) - output = q.sandwich(v) + op = ga.operator.sandwich(Q, V) + op_partial = op.partial({0: q, 2: q}) + op_partial2 = q.sandwich_map(V) # same as above + assert np.allclose(op_partial.kernel, op_partial2.kernel) + + r0 = (q * v * ~q).select_subspace(v.subspace) # direct sandwich product + r1 = q.sandwich(v) # using optimized sandwich operator + r2 = op_partial(v) # using pre-bound q argument + r3 = op_partial2(v) - op = q.sandwich_map() - print('our projective matrix') - print(op.operator.kernel) - print(op.operator.axes) - print(op(v)) + assert np.allclose(r0.values, r1.values) + assert np.allclose(r0.values, r2.values) + assert np.allclose(r0.values, r3.values) def print_op(op): print(op.operator.axes) @@ -108,11 +116,29 @@ def test_operator_composition(): print() -def test_inverse(): - def check_inverse(x, i): - assert np.allclose((x * i - 1).values, 0, atol=1e-9) - assert np.allclose((i * x - 1).values, 0, atol=1e-9) +def check_inverse(x, i): + assert np.allclose((x * i - 1).values, 0, atol=1e-9) + assert np.allclose((i * x - 1).values, 0, atol=1e-9) + + +def test_inverse_counterexample(): + """test counterexample that inverse_factor method fails to solve + afaik, the inverse_factor method should work for all multivectors < 6d + """ + ga = NumpyContext('x+y+z+a+b+c+') + mv = ga.multivector + x = 2 + mv.xy + mv.ab + mv.xyzabc + + i = x.la_inverse() + check_inverse(x, i) + with pytest.raises(Exception): + i = x.inverse() + check_inverse(x, i) + + +def test_inverse(): + """test some general inversion cases""" ga = NumpyContext('x+y+z+w+') V = ga.subspace.vector() x = ga.multivector.vector(values=np.random.normal(size=(2, len(V)))) @@ -139,6 +165,8 @@ def check_inverse(x, i): q = x.inverse_factor() check_inverse(x, q / x.scalar_product(q)) + +def test_inverse_degenerate(): with pytest.raises(Exception): ga = NumpyContext('x+w0') x = ga.multivector.w diff --git a/numga/operator/test/test_operator.py b/numga/operator/test/test_operator.py index fcff742..46a7cfc 100644 --- a/numga/operator/test/test_operator.py +++ b/numga/operator/test/test_operator.py @@ -78,18 +78,20 @@ def test_basic(): def test_commutator(): """Visualize the output grade of the commutator of pairs of i-j vectors""" print() - algebra = Algebra('x+y+z+w+') + algebra = Algebra.from_pqr(8, 0, 0) n = algebra.n_dimensions + 1 dims = range(n) v = [algebra.subspace.k_vector(i) for i in dims] - r = -np.ones((n, n)) + r = np.empty((n, n), dtype=object) + for i in dims: for j in dims: try: - r[i, j] = algebra.operator.anti_commutator(v[i], v[j]).output.grade() + s = str(np.unique(algebra.operator.commutator(v[i], v[j]).output.grades())) + r[i, j] = s.rjust(5, ' ') except: pass - print(r) + print('\n'.join(' '.join(q) for q in r)) def test_square_signs(): @@ -103,6 +105,7 @@ def test_square_signs(): from numga.multivector.test.util import random_subspace import numpy.testing as npt + def test_inertia(): """Test equivalence of composed ternary operators to their direct expression form""" from numga.backend.numpy.context import NumpyContext diff --git a/numga/subspace/base.py b/numga/subspace/base.py index d5609f5..ba12964 100644 --- a/numga/subspace/base.py +++ b/numga/subspace/base.py @@ -39,7 +39,7 @@ def n_blades(self): def __hash__(self): return hash((id(self.algebra), hash(self.blades.tostring()))) def __eq__(self, other: "SubSpaceInterface"): - return (self is other) or ((self.algebra is other.algebra) and np.alltrue(self.blades == other.blades)) + return (self is other) or ((self.algebra is other.algebra) and np.all(self.blades == other.blades)) def __contains__(self, other: "SubSpaceInterface"): return set(other.blades) in set(self.blades)