PolyhedralGeometry.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction to Combinatorial and Polyhedral Geometry in Sage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Important pages available online:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- General Page: https://doc.sagemath.org/html/en/reference/discrete_geometry/index.html\n",
    "- Thematic Tutorials: https://doc.sagemath.org/html/en/thematic_tutorials/geometry.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wiki Pages on Development:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Development Wiki page: https://wiki.sagemath.org/OptiPolyGeom    (Older but still used)\n",
    "- Trac Wiki page: https://trac.sagemath.org/wiki/SagePolyhedralGeometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Basics: $V$-representation\n",
    "\n",
    "First, let’s define a polyhedron object as the convex hull of a set of points and some rays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P1 = Polyhedron(vertices=[[2, 0], [0, 2]], rays=[[1, 1]], base_ring=ZZ); P1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The string representation already gives a lot of information:\n",
    "\n",
    "* the **dimension** of the polyhedron (the smallest affine space containing it)\n",
    "* the **dimension of the space** in which it is defined\n",
    "* the **base ring** ($\\mathbb{Z}$) over which the polyhedron lives (this specifies the parent class)\n",
    "* the **number of vertices**\n",
    "* the **number of rays**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The usual tab completion shows us everything we can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### We can also add a lineality space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P2 = Polyhedron(vertices=[[0.5, 0, 0], [0, 0.5, 0]], rays=[[1, 1, 0]], lines=[[0, 0, 1]]);P2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "P2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The objects have a specific parent (base ring) and a specific backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P1.parent(), P1.backend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P2.parent(), P2.backend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### These two (base ring & backend) determine which methods and algorithms are available/used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## $H$-representation\n",
    "\n",
    "If a polyhedron object was constructed via a $V$-representation, Sage can provide the $H$-representation of the object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for h in P1.Hrepresentation():\n",
    "    print(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newP1 = Polyhedron(ieqs=[[2,1,-1],[-2,1,1],[2,-1,1]]); newP1 == P1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Polyhedron?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### A large library of precomputed examples:\n",
    "\n",
    "http://doc.sagemath.org/html/en/reference/discrete_geometry/sage/geometry/polyhedron/library.html\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oht = polytopes.one_hundred_twenty_cell(exact=True,backend='normaliz');oht.show(frame=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oht.f_vector(), oht.base_ring()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Representation of objects\n",
    "\n",
    "Sage has classes implemented for $H$- and $V$-representations and faces of polyhedron."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H$-representation\n",
    "\n",
    "You can store the $H$-representation in a variable and use the inequalities and equalities as objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HRep = P2.Hrepresentation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H1 = HRep[0]; H1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H1.A(), H1.b()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H1.contains(vector([0,0,0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### $V$-representation\n",
    "\n",
    "Similarly, you can access vertices, rays and lines of the polyhedron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VRep = P2.Vrepresentation(); VRep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = VRep[0]; L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(L.neighbors())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Faces:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aface = P1.faces(1)[1]; aface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# The algebra of polytopes: classical operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polar and Dilation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cube = polytopes.cube()             # [-1,1]-cube\n",
    "Octahedron = 3/2 * Cube.polar()     # Dilation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minkowski sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cuboctahedron = Cube + Octahedron   # Minkowski sum\n",
    "Cuboctahedron.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Minkowski difference (when it makes sense)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mdiff = Cuboctahedron - Cube        # Minkowski difference\n",
    "Mdiff.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Intersection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Permu = Cube & Octahedron           # Intersection, same as   Cube.intersection(Octahedron)\n",
    "Permu.show(opacity=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question to Sage developers listening: is this a bug ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Cartesian product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "CP = Cube * Octahedron; CP          # Cartesian product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CP.f_vector()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Affine transformation through coersion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zo_cube = polytopes.cube(intervals='zero_one').change_ring(QQ) # (0,1) 3-d cube\n",
    "transformation = matrix(AA,[[1,sqrt(2),sqrt(3)],[0,1,sqrt(5)],[0,0,1]])\n",
    "skew_cube = (transformation * zo_cube) + vector([0,0,1/2])   # affinely transformed\n",
    "skew_cube.plot(color='blue',opacity=0.75) + zo_cube.plot(color='red',opacity=0.75)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(\"The cube base ring is: {}\".format(zo_cube.base_ring()))\n",
    "print(\"The skew cube base ring is: {}\\n\".format(skew_cube.base_ring()))\n",
    "print(\"The cube backend is: {}\".format(zo_cube.backend()))\n",
    "print(\"The skew cube backend is: {}\".format(skew_cube.backend()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Combinatorial equivalence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "co = polytopes.cuboctahedron()\n",
    "another_co = Cube & (2 * Cube.polar())\n",
    "yet_another_co = Cube.truncation(1/2)\n",
    "\n",
    "co.is_combinatorially_isomorphic(another_co), co.is_combinatorially_isomorphic(yet_another_co)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### There are many more operations that are possible: http://doc.sagemath.org/html/en/thematic_tutorials/geometry/polyhedra_quickref.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A cool trick to get back the input of an object:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "An_important_polytope = skew_cube & zo_cube  # This was hard to compute... let's say!\n",
    "sage_input(An_important_polytope)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## What if I want to include a polytope to my article in latex?\n",
    "### First, you may want to have a look at the inequalities, just for sanity check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Permu.Hrepresentation_str(style='<='))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Then, you can create the latex code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "latex_code = LatexExpr(Permu.Hrepresentation_str(latex=True))\n",
    "print(latex_code)\n",
    "view(latex_code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Then, maybe you would like to add a tikz picture of your polytope..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Permu.plot(opacity=0.75,frame=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latex.add_to_preamble(\"\\\\usepackage{tikz}\")     # internal pdf compiler\n",
    "tikz_code = Permu.projection().tikz([-0.3335,-0.5165,-0.7887],126.02,scale=4,opacity=0.25,facet_color='green')\n",
    "view(tikz_code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The code is easy to read:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print('\\n'.join(tikz_code.splitlines()[:20]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## See more details here: http://page.mi.fu-berlin.de/labbe/polytope.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### You may have computed the f-vector, and would like to save the object for later use or to send to a colleague..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Permu.f_vector()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Saves in the current folder\n",
    "Permu.save('somecoolpolytope.sobj')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loadedpolytope = load(\"somecoolpolytope\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The f-vector is preserved after saving:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loadedpolytope.f_vector.is_in_cache()      # Since Sage 9.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Internal backends for computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### There are several possible backends to do computations with polyhedral objects. The performance and possible base rings differ depending on the chosen backend/engine:\n",
    "\n",
    "- ppl ($\\mathbb{Z},\\mathbb{Q}$)\n",
    "- cdd ($\\mathbb{Q}$,RDF)\n",
    "- normaliz ($\\mathbb{Z},\\mathbb{Q}$, Algebraic numbers)\n",
    "- native python \"field\" ($\\mathbb{Z},\\mathbb{Q}$, Algebraic numbers)\n",
    "- polymake ($\\mathbb{Z},\\mathbb{Q}$, quadratic extensions)\n",
    "\n",
    "\n",
    "### Some engines: \n",
    "- topcom (triangulations)\n",
    "- lrs (triangulations)\n",
    "- latte_int (triangulations, volume, integer counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The features of each backend are steadily under development and constantly improving. (We need feedback to know which features of which backends are desired!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### A good example is to compute the volume of a polytope:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Permu.volume()                    # default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Permu.volume(engine='latte')      # uses LattE Integrale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uses normaliz, which requires to change the backend.\n",
    "Permu.base_extend(base_ring=QQ,backend='normaliz').volume(engine='normaliz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Permu.volume(engine='lrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Certain backends have some more features that are interfaced. $h^*$-vector is an example with normaliz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cube = polytopes.cube(intervals='zero_one',backend='normaliz')   # [0,1]-cube\n",
    "Cube.h_star_vector()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ehrhart polynomials are also available through latte or normaliz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cube.ehrhart_polynomial(engine='latte').factor()  # Computed via latte"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cube.ehrhart_polynomial(engine='normaliz').factor()  # Computed via normaliz, the backend should be normaliz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Algebraic polyhedra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In https://arxiv.org/abs/1910.05241 [Doolittle-L-Lange-Sinn-Spreer-Ziegler] we used Sage to certify that the cyclic polytope $C_4(7)$ has realizations with all its vertices on a sphere."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\Rightarrow$ This required fast, exact arithmetic computations of convex hulls over algebraic number fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thanks to the e-antic, (py)normaliz, and Sage development teams, this was made possible!\n",
    "\n",
    "### (An article is in preparation with Bruns, Delecroix, Gutsche, Köppe, L)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.1",
   "language": "sage",
   "name": "sagemath"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}