This page was generated from docs/tutorials/apollonius-cga-augmented.ipynb. Interactive online version: .

Working with custom algebras

This notebook explores the algebra defined in The Lie Model for Euclidean Geometry (Hongbo Li), and its application to solving Apollonius’ Problem. It also shows

The algebra is constructed with basis elements \(e_{-2}, e_{-1}, e_1, \cdots, e_n, e_{n+1}\), where \(e_{-2}^2 = -e_{-1}^2 = -e_{n+1}^2 = 1\). This is an extension of a standard conformal algebra, with an extra \(e_{n+1}\) basis vector.

Note that we permuted the order in the source code below to make ConformalLayout happy.

from clifford import ConformalLayout, BasisVectorIds, MultiVector, transformations

class OurCustomLayout(ConformalLayout):
    def __init__(self, ndims):
        self.ndims = ndims

        euclidean_vectors = [str(i + 1) for i in range(ndims)]
        conformal_vectors = ['m2', 'm1']

        # Construct our custom algebra. Note that ConformalLayout requires the e- and e+ basis vectors to be last.
        ConformalLayout.__init__(
            self,
            [1]*ndims + [-1] + [1, -1],
            ids=BasisVectorIds(euclidean_vectors + ['np1'] + conformal_vectors)
        )
        self.enp1 = self.basis_vectors_lst[ndims]

        # Construct a base algebra without the extra `enp1`, which would not be understood by pyganja.
        self.conformal_base = ConformalLayout(
            [1]*ndims + [1, -1],
            ids=BasisVectorIds(euclidean_vectors + conformal_vectors)
        )

        # this lets us convert between the two layouts
        self.to_conformal = transformations.between_basis_vectors(self, self.conformal_base)

The code above also defines a stardard conformal \(\mathbb{R}^{N+1,1}\) layout without this new basis vector. This is primarily to support rendering with pyganja, which doesn’t support the presence of this extra vector. BasisVectorMap defaults to preserving vectors by name between one algebra and another, while throwing away blades containing vectors missing from the destination algebra.

We define an ups function which maps conformal dual-spheres into this algebra, as \(s^\prime = s + \left|s\right|e_{n+1}\), and a downs that applies the correct sign. The s suffix here is chosen to mean sphere.

def ups(self, s):
    return s + self.enp1*abs(s)

OurCustomLayout.ups = ups; del ups

def downs(self, mv):
    if (mv | self.enp1)[()] > 0:
        mv = -mv
    return mv

OurCustomLayout.downs = downs; del downs

Before we start looking at specified dimensions of euclidean space, we build a helper to construct conformal dual circles and spheres, with the word round being a general term intended to cover both circles and spheres.

def dual_round(at, r):
    l = at.layout
    return l.up(at) - 0.5*l.einf*r*r

In order to render with pyganja, we’ll need a helper to convert from our custom \(\mathbb{R}^{N+1,2}\) layout into a standard conformal \(\mathbb{R}^{N+1,1}\) layout. clifford maps indices in .value to basis blades via layout._basis_blade_order.index_to_bitmap, which we can use to convert the indices in one layout to the indices in another.

Visualization

Finally, we’ll define a plotting function, which plots the problem and solution circles in suitable colors via pyganja. Note that this all works because of our definition of the to_conformal BasisVectorMap.

import itertools
from pyganja import GanjaScene, draw

def plot_rounds(in_rounds, out_rounds, scale=1):
    colors = itertools.cycle([
        (255, 0, 0),
        (0, 255, 0),
        (0, 0, 255),
        (0, 255, 255),
    ])
    # note: .dual() neede here because we're passing in dual rounds, but ganja expects direct rounds
    s = GanjaScene()
    for r, color in zip(in_rounds, colors):
        s.add_object(r.layout.to_conformal(r).dual(), color=color)
    for r in out_rounds:
        s.add_object(r.layout.to_conformal(r).dual(), color=(64, 64, 64))
    draw(s, sig=r.layout.conformal_base.sig, scale=scale)

/home/docs/checkouts/readthedocs.org/user_builds/clifford/envs/latest/lib/python3.8/site-packages/pyganja/__init__.py:2: UserWarning: Failed to import cef_gui, cef functions will be unavailable
  from .script_api import *

Apollonius’ problem in \(\mathbb{R}^2\) with circles

l2 = OurCustomLayout(ndims=2)
e1, e2 = l2.basis_vectors_lst[:2]

This gives us the Layout l2 with the desired metric,

import pandas as pd  # convenient but somewhat slow trick for showing tables
pd.DataFrame(l2.metric, index=l2.basis_names, columns=l2.basis_names)

	e1	e2	enp1	em2	em1
e1	1.0	0.0	0.0	0.0	0.0
e2	0.0	1.0	0.0	0.0	0.0
enp1	0.0	0.0	-1.0	0.0	0.0
em2	0.0	0.0	0.0	1.0	0.0
em1	0.0	0.0	0.0	0.0	-1.0

Now we can build some dual circles:

# add minus signs before `dual_round` to flip circle directions
c1 = dual_round(-e1-e2, 1)
c2 = dual_round(e1-e2, 0.75)
c3 = dual_round(e2, 0.5)

Compute the space orthogonal to all of them, which is an object of grade 2:

pp = (l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual()
pp.grades()

{2}

We hypothesize that this object is of the form l2.ups(c4) ^ l2.ups(c5). Taking a step not mentioned in the original paper, we decide to treat this as a regular conformal point pair, which allows us to project out the two factors with the approach taken in [LLW04]. Here, we normalize with \(e_{n+1}\) instead of the usual \(n_\infty\):

def pp_ends(pp):
    P = (1 + pp.normal()) / 2
    return P * (pp | pp.layout.enp1), ~P * (pp | pp.layout.enp1)

c4u, c5u = pp_ends(pp)

And finally, plot our circles:

plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)], scale=0.75)

This works for colinear circles too:

c1 = dual_round(-1.5*e1, 0.5)
c2 = dual_round(e1*0, 0.5)
c3 = dual_round(1.5*e1, 0.5)
c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())

plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])

c1 = dual_round(-3*e1, 1.5)
c2 = dual_round(-2*e1, 1)
c3 = -dual_round(2*e1, 1)
c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())

plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])

Apollonius’ problem in \(\mathbb{R}^3\) with spheres

l3 = OurCustomLayout(ndims=3)
e1, e2, e3 = l3.basis_vectors_lst[:3]

Again, we can check the metric:

pd.DataFrame(l3.metric, index=l3.basis_names, columns=l3.basis_names)

	e1	e2	e3	enp1	em2	em1
e1	1.0	0.0	0.0	0.0	0.0	0.0
e2	0.0	1.0	0.0	0.0	0.0	0.0
e3	0.0	0.0	1.0	0.0	0.0	0.0
enp1	0.0	0.0	0.0	-1.0	0.0	0.0
em2	0.0	0.0	0.0	0.0	1.0	0.0
em1	0.0	0.0	0.0	0.0	0.0	-1.0

And apply the solution to some spheres, noting that we now need 4 in order to constrain our solution

c1 = dual_round(e1+e2+e3, 1)
c2 = dual_round(-e1+e2-e3, 0.25)
c3 = dual_round(e1-e2-e3, 0.5)
c4 = dual_round(-e1-e2+e3, 1)
c5u, c6u = pp_ends((l3.ups(c1) ^ l3.ups(c2) ^ l3.ups(c3) ^ l3.ups(c4)).dual())

plot_rounds([c1, c2, c3, c4], [l3.downs(c6u), l3.downs(c5u)], scale=0.25)

Note that the figure above can be rotated!