PyNormaliz_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introducing PyNormaliz \n",
    "\n",
    "(Winfried Bruns, Sebastian Gutsche, Justhin Shenk, Richard Sieg)\n",
    "\n",
    "This small tutorial introduces the main syntax of the Python interface PyNormaliz for Normaliz.\n",
    "\n",
    "To install PyNormaliz on your machine, download the PyNormaliz and run \n",
    "<pre>\n",
    "python setup.py install\n",
    "</pre>\n",
    "If you are using python3 you might need to enter python3 instead of python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example 1: A cone in dimension 2\n",
    "We want to investigate the cone $C=\\mathbb{R}_{+}(2,1)+\\mathbb{R}_{+}(1,3)\\subset\\mathbb{R}^2$.\n",
    "\n",
    "In the cone constructor of PyNormaliz, we specify the input type (`cone`) and an input matrix for this input type. \n",
    "\n",
    "To compute its Hilbert basis we use the point operator and call `HilbertBasis()` on our cone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 1], [1, 2], [1, 3], [2, 1]]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from PyNormaliz import * # or import PyNormaliz\n",
    "\n",
    "gen = [[1,3],[2,1]]\n",
    "C = Cone(cone=gen)\n",
    "C.HilbertBasis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also use inequalities as an input. We also check that the extreme rays of the new cone are really the generators of the previous one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "[[1, 1], [1, 2], [1, 3], [2, 1]]\n"
     ]
    }
   ],
   "source": [
    "ineq = [[-1,2],[3,-1]]\n",
    "C2 = Cone(inequalities=ineq)\n",
    "gen2 = C2.ExtremeRays()\n",
    "print(gen==gen2)\n",
    "HB2 = C2.HilbertBasis()\n",
    "print(HB2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general construction of objects in PyNormaliz follows the same fashion: You can generate a cone via `Cone` and then specify the usual input type(s) as parameters with the input data as `type=input data`. Once the cone is created, all its possible data is accessible using the `.` operator, like `C.HilbertBasis()`. \n",
    "\n",
    "Basically all input types and output data is available. An overview is given in the usual manual."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Example 2: A lattice polytope\n",
    "Now we investigate a lattice simplex, which is not normal/integrally closed.\n",
    "\n",
    "The input type is `polytope`. Note that we do not need a homogenizing variable but really specify the points in their actual dimension.\n",
    "\n",
    "We compute its Hilbert series and Hilbert polynomial, which are also sometimes referred to as Ehrhart series and polynomial:\n",
    "\n",
    "Hilbert series: $H(t)=\\frac{1+14t+15t^2}{(1-t)^4}$\n",
    "\n",
    "Hilbert polynomial: $p(k)=1+4k+8k^2+5k^3$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19\n",
      "False\n",
      "[[1, 14, 15], [1, 1, 1, 1], 0]\n",
      "[[1, 4, 8, 5], 1]\n"
     ]
    }
   ],
   "source": [
    "vertices = [[0,0,0],[2,0,0],[0,3,0],[0,0,5]]\n",
    "poly = Cone(polytope=vertices)\n",
    "HB = poly.HilbertBasis()\n",
    "print(len(HB))\n",
    "print(poly.IsDeg1HilbertBasis())\n",
    "print(poly.HilbertSeries()) \n",
    "print(poly.HilbertQuasiPolynomial())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output for the Hilbert series and quasipolynomial has to be read as follows: \n",
    "\n",
    "For the series, the first vector consists of the coefficients of the numerator polynomial, e.g., `[1,14,15]` is $1+14t+15t^2$. The second vector collects the exponents in the denominator. For example `[1,1,1,1]` represents $(1-t)(1-t)(1-t)(1-t)=(1-t)^4$ and `[1,2]` is $(1-t)(1-t^2)$. The last entry is the possible shift of the Hilbert series.\n",
    "\n",
    "The output of the quasipolynomial usually consists of a list of vectors which show the coefficients of the respective polynomial associated to the congruence class modulo the period of the quasipolynomial. In the above example it is a proper polynomial, so only one vector is listed. The last entry is the greatest common divisor of all denominators in these polynomials. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If you are only interested in the lattice points, you can compute them via the command `Deg1Elements()`. Note that the output now has the homogenizing variable as last coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18\n",
      "[[0, 0, 0, 1], [0, 0, 1, 1], [0, 0, 2, 1], [0, 0, 3, 1], [0, 0, 4, 1], [0, 0, 5, 1], [0, 1, 0, 1], [0, 1, 1, 1], [0, 1, 2, 1], [0, 1, 3, 1], [0, 2, 0, 1], [0, 2, 1, 1], [0, 3, 0, 1], [1, 0, 0, 1], [1, 0, 1, 1], [1, 0, 2, 1], [1, 1, 0, 1], [2, 0, 0, 1]]\n"
     ]
    }
   ],
   "source": [
    "LP = poly.Deg1Elements()\n",
    "print(len(LP))\n",
    "print(LP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "If you would like to have a list of all already computed properties and data, you can use the `print_properties()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generators:\n",
      "[[0, 0, 0, 1], [0, 0, 5, 1], [0, 3, 0, 1], [2, 0, 0, 1]]\n",
      "\n",
      "\n",
      "ExtremeRays:\n",
      "[[0, 0, 0, 1], [0, 0, 5, 1], [0, 3, 0, 1], [2, 0, 0, 1]]\n",
      "\n",
      "\n",
      "SupportHyperplanes:\n",
      "[[-15, -10, -6, 30], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]]\n",
      "\n",
      "\n",
      "HilbertBasis:\n",
      "[[0, 0, 0, 1], [0, 0, 1, 1], [0, 0, 2, 1], [0, 0, 3, 1], [0, 0, 4, 1], [0, 0, 5, 1], [0, 1, 0, 1], [0, 1, 1, 1], [0, 1, 2, 1], [0, 1, 3, 1], [0, 2, 0, 1], [0, 2, 1, 1], [0, 3, 0, 1], [1, 0, 0, 1], [1, 0, 1, 1], [1, 0, 2, 1], [1, 1, 0, 1], [2, 0, 0, 1], [1, 2, 4, 2]]\n",
      "\n",
      "\n",
      "Deg1Elements:\n",
      "[[0, 0, 0, 1], [0, 0, 1, 1], [0, 0, 2, 1], [0, 0, 3, 1], [0, 0, 4, 1], [0, 0, 5, 1], [0, 1, 0, 1], [0, 1, 1, 1], [0, 1, 2, 1], [0, 1, 3, 1], [0, 2, 0, 1], [0, 2, 1, 1], [0, 3, 0, 1], [1, 0, 0, 1], [1, 0, 1, 1], [1, 0, 2, 1], [1, 1, 0, 1], [2, 0, 0, 1]]\n",
      "\n",
      "\n",
      "Sublattice:\n",
      "[[[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]], [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]], 1]\n",
      "\n",
      "\n",
      "OriginalMonoidGenerators:\n",
      "[[0, 0, 0, 1], [2, 0, 0, 1], [0, 3, 0, 1], [0, 0, 5, 1]]\n",
      "\n",
      "\n",
      "MaximalSubspace:\n",
      "[]\n",
      "\n",
      "\n",
      "Grading:\n",
      "[0, 0, 0, 1, 1]\n",
      "\n",
      "\n",
      "TriangulationSize:\n",
      "1\n",
      "\n",
      "\n",
      "TriangulationDetSum:\n",
      "30\n",
      "\n",
      "\n",
      "GradingDenom:\n",
      "1\n",
      "\n",
      "\n",
      "UnitGroupIndex:\n",
      "1\n",
      "\n",
      "\n",
      "InternalIndex:\n",
      "30\n",
      "\n",
      "\n",
      "Multiplicity:\n",
      "[30, 1]\n",
      "\n",
      "\n",
      "Rank:\n",
      "4\n",
      "\n",
      "\n",
      "EmbeddingDim:\n",
      "4\n",
      "\n",
      "\n",
      "IsPointed:\n",
      "True\n",
      "\n",
      "\n",
      "IsDeg1ExtremeRays:\n",
      "True\n",
      "\n",
      "\n",
      "IsDeg1HilbertBasis:\n",
      "False\n",
      "\n",
      "\n",
      "IsIntegrallyClosed:\n",
      "False\n",
      "\n",
      "\n",
      "IsInhomogeneous:\n",
      "False\n",
      "\n",
      "\n",
      "HilbertSeries:\n",
      "[[1, 14, 15], [1, 1, 1, 1], 0]\n",
      "\n",
      "\n",
      "HilbertQuasiPolynomial:\n",
      "[[1, 4, 8, 5], 1]\n",
      "\n",
      "\n",
      "IsTriangulationNested:\n",
      "False\n",
      "\n",
      "\n",
      "IsTriangulationPartial:\n",
      "False\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "poly.print_properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Example 3: A rational polytope\n",
    "\n",
    "We construct a polytope with vertices $(5/2,3/2),(-2/3,-4/3),(1/4,-7/4)$.\n",
    "\n",
    "This is the first time we enter two different input types into the constructor. With the grading we declare the last coordinate to be the denominator of the vertices.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1, 1, 4, 4, 5, 2, 4, 6, 3, 2, 6, 4, 2, 3], [1, 1, 12], 0]\n",
      "[[24, 6, 47], [19, 6, 47], [16, 6, 47], [15, 6, 47], [40, 6, 47], [19, 6, 47], [0, 6, 47], [31, 6, 47], [40, 6, 47], [3, 6, 47], [16, 6, 47], [31, 6, 47], 24]\n"
     ]
    }
   ],
   "source": [
    "rat_vert = [[5,3,2],[-2,-4,3],[1,-7,4]]\n",
    "g = [0,0,1]\n",
    "rat_poly = Cone(cone=rat_vert,grading=g)\n",
    "print(rat_poly.HilbertSeries())\n",
    "print(rat_poly.HilbertQuasiPolynomial())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We compute the integer hull of this polytope which is again a cone object and thus we can access its vertices or support hyperplanes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0, -1, 1], [2, 1, 1]]\n",
      "[[1, -2, 0], [1, 0, 0]]\n",
      "[[1, -1, -1]]\n"
     ]
    }
   ],
   "source": [
    "int_hull = rat_poly.IntegerHull()\n",
    "print(int_hull.VerticesOfPolyhedron())\n",
    "print(int_hull.SupportHyperplanes())\n",
    "print(int_hull.Equations())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Example 4: A polyhedron\n",
    "\n",
    "We define a polyhedron by inequalties:\n",
    "$2x_2 \\geq -1$,\n",
    "$2x_2 \\leq 3$,\n",
    "$-2x_1+2x_2 \\leq 3$\n",
    "\n",
    "The input type is `inhom_inequalities` where a row $(a_1,\\ldots,a_n,b)$ in the input matrix defines the inequality\n",
    "\n",
    "\\begin{equation*}a_1x_1 +\\dots+a_nx_n+b\\geq0.\\end{equation*}\n",
    "\n",
    "In this case, `HilbertBasis()` returns the Hilbert basis of the recession cone associated to the polyhedron.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1, 0, 0]]\n",
      "[[-1, 0, 1], [0, 1, 1]]\n",
      "[[-4, -1, 2], [0, 3, 2]]\n"
     ]
    }
   ],
   "source": [
    "ineq2 = [[0,2,1],[0,-2,3],[2,-2,3]]\n",
    "polyhedron = Cone(inhom_inequalities=ineq2)\n",
    "HB_rec = polyhedron.HilbertBasis()\n",
    "print(HB_rec)\n",
    "\n",
    "module_gen = polyhedron.ModuleGenerators()\n",
    "print(module_gen)\n",
    "\n",
    "vert_poly = polyhedron.VerticesOfPolyhedron()\n",
    "print(vert_poly)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also define the polyhedron via its generators, i.e. its vertices and generators of its recession cone.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1, 0, 1], [0, 1, 1]]\n",
      "[[1, 0, 0]]\n"
     ]
    }
   ],
   "source": [
    "poly2 = Cone(vertices=vert_poly,cone=[[1,0]])\n",
    "\n",
    "int_hull2 = poly2.IntegerHull()\n",
    "print(int_hull2.VerticesOfPolyhedron())\n",
    "print(int_hull2.ExtremeRays())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Creating new functions\n",
    "As an example for the usefulness of PyNormaliz we write a small function that creates the intersection of two cones - something which would be quite messy if we work with input files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1, 1], [1, 2]]\n"
     ]
    }
   ],
   "source": [
    "def intersection(cone1, cone2):\n",
    "    intersection_ineq = cone1.SupportHyperplanes()+cone2.SupportHyperplanes()\n",
    "    C = Cone(inequalities = intersection_ineq)\n",
    "    return C\n",
    "\n",
    "C1 = Cone(cone=[[1,2],[2,1]])\n",
    "C2 = Cone(cone=[[1,1],[1,3]])\n",
    "print(intersection(C1,C2).ExtremeRays())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "hideCode": true,
    "hideOutput": true,
    "hidePrompt": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 1, 1, -1, -1, -1, 0, 0, 0]\n",
      "[1, 1, 1, 0, 0, 0, -1, -1, -1]\n",
      "[0, 1, 1, -1, 0, 0, -1, 0, 0]\n",
      "[1, 0, 1, 0, -1, 0, 0, -1, 0]\n",
      "[1, 1, 0, 0, 0, -1, 0, 0, -1]\n",
      "[0, 1, 1, 0, -1, 0, 0, 0, -1]\n",
      "[1, 1, 0, 0, -1, 0, -1, 0, 0]\n"
     ]
    }
   ],
   "source": [
    "def magic_square(n):\n",
    "    equation = [1 if x<n else 0 for x in range(n**2)]\n",
    "    eqs = []\n",
    "\n",
    "    for i in range(1,n):\n",
    "        equation = [1 if x<n else 0 for x in range(n**2)]\n",
    "        for j in range(n):\n",
    "            equation[i * n + j] = -1\n",
    "        eqs.append(equation)\n",
    "\n",
    "    for i in range(n):\n",
    "        equation = [1 if x<n else 0 for x in range(n**2)]\n",
    "        equation[i] = 0\n",
    "        for j in range(1,n):\n",
    "            equation[n*j + i] = -1\n",
    "        eqs.append(equation)\n",
    "\n",
    "    equation = [1 if x<n else 0 for x in range(n**2)]\n",
    "    diagonal = list(equation)\n",
    "    anti_diagonal = list(equation)\n",
    "    diagonal[0] = 0\n",
    "    anti_diagonal[n-1]=0\n",
    "\n",
    "    for i in range(1, n):\n",
    "        diagonal[n*i+i] = -1\n",
    "        anti_diagonal[n*(i+1)-(i+1)] = -1\n",
    "    eqs.append(diagonal)\n",
    "    eqs.append(anti_diagonal)\n",
    "    return eqs\n",
    "\n",
    "eqs = magic_square(3)\n",
    "\n",
    "for line in eqs:\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "A magic square is an $n\\times n$ table in which each row and column and the two diagonals sum up to the same magic constant $\\mathcal{M}$. The function\n",
    "<pre>magic_square(n)</pre>\n",
    "generates all defining equations of the $n\\times n$ magic squares. \n",
    "\n",
    "So only a few lines of code are necessary to create the respective cone and compute its Hilbert basis (in the dual mode) and Hilbert series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n",
      "[[1, 4, 18, 36, 50, 36, 18, 4, 1], [1, 1, 1, 1, 2, 2, 2, 2], 0]\n"
     ]
    }
   ],
   "source": [
    "n=4\n",
    "grading = [1 if x<n else 0 for x in range(n**2)]\n",
    "M = Cone(equations=magic_square(n),grading=grading)\n",
    "print(len(M.HilbertBasis(DualMode=True)))\n",
    "print(M.HilbertSeries())"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Hide code",
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}