{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Here are some preliminary functions. If you just want to read the document, then you can ignore this first cell." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#preliminary functions\n", "\n", "def hyperplane_arrangement_to_fan(hyperplane_arrangement):\n", " '''\n", " Input: Linear hyperplane arrangement\n", " Output: Maximal cones of the induced fan. Each maximal cone is the closure of a connected component \n", " of R^d after removal of the hyperplane arrangement\n", " '''\n", " n = len(hyperplane_arrangement)\n", " p = list(powerset(range(n)))\n", " d = hyperplane_arrangement[0].ambient_dim()\n", " cones = []\n", " for entry in p:\n", " signs = [1 if i in entry else -1 for i in range(n)]\n", " inequalities = []\n", " for i in range(n):\n", " hyperplane = hyperplane_arrangement[i]\n", " A = hyperplane.Hrepresentation()[0].A()\n", " ieq = [0]+[signs[i]*a for a in A]\n", " inequalities.append(ieq)\n", " C = Polyhedron(ieqs = inequalities)\n", " if C.dimension() == d:\n", " cones.append(C)\n", " return cones\n", "\n", "def intersecting_edges(polytope, u):\n", " '''\n", " Input: A polytope and a normal vector u\n", " Output: Set of edges of polytope that has nonempty intersection with the hyperplane orthogonal to u\n", " '''\n", " u_perp = Polyhedron(eqns= [ [0]+list(u) ])\n", " intersecting_edges = []\n", " for edge in polytope.faces(1):\n", " if edge.as_polyhedron().intersection(u_perp).dimension()>-1:\n", " intersecting_edges.append(edge)\n", " return intersecting_edges" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Computing the intersection body of the cube." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider the $3$-dimensional cube $[-1,1]^3$ and compute the radial function of its intersection body step by step." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "P = polytopes.cube()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The combinatorial type of the intersection of $P$ with a hyperplane through the origin is fixed for each region of the normal fan of the following zonotope Z. You can visualize this fan by looking at the dual polytope." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "Graphics3d Object" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Z = sum([Polyhedron(vertices = [-vector(v),vector(v)]) for v in P.vertices()])\n", "dualZ = Z.polar()\n", "dualZ.show(frame=False, point = False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now construct the associated hyperplane arrangement $H$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "H = list(set([Polyhedron(eqns = [[0]+list(vertex)]) for vertex in P.vertices()]))\n", "Fan = hyperplane_arrangement_to_fan(H)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us focus now on one of the open chambers of this hyperplane arrangement, and compute the radial function in here. \n", "We sample a vector in the interior of this region." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-2, 2, 2)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "C = Fan[1]\n", "sample_vector = sum([vector(ray) for ray in C.rays()])\n", "sample_vector" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now collect all facets of $P$ that intersect the hyperplane `perp` orthogonal to `sample_vector`. Note that in this case $Q = P\\cap$`perp` is a hexagon: indeed we are intersecting $6$ facets of $P$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,\n", " A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,\n", " A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,\n", " A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,\n", " A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,\n", " A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "perp = Polyhedron(eqns= [ [0]+list(sample_vector) ])\n", "\n", "intersecting_facets = []\n", "if P.intersection(perp).dimension()==2:\n", " for facet in P.facets():\n", " if facet.as_polyhedron().intersection(perp).dimension() == 1:\n", " intersecting_facets.append(facet)\n", "intersecting_facets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The vertices of such a hexagon arise as the intersection of `perp` with some edges of the above facets of $P$. We now compute these edges." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,\n", " A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],\n", " [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,\n", " A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],\n", " [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,\n", " A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],\n", " [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,\n", " A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],\n", " [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,\n", " A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],\n", " [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,\n", " A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices]]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "intersecting_edges_of_facets = []\n", "for facet in intersecting_facets:\n", " edges = intersecting_edges(facet.as_polyhedron(),sample_vector)\n", " intersecting_edges_of_facets.append(edges)\n", "intersecting_edges_of_facets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let $x \\in C$ arbitrary. By the theory we know that $x^\\perp$ intersects $P$ in the set of edges in `intersecting_edges_of_facets`. Therefore, for each such edge $\\hbox{conv}(a,b)$ we write the corresponding vertex $v$ of $P\\cap x^\\perp$ as:\n", "$$\n", "v = \\frac{(b\\cdot x)\\,a - (a\\cdot x)\\,b}{(b-a)\\cdot x} \\,.\n", "$$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[((x2 - x3)/x1, -1, 1), (-1, -1, (x1 + x2)/x3)],\n", " [(-(x2 - x3)/x1, 1, -1), (-1, (x1 + x3)/x2, -1)],\n", " [(-1, (x1 + x3)/x2, -1), (-1, -1, (x1 + x2)/x3)],\n", " [(1, -(x1 + x3)/x2, 1), ((x2 - x3)/x1, -1, 1)],\n", " [(1, 1, -(x1 + x2)/x3), (-(x2 - x3)/x1, 1, -1)],\n", " [(1, -(x1 + x3)/x2, 1), (1, 1, -(x1 + x2)/x3)]]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = vector(var('x1,x2,x3'))\n", "x1, x2, x3 = x\n", "v = []\n", "for edges_of_facet in intersecting_edges_of_facets:\n", " intersecting_points_of_facets = []\n", " for edge in edges_of_facet:\n", " [a,b] = list(edge.vertices())\n", " a = vector(a)\n", " b = vector(b)\n", " intersection_point = ((b.inner_product(x)*a - a.inner_product(x)*b )/ (b-a).inner_product(x))\n", " intersecting_points_of_facets.append(intersection_point)\n", " v.append(intersecting_points_of_facets)\n", "v" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are now ready to compute the volume of the hexagon. This will be done by subdividing the polygon into simplices, and by computing the volume of each simplex as a determinant." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[-1, -x1*((x1 + x2)/x3 - 1) - x2 - (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 - x3],\n", " [1, x1*((x1 + x3)/x2 - 1) - (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 + x2 + x3],\n", " [1,\n", " x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3],\n", " [-1, -x1*((x1 + x3)/x2 - 1) + (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 - x2 - x3],\n", " [1, x1*((x1 + x2)/x3 - 1) + x2 + (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 + x3],\n", " [1,\n", " x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3]]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "determinants = []\n", "for simplex in v:\n", " M = matrix(\n", " simplex + [x]\n", " )\n", " det_M = M.determinant()\n", " if det:\n", " sgn = sign(det_M(x1 = sample_vector[0], x2 = sample_vector[1], x3 = sample_vector[2]))\n", " determinants.append([sgn,M.determinant()])\n", "determinants\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, by adding all the determinants with a correction factor, we obtain the radial function of $IP$ restricted to the region $C$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x1^2 + 2*x1*x2 + x2^2 + 2*x1*x3 - 2*x2*x3 + x3^2)/(x1*x2*x3)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "correcting_factor = 1/(factorial(P.dimension()-1) * sum([x_i^2 for x_i in x]))\n", "rho = correcting_factor*sum([sng*det for [sng,det] in determinants])\n", "rho.factor()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the radial function we can compute the algebraic boundary of $IP$. In $C$ this is the zero set of the following polynomial." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "x1*x2*x3 - x1^2 - 2*x1*(x2 + x3) - x2^2 + 2*x2*x3 - x3^2" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p, q = rho.numerator_denominator()\n", "alg_boundary = q-p\n", "alg_boundary" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.2", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }