Skip to content

Commit

Permalink
Add a tutorial notebook showing how to hack together custom algebras
Browse files Browse the repository at this point in the history
  • Loading branch information
@eric-wieser
eric-wieser committed Apr 1, 2020
1 parent 3c70a07 commit 97a04bd
Show file tree
Hide file tree
Showing 3 changed files with 378 additions and 0 deletions.
1 change: 1 addition & 0 deletions docs/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ Scalars, vectors, and higher-grade entities can be mixed freely and consistently
tutorials/InterfacingOtherMathSystems
tutorials/PerformanceCliffordTutorial
tutorials/cga/index
tutorials/apollonius-cga-augmented

.. toctree::
:maxdepth: 1
Expand Down
1 change: 1 addition & 0 deletions docs/requirements.txt
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
# extra requiremenets needed by the docs themselves
# numpydoc==0.4
pandas
ipykernel
nbsphinx
ipywidgets
Expand Down
376 changes: 376 additions & 0 deletions docs/tutorials/apollonius-cga-augmented.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,376 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Apollonius' problem via augmented CGA\n",
"\n",
"This notebook explores the algebra defined in [The Lie Model for Euclidean Geometry (Hongbo Li)](https://link.springer.com/chapter/10.1007/10722492_7), and its application to solving [Apollonius' Problem](https://mathworld.wolfram.com/ApolloniusProblem.html)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from clifford import Layout, BasisVectorIds, ConformalLayout, MultiVector\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The algebra is constructed with basis elements $e_1, \\cdots, e_n, e_{n+1}, e_{-2}, e_{-1}$, where $e_{-2}^2 = -e_{-1}^2 = -e_{n+1}^2 = 1$.\n",
"Note that we permuted the order in the source code below to make `ConformalLayout` happy."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"class OurCustomLayout(ConformalLayout):\n",
" def __init__(self, ndims):\n",
" # construct our custom algebra\n",
" ConformalLayout.__init__(self,\n",
" [1]*ndims + [-1] + [1, -1],\n",
" ids=BasisVectorIds([str(i + 1) for i in range(ndims)] + ['np1', 'm2', 'm1']))\n",
"\n",
" # and import a matching regular conformal one, which ganja.js supports\n",
" if ndims == 3:\n",
" from clifford.g3c import layout as l_base\n",
" elif ndims == 2:\n",
" from clifford.g2c import layout as l_base\n",
" else:\n",
" raise NotImplementedError\n",
" self.conformal_base = l_base\n",
" self.enp1 = self.basis_vectors_lst[ndims]\n",
" self.ndims = ndims"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## In 2d space with circles"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"l2 = OurCustomLayout(ndims=2)\n",
"e1, e2 = l2.basis_vectors_lst[:2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This gives us the `Layout` `l2` with the desired metric,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd # convenient trick for showing tables\n",
"pd.DataFrame(l2.metric, index=l2.basis_names, columns=l2.basis_names)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We define an `ups` function which maps conformal dual-spheres into this algebra, as $s + \\left|s\\right|e_{n+1}$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# s for spherical\n",
"def ups(self, s):\n",
" return s + self.enp1*abs(s)\n",
"OurCustomLayout.ups = ups\n",
"del ups\n",
"\n",
"def downs(self, mv):\n",
" if (mv | self.enp1)[()] > 0:\n",
" mv = -mv\n",
" return mv\n",
"OurCustomLayout.downs = downs\n",
"del downs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And a helper to construct conformal dual spheres"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def dual_circ(at, r):\n",
" l = at.layout\n",
" return l.up(at) - 0.5*l.einf*r*r"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some ugly code to convert from $\\mathbb{R}^{N+1,2}$ to $\\mathbb{R}^{N+1,1}$ in a way that pyganja will understand, by throwing away unwanted coefficients."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from typing import List\n",
"\n",
"class BitTransformer:\n",
" \"\"\" Convert bitmap blade representations from one algebra to another \"\"\"\n",
" def __init__(self, from_list: List[int]) -> None:\n",
" self.from_list = from_list\n",
"\n",
" def __call__(self, bits: int) -> int:\n",
" ret = 0\n",
" for t, f in enumerate(self.from_list):\n",
" if f is None:\n",
" continue\n",
" ret |= ((bits >> f) & 1) << t\n",
" bits &= ~(1 << f)\n",
" if bits:\n",
" raise ValueError(\"Unsupported bits set: 0b{:b}\".format(bits))\n",
" return ret\n",
"\n",
" @property\n",
" def inverse(self) -> 'BitTransformer':\n",
" to_list = [None]*(max(self.from_list, default=-1) + 1)\n",
" for t, f in enumerate(self.from_list):\n",
" to_list[f] = t\n",
" return BitTransformer(to_list)\n",
"\n",
"\n",
"def to_conformal(self, mv: MultiVector) -> MultiVector:\n",
" \"\"\" Convert a multivector in the custom algebra to the standard conformal one\"\"\"\n",
" \n",
" # remove the np1 dimension\n",
" bitmap_to_conformal = BitTransformer([*range(self.ndims), self.ndims + 1, self.ndims + 2]).inverse\n",
" \n",
" mv_base_vals = np.zeros(self.conformal_base.gaDims, dtype=mv.value.dtype)\n",
" \n",
" for i_base, bits_base in enumerate(self.conformal_base._basis_blade_order.index_to_bitmap):\n",
" bits = bitmap_to_conformal(bits_base)\n",
" i = self._basis_blade_order.bitmap_to_index[bits]\n",
" mv_base_vals[i_base] = mv.value[i]\n",
" return self.conformal_base.MultiVector(mv_base_vals)\n",
"\n",
"OurCustomLayout.to_conformal = to_conformal\n",
"del to_conformal"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now build some dual circles in CGA:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# adjust minus signs to get all solutions\n",
"c1 = dual_circ(-e1-e2, 1)\n",
"c2 = dual_circ(e1-e2, 0.75)\n",
"c3 = dual_circ(e2, 0.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the space orthogonal to all of them, which is an object of grade 2"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pp = (l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual()\n",
"pp.grades()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Project out the point-pair ends, using the approach taken in [A Covariant Approach to Geometry using Geometric Algebra](https://pdfs.semanticscholar.org/baba/976fd7f6577eeaa1d3ef488c1db13ec24652.pdf), noting that we now normalize with $e_{n+1}$ instead of $n_\\infty$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def pp_ends(pp):\n",
" P = (1 + pp.normal()) / 2\n",
" return P * (pp | pp.layout.enp1), ~P * (pp | pp.layout.enp1)\n",
"\n",
"c4u, c5u = pp_ends(pp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And we can now plot all our circles, noting that ganja requires us to work in the direct form:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import itertools\n",
"from pyganja import GanjaScene, draw\n",
"\n",
"def plot_circles(in_circles, out_circles, scale=1):\n",
" colors = itertools.cycle([\n",
" (255, 0, 0),\n",
" (0, 255, 0),\n",
" (0, 0, 255),\n",
" (0, 255, 255),\n",
" ])\n",
" s = GanjaScene()\n",
" for c, color in zip(in_circles, colors):\n",
" s.add_object(c.layout.to_conformal(c).dual(), color=color)\n",
" for c in out_circles:\n",
" s.add_object(c.layout.to_conformal(c).dual(), color=(64, 64, 64))\n",
" draw(s, sig=c.layout.conformal_base.sig, scale=scale)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plot_circles([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)], scale=0.75)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This works for colinear circles too:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"c1 = dual_circ(-1.5*e1, 0.5)\n",
"c2 = dual_circ(e1*0, 0.5)\n",
"c3 = dual_circ(1.5*e1, 0.5)\n",
"c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())\n",
"\n",
"plot_circles([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"c1 = dual_circ(-3*e1, 1.5)\n",
"c2 = dual_circ(-2*e1, 1)\n",
"c3 = -dual_circ(2*e1, 1)\n",
"c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())\n",
"\n",
"plot_circles([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## In 3d space with spheres"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"l3 = OurCustomLayout(ndims=3)\n",
"e1, e2, e3 = l3.basis_vectors_lst[:3]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"c1 = dual_circ(e1+e2+e3, 1)\n",
"c2 = dual_circ(-e1+e2-e3, 0.25)\n",
"c3 = dual_circ(e1-e2-e3, 0.5)\n",
"c4 = dual_circ(-e1-e2+e3, 1)\n",
"c5u, c6u = pp_ends((l3.ups(c1) ^ l3.ups(c2) ^ l3.ups(c3) ^ l3.ups(c4)).dual())\n",
"\n",
"plot_circles([c1, c2, c3, c4], [l3.downs(c6u), l3.downs(c5u)], scale=0.25)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.0"
}
},
"nbformat": 4,
"nbformat_minor": 1
}

0 comments on commit 97a04bd

Please sign in to comment.