{
 "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, transformations\n",
    "\n",
    "class OurCustomLayout(ConformalLayout):\n",
    "    def __init__(self, ndims):\n",
    "        self.ndims = ndims\n",
    "        \n",
    "        euclidean_vectors = [str(i + 1) for i in range(ndims)]\n",
    "        conformal_vectors = ['m2', 'm1']\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(euclidean_vectors + ['np1'] + conformal_vectors)\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(euclidean_vectors + conformal_vectors)\n",
    "        )\n",
    "        \n",
    "        # this lets us convert between the two layouts\n",
    "        self.to_conformal = transformations.between_basis_vectors(self, self.conformal_base)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "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": [
    "\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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "\n",
    "Finally, we'll define a plotting function, which plots the problem and solution circles in suitable colors via `pyganja`.\n",
    "Note that this all works because of our definition of the `to_conformal` `BasisVectorMap`."
   ]
  },
  {
   "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 <cite data-cite=\"lasenby-covariant-approach\">A Covariant Approach to Geometry using Geometric Algebra</cite>.\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
}