Adding Double Conformal Geometric Algebra stub and tests (#312)

* Added dg3c stub and tests

* Add dg3c to the docs and to flake8 ignores

Co-authored-by: Eric Wieser <wieser.eric@gmail.com>

clifford/dg3c.py

@@ -0,0 +1,123 @@
+"""
+The Cl(8,2) Double Conformal Geometric Algebra
+
+Easter, R.B., Hitzer, E. Double Conformal Geometric Algebra.
+Adv. Appl. Clifford Algebras 27, 2175–2199 (2017).
+https://doi.org/10.1007/s00006-017-0784-0
+"""
+
+from . import Layout
+
+# The metric of the algebra is two CGAs glued together
+layout = Layout([1]*4+[-1] + [1]*4+[-1])
+blades = layout.bases()
+
+locals().update(blades)
+
+
+# for shorter reprs
+layout.__name__ = 'layout'
+layout.__module__ = __name__
+
+
+# The two pseudo-scalars, infinities and origins
+IC1 = e12345
+IC2 = e678910
+einf1 = e4 + e5
+eo1 = 0.5*(e5 - e4)
+einf2 = e9 + e10
+eo2 = 0.5*(e10 - e9)
+
+# The joint infinity and origin
+eo = eo1^eo2
+einf = einf1^einf2
+
+
+def up_cga1(pnt_vector):
+    """
+    Take a vector and embed it as a point in the first
+    copy of cga
+    """
+    euc_point = pnt_vector[0]*e1 + pnt_vector[1]*e2 + pnt_vector[2]*e3
+    return euc_point + 0.5*(euc_point|euc_point)*einf1 + eo1
+
+
+def down_cga1(point_cga1):
+    """
+    Take a point in CGA
+    """
+    return (point_cga1/-(point_cga1|einf1)[0]).value[1:4]
+
+
+def up_cga2(pnt_vector):
+    """
+    Take a vector and embed it as a point in the second
+    copy of cga
+    """
+    euc_point = pnt_vector[0]*e6 + pnt_vector[1]*e7 + pnt_vector[2]*e8
+    return euc_point + 0.5*euc_point**2*einf2 + eo2
+
+
+def up(pnt_vector):
+    """
+    Take a vector and embed it as a dcga point
+    """
+    return up_cga1(pnt_vector)^up_cga2(pnt_vector)
+
+
+def down(dcga_point):
+    """
+    Take a dcga_point and return the 3d vector it represents
+    """
+    cga_pnt = ((dcga_point|einf2)|IC1)*IC1
+    return down_cga1(cga_pnt)
+
+
+"""
+These cyclide_ops are the elements that make up a general Darboux cyclide
+See Table 1 and Table 2 from Easter, Hitzer, Double Conformal Geometric Algebra (2017)
+"""
+cyclide_ops = {
+    "Tx": 0.5 * (e1 * einf2 + einf1 * e6),
+    "Ty":  0.5 * (e2 * einf2 + einf1 * e7),
+    "Tz":  0.5 * (e3 * einf2 + einf1 * e8),
+
+    "Tx2":  e6 * e1,
+    "Ty2":  e7 * e2,
+    "Tz2":  e8 * e3,
+
+    "Txy":  0.5 * (e7 * e1 + e6 * e2),
+    "Tyz":  0.5 * (e7 * e3 + e8 * e2),
+    "Tzx":  0.5 * (e8 * e1 + e6 * e3),
+
+    "Txt2":  e1 * eo2 + eo1 * e6,
+    "Tyt2":  e2 * eo2 + eo1 * e7,
+    "Tzt2":  e3 * eo2 + eo1 * e8,
+
+    "T1":  -einf,
+    "Tt2":  eo2*einf1 + einf2*eo1,
+    "Tt4":  -4*eo
+}
+
+cyclide_ops_reciprocal = {
+    "Tx":  cyclide_ops['Txt2'],
+    "Ty":  cyclide_ops['Tyt2'],
+    "Tz":  cyclide_ops['Tzt2'],
+
+    "Tx2":  -cyclide_ops['Tx2'],
+    "Ty2":  -cyclide_ops['Ty2'],
+    "Tz2":  -cyclide_ops['Tz2'],
+
+    "Txy":  -2*cyclide_ops['Txy'],
+    "Tyz":  -2*cyclide_ops['Tyz'],
+    "Tzx":  -2*cyclide_ops['Tzx'],
+
+    "Txt2":  cyclide_ops['Tx'],
+    "Tyt2":  cyclide_ops['Ty'],
+    "Tzt2":  cyclide_ops['Tz'],
+
+    "T1":  -cyclide_ops['Tt4']/4,
+    "Tt2":  -cyclide_ops['Tt2']/2,
+    "Tt4":  -cyclide_ops['T1']/4
+}
+

clifford/test/test_dg3c.py

@@ -0,0 +1,243 @@
+import unittest
+import numpy as np
+from ..dg3c import *
+
+"""
+All basic test identities come from:
+
+Easter, R.B., Hitzer, E. Double Conformal Geometric Algebra.
+Adv. Appl. Clifford Algebras 27, 2175–2199 (2017).
+https://doi.org/10.1007/s00006-017-0784-0
+"""
+
+
+class BasicTests(unittest.TestCase):
+    def test_metric(self):
+        """
+        Ensure that the metric comes out with a double copy of
+        the CGA metric
+        """
+        assert np.all(layout.metric == np.array([
+            [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
+            [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
+            [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
+            [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
+            [0., 0., 0., 0., -1., 0., 0., 0., 0., 0.],
+            [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
+            [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
+            [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
+            [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
+            [0., 0., 0., 0., 0., 0., 0., 0., 0., -1.]]))
+
+    def test_up_down(self):
+        """
+        Test that we can map points up and down into the dpga
+        """
+        rng = np.random.RandomState()
+        for i in range(100):
+            pnt_vector = rng.randn(3)
+            pnt = up(pnt_vector)
+            res = down(100*pnt)
+            np.testing.assert_allclose(res, pnt_vector)
+
+    def test_up_down_cga1(self):
+        """
+        Test that we can map points up and down from cga1
+        """
+        rng = np.random.RandomState()
+        pnt_vector = rng.randn(3)
+        for i in range(100):
+            pnt = up_cga1(pnt_vector)
+            res = down_cga1(100*pnt)
+            np.testing.assert_allclose(res, pnt_vector)
+
+
+class GeometricPrimitiveTests(unittest.TestCase):
+    def test_reciprocality(self):
+        """
+        Ensure that the cyclide ops and the reciprocal frame are
+        actually reciprocal...
+        """
+        for key, cyc_op in cyclide_ops.items():
+            for key2, cyc_op_recip in cyclide_ops_reciprocal.items():
+                assert cyc_op|cyc_op_recip == layout.scalar*(key == key2)
+
+    def test_general_elipsoid(self):
+        """
+        Test the construction of a general elipsoid as per
+        Appendix A.1
+        """
+        px = 0
+        py = 1
+        pz = 0
+        rx = 3
+        ry = 1
+        rz = 2.5
+        E = sum([
+            (-2 * px / rx ** 2) * cyclide_ops['Tx'],
+            (-2 * py / ry ** 2) * cyclide_ops['Ty'],
+            (-2 * pz / rz ** 2) * cyclide_ops['Tz'],
+            (1 / rx ** 2) * cyclide_ops['Tx2'],
+            (1 / ry ** 2) * cyclide_ops['Ty2'],
+            (1 / rz ** 2) * cyclide_ops['Tz2'],
+            (px**2 / rx ** 2 + py**2 / ry ** 2 + pz**2 / rz ** 2 - 1) * cyclide_ops['T1']
+        ])
+        # The cyclides are an IPNS
+        assert E|eo == 0*e1
+
+    def test_line(self):
+        rng = np.random.RandomState()
+        # Make a dcga line
+        pnt_vec_a = rng.randn(3)
+        pnt_vec_b = rng.randn(3)
+        Lcga1 = IC1*(up_cga1(pnt_vec_a) ^ up_cga1(pnt_vec_b) ^ einf1)
+        Lcga2 = IC2*(up_cga2(pnt_vec_a) ^ up_cga2(pnt_vec_b) ^ einf2)
+        Ldcga = Lcga1 ^ Lcga2
+
+        # Assert that it is an IPNS
+        assert Ldcga | up(pnt_vec_a) == 0*eo
+        assert Ldcga | up(pnt_vec_b) == 0 * eo
+        assert Ldcga | up(0.5*pnt_vec_a + 0.5*pnt_vec_b) == 0 * eo
+
+    def test_translation(self):
+        rng = np.random.RandomState()
+        # Make a dcga line
+        pnt_vec = rng.randn(3)
+        direction_vec = rng.randn(3)
+        Lcga1 = IC1 * (up_cga1(pnt_vec) ^ up_cga1(pnt_vec + direction_vec) ^ einf1)
+        Lcga2 = IC2 * (up_cga2(pnt_vec) ^ up_cga2(pnt_vec + direction_vec) ^ einf2)
+        Ldcga = Lcga1 ^ Lcga2
+
+        # Make a dcga translation rotor in direction of the line
+        Tc1 = 1 - (direction_vec[0] * e1 + direction_vec[1] * e2 + direction_vec[2] * e3) * einf1
+        Tc2 = 1 - (direction_vec[0] * e6 + direction_vec[1] * e7 + direction_vec[2] * e8) * einf2
+        Tdcga = (Tc1*Tc2).normal()
+
+        # Assert the rotor is normalised
+        assert Tdcga*~Tdcga == layout.scalar
+
+        # Apply the rotor to the line
+        np.testing.assert_allclose((Tdcga*Ldcga*~Tdcga).value, Ldcga.value, rtol=1E-4, atol=1E-6)
+
+        # Apply the rotor to a point on the line
+        np.testing.assert_allclose(((Tdcga * up(pnt_vec) * ~Tdcga)|Ldcga).value, 0, rtol=1E-4, atol=1E-6)
+
+        # Construct and ellipsoid at the origin
+        px = 0
+        py = 0
+        pz = 0
+        rx = 3
+        ry = 1
+        rz = 2.5
+        E = sum([
+            (-2 * px / rx ** 2) * cyclide_ops['Tx'],
+            (-2 * py / ry ** 2) * cyclide_ops['Ty'],
+            (-2 * pz / rz ** 2) * cyclide_ops['Tz'],
+            (1 / rx ** 2) * cyclide_ops['Tx2'],
+            (1 / ry ** 2) * cyclide_ops['Ty2'],
+            (1 / rz ** 2) * cyclide_ops['Tz2'],
+            (px ** 2 / rx ** 2 + py ** 2 / ry ** 2 + pz ** 2 / rz ** 2 - 1) * cyclide_ops['T1']
+        ])
+        # Before moving the elipsoid surface is not touching the origin
+        assert E|eo != 0*e1
+
+        # Make a dcga translation rotor to move the ellipsoid
+        Tc1 = 1 - 0.5 * rx * e1 * einf1
+        Tc2 = 1 - 0.5 * rx * e6 * einf2
+        Tdcga = (Tc1 * Tc2).normal()
+        assert (Tdcga*E*~Tdcga) | eo == 0 * e1
+
+        # Make a dcga translation rotor to move the ellipsoid
+        Tc1 = 1 - 0.5 * ry * e2 * einf1
+        Tc2 = 1 - 0.5 * ry * e7 * einf2
+        Tdcga = (Tc1 * Tc2).normal()
+        assert (Tdcga * E * ~Tdcga) | eo == 0 * e1
+
+        # Make a dcga translation rotor to move the ellipsoid
+        Tc1 = 1 - 0.5 * rz * e3 * einf1
+        Tc2 = 1 - 0.5 * rz * e8 * einf2
+        Tdcga = (Tc1 * Tc2).normal()
+        assert (Tdcga * E * ~Tdcga) | eo == 0 * e1
+
+    def test_rotation(self):
+        theta = np.pi/2
+        RC1 = np.e ** (-0.5*theta*e12)
+        RC2 = np.e ** (-0.5*theta*e67)
+        Rdcga = (RC1 * RC2).normal()
+        assert Rdcga * ~Rdcga == 1.0 + 0 * eo
+
+        # Construct a line
+        pnt_vec = np.array([1, 0, 0])
+        direction_vec = np.array([0, 0, 1])
+        Lcga1 = IC1 * (up_cga1(pnt_vec) ^ up_cga1(pnt_vec + direction_vec) ^ einf1)
+        Lcga2 = IC2 * (up_cga2(pnt_vec) ^ up_cga2(pnt_vec + direction_vec) ^ einf2)
+        Ldcga = Lcga1 ^ Lcga2
+
+        # Construct a second line
+        pnt_vec_rotated = np.array([0, 1, 0])
+        Lcga1_rotated = IC1 * (up_cga1(pnt_vec_rotated) ^ up_cga1(pnt_vec_rotated + direction_vec) ^ einf1)
+        Lcga2_rotated = IC2 * (up_cga2(pnt_vec_rotated) ^ up_cga2(pnt_vec_rotated + direction_vec) ^ einf2)
+        Ldcga_rotated = Lcga1_rotated ^ Lcga2_rotated
+
+        # Assert the rotor rotates it
+        assert (Rdcga * Ldcga * ~Rdcga)|up(pnt_vec_rotated) == 0*e1
+        np.testing.assert_allclose((Rdcga * Ldcga * ~Rdcga).value, Ldcga_rotated.value, rtol=1E-4, atol=1E-6)
+
+        # Construct and ellipsoid
+        px = 0
+        py = 2.5
+        pz = 0
+        rx = 3
+        ry = 1
+        rz = 2.5
+        E = sum([
+            (-2 * px / rx ** 2) * cyclide_ops['Tx'],
+            (-2 * py / ry ** 2) * cyclide_ops['Ty'],
+            (-2 * pz / rz ** 2) * cyclide_ops['Tz'],
+            (1 / rx ** 2) * cyclide_ops['Tx2'],
+            (1 / ry ** 2) * cyclide_ops['Ty2'],
+            (1 / rz ** 2) * cyclide_ops['Tz2'],
+            (px ** 2 / rx ** 2 + py ** 2 / ry ** 2 + pz ** 2 / rz ** 2 - 1) * cyclide_ops['T1']
+        ])
+        assert E | eo != 0 * eo
+
+        # Make a dcga translation rotor to move the ellipsoid
+        Tc1 = 1 - 0.5 * py * e2 * einf1
+        Tc2 = 1 - 0.5 * py * e7 * einf2
+        Tdcga = (Tc1 * Tc2).normal()
+
+        # Construct a rotation rotor
+        theta = np.pi / 2
+        RC1 = np.e ** (-0.5 * theta * e23)
+        RC2 = np.e ** (-0.5 * theta * e78)
+        Rdcga = (RC1 * RC2).normal()
+
+        Comborotor = (Tdcga*Rdcga*~Tdcga).normal()
+
+        Erot = Comborotor*E*~Comborotor
+
+        assert Erot|eo == 0*eo
+
+    def test_bivector_orthogonality(self):
+        """
+        Rotors in each algebra should be orthogonal
+        """
+        theta = np.pi / 2
+        RC1 = np.e ** (-0.5 * theta * e12)
+        RC2 = np.e ** (-0.5 * theta * e67)
+        Rdcga = (RC1 * RC2).normal()
+        Rexp = np.e**(-0.5*theta*(e12 + e67))
+        np.testing.assert_allclose(Rexp.value, Rdcga.value,
+                                   rtol=1E-4, atol=1E-6)
+
+        mag = 5.0
+        Tc1 = 1 - 0.5 * mag * e2 * einf1
+        Tc2 = 1 - 0.5 * mag * e7 * einf2
+        Tdcga = (Tc1 * Tc2).normal()
+        Texp = np.e ** (-0.5 * mag * (e2 * einf1 + e7 * einf2))
+        np.testing.assert_allclose(Texp.value, Tdcga.value,
+                                   rtol=1E-4, atol=1E-6)
+
+
+if __name__ == '__main__':
+    unittest.main()

docs/predefined-algebras.rst

@@ -27,6 +27,9 @@ The easiest way to get started with ``clifford`` is to use one of several predef
 * .. module:: clifford.dpga
 
   ``dpga``: Double PGA also referred to as the Mother Algebra, ``Cl(4, 4)``.
+* .. module:: clifford.dg3c
+
+  ``dg3c``: Double Conformal Geometric Algebra, effectively two g3c algebras glued together ``Cl(8, 2)``.