fan_rtree_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sage.all\n",
    "from collections.abc import Callable, Container\n",
    "from copy import copy\n",
    "from warnings import warn\n",
    "\n",
    "from sage.structure.richcmp import richcmp_method, richcmp\n",
    "from sage.combinat.combination import Combinations\n",
    "from sage.combinat.posets.posets import FinitePoset\n",
    "from sage.geometry.cone import (_ambient_space_point,\n",
    "                                Cone,\n",
    "                                ConvexRationalPolyhedralCone,\n",
    "                                IntegralRayCollection,\n",
    "                                is_Cone,\n",
    "                                normalize_rays)\n",
    "from sage.geometry.hasse_diagram import lattice_from_incidences\n",
    "from sage.geometry.point_collection import PointCollection\n",
    "from sage.geometry.toric_lattice import ToricLattice, is_ToricLattice\n",
    "from sage.geometry.toric_plotter import ToricPlotter\n",
    "from sage.graphs.digraph import DiGraph\n",
    "from sage.matrix.all import matrix\n",
    "from sage.misc.cachefunc import cached_method\n",
    "from sage.misc.all import walltime, prod\n",
    "from sage.modules.all import vector, span\n",
    "from sage.rings.all import QQ, ZZ\n",
    "from sage.functions.other import floor, ceil\n",
    "\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "def timeit(func):\n",
    "    def wrapper(*args, **kw):\n",
    "        tik = time.time()\n",
    "        func(*args, **kw)\n",
    "        tok = time.time()\n",
    "        print(\"time is {}\".format(tok-tik))\n",
    "    return wrapper\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sage.geometry.fan import RationalPolyhedralFan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "@richcmp_method\n",
    "class RationalPolyhedralFan_rtree(RationalPolyhedralFan):\n",
    "    r\"\"\"\n",
    "    Create a pointed rational polyhedral fan with rtree data structure, \n",
    "    which reduces the time for determining whether a point is in the support \n",
    "    of the fan or not.\n",
    "\n",
    "    .. WARNING::\n",
    "\n",
    "        This class does not perform any checks of correctness of input nor\n",
    "        does it convert input into the standard representation. Use\n",
    "        :func:`Fan` to construct fans from \"raw data\" or :func:`FaceFan` and\n",
    "        :func:`NormalFan` to get fans associated to polytopes.\n",
    "\n",
    "    Fans are immutable, but they cache most of the returned values.\n",
    "\n",
    "    INPUT:\n",
    "\n",
    "    - ``cones`` -- list of generating cones of the fan, each cone given as a\n",
    "      list of indices of its generating rays in ``rays``;\n",
    "\n",
    "    - ``rays`` -- list of immutable primitive vectors in ``lattice``\n",
    "      consisting of exactly the rays of the fan (i.e. no \"extra\" ones);\n",
    "\n",
    "    - ``lattice`` -- :class:`ToricLattice\n",
    "      <sage.geometry.toric_lattice.ToricLatticeFactory>`, `\\ZZ^n`, or any\n",
    "      other object that behaves like these. If ``None``, it will be determined\n",
    "      as :func:`parent` of the first ray. Of course, this cannot be done if\n",
    "      there are no rays, so in this case you must give an appropriate\n",
    "      ``lattice`` directly;\n",
    "\n",
    "    - ``is_complete`` -- if given, must be ``True`` or ``False`` depending on\n",
    "      whether this fan is complete or not. By default, it will be determined\n",
    "      automatically if necessary;\n",
    "\n",
    "    - ``virtual_rays`` -- if given, must be a list of immutable primitive\n",
    "      vectors in ``lattice``, see :meth:`virtual_rays` for details. By default,\n",
    "      it will be determined automatically if necessary.\n",
    "\n",
    "    OUTPUT:\n",
    "\n",
    "    - a rational polyhedral fan with rtree if it is pointed, or a rational \n",
    "        polyhedral fan without rtree if it is NOT pointed.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, cones, rays, lattice,\n",
    "                 is_complete=None, virtual_rays=None):\n",
    "        r\"\"\"\n",
    "        See :class:`RationalPolyhedralFan_rtree` for documentation.\n",
    "\n",
    "        TESTS::\n",
    "\n",
    "            sage: v = vector([0,1])\n",
    "            sage: v.set_immutable()\n",
    "            sage: f = sage.geometry.fan.RationalPolyhedralFan(\n",
    "            ....:                       [(0,)], [v], None)\n",
    "            sage: f.rays()\n",
    "            (0, 1)\n",
    "            in Ambient free module of rank 2\n",
    "            over the principal ideal domain Integer Ring\n",
    "            sage: TestSuite(f).run()\n",
    "            sage: f = Fan([(0,)], [(0,1)])\n",
    "            sage: TestSuite(f).run()\n",
    "        \"\"\"\n",
    "        #copy generator\n",
    "        tik = time.time()\n",
    "        fan_indices = list(cones)\n",
    "        cones = (c_idx for c_idx in fan_indices)\n",
    "        super().__init__(cones, rays, lattice, is_complete=None, virtual_rays=None)\n",
    "        self._is_pointed = None\n",
    "        \n",
    "        # Determining whether the fan is pointed automatically. \n",
    "        \n",
    "        if len(self._generating_cones) != 0:\n",
    "            #Trying LP solver\n",
    "            p = MixedIntegerLinearProgram(solver = 'GLPK', maximization = True)\n",
    "            x = p.new_variable()\n",
    "            b = p.new_variable()\n",
    "            B = p.new_variable()\n",
    "            p.solver_parameter(\"simplex_or_intopt\", \"simplex_only\") # use simplex only\n",
    "            p.set_objective(B[0])\n",
    "            for i in range(len(rays)):\n",
    "                p.add_constraint(p.sum(x[_]*rays[i][_] for _ in range(self.dim()))<= -b[i])\n",
    "            p.add_constraint(1>= b[i] >= B[0])\n",
    "            for _ in range(self.dim()):\n",
    "                p.add_constraint( -1 <= x[_] <= 1 )\n",
    "            p.add_constraint(1>=B[0]>=0)\n",
    "\n",
    "            try:\n",
    "                # Verifying whether it is a normal vector or not\n",
    "                if p.solve() > 10e-15:\n",
    "                    temp_x = vector(p.get_values(x))\n",
    "                    if all(temp_x.dot_product(ray)< 0 for ray in rays):\n",
    "                        # Been verified, thus, the fan is pointed.\n",
    "                        self._normal_vector = temp_x     \n",
    "                        self._is_pointed = True\n",
    "            except sage.numerical.mip.MIPSolverException:\n",
    "                pass\n",
    "            # If the LP-solver fails, we may use the following facets-based method.\n",
    "            if not self._is_pointed:\n",
    "                # Trying Facets-based method.\n",
    "                # Facets-based method, computing the common intersection of all polar \n",
    "                # cones of the generating cones in the fan. if the common intersection \n",
    "                # is empty, then it is not pointed; otherwise, it is pointed.\n",
    "                empty_cone = Cone([],lattice = lattice)\n",
    "                common_intersection = empty_cone.dual()\n",
    "                for c in self._generating_cones:\n",
    "                    common_intersection = common_intersection.intersection(c.polar())\n",
    "                    if common_intersection.dim() == 0:\n",
    "                        # If common_intersection is empty, then the fan is not pointed.\n",
    "                        self._is_pointed = False\n",
    "                        break\n",
    "                if self._is_pointed == None:\n",
    "                    temp_x = sum(vector(ray, QQ)/vector(ray, QQ).norm(infinity) for ray in common_intersection.rays())  \n",
    "                    if all(temp_x.dot_product(ray)< 0 for ray in rays):\n",
    "                        # Been verified, thus, the fan is pointed.\n",
    "                        self._normal_vector = temp_x      \n",
    "                        self._is_pointed = True \n",
    "        else:\n",
    "            # An empty fan is considered as pointed, but without a normal vector.\n",
    "            self._is_pointed = True\n",
    "\n",
    "        # If the fan is pointed, then we build the rtree for it. ``self._has_rtree`` is a \n",
    "        # trigger for the usage of rtree. Since rtree package can only store 2 to high \n",
    "        # dimensional boxes, we give restriction below.\n",
    "        self._has_rtree = False\n",
    "        if self._is_pointed and self.dim()>=2 and len(self._generating_cones)>0:\n",
    "            # Generating a random family of directions for projections, by \n",
    "            # neuralization of inner product of a direction with the normal vector.\n",
    "            self._has_rtree = True\n",
    "            self._normal_vector = vector(self._normal_vector)\n",
    "            self._directions = []\n",
    "            ndirections = self.dim()\n",
    "            # Finding a coordinate that is non-zero of the normal vector.\n",
    "            for i in range(self.dim()):\n",
    "                if self._normal_vector[i] != 0:\n",
    "                    coord = i\n",
    "                    break\n",
    "\n",
    "            for i in range(ndirections):\n",
    "                direction = vector([choice([-1,1])*randint(1,127) for i in range(self.dim())], QQ)\n",
    "                diff = direction.dot_product(self._normal_vector)\n",
    "                direction[coord] = direction[coord] - diff/self._normal_vector[coord]\n",
    "                self._directions.append(direction)\n",
    "            table = []\n",
    "            for (i,d) in enumerate(self._directions):\n",
    "                t_i = []\n",
    "                abs_max = 0\n",
    "                for r in rays:\n",
    "                    val = r.dot_product(d)/r.dot_product(self._normal_vector)\n",
    "                    if abs_max < abs(val):\n",
    "                        abs_max = abs(val)\n",
    "                    t_i.append(val)\n",
    "                if abs_max > 10e-15:\n",
    "                    ratio = 2**26/ceil(abs_max)\n",
    "                    self._directions[i] = ratio * self._directions[i]\n",
    "                    t_i = [val * ratio for val in t_i]\n",
    "                table.append(t_i)\n",
    "\n",
    "            # The intersection between each cone and the separating hyperplane, which is -1 distance away\n",
    "            # from the origin, is defined as the polytopical token for each cone.\n",
    "            # Let r be the point that is the intersection between a ray and the separating hyperplane, x\n",
    "            # be the unit vector of the normal vector, u be the unit vector of the ray, r = k*u, k>0.\n",
    "            # We have <-x,u>= cos t = 1 / k. Thus, k = -1/<x,u>, r = k*u = -u/<x,u>. For any vector v \n",
    "            # whose direction is along with the ray r, v = j * u, because u = v/j, we have \n",
    "            # r = -(v/j)/<x,v/j>, = -v/<x,v>.\n",
    "            tok = time.time()\n",
    "            print(\"Before rtree{}s.\".format(tok-tik))\n",
    "            from rtree import index\n",
    "            p = index.Property()\n",
    "            p.dimension = ndirections\n",
    "            self._rtree = index.Index(properties = p, interleaved = False)\n",
    "            time_sum = 0\n",
    "            for (i, indices) in enumerate(fan_indices):\n",
    "                box = []\n",
    "                for _ in range(ndirections):\n",
    "                    proj_min, proj_max = infinity, -infinity\n",
    "                    for idx in indices:\n",
    "                        if table[_][idx] < proj_min:\n",
    "                            proj_min = table[_][idx]\n",
    "                        if table[_][idx] > proj_max:\n",
    "                            proj_max = table[_][idx]\n",
    "                    box = box + self._adjusting_gauge(proj_min, proj_max)\n",
    "                self._rtree.insert(0, box, obj = i)\n",
    "\n",
    "            \n",
    "    def _adjusting_gauge(self, lower, upper):\n",
    "        # This function is used for inserting boxes in rtree. Since rtree has limited precision, the rational \n",
    "        # interval has to been converted to a compatible integer interval.\n",
    "        # The domain of rtree only consists of -\\infty, integers in [-2^52,2^52], +\\infty.\n",
    "        prec_range = (-2**52, 2**52)\n",
    "        if lower<=prec_range[0]:\n",
    "            x_k = -infinity\n",
    "        elif lower<prec_range[1]+1:\n",
    "            x_k = floor(lower)\n",
    "        else:\n",
    "            x_k = infinity\n",
    "        if upper<prec_range[0]-1:\n",
    "            y_k = -infinity\n",
    "        elif upper<=prec_range[1]:\n",
    "            y_k = ceil(upper)\n",
    "        else:\n",
    "            y_k = infinity\n",
    "        return [x_k, y_k]\n",
    "\n",
    "\n",
    "    def is_pointed(self):\n",
    "        return self._is_pointed\n",
    "    \n",
    "    def support_contains(self, *args):\n",
    "        r\"\"\"\n",
    "        Check if a point is contained in the support of the fan.\n",
    "\n",
    "        The support of a fan is the union of all cones of the fan. If\n",
    "        you want to know whether the fan contains a given cone, you\n",
    "        should use :meth:`contains` instead.\n",
    "\n",
    "        INPUT:\n",
    "\n",
    "        - ``*args`` -- an element of ``self.lattice()`` or something\n",
    "          that can be converted to it (for example, a list of\n",
    "          coordinates).\n",
    "\n",
    "        OUTPUT:\n",
    "\n",
    "        - ``True`` if ``point`` is contained in the support of the\n",
    "          fan, ``False`` otherwise.\n",
    "\n",
    "        TESTS::\n",
    "\n",
    "            sage: cone1 = Cone([(0,-1), (1,0)])\n",
    "            sage: cone2 = Cone([(1,0), (0,1)])\n",
    "            sage: f = Fan([cone1, cone2])\n",
    "\n",
    "        We check if some points are in this fan::\n",
    "\n",
    "            sage: f.support_contains(f.lattice()(1,0))\n",
    "            True\n",
    "            sage: f.support_contains(cone1)    # a cone is not a point of the lattice\n",
    "            False\n",
    "            sage: f.support_contains((1,0))\n",
    "            True\n",
    "            sage: f.support_contains(1,1)\n",
    "            True\n",
    "            sage: f.support_contains((-1,0))\n",
    "            False\n",
    "            sage: f.support_contains(f.lattice().dual()(1,0)) #random output (warning)\n",
    "            False\n",
    "            sage: f.support_contains(f.lattice().dual()(1,0))\n",
    "            False\n",
    "            sage: f.support_contains(1)\n",
    "            False\n",
    "            sage: f.support_contains(0)   # 0 converts to the origin in the lattice\n",
    "            True\n",
    "            sage: f.support_contains(1/2, sqrt(3))\n",
    "            True\n",
    "            sage: f.support_contains(-1/2, sqrt(3))\n",
    "            False\n",
    "        \"\"\"\n",
    "\n",
    "        if len(args) == 1:\n",
    "            point = args[0]\n",
    "        else:\n",
    "            point = args\n",
    "        try:\n",
    "            point = _ambient_space_point(self, point)\n",
    "        except TypeError as ex:\n",
    "            if str(ex).endswith(\"have incompatible lattices!\"):\n",
    "                warn(\"you have checked if a fan contains a point \"\n",
    "                     \"from an incompatible lattice, this is False!\",\n",
    "                     stacklevel=3)\n",
    "            return False\n",
    "        if self.is_complete():\n",
    "            return True\n",
    "        \n",
    "        # Using rtree for overestimating.\n",
    "        if self._has_rtree:\n",
    "            r = vector(point)\n",
    "            # See whether the point has a positive direction with the normal vector.\n",
    "            if r.dot_product(self._normal_vector)>=0: \n",
    "                return False\n",
    "            box = []\n",
    "            for d in self._directions:\n",
    "                val = r.dot_product(d)/r.dot_product(self._normal_vector)\n",
    "                box = box + self._adjusting_gauge(val, val)\n",
    "            if self._rtree.count(box) == 0:\n",
    "                return False\n",
    "\n",
    "        return any(point in cone for cone in self)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": [
     378
    ]
   },
   "outputs": [],
   "source": [
    "def Fan_rtree(cones, rays=None, lattice=None, check=True, normalize=True,\n",
    "        is_complete=None, virtual_rays=None, discard_faces=False,\n",
    "        allow_arrangement=False):\n",
    "    r\"\"\"\n",
    "    Construct a rational polyhedral fan.\n",
    "\n",
    "    .. NOTE::\n",
    "\n",
    "        Approximate time to construct a fan consisting of `n` cones is `n^2/5`\n",
    "        seconds. That is half an hour for 100 cones. This time can be\n",
    "        significantly reduced in the future, but it is still likely to be\n",
    "        `\\sim n^2` (with, say, `/500` instead of `/5`). If you know that your\n",
    "        input does form a valid fan, use ``check=False`` option to skip\n",
    "        consistency checks.\n",
    "\n",
    "    INPUT:\n",
    "\n",
    "    - ``cones`` -- list of either\n",
    "      :class:`Cone<sage.geometry.cone.ConvexRationalPolyhedralCone>` objects\n",
    "      or lists of integers interpreted as indices of generating rays in\n",
    "      ``rays``. These must be only **maximal** cones of the fan, unless\n",
    "      ``discard_faces=True`` or ``allow_arrangement=True`` option is specified;\n",
    "\n",
    "    - ``rays`` -- list of rays given as list or vectors convertible to the\n",
    "      rational extension of ``lattice``. If ``cones`` are given by\n",
    "      :class:`Cone<sage.geometry.cone.ConvexRationalPolyhedralCone>` objects\n",
    "      ``rays`` may be determined automatically. You still may give them\n",
    "      explicitly to ensure a particular order of rays in the fan. In this case\n",
    "      you must list all rays that appear in ``cones``. You can give \"extra\"\n",
    "      ones if it is convenient (e.g. if you have a big list of rays for\n",
    "      several fans), but all \"extra\" rays will be discarded;\n",
    "\n",
    "    - ``lattice`` -- :class:`ToricLattice\n",
    "      <sage.geometry.toric_lattice.ToricLatticeFactory>`, `\\ZZ^n`, or any\n",
    "      other object that behaves like these. If not specified, an attempt will\n",
    "      be made to determine an appropriate toric lattice automatically;\n",
    "\n",
    "    - ``check`` -- by default the input data will be checked for correctness\n",
    "      (e.g. that intersection of any two given cones is a face of each),\n",
    "      unless ``allow_arrangement=True`` option is specified. If you\n",
    "      know for sure that the input is correct, you may significantly decrease\n",
    "      construction time using ``check=False`` option;\n",
    "\n",
    "    - ``normalize`` -- you can further speed up construction using\n",
    "      ``normalize=False`` option. In this case ``cones`` must be a list of\n",
    "      **sorted** :class:`tuples` and ``rays`` must be immutable primitive\n",
    "      vectors in ``lattice``. In general, you should not use this option, it\n",
    "      is designed for code optimization and does not give as drastic\n",
    "      improvement in speed as the previous one;\n",
    "\n",
    "    - ``is_complete`` -- every fan can determine on its own if it is complete\n",
    "      or not, however it can take quite a bit of time for \"big\" fans with many\n",
    "      generating cones. On the other hand, in some situations it is known in\n",
    "      advance that a certain fan is complete. In this case you can pass\n",
    "      ``is_complete=True`` option to speed up some computations. You may also\n",
    "      pass ``is_complete=False`` option, although it is less likely to be\n",
    "      beneficial. Of course, passing a wrong value can compromise the\n",
    "      integrity of data structures of the fan and lead to wrong results, so\n",
    "      you should be very careful if you decide to use this option;\n",
    "\n",
    "    - ``virtual_rays`` -- (optional, computed automatically if needed) a list of\n",
    "      ray generators to be used for :meth:`virtual_rays`;\n",
    "\n",
    "    - ``discard_faces`` -- by default, the fan constructor expects the list of\n",
    "      **maximal** cones, unless ``allow_arrangement=True`` option is specified.\n",
    "      If you provide \"extra\" ones and leave ``allow_arrangement=False`` (default)\n",
    "      and ``check=True`` (default), an exception will be raised.\n",
    "      If you provide \"extra\" cones and set ``allow_arrangement=False`` (default)\n",
    "      and ``check=False``, you may get wrong results as assumptions on internal\n",
    "      data structures will be invalid. If you want the fan constructor to\n",
    "      select the maximal cones from the given input, you may provide\n",
    "      ``discard_faces=True`` option (it works both for ``check=True`` and\n",
    "      ``check=False``).\n",
    "\n",
    "    - ``allow_arrangement`` -- by default (``allow_arrangement=False``),\n",
    "      the fan constructor expects that the intersection of any two given cones is\n",
    "      a face of each. If ``allow_arrangement=True`` option is specified, then\n",
    "      construct a rational polyhedralfan from the cone arrangement, so that the\n",
    "      union of the cones in the polyhedral fan equals to the union of the given\n",
    "      cones, and each given cone is the union of some cones in the polyhedral fan.\n",
    "\n",
    "    OUTPUT:\n",
    "\n",
    "    - a :class:`fan <RationalPolyhedralFan>`.\n",
    "\n",
    "    .. SEEALSO::\n",
    "\n",
    "        In 2 dimensions you can cyclically order the rays. Hence the\n",
    "        rays determine a unique maximal fan without having to specify\n",
    "        the cones, and you can use :func:`Fan2d` to construct this\n",
    "        fan from just the rays.\n",
    "\n",
    "    EXAMPLES:\n",
    "\n",
    "    Let's construct a fan corresponding to the projective plane in several\n",
    "    ways::\n",
    "\n",
    "        sage: cone1 = Cone([(1,0), (0,1)])\n",
    "        sage: cone2 = Cone([(0,1), (-1,-1)])\n",
    "        sage: cone3 = Cone([(-1,-1), (1,0)])\n",
    "        sage: P2 = Fan([cone1, cone2, cone2])\n",
    "        Traceback (most recent call last):\n",
    "        ...\n",
    "        ValueError: you have provided 3 cones, but only 2 of them are maximal!\n",
    "        Use discard_faces=True if you indeed need to construct a fan from\n",
    "        these cones.\n",
    "\n",
    "    Oops! There was a typo and ``cone2`` was listed twice as a generating cone\n",
    "    of the fan. If it was intentional (e.g. the list of cones was generated\n",
    "    automatically and it is possible that it contains repetitions or faces of\n",
    "    other cones), use ``discard_faces=True`` option::\n",
    "\n",
    "        sage: P2 = Fan([cone1, cone2, cone2], discard_faces=True)\n",
    "        sage: P2.ngenerating_cones()\n",
    "        2\n",
    "\n",
    "    However, in this case it was definitely a typo, since the fan of\n",
    "    `\\mathbb{P}^2` has 3 maximal cones::\n",
    "\n",
    "        sage: P2 = Fan([cone1, cone2, cone3])\n",
    "        sage: P2.ngenerating_cones()\n",
    "        3\n",
    "\n",
    "    Looks better. An alternative way is ::\n",
    "\n",
    "        sage: rays = [(1,0), (0,1), (-1,-1)]\n",
    "        sage: cones = [(0,1), (1,2), (2,0)]\n",
    "        sage: P2a = Fan(cones, rays)\n",
    "        sage: P2a.ngenerating_cones()\n",
    "        3\n",
    "        sage: P2 == P2a\n",
    "        False\n",
    "\n",
    "    That may seem wrong, but it is not::\n",
    "\n",
    "        sage: P2.is_equivalent(P2a)\n",
    "        True\n",
    "\n",
    "    See :meth:`~RationalPolyhedralFan.is_equivalent` for details.\n",
    "\n",
    "    Yet another way to construct this fan is ::\n",
    "\n",
    "        sage: P2b = Fan(cones, rays, check=False)\n",
    "        sage: P2b.ngenerating_cones()\n",
    "        3\n",
    "        sage: P2a == P2b\n",
    "        True\n",
    "\n",
    "    If you try the above examples, you are likely to notice the difference in\n",
    "    speed, so when you are sure that everything is correct, it is a good idea\n",
    "    to use ``check=False`` option. On the other hand, it is usually **NOT** a\n",
    "    good idea to use ``normalize=False`` option::\n",
    "\n",
    "        sage: P2c = Fan(cones, rays, check=False, normalize=False)\n",
    "        Traceback (most recent call last):\n",
    "        ...\n",
    "        AttributeError: 'tuple' object has no attribute 'parent'\n",
    "\n",
    "    Yet another way is to use functions :func:`FaceFan` and :func:`NormalFan`\n",
    "    to construct fans from :class:`lattice polytopes\n",
    "    <sage.geometry.lattice_polytope.LatticePolytopeClass>`.\n",
    "\n",
    "    We have not yet used ``lattice`` argument, since if was determined\n",
    "    automatically::\n",
    "\n",
    "        sage: P2.lattice()\n",
    "        2-d lattice N\n",
    "        sage: P2b.lattice()\n",
    "        2-d lattice N\n",
    "\n",
    "    However, it is necessary to specify it explicitly if you want to construct\n",
    "    a fan without rays or cones::\n",
    "\n",
    "        sage: Fan([], [])\n",
    "        Traceback (most recent call last):\n",
    "        ...\n",
    "        ValueError: you must specify the lattice\n",
    "        when you construct a fan without rays and cones!\n",
    "        sage: F = Fan([], [], lattice=ToricLattice(2, \"L\"))\n",
    "        sage: F\n",
    "        Rational polyhedral fan in 2-d lattice L\n",
    "        sage: F.lattice_dim()\n",
    "        2\n",
    "        sage: F.dim()\n",
    "        0\n",
    "\n",
    "    In the following examples, we test the ``allow_arrangement=True`` option.\n",
    "    See :trac:`25122`.\n",
    "\n",
    "    The intersection of the two cones is not a face of each. Therefore,\n",
    "    they do not belong to the same rational polyhedral fan::\n",
    "\n",
    "        sage: c1 = Cone([(-2,-1,1), (-2,1,1), (2,1,1), (2,-1,1)])\n",
    "        sage: c2 = Cone([(-1,-2,1), (-1,2,1), (1,2,1), (1,-2,1)])\n",
    "        sage: c1.intersection(c2).is_face_of(c1)\n",
    "        False\n",
    "        sage: c1.intersection(c2).is_face_of(c2)\n",
    "        False\n",
    "        sage: Fan([c1, c2])\n",
    "        Traceback (most recent call last):\n",
    "        ...\n",
    "        ValueError: these cones cannot belong to the same fan!\n",
    "        ...\n",
    "\n",
    "    Let's construct the fan using ``allow_arrangement=True`` option::\n",
    "\n",
    "        sage: fan = Fan([c1, c2], allow_arrangement=True)\n",
    "        sage: fan.ngenerating_cones()\n",
    "        5\n",
    "\n",
    "    Another example where cone c2 is inside cone c1::\n",
    "\n",
    "        sage: c1 = Cone([(4, 0, 0), (0, 4, 0), (0, 0, 4)])\n",
    "        sage: c2 = Cone([(2, 1, 1), (1, 2, 1), (1, 1, 2)])\n",
    "        sage: fan = Fan([c1, c2], allow_arrangement=True)\n",
    "        sage: fan.ngenerating_cones()\n",
    "        7\n",
    "        sage: fan.plot()\n",
    "        Graphics3d Object\n",
    "\n",
    "    Cones of different dimension::\n",
    "\n",
    "        sage: c1 = Cone([(1,0),(0,1)])\n",
    "        sage: c2 = Cone([(2,1)])\n",
    "        sage: c3 = Cone([(-1,-2)])\n",
    "        sage: fan = Fan([c1, c2, c3], allow_arrangement=True)\n",
    "        sage: for cone in sorted(fan.generating_cones()): print(sorted(cone.rays()))\n",
    "        [N(-1, -2)]\n",
    "        [N(0, 1), N(1, 2)]\n",
    "        [N(1, 0), N(2, 1)]\n",
    "        [N(1, 2), N(2, 1)]\n",
    "\n",
    "    A 3-d cone and a 1-d cone::\n",
    "\n",
    "        sage: c3 = Cone([[0, 1, 1], [1, 0, 1], [0, -1, 1], [-1, 0, 1]])\n",
    "        sage: c1 = Cone([[0, 0, 1]])\n",
    "        sage: fan1 = Fan([c1, c3], allow_arrangement=True)\n",
    "        sage: fan1.plot()\n",
    "        Graphics3d Object\n",
    "\n",
    "    A 3-d cone and two 2-d cones::\n",
    "\n",
    "        sage: c2v = Cone([[0, 1, 1], [0, -1, 1]])\n",
    "        sage: c2h = Cone([[1, 0, 1], [-1, 0, 1]])\n",
    "        sage: fan2 = Fan([c2v, c2h, c3], allow_arrangement=True)\n",
    "        sage: fan2.is_simplicial()\n",
    "        True\n",
    "        sage: fan2.is_equivalent(fan1)\n",
    "        True\n",
    "    \"\"\"\n",
    "    def result():\n",
    "        # \"global\" does not work here...\n",
    "        R, V = rays, virtual_rays\n",
    "        if V is not None:\n",
    "            if normalize:\n",
    "                V = normalize_rays(V, lattice)\n",
    "            if check:\n",
    "                R = PointCollection(V, lattice)\n",
    "                V = PointCollection(V, lattice)\n",
    "                d = lattice.dimension()\n",
    "                if len(V) != d - R.dim() or (R + V).dim() != d:\n",
    "                    raise ValueError(\"virtual rays must be linearly \"\n",
    "                    \"independent and with other rays span the ambient space.\")\n",
    "        return RationalPolyhedralFan_rtree(cones, R, lattice, is_complete, V)\n",
    "\n",
    "    if not check and not normalize and not discard_faces and not allow_arrangement:\n",
    "        return result()\n",
    "    if not isinstance(cones, list):\n",
    "        try:\n",
    "            cones = list(cones)\n",
    "        except TypeError:\n",
    "            raise TypeError(\n",
    "                \"cones must be given as an iterable!\"\n",
    "                \"\\nGot: %s\" % cones)\n",
    "    if not cones:\n",
    "        if lattice is None:\n",
    "            if rays is not None and rays:\n",
    "                lattice = normalize_rays(rays, lattice)[0].parent()\n",
    "            else:\n",
    "                raise ValueError(\"you must specify the lattice when you \"\n",
    "                                 \"construct a fan without rays and cones!\")\n",
    "        cones = ((), )\n",
    "        rays = ()\n",
    "        return result()\n",
    "    if is_Cone(cones[0]):\n",
    "        # Construct the fan from Cone objects\n",
    "        if lattice is None:\n",
    "            lattice = cones[0].lattice()\n",
    "            # If we determine the lattice automatically, we do not\n",
    "            # want to force any conversion. TODO: take into account\n",
    "            # coercions?\n",
    "            if check:\n",
    "                for cone in cones:\n",
    "                    if cone.lattice() != lattice:\n",
    "                        raise ValueError(\"cones belong to different lattices \"\n",
    "                            \"(%s and %s), cannot determine the lattice of the \"\n",
    "                            \"fan!\" % (lattice, cone.lattice()))\n",
    "        for i, cone in enumerate(cones):\n",
    "            if cone.lattice() != lattice:\n",
    "                cones[i] = Cone(cone.rays(), lattice, check=False)\n",
    "        if check:\n",
    "            for cone in cones:\n",
    "                if not cone.is_strictly_convex():\n",
    "                    raise ValueError(\n",
    "                                    \"cones of a fan must be strictly convex!\")\n",
    "        # Optimization for fans generated by a single cone\n",
    "        if len(cones) == 1 and rays is None:\n",
    "            cone = cones[0]\n",
    "            cones = (tuple(range(cone.nrays())), )\n",
    "            rays = cone.rays()\n",
    "            is_complete = lattice.dimension() == 0\n",
    "            return result()\n",
    "        if allow_arrangement:\n",
    "            cones = _refine_arrangement_to_fan(cones)\n",
    "            cones = _discard_faces(cones)\n",
    "        elif check:\n",
    "            # Maybe we should compute all faces of all cones and save them for\n",
    "            # later if we are doing this check?\n",
    "            generating_cones = []\n",
    "            for cone in sorted(cones, key=lambda cone: cone.dim(),\n",
    "                               reverse=True):\n",
    "                is_generating = True\n",
    "                for g_cone in generating_cones:\n",
    "                    i_cone = cone.intersection(g_cone)\n",
    "                    if i_cone.is_face_of(cone) and i_cone.is_face_of(g_cone):\n",
    "                        if i_cone.dim() == cone.dim():\n",
    "                            is_generating = False  # cone is a face of g_cone\n",
    "                            break\n",
    "                    else:\n",
    "                        raise ValueError(\n",
    "                                \"these cones cannot belong to the same fan!\"\n",
    "                                \"\\nCone 1 rays: %s\\nCone 2 rays: %s\"\n",
    "                                % (g_cone.rays(), cone.rays()))\n",
    "                if is_generating:\n",
    "                    generating_cones.append(cone)\n",
    "            if len(cones) > len(generating_cones):\n",
    "                if discard_faces:\n",
    "                    cones = generating_cones\n",
    "                else:\n",
    "                    raise ValueError(\"you have provided %d cones, but only %d \"\n",
    "                        \"of them are maximal! Use discard_faces=True if you \"\n",
    "                        \"indeed need to construct a fan from these cones.\" %\n",
    "                        (len(cones), len(generating_cones)))\n",
    "        elif discard_faces:\n",
    "            cones = _discard_faces(cones)\n",
    "        ray_set = set([])\n",
    "        for cone in cones:\n",
    "            ray_set.update(cone.rays())\n",
    "        if rays:    # Preserve the initial order of rays, if they were given\n",
    "            rays = normalize_rays(rays, lattice)\n",
    "            new_rays = []\n",
    "            for ray in rays:\n",
    "                if ray in ray_set and ray not in new_rays:\n",
    "                    new_rays.append(ray)\n",
    "            if len(new_rays) != len(ray_set):\n",
    "                raise ValueError(\n",
    "                  \"if rays are given, they must include all rays of the fan!\")\n",
    "            rays = new_rays\n",
    "        else:\n",
    "            rays = tuple(sorted(ray_set))\n",
    "        cones = (tuple(sorted(rays.index(ray) for ray in cone.rays()))\n",
    "                 for cone in cones)\n",
    "        return result()\n",
    "    # Construct the fan from rays and \"tuple cones\"\n",
    "    rays = normalize_rays(rays, lattice)\n",
    "    for n, cone in enumerate(cones):\n",
    "        try:\n",
    "            cones[n] = sorted(cone)\n",
    "        except TypeError:\n",
    "            raise TypeError(\"cannot interpret %s as a cone!\" % cone)\n",
    "    if not check and not discard_faces and not allow_arrangement:\n",
    "        return result()\n",
    "    # If we do need to make all the check, build explicit cone objects first\n",
    "    return Fan_rtree((Cone([rays[n] for n in cone], lattice) for cone in cones),\n",
    "               rays, lattice, is_complete=is_complete,\n",
    "               virtual_rays=virtual_rays, discard_faces=discard_faces,\n",
    "               allow_arrangement=allow_arrangement)\n",
    "\n",
    "def discard_faces(cones):\n",
    "    r\"\"\"\n",
    "    Return the cones of the given list which are not faces of each other.\n",
    "\n",
    "    INPUT:\n",
    "\n",
    "    - ``cones`` -- a list of\n",
    "      :class:`cones <sage.geometry.cone.ConvexRationalPolyhedralCone>`.\n",
    "\n",
    "    OUTPUT:\n",
    "\n",
    "    - a list of\n",
    "      :class:`cones <sage.geometry.cone.ConvexRationalPolyhedralCone>`,\n",
    "      sorted by dimension in decreasing order.\n",
    "\n",
    "    EXAMPLES:\n",
    "\n",
    "    Consider all cones of a fan::\n",
    "\n",
    "        sage: Sigma = toric_varieties.P2().fan()\n",
    "        sage: cones = flatten(Sigma.cones())\n",
    "        sage: len(cones)\n",
    "        7\n",
    "\n",
    "    Most of them are not necessary to generate this fan::\n",
    "\n",
    "        sage: from sage.geometry.fan import discard_faces\n",
    "        sage: len(discard_faces(cones))\n",
    "        3\n",
    "        sage: Sigma.ngenerating_cones()\n",
    "        3\n",
    "    \"\"\"\n",
    "    # Convert to a list or make a copy, so that the input is unchanged.\n",
    "    cones = list(cones)\n",
    "    cones.sort(key=lambda cone: cone.dim(), reverse=True)\n",
    "    generators = []\n",
    "    for cone in cones:\n",
    "        if not any(cone.is_face_of(other) for other in generators):\n",
    "            generators.append(cone)\n",
    "    return generators\n",
    "\n",
    "\n",
    "_discard_faces = discard_faces  # Due to a name conflict in Fan constructor\n",
    "\n",
    "\n",
    "def _refine_arrangement_to_fan(cones):\n",
    "    \"\"\"\n",
    "    Refine the cones of the given list so that they can belong to the same fan.\n",
    "\n",
    "    INPUT:\n",
    "\n",
    "    - ``cones`` -- a list of rational cones that are possibly overlapping.\n",
    "\n",
    "    OUTPUT:\n",
    "\n",
    "    - a list of refined cones.\n",
    "\n",
    "    EXAMPLES::\n",
    "\n",
    "        sage: from sage.geometry.fan import _refine_arrangement_to_fan\n",
    "        sage: c1 = Cone([(-2,-1,1), (-2,1,1), (2,1,1), (2,-1,1)])\n",
    "        sage: c2 = Cone([(-1,-2,1), (-1,2,1), (1,2,1), (1,-2,1)])\n",
    "        sage: refined_cones = _refine_arrangement_to_fan([c1, c2])\n",
    "        sage: for cone in refined_cones: print(cone.rays())\n",
    "        N(-1,  1, 1),\n",
    "        N(-1, -1, 1),\n",
    "        N( 1, -1, 1),\n",
    "        N( 1,  1, 1)\n",
    "        in 3-d lattice N\n",
    "        N(1, -1, 1),\n",
    "        N(1,  1, 1),\n",
    "        N(2, -1, 1),\n",
    "        N(2,  1, 1)\n",
    "        in 3-d lattice N\n",
    "        N(-2,  1, 1),\n",
    "        N(-1, -1, 1),\n",
    "        N(-1,  1, 1),\n",
    "        N(-2, -1, 1)\n",
    "        in 3-d lattice N\n",
    "        N(-1, 1, 1),\n",
    "        N(-1, 2, 1),\n",
    "        N( 1, 1, 1),\n",
    "        N( 1, 2, 1)\n",
    "        in 3-d lattice N\n",
    "        N(-1, -1, 1),\n",
    "        N(-1, -2, 1),\n",
    "        N( 1, -2, 1),\n",
    "        N( 1, -1, 1)\n",
    "        in 3-d lattice N\n",
    "    \"\"\"\n",
    "    dual_lattice = cones[0].dual_lattice()\n",
    "    is_face_to_face = True\n",
    "    for i in range(len(cones)):\n",
    "        ci = cones[i]\n",
    "        for j in range(i):\n",
    "            cj = cones[j]\n",
    "            c = ci.intersection(cj)\n",
    "            if not (c.is_face_of(ci)) or not (c.is_face_of(cj)):\n",
    "                is_face_to_face = False\n",
    "                break\n",
    "        if not is_face_to_face:\n",
    "            break\n",
    "    if is_face_to_face:\n",
    "        return cones\n",
    "    facet_normal_vectors = []\n",
    "    for c in cones:\n",
    "        for l in c.polyhedron().Hrepresentation():\n",
    "            v = l[1::]\n",
    "            is_new = True\n",
    "            for fnv in facet_normal_vectors:\n",
    "                if span([v, fnv]).dimension() < 2:\n",
    "                    is_new = False\n",
    "                    break\n",
    "            if is_new:\n",
    "                facet_normal_vectors.append(v)\n",
    "    for v in facet_normal_vectors:\n",
    "        halfspace1 = Cone([v], lattice=dual_lattice).dual()\n",
    "        halfspace2 = Cone([-v], lattice=dual_lattice).dual()\n",
    "        subcones = []\n",
    "        for c in cones:\n",
    "            subc1 = c.intersection(halfspace1)\n",
    "            subc2 = c.intersection(halfspace2)\n",
    "            for subc in [subc1, subc2]:\n",
    "                if subc.dim() == c.dim():\n",
    "                    is_new = True\n",
    "                    for subcone in subcones:\n",
    "                        if subc.dim() == subcone.dim() and subc.is_equivalent(subcone):\n",
    "                            is_new = False\n",
    "                            break\n",
    "                    if is_new:\n",
    "                        subcones.append(subc)\n",
    "        cones = subcones\n",
    "    return cones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sage.misc.prandom import random, choice\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 5\n",
    "num = 5\n",
    "rays_in_cone = 5\n",
    "nrays = 50\n",
    "ntests = 1000\n",
    "\n",
    "def inner_product(x,y):\n",
    "    return sum([x[_]*y[_] for _ in range(len(x))])\n",
    "\n",
    "direction_ray = tuple(choice([-1,1])*randint(1,100) for i in range(dim))\n",
    "generating_cones = []\n",
    "i = 0\n",
    "rays = []\n",
    "while i<nrays:\n",
    "    ray = tuple(choice([-1,1])*randint(1,100) for i in range(dim))\n",
    "    if inner_product(direction_ray, ray) <= 0:\n",
    "        continue\n",
    "    else:\n",
    "        rays.append(ray)\n",
    "        i+=1\n",
    "\n",
    "generating_cones = []\n",
    "for i in range(num):\n",
    "    c_nrays = choice([1 .. rays_in_cone])\n",
    "    c_rays = []\n",
    "    for j in [1 .. c_nrays]:\n",
    "        c_rays.append(rays[choice([0..nrays-1])])\n",
    "    c = Cone(rays = c_rays)\n",
    "    generating_cones.append(c)\n",
    "        \n",
    "testing_cone = []\n",
    "for i in range(ntests):\n",
    "    c_nrays = choice([1 .. rays_in_cone])\n",
    "    c_rays = []\n",
    "    for j in [1 .. c_nrays]:\n",
    "        c_rays.append(rays[choice([0..nrays-1])])\n",
    "    c = Cone(rays = c_rays)\n",
    "    testing_cone.append(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_points = []\n",
    "for i in range(1000):\n",
    "    ray = tuple(randint(-1000,1000) for i in range(dim))\n",
    "    testing_points.append(ray) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time is 102.33762097358704\n",
      "Before rtree60.47275185585022s.\n",
      "time is 113.38309407234192\n",
      "2728\n"
     ]
    }
   ],
   "source": [
    "tik = time.time()\n",
    "D = Fan(generating_cones,discard_faces=True,allow_arrangement=True)\n",
    "tok = time.time()\n",
    "print(\"time is {}\".format(tok-tik))\n",
    "\n",
    "tik = time.time()\n",
    "C = Fan_rtree(generating_cones,discard_faces=True,allow_arrangement=True)\n",
    "tok = time.time()\n",
    "print(\"time is {}\".format(tok-tik))\n",
    "\n",
    "print(len(C.generating_cones()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "@timeit\n",
    "def D_result(L):\n",
    "    for c in testing_points:\n",
    "        L.append(D.support_contains(c))\n",
    "        \n",
    "@timeit\n",
    "def C_result(L):\n",
    "    for c in testing_points:\n",
    "        L.append(C.support_contains(c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time is 140.72904682159424\n",
      "time is 55.95819425582886\n"
     ]
    }
   ],
   "source": [
    "E = []\n",
    "F = []\n",
    "D_result(E)\n",
    "C_result(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "good!\n"
     ]
    }
   ],
   "source": [
    "if all(E[i]==F[i] for i,x in enumerate(E)):\n",
    "    print(\"good!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rtree.index.Index(bounds=[-67066469.0, 22895649.0, -44152650.0, 66957976.0, -67107742.0, 24429166.0, -33160376.0, 67059953.0, -26333891.0, 67082013.0], size=2728)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C._rtree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.5.beta0",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}