diff --git a/docs/index.rst b/docs/index.rst index 2086d5f7..49d38dc6 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -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 diff --git a/docs/requirements.txt b/docs/requirements.txt index a6052a50..946d1caf 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,5 +1,6 @@ # extra requiremenets needed by the docs themselves # numpydoc==0.4 +pandas ipykernel nbsphinx ipywidgets diff --git a/docs/tutorials/apollonius-cga-augmented.ipynb b/docs/tutorials/apollonius-cga-augmented.ipynb new file mode 100644 index 00000000..d14093e9 --- /dev/null +++ b/docs/tutorials/apollonius-cga-augmented.ipynb @@ -0,0 +1,393 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Working with custom algebras\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). It also shows \n", + "\n", + "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$.\n", + "This is an extension of a standard conformal algebra, with an extra $e_{n+1}$ basis vector.\n", + "\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": [ + "from clifford import ConformalLayout, BasisVectorIds, MultiVector\n", + "\n", + "class OurCustomLayout(ConformalLayout):\n", + " def __init__(self, ndims):\n", + " self.ndims = ndims\n", + "\n", + " # Construct our custom algebra. Note that ConformalLayout requires the e- and e+ basis vectors to be last.\n", + " ConformalLayout.__init__(\n", + " self,\n", + " [1]*ndims + [-1] + [1, -1],\n", + " ids=BasisVectorIds([str(i + 1) for i in range(ndims)] + ['np1', 'm2', 'm1'])\n", + " )\n", + " self.enp1 = self.basis_vectors_lst[ndims]\n", + "\n", + " # Construct a base algebra without the extra `enp1`, which would not be understood by pyganja.\n", + " self.conformal_base = ConformalLayout(\n", + " [1]*ndims + [1, -1],\n", + " ids=BasisVectorIds([str(i + 1) for i in range(ndims)] + ['m2', 'm1'])\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def ups(self, s):\n", + " return s + self.enp1*abs(s)\n", + "\n", + "OurCustomLayout.ups = ups; del ups\n", + "\n", + "def downs(self, mv):\n", + " if (mv | self.enp1)[()] > 0:\n", + " mv = -mv\n", + " return mv\n", + "\n", + "OurCustomLayout.downs = downs; del downs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def dual_round(at, r):\n", + " l = at.layout\n", + " return l.up(at) - 0.5*l.einf*r*r" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization of custom algebras\n", + "\n", + "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.\n", + "\n", + "
\n", + "\n", + "Warning\n", + "\n", + "This function uses the private attribute `Layout._basis_blade_order`, which may change or be removed in future.\n", + "Future releases of `clifford` may contain tools to make switching between algebras easier.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def to_conformal(self, mv: MultiVector) -> MultiVector:\n", + " \"\"\" Convert a mv with `mv.layout == self` to one with `mv_new.layout == self.conformal_base`\n", + " \n", + " This is done by discarding coefficients of blades containing enp1\n", + " \"\"\"\n", + " bits_base = self.conformal_base._basis_blade_order.index_to_bitmap\n", + " # insert the np1 bit into each basis blade bitmask, by...\n", + " bits = (\n", + " # ... leaving e1...en in place ...\n", + " bits_base & ~(~0 << self.ndims) |\n", + " # and moving em1 and em2 one bit to the left, leaving a 0 in the np1 bit\n", + " (bits_base & (0b11 << self.ndims)) << 1\n", + " )\n", + " # convert the new bitmaps into indices into the augmented value\n", + " inds = self._basis_blade_order.bitmap_to_index[bits]\n", + "\n", + " # and pick only those indicess\n", + " return self.conformal_base.MultiVector(mv.value[inds])\n", + "\n", + "OurCustomLayout.to_conformal = to_conformal; del to_conformal" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we'll define a plotting function, which plots the problem and solution circles in suitable colors via `pyganja`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import itertools\n", + "from pyganja import GanjaScene, draw\n", + "\n", + "def plot_rounds(in_rounds, out_rounds, scale=1):\n", + " colors = itertools.cycle([\n", + " (255, 0, 0),\n", + " (0, 255, 0),\n", + " (0, 0, 255),\n", + " (0, 255, 255),\n", + " ])\n", + " # note: .dual() neede here because we're passing in dual rounds, but ganja expects direct rounds\n", + " s = GanjaScene()\n", + " for r, color in zip(in_rounds, colors):\n", + " s.add_object(r.layout.to_conformal(r).dual(), color=color)\n", + " for r in out_rounds:\n", + " s.add_object(r.layout.to_conformal(r).dual(), color=(64, 64, 64))\n", + " draw(s, sig=r.layout.conformal_base.sig, scale=scale)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Apollonius' problem in $\\mathbb{R}^2$ 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 but somewhat slow trick for showing tables\n", + "pd.DataFrame(l2.metric, index=l2.basis_names, columns=l2.basis_names)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can build some dual circles:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# add minus signs before `dual_round` to flip circle directions\n", + "c1 = dual_round(-e1-e2, 1)\n", + "c2 = dual_round(e1-e2, 0.75)\n", + "c3 = dual_round(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": [ + "We hypothesize that this object is of the form `l2.ups(c4) ^ l2.ups(c5)`.\n", + "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 [A Covariant Approach to Geometry using Geometric Algebra](https://pdfs.semanticscholar.org/baba/976fd7f6577eeaa1d3ef488c1db13ec24652.pdf).\n", + "Here, we normalize with $e_{n+1}$ instead of the usual $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 finally, plot our circles:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_rounds([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_round(-1.5*e1, 0.5)\n", + "c2 = dual_round(e1*0, 0.5)\n", + "c3 = dual_round(1.5*e1, 0.5)\n", + "c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())\n", + "\n", + "plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "c1 = dual_round(-3*e1, 1.5)\n", + "c2 = dual_round(-2*e1, 1)\n", + "c3 = -dual_round(2*e1, 1)\n", + "c4u, c5u = pp_ends((l2.ups(c1) ^ l2.ups(c2) ^ l2.ups(c3)).dual())\n", + "\n", + "plot_rounds([c1, c2, c3], [l2.downs(c4u), l2.downs(c5u)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Apollonius' problem in $\\mathbb{R}^3$ 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": "markdown", + "metadata": {}, + "source": [ + "Again, we can check the metric:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pd.DataFrame(l3.metric, index=l3.basis_names, columns=l3.basis_names)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And apply the solution to some spheres, noting that we now need 4 in order to constrain our solution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "c1 = dual_round(e1+e2+e3, 1)\n", + "c2 = dual_round(-e1+e2-e3, 0.25)\n", + "c3 = dual_round(e1-e2-e3, 0.5)\n", + "c4 = dual_round(-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_rounds([c1, c2, c3, c4], [l3.downs(c6u), l3.downs(c5u)], scale=0.25)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the figure above can be rotated!" + ] + } + ], + "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 +}