fan_morphism.py

r"""
Morphisms between toric lattices compatible with fans

This module is a part of the framework for toric varieties
(:mod:`~sage.schemes.toric.variety`,
:mod:`~sage.schemes.toric.fano_variety`). Its main purpose is to
provide support for working with lattice morphisms compatible with fans via
:class:`FanMorphism` class.

AUTHORS:

- Andrey Novoseltsev (2010-10-17): initial version.
- Andrey Novoseltsev (2011-04-11): added tests for injectivity/surjectivity,
    fibration, bundle, as well as some related methods.

EXAMPLES:

Let's consider the face and normal fans of the "diamond" and the projection
to the `x`-axis::

    sage: diamond = lattice_polytope.octahedron(2)
    sage: face = FaceFan(diamond)
    sage: normal = NormalFan(diamond)
    sage: N = face.lattice()
    sage: H = End(N)
    sage: phi = H([N.0, 0])
    sage: phi
    Free module morphism defined by the matrix
    [1 0]
    [0 0]
    Domain: 2-d lattice N
    Codomain: 2-d lattice N
    sage: FanMorphism(phi, normal, face)
    Traceback (most recent call last):
    ...
    ValueError: the image of generating cone #1 of the domain fan
    is not contained in a single cone of the codomain fan!

Some of the cones of the normal fan fail to be mapped to a single cone of the
face fan. We can rectify the situation in the following way::

    sage: fm = FanMorphism(phi, normal, face, subdivide=True)
    sage: fm
    Fan morphism defined by the matrix
    [1 0]
    [0 0]
    Domain fan: Rational polyhedral fan in 2-d lattice N
    Codomain fan: Rational polyhedral fan in 2-d lattice N
    sage: fm.domain_fan().rays()
    N(-1,  1),
    N( 1,  1),
    N(-1, -1),
    N( 1, -1),
    N( 0, -1),
    N( 0,  1)
    in 2-d lattice N
    sage: normal.rays()
    N(-1,  1),
    N( 1,  1),
    N(-1, -1),
    N( 1, -1)
    in 2-d lattice N

As you see, it was necessary to insert two new rays (to prevent "upper" and
"lower" cones of the normal fan from being mapped to the whole `x`-axis).
"""


#*****************************************************************************
#       Copyright (C) 2010 Andrey Novoseltsev <[email protected]>
#       Copyright (C) 2010 William Stein <[email protected]>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#
#                  http://www.gnu.org/licenses/
#*****************************************************************************


import operator

from sage.categories.all import Hom
from sage.geometry.cone import Cone
from sage.geometry.fan import Fan, is_Fan, discard_faces
from sage.matrix.all import identity_matrix, matrix
from sage.matrix.matrix import is_Matrix
from sage.misc.all import cached_method, latex, prod, walltime
from sage.modules.free_module_morphism import (FreeModuleMorphism,
                                               is_FreeModuleMorphism)
from sage.rings.all import Infinity, ZZ
from sage.rings.infinity import is_Infinite

class FanMorphism(FreeModuleMorphism):
    r"""
    Create a fan morphism.

    Let `\Sigma_1` and `\Sigma_2` be two fans in lattices `N_1` and `N_2`
    respectively. Let `\phi` be a morphism (i.e. a linear map) from `N_1` to
    `N_2`. We say that `\phi` is *compatible* with `\Sigma_1` and `\Sigma_2`
    if every cone `\sigma_1\in\Sigma_1` is mapped by `\phi` into a single cone
    `\sigma_2\in\Sigma_2`, i.e. `\phi(\sigma_1)\subset\sigma_2` (`\sigma_2`
    may be different for different `\sigma_1`).

    By a **fan morphism** we understand a morphism between two lattices
    compatible with specified fans in these lattices. Such morphisms behave in
    exactly the same way as "regular" morphisms between lattices, but:

    * fan morphisms have a special constructor allowing some automatic
      adjustments to the initial fans (see below);
    * fan morphisms are aware of the associated fans and they can be
      accessed via :meth:`codomain_fan` and :meth:`domain_fan`;
    * fan morphisms can efficiently compute :meth:`image_cone` of a given
      cone of the domain fan and :meth:`preimage_cones` of a given cone of
      the codomain fan.

    INPUT:

    - ``morphism`` -- either a morphism between domain and codomain, or an
      integral matrix defining such a morphism;

    - ``domain_fan`` -- a :class:`fan
      <sage.geometry.fan.RationalPolyhedralFan>` in the domain;

    - ``codomain`` -- (default: ``None``) either a codomain lattice or a fan in
      the codomain. If the codomain fan is not given, the image fan (fan
      generated by images of generating cones) of ``domain_fan`` will be used,
      if possible;

    - ``subdivide`` -- (default: ``False``) if ``True`` and ``domain_fan`` is
      not compatible with the codomain fan because it is too coarse, it will be
      automatically refined to become compatible (the minimal refinement is
      canonical, so there are no choices involved);

    - ``check`` -- (default: ``True``) if ``False``, given fans and morphism
      will be assumed to be compatible. Be careful when using this option,
      since wrong assumptions can lead to wrong and hard-to-detect errors. On
      the other hand, this option may save you some time;

    - ``verbose`` -- (default: ``False``) if ``True``, some information may be
      printed during construction of the fan morphism.

    OUTPUT:

    - a fan morphism.

    EXAMPLES:

    Here we consider the face and normal fans of the "diamond" and the
    projection to the `x`-axis::

        sage: diamond = lattice_polytope.octahedron(2)
        sage: face = FaceFan(diamond)
        sage: normal = NormalFan(diamond)
        sage: N = face.lattice()
        sage: H = End(N)
        sage: phi = H([N.0, 0])
        sage: phi
        Free module morphism defined by the matrix
        [1 0]
        [0 0]
        Domain: 2-d lattice N
        Codomain: 2-d lattice N
        sage: fm = FanMorphism(phi, face, normal)
        sage: fm.domain_fan() is face
        True

    Note, that since ``phi`` is compatible with these fans, the returned
    fan is exactly the same object as the initial ``domain_fan``. ::

        sage: FanMorphism(phi, normal, face)
        Traceback (most recent call last):
        ...
        ValueError: the image of generating cone #1 of the domain fan
        is not contained in a single cone of the codomain fan!
        sage: fm = FanMorphism(phi, normal, face, subdivide=True)
        sage: fm.domain_fan() is normal
        False
        sage: fm.domain_fan().ngenerating_cones()
        6

    We had to subdivide two of the four cones of the normal fan, since
    they were mapped by ``phi`` into non-strictly convex cones.

    It is possible to omit the codomain fan, in which case the image fan will
    be used instead of it::

        sage: fm = FanMorphism(phi, face)
        sage: fm.codomain_fan()
        Rational polyhedral fan in 2-d lattice N
        sage: fm.codomain_fan().rays()
        N( 1, 0),
        N(-1, 0)
        in 2-d lattice N

    Now we demonstrate a more subtle example. We take the first quadrant as our
    domain fan. Then we divide the first quadrant into three cones, throw away
    the middle one and take the other two as our codomain fan. These fans are
    incompatible with the identity lattice morphism since the image of the
    domain fan is out of the support of the codomain fan::

        sage: N = ToricLattice(2)
        sage: phi = End(N).identity()
        sage: F1 = Fan(cones=[(0,1)], rays=[(1,0), (0,1)])
        sage: F2 = Fan(cones=[(0,1), (2,3)],
        ...            rays=[(1,0), (2,1), (1,2), (0,1)])
        sage: FanMorphism(phi, F1, F2)
        Traceback (most recent call last):
        ...
        ValueError: the image of generating cone #0 of the domain fan
        is not contained in a single cone of the codomain fan!
        sage: FanMorphism(phi, F1, F2, subdivide=True)
        Traceback (most recent call last):
        ...
        ValueError: morphism defined by
        [1 0]
        [0 1]
        does not map
        Rational polyhedral fan in 2-d lattice N
        into the support of
        Rational polyhedral fan in 2-d lattice N!

    The problem was detected and handled correctly (i.e. an exception was
    raised). However, the used algorithm requires extra checks for this
    situation after constructing a potential subdivision and this can take
    significant time. You can save about half the time using
    ``check=False`` option, if you know in advance that it is possible to
    make fans compatible with the morphism by subdividing the domain fan.
    Of course, if your assumption was incorrect, the result will be wrong
    and you will get a fan which *does* map into the support of the
    codomain fan, but is **not** a subdivision of the domain fan. You
    can test it on the example above::

        sage: fm = FanMorphism(phi, F1, F2, subdivide=True,
        ...                    check=False, verbose=True)
        Placing ray images (... ms)
        Computing chambers (... ms)
        Number of domain cones: 1.
        Number of chambers: 2.
        Cone 0 sits in chambers 0 1 (... ms)
        sage: fm.domain_fan().is_equivalent(F2)
        True
    """

    def __init__(self, morphism, domain_fan,
                 codomain=None,
                 subdivide=False,
                 check=True,
                 verbose=False):
        r"""
        Create a fan morphism.

        See :class:`FanMorphism` for documentation.

        TESTS::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, quadrant)
            sage: fm
            Fan morphism defined by the matrix
            [1 0]
            [0 1]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in 2-d lattice N
            sage: TestSuite(fm).run(skip="_test_category")
        """
        assert is_Fan(domain_fan)
        if is_Fan(codomain):
            codomain, codomain_fan = codomain.lattice(), codomain
        else:
            codomain_fan = None
        if is_FreeModuleMorphism(morphism):
            parent = morphism.parent()
            A = morphism.matrix()
        elif is_Matrix(morphism):
            A = morphism
            if codomain is None:
                raise ValueError("codomain (fan) must be given explicitly if "
                                 "morphism is given by a matrix!")
            parent = Hom(domain_fan.lattice(), codomain)
        else:
            raise TypeError("morphism must be either a FreeModuleMorphism "
                            "or a matrix!\nGot: %s" % morphism)
        super(FanMorphism, self).__init__(parent, A)
        self._domain_fan = domain_fan
        self._image_cone = dict()
        self._preimage_cones = dict()
        self._preimage_fans = dict()
        self._primitive_preimage_cones = dict()
        if codomain_fan is None:
            self._construct_codomain_fan(check)
        else:
            self._codomain_fan = codomain_fan
            if subdivide:
                self._subdivide_domain_fan(check, verbose)
            elif check:
                self._validate()

    def __mul__(self, right):
        """
        Return the composition of ``self`` and ``right``.

        INPUT:

        - ``right`` -- a :class:`FanMorphism` whose :meth:`codomain_fan` can be
          mapped to :meth:`domain_fan` of ``self`` via identity map of lattices.

        OUTPUT:

        - a :class:`FanMorphism`.

        EXAMPLES::

            sage: A2 = toric_varieties.A2()
            sage: P3 = toric_varieties.P(3)
            sage: m = matrix([(2,0,0), (1,1,0)])
            sage: phi = A2.hom(m, P3).fan_morphism()
            sage: phi
            Fan morphism defined by the matrix
            [2 0 0]
            [1 1 0]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in 3-d lattice N
            sage: prod(phi.factor()) # indirect test
            Fan morphism defined by the matrix
            [2 0 0]
            [1 1 0]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in 3-d lattice N
        """
        if not isinstance(right, FanMorphism):
            raise TypeError(
                "fan morphisms should be composed with fan morphisms")
        # We don't need it, we just check compatibility of fans:
        FanMorphism(identity_matrix(self.domain().dimension()),
                    right.codomain_fan(), self.domain_fan())
        m = right.matrix() * self.matrix()
        return FanMorphism(m, right.domain_fan(), self.codomain_fan())

    def _RISGIS(self):
        r"""
        Return Ray Images Star Generator Indices Sets.

        OUTPUT:

        - a :class:`tuple` of :class:`frozenset`s of integers, the `i`-th set
          is the set of indices of star generators for the minimal cone of the
          :meth:`codomain_fan` containing the image of the `i`-th ray of the
          :meth:`domain_fan`.

        TESTS::

            sage: diamond = lattice_polytope.octahedron(2)
            sage: face = FaceFan(diamond)
            sage: normal = NormalFan(diamond)
            sage: N = face.lattice()
            sage: fm = FanMorphism(identity_matrix(2),
            ...           normal, face, subdivide=True)
            sage: fm._RISGIS()
            (frozenset([3]), frozenset([2]),
             frozenset([1]), frozenset([0]),
             frozenset([1, 3]), frozenset([0, 1]),
             frozenset([0, 2]), frozenset([2, 3]))
        """
        if "_RISGIS_" not in self.__dict__:
            try:
                cones = [self._codomain_fan.cone_containing(self(ray))
                         for ray in self._domain_fan.rays()]
            except ValueError:
                self._support_error()
            self._RISGIS_ = tuple(frozenset(cone.star_generator_indices())
                                  for cone in cones)
        return self._RISGIS_

    def _chambers(self):
        r"""
        Return chambers in the domain corresponding to the codomain fan.

        This function is used during automatic refinement of the domain fans,
        see :meth:`_subdivide_domain_fan`.

        OUTPUT:

        - a :class:`tuple` ``(chambers, cone_to_chamber)``, where

        - ``chambers`` is a :class:`list` of :class:`cones
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` in the domain of
          ``self``;

        - ``cone_to_chamber`` is a :class:`list` of integers, if its `i`-th
          element is `j`, then the `j`-th element of ``chambers`` is the
          inverse image of the `i`-th generating cone of the codomain fan.

        TESTS::

            sage: F = NormalFan(lattice_polytope.octahedron(2))
            sage: N = F.lattice()
            sage: H = End(N)
            sage: phi = H([N.0, 0])
            sage: fm = FanMorphism(phi, F, F, subdivide=True)
            sage: fm._chambers()
            ([2-d cone in 2-d lattice N,
              1-d cone in 2-d lattice N,
              2-d cone in 2-d lattice N],
             [0, 1, 2, 1])
        """
        kernel_rays = []
        for ray in self.kernel().basis():
            kernel_rays.append(ray)
            kernel_rays.append(-ray)
        image_rays = []
        for ray in self.image().basis():
            image_rays.append(ray)
            image_rays.append(-ray)
        image = Cone(image_rays)
        chambers = []
        cone_to_chamber = []
        for cone in self._codomain_fan:
            chamber = Cone([self.lift(ray) for ray in cone.intersection(image)]
                           + kernel_rays, lattice=self.domain())
            cone_to_chamber.append(len(chambers))
            for i, old_chamber in enumerate(chambers):
                if old_chamber.is_equivalent(chamber):
                    cone_to_chamber[-1] = i
                    break
            if cone_to_chamber[-1] == len(chambers):
                chambers.append(chamber)
        return (chambers, cone_to_chamber)

    def _construct_codomain_fan(self, check):
        r"""
        Construct the codomain fan as the image of the domain one.

        .. WARNING::

            This method should be called only during initialization.

        INPUT:

        - ``check`` -- passed on to the fan constructor.

        OUTPUT:

        - none, but the codomain fan of ``self`` is set to the constructed fan.

        TESTS::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl)
            Traceback (most recent call last):
            ...
            ValueError: codomain (fan) must be given explicitly
            if morphism is given by a matrix!
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, ZZ^2)
            sage: fm.codomain_fan().rays()  # indirect doctest
            (1, 0),
            (0, 1),
            (1, 1)
            in Ambient free module of rank 2
            over the principal ideal domain Integer Ring
        """
        # We literally try to construct the image fan and hope that it works.
        # If it does not, the fan constructor will raise an exception.
        domain_fan = self._domain_fan
        self._codomain_fan = Fan(cones=(domain_cone.ambient_ray_indices()
                                        for domain_cone in domain_fan),
                                 rays=(self(ray) for ray in domain_fan.rays()),
                                 lattice=self.codomain(),
                                 discard_faces=True, check=check)

    def _latex_(self):
        r"""
        Return the `\LaTeX` representation of ``self``.

        OUTPUT:

        - a :class:`string`.

        EXAMPLES::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, quadrant)
            sage: fm
            Fan morphism defined by the matrix
            [1 0]
            [0 1]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in 2-d lattice N
            sage: print fm._latex_()
            \left(\begin{array}{rr}
            1 & 0 \\
            0 & 1
            \end{array}\right) : \Sigma^{2} \to \Sigma^{2}
        """
        return (r"%s : %s \to %s" % (latex(self.matrix()),
                        latex(self.domain_fan()), latex(self.codomain_fan())))

    @cached_method
    def _ray_index_map(self):
        r"""
        Return the map between indices of rays in domain and codomain fans.

        OUTPUT:

        - a tuple of integers. If the `i`-th entry is -1, the `i`-th ray of the
          domain fan is mapped to the origin. If it is `j`, then the `i`-th ray
          of the domain fan is mapped onto the `j`-th ray of the codomain fan.
          If there is a ray in the domain fan which is mapped into the relative
          interior of a higher dimensional cone, a ``ValueError`` exception is
          raised.

        .. NOTE::

            This method is used by :meth:`is_bundle` and :meth:`is_fibration`.

        TESTS::

            sage: Sigma = toric_varieties.dP8().fan()
            sage: Sigma_p = toric_varieties.P1().fan()
            sage: phi = FanMorphism(matrix([[1], [-1]]), Sigma, Sigma_p)
            sage: phi._ray_index_map()
            (-1, 1, -1, 0)
            sage: xi = FanMorphism(matrix([[1, 0]]), Sigma_p, Sigma)
            sage: xi._ray_index_map()
            Traceback (most recent call last):
            ...
            ValueError: ray #1 is mapped into a 2-d cone!
        """
        Sigma = self.domain_fan()
        ray_index_map = [-1] * Sigma.nrays()
        for i, rho in enumerate(Sigma(1)):
            sigma_p = self.image_cone(rho)
            if sigma_p.nrays() > 1:
                raise ValueError("ray #%d is mapped into a %d-d cone!" %
                                 (i, sigma_p.dim()))
            elif sigma_p.nrays() == 1:
                ray_index_map[i] = sigma_p.ambient_ray_indices()[0]
        return tuple(ray_index_map)

    def _repr_(self):
        r"""
        Return the string representation of ``self``.

        OUTPUT:

        - a :class:`string`.

        EXAMPLES::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, quadrant)
            sage: print fm._repr_()
            Fan morphism defined by the matrix
            [1 0]
            [0 1]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in 2-d lattice N
        """
        return ("Fan morphism defined by the matrix\n"
                "%s\n"
                "Domain fan: %s\n"
                "Codomain fan: %s"
                % (self.matrix(), self.domain_fan(), self.codomain_fan()))

    def _subdivide_domain_fan(self, check, verbose):
        r"""
        Subdivide the domain fan to make it compatible with the codomain fan.

        .. WARNING::

            This method should be called only during initialization.

        INPUT:

        - ``check`` -- (default: ``True``) if ``False``, some of the
          consistency checks will be omitted, which saves time but can
          potentially lead to wrong results. Currently, with
          ``check=False`` option there will be no check that ``domain_fan``
          maps to ``codomain_fan``;

        - ``verbose`` -- (default: ``False``) if ``True``, some timing
          information will be printed in the process.

        OUTPUT:

        - none, but the domain fan of self is replaced with its minimal
          refinement, if possible. Otherwise a ``ValueError`` exception is
          raised.

        TESTS::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant, quadrant_bl)
            Traceback (most recent call last):
            ...
            ValueError: the image of generating cone #0 of the domain fan
            is not contained in a single cone of the codomain fan!
            sage: fm = FanMorphism(identity_matrix(2), quadrant,
            ...                    quadrant_bl, subdivide=True)
            sage: fm.domain_fan().rays()  # indirect doctest
            N(1, 0),
            N(0, 1),
            N(1, 1)
            in 2-d lattice N

        Now we demonstrate a more subtle example. We take the first quadrant
        as our ``domain_fan``. Then we divide the first quadrant into three
        cones, throw away the middle one and take the other two as our
        ``codomain_fan``. These fans are incompatible with the identity
        lattice morphism since the image of ``domain_fan`` is out of the
        support of ``codomain_fan``::

            sage: N = ToricLattice(2)
            sage: phi = End(N).identity()
            sage: F1 = Fan(cones=[(0,1)], rays=[(1,0), (0,1)])
            sage: F2 = Fan(cones=[(0,1), (2,3)],
            ...            rays=[(1,0), (2,1), (1,2), (0,1)])
            sage: FanMorphism(phi, F1, F2)
            Traceback (most recent call last):
            ...
            ValueError: the image of generating cone #0 of the domain fan
            is not contained in a single cone of the codomain fan!
            sage: FanMorphism(phi, F1, F2, subdivide=True)
            Traceback (most recent call last):
            ...
            ValueError: morphism defined by
            [1 0]
            [0 1]
            does not map
            Rational polyhedral fan in 2-d lattice N
            into the support of
            Rational polyhedral fan in 2-d lattice N!

        We check that Trac #10943 is fixed::

            sage: Sigma = Fan(rays=[(1,1,0), (1,-1,0)], cones=[(0,1)])
            sage: Sigma_prime = FaceFan(lattice_polytope.octahedron(3))
            sage: fm = FanMorphism(identity_matrix(3),
            ...                    Sigma, Sigma_prime, subdivide=True)
            sage: fm.domain_fan().rays()
            N(1,  1, 0),
            N(1, -1, 0),
            N(1,  0, 0)
            in 3-d lattice N
            sage: [cone.ambient_ray_indices() for cone in fm.domain_fan()]
            [(0, 2), (1, 2)]

            sage: sigma = Cone([(0,1), (3,1)])
            sage: Sigma = Fan([sigma])
            sage: Sigma_prime = Sigma.subdivide([(1,1), (2,1)])
            sage: FanMorphism(identity_matrix(2),
            ...               Sigma, Sigma_prime, subdivide=True)
            Fan morphism defined by the matrix
            [1 0]
            [0 1]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in 2-d lattice N
        """
        domain_fan = self._domain_fan
        codomain_fan = self._codomain_fan
        lattice_dim = self.domain().dimension()
        if verbose:
            start = walltime()
            print "Placing ray images",
        # Figure out where 1-dimensional cones (i.e. rays) are mapped.
        RISGIS = self._RISGIS()
        if verbose:
            print "(%.3f ms)" % walltime(start)
        # Subdivide cones that require it.
        chambers = None # preimages of codomain cones, computed if necessary
        new_cones = []
        for cone_index, domain_cone in enumerate(domain_fan):
            if reduce(operator.and_, (RISGIS[i]
                                for i in domain_cone.ambient_ray_indices())):
                # There is a codomain cone containing all rays of this domain
                # cone, no need to subdivide it.
                new_cones.append(domain_cone)
                continue
            dim = domain_cone.dim()
            if chambers is None:
                if verbose:
                    start = walltime()
                    print "Computing chambers",
                chambers, cone_to_chamber = self._chambers()
                if verbose:
                    print "(%.3f ms)" % walltime(start)
                    print ("Number of domain cones: %d.\n"
                           "Number of chambers: %d." %
                           (domain_fan.ngenerating_cones(), len(chambers)))
            # Subdivide domain_cone.
            if verbose:
                start = walltime()
                print "Cone %d sits in chambers" % cone_index,
            # It seems that there is no quick way to determine which chambers
            # do NOT intersect this domain cone, so we go over all of them.
            parts = []
            for chamber_index, chamber in enumerate(chambers):
                # If chamber is small, intersection will be small.
                if chamber.dim() < dim:
                    continue
                new_part = domain_cone.intersection(chamber)
                # If intersection is small, it is not a generating cone.
                if new_part.dim() < dim:
                    continue
                # Small cones may have repetitive intersections with chambers.
                if (dim == lattice_dim or
                    not any(part.is_equivalent(new_part) for part in parts)):
                    parts.append(new_part)
                    if verbose:
                        print chamber_index,
            if check:
                # Check if the subdivision is complete, i.e. there are no
                # missing pieces of domain_cone. To do this, we construct a
                # fan from the obtained parts and check that interior points
                # of boundary cones of this fan are in the interior of the
                # original cone. In any case we know that we are constructing
                # a valid fan, so passing check=False to Fan(...) is OK.
                if verbose:
                    print "(%.3f ms)" % walltime(start)
                    start = walltime()
                    print "Checking for missing pieces",
                cone_subdivision = Fan(parts, check=False)
                for cone in cone_subdivision(dim - 1):
                    if len(cone.star_generators()) == 1:
                        if domain_cone.relative_interior_contains(
                                                            sum(cone.rays())):
                            self._support_error()
            new_cones.extend(parts)
            if verbose:
                print "(%.3f ms)" % walltime(start)
        if len(new_cones) > domain_fan.ngenerating_cones():
            # Construct a new fan keeping old rays in the same order
            new_rays = list(domain_fan.rays())
            for cone in new_cones:
                for ray in cone:
                    if ray not in new_rays:
                        new_rays.append(ray)
            # Replace domain_fan, this is OK since this method must be called
            # only during initialization of the FanMorphism.
            self._domain_fan = Fan(new_cones, new_rays, domain_fan.lattice(),
                                   check=False)
            # Also remove RISGIS for the old fan
            del self._RISGIS_

    def _support_error(self):
        r"""
        Raise a ``ValueError`` exception due to support incompatibility.

        OUTPUT:

        - none, a ``ValueError`` exception is raised.

        TESTS:

        We deliberately construct an invalid morphism for this demonstration::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2),
            ...             quadrant, quadrant_bl, check=False)

        Now we report that the morphism is invalid::

            sage: fm._support_error()
            Traceback (most recent call last):
            ...
            ValueError: morphism defined by
            [1 0]
            [0 1]
            does not map
            Rational polyhedral fan in 2-d lattice N
            into the support of
            Rational polyhedral fan in 2-d lattice N!
        """
        raise ValueError("morphism defined by\n"
                         "%s\n"
                         "does not map\n"
                         "%s\n"
                         "into the support of\n"
                         "%s!"
                    % (self.matrix(), self.domain_fan(), self.codomain_fan()))

    def _validate(self):
        r"""
        Check if ``self`` is indeed compatible with domain and codomain fans.

        OUTPUT:

        - none, but a ``ValueError`` exception is raised if there is a cone of
          the domain fan of ``self`` which is not completely contained in a
          single cone of the codomain fan of ``self``, or if one of these fans
          does not sit in the appropriate lattice.

        EXAMPLES::

            sage: N3 = ToricLattice(3, "N3")
            sage: N2 = ToricLattice(2, "N2")
            sage: H = Hom(N3, N2)
            sage: phi = H([N2.0, N2.1, N2.0])
            sage: F1 = Fan(cones=[(0,1,2), (1,2,3)],
            ...            rays=[(1,1,1), (1,1,-1), (1,-1,1), (1,-1,-1)],
            ...            lattice=N3)
            sage: F1.rays()
            N3(1,  1,  1),
            N3(1,  1, -1),
            N3(1, -1,  1),
            N3(1, -1, -1)
            in 3-d lattice N3
            sage: [phi(ray) for ray in F1.rays()]
            [N2(2, 1), N2(0, 1), N2(2, -1), N2(0, -1)]
            sage: F2 = Fan(cones=[(0,1,2), (1,2,3)],
            ...            rays=[(1,1,1), (1,1,-1), (1,2,1), (1,2,-1)],
            ...            lattice=N3)
            sage: F2.rays()
            N3(1, 1,  1),
            N3(1, 1, -1),
            N3(1, 2,  1),
            N3(1, 2, -1)
            in 3-d lattice N3
            sage: [phi(ray) for ray in F2.rays()]
            [N2(2, 1), N2(0, 1), N2(2, 2), N2(0, 2)]
            sage: F3 = Fan(cones=[(0,1), (1,2)],
            ...            rays=[(1,0), (2,1), (0,1)],
            ...            lattice=N2)
            sage: FanMorphism(phi, F2, F3)
            Fan morphism defined by the matrix
            [1 0]
            [0 1]
            [1 0]
            Domain fan: Rational polyhedral fan in 3-d lattice N3
            Codomain fan: Rational polyhedral fan in 2-d lattice N2
            sage: FanMorphism(phi, F1, F3)  # indirect doctest
            Traceback (most recent call last):
            ...
            ValueError: morphism defined by
            [1 0]
            [0 1]
            [1 0]
            does not map
            Rational polyhedral fan in 3-d lattice N3
            into the support of
            Rational polyhedral fan in 2-d lattice N2!
        """
        domain_fan = self._domain_fan
        if domain_fan.lattice() is not self.domain():
            raise ValueError("%s does not sit in %s!"
                             % (domain_fan, self.domain()))
        codomain_fan = self._codomain_fan
        if codomain_fan.lattice() is not self.codomain():
            raise ValueError("%s does not sit in %s!"
                             % (codomain_fan, self.codomain()))
        RISGIS = self._RISGIS()
        for n, domain_cone in enumerate(domain_fan):
            if not reduce(operator.and_,
                (RISGIS[i] for i in domain_cone.ambient_ray_indices())):
                raise ValueError("the image of generating cone #%d of the "
                                 "domain fan is not contained in a single "
                                 "cone of the codomain fan!" % n)

    def codomain_fan(self, dim=None, codim=None):
        r"""
        Return the codomain fan of ``self``.

        INPUT:

        - ``dim`` -- dimension of the requested cones;

        - ``codim`` -- codimension of the requested cones.

        OUTPUT:

        - :class:`rational polyhedral fan
          <sage.geometry.fan.RationalPolyhedralFan>` if no parameters were
          given, :class:`tuple` of :class:`cones
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` otherwise.

        EXAMPLES::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, quadrant)
            sage: fm.codomain_fan()
            Rational polyhedral fan in 2-d lattice N
            sage: fm.codomain_fan() is quadrant
            True
        """
        return self._codomain_fan(dim=dim, codim=codim)

    def domain_fan(self, dim=None, codim=None):
        r"""
        Return the codomain fan of ``self``.

        INPUT:

        - ``dim`` -- dimension of the requested cones;

        - ``codim`` -- codimension of the requested cones.

        OUTPUT:

        - :class:`rational polyhedral fan
          <sage.geometry.fan.RationalPolyhedralFan>` if no parameters were
          given, :class:`tuple` of :class:`cones
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` otherwise.

        EXAMPLES::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, quadrant)
            sage: fm.domain_fan()
            Rational polyhedral fan in 2-d lattice N
            sage: fm.domain_fan() is quadrant_bl
            True
        """
        return self._domain_fan(dim=dim, codim=codim)

    def image_cone(self, cone):
        r"""
        Return the cone of the codomain fan containing the image of ``cone``.

        INPUT:

        - ``cone`` -- a :class:`cone
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` equivalent to a
          cone of the :meth:`domain_fan` of ``self``.

        OUTPUT:

        - a :class:`cone <sage.geometry.cone.ConvexRationalPolyhedralCone>`
          of the :meth:`codomain_fan` of ``self``.

        EXAMPLES::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, quadrant)
            sage: fm.image_cone(Cone([(1,0)]))
            1-d cone of Rational polyhedral fan in 2-d lattice N
            sage: fm.image_cone(Cone([(1,1)]))
            2-d cone of Rational polyhedral fan in 2-d lattice N

        TESTS:

        We check that complete codomain fans are handled correctly, since a
        different algorithm is used in this case::

            sage: diamond = lattice_polytope.octahedron(2)
            sage: face = FaceFan(diamond)
            sage: normal = NormalFan(diamond)
            sage: N = face.lattice()
            sage: fm = FanMorphism(identity_matrix(2),
            ...           normal, face, subdivide=True)
            sage: fm.image_cone(Cone([(1,0)]))
            1-d cone of Rational polyhedral fan in 2-d lattice N
            sage: fm.image_cone(Cone([(1,1)]))
            2-d cone of Rational polyhedral fan in 2-d lattice N
        """
        cone = self._domain_fan.embed(cone)
        if cone not in self._image_cone:
            codomain_fan = self._codomain_fan()
            if cone.is_trivial():
                self._image_cone[cone] = codomain_fan(0)[0]
            elif codomain_fan.is_complete():
                # Optimization for a common case
                RISGIS = self._RISGIS()
                CSGIS = set(reduce(operator.and_,
                            (RISGIS[i] for i in cone.ambient_ray_indices())))
                image_cone = codomain_fan.generating_cone(CSGIS.pop())
                for i in CSGIS:
                    image_cone = image_cone.intersection(
                                            codomain_fan.generating_cone(i))
                self._image_cone[cone] = image_cone
            else:
                self._image_cone[cone] = codomain_fan.cone_containing(
                                                    self(ray) for ray in cone)
        return self._image_cone[cone]

    def index(self, cone=None):
        r"""
        Return the index of ``self`` as a map between lattices.

        INPUT:

        - ``cone`` -- (default: ``None``) a :class:`cone
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` of the
          :meth:`codomain_fan` of ``self``.

        OUTPUT:

        - an integer, infinity, or ``None``.

        If no cone was specified, this function computes the index of the
        image of ``self`` in the codomain. If a cone `\sigma` was given, the
        index of ``self`` over `\sigma` is computed in the sense of
        Definition 2.1.7 of [HLY]: if `\sigma'` is any cone of the
        :meth:`domain_fan` of ``self`` whose relative interior is mapped to the
        relative interior of `\sigma`, it is the index of the image of
        `N'(\sigma')` in `N(\sigma)`, where `N'` and `N` are domain and codomain
        lattices respectively. While that definition was formulated for the case
        of the finite index only, we extend it to the infinite one as well and
        return ``None`` if there is no `\sigma'` at all. See examples below for
        situations when such things happen. Note also that the index of ``self``
        is the same as index over the trivial cone.

        EXAMPLES::

            sage: Sigma = toric_varieties.dP8().fan()
            sage: Sigma_p = toric_varieties.P1().fan()
            sage: phi = FanMorphism(matrix([[1], [-1]]), Sigma, Sigma_p)
            sage: phi.index()
            1
            sage: psi = FanMorphism(matrix([[2], [-2]]), Sigma, Sigma_p)
            sage: psi.index()
            2
            sage: xi = FanMorphism(matrix([[1, 0]]), Sigma_p, Sigma)
            sage: xi.index()
            +Infinity

        Infinite index in the last example indicates that the image has positive
        codimension in the codomain. Let's look at the rays of our fans::

            sage: Sigma_p.rays()
            N( 1),
            N(-1)
            in 1-d lattice N
            sage: Sigma.rays()
            N( 1,  1),
            N( 0,  1),
            N(-1, -1),
            N( 1,  0)
            in 2-d lattice N
            sage: xi.factor()[0].domain_fan().rays()
            N( 1, 0),
            N(-1, 0)
            in Sublattice <N(1, 0)>

        We see that one of the rays of the fan of ``P1`` is mapped to a ray,
        while the other one to the interior of some 2-d cone. Both rays
        correspond to single points on ``P1``, yet one is mapped to the
        distinguished point of a torus invariant curve of ``dP8`` (with the
        rest of this curve being uncovered) and the other to a fixed point
        of ``dP8`` (thus completely covering this torus orbit in ``dP8``).

        We should therefore expect the following behaviour: all indices over
        1-d cones are ``None``, except for one which is infinite, and all
        indices over 2-d cones are ``None``, except for one which is 1::

            sage: [xi.index(cone) for cone in Sigma(1)]
            [None, None, None, +Infinity]
            sage: [xi.index(cone) for cone in Sigma(2)]
            [None, 1, None, None]

        TESTS::

            sage: Sigma = toric_varieties.dP8().fan()
            sage: Sigma_p = toric_varieties.Cube_nonpolyhedral().fan()
            sage: m = matrix([[2,6,10], [7,11,13]])
            sage: zeta = FanMorphism(m, Sigma, Sigma_p, subdivide=True)
            sage: [zeta.index(cone) for cone in flatten(Sigma_p.cones())]
            [+Infinity, None, None, None, None, None, None, None, None, None,
             4, 4, None, 4, None, None, 2, None, 4, None, 4, 1, 1, 1, 1, 1, 1]
            sage: zeta = prod(zeta.factor()[1:])
            sage: Sigma_p = zeta.codomain_fan()
            sage: [zeta.index(cone) for cone in flatten(Sigma_p.cones())]
            [4, 4, 1, 4, 4, 4, 4, 1, 1, 1, 1, 1, 1]
            sage: zeta.index() == zeta.index(Sigma_p(0)[0])
            True
        """
        if cone is None:
            try:
                return self.matrix().image().index_in(self.codomain())
            except ArithmeticError:
                cone = Cone([], lattice=self.codomain())
        cone = self._codomain_fan.embed(cone)
        PPCs = self.primitive_preimage_cones(cone)
        if not PPCs:
            return None
        Q = cone.sublattice_quotient()
        S = Q.submodule([self(g)
                         for g in PPCs[0].sublattice_complement().gens()])
        i = prod((Q/S).invariants())
        return i if i > 0 else Infinity

    def is_birational(self):
        r"""
        Check if ``self`` is birational.

        OUTPUT:

        - ``True`` if ``self`` is birational, ``False`` otherwise.

        For fan morphisms this check is equivalent to ``self.index() == 1`` and
        means that the corresponding map between toric varieties is birational.

        EXAMPLES::

            sage: Sigma = toric_varieties.dP8().fan()
            sage: Sigma_p = toric_varieties.P1().fan()
            sage: phi = FanMorphism(matrix([[1], [-1]]), Sigma, Sigma_p)
            sage: psi = FanMorphism(matrix([[2], [-2]]), Sigma, Sigma_p)
            sage: xi = FanMorphism(matrix([[1, 0]]), Sigma_p, Sigma)
            sage: phi.index(), psi.index(), xi.index()
            (1, 2, +Infinity)
            sage: phi.is_birational(), psi.is_birational(), xi.is_birational()
            (True, False, False)
        """
        return self.index() == 1

    @cached_method
    def is_bundle(self):
        r"""
        Check if ``self`` is a bundle.

        OUTPUT:

        - ``True`` if ``self`` is a bundle, ``False`` otherwise.

        Let `\phi: \Sigma \to \Sigma'` be a fan morphism such that the
        underlying lattice morphism `\phi: N \to N'` is surjective. Let
        `\Sigma_0` be the kernel fan of `\phi`. Then `\phi` is a **bundle** (or
        splitting) if there is a subfan `\widehat{\Sigma}` of `\Sigma` such
        that the following two conditions are satisfied:

            #. Cones of `\Sigma` are precisely the cones of the form
               `\sigma_0 + \widehat{\sigma}`, where `\sigma_0 \in \Sigma_0` and
               `\widehat{\sigma} \in \widehat{\Sigma}`.

            #. Cones of `\widehat{\Sigma}` are in bijection with cones of
               `\Sigma'` induced by `\phi` and `\phi` maps lattice points in
               every cone `\widehat{\sigma} \in \widehat{\Sigma}` bijectively
               onto lattice points in `\phi(\widehat{\sigma})`.

        If a fan morphism `\phi: \Sigma \to \Sigma'` is a bundle, then
        `X_\Sigma` is a fiber bundle over `X_{\Sigma'}` with fibers
        `X_{\Sigma_0, N_0}`, where `N_0` is the kernel lattice of `\phi`. See
        [CLS11]_ for more details.

        .. seealso:: :meth:`is_fibration`, :meth:`kernel_fan`.

        REFERENCES:

        ..  [CLS11]
            David A. Cox, John Little, and Hal Schenck.
            *Toric Varieties*.
            Volume 124 of *Graduate Studies in Mathematics*.
            American Mathematical Society, Providence, RI, 2011.

        EXAMPLES:

        We consider several maps between fans of a del Pezzo surface and the
        projective line::

            sage: Sigma = toric_varieties.dP8().fan()
            sage: Sigma_p = toric_varieties.P1().fan()
            sage: phi = FanMorphism(matrix([[1], [-1]]), Sigma, Sigma_p)
            sage: psi = FanMorphism(matrix([[2], [-2]]), Sigma, Sigma_p)
            sage: xi = FanMorphism(matrix([[1, 0]]), Sigma_p, Sigma)
            sage: phi.is_bundle()
            True
            sage: phi.is_fibration()
            True
            sage: phi.index()
            1
            sage: psi.is_bundle()
            False
            sage: psi.is_fibration()
            True
            sage: psi.index()
            2
            sage: xi.is_fibration()
            False
            sage: xi.index()
            +Infinity

        The first of these maps induces not only a fibration, but a fiber
        bundle structure. The second map is very similar, yet it fails to be
        a bundle, as its index is 2. The last map is not even a fibration.
        """
        if self.index() != 1:
            return False    # Not surjective between lattices.
        Sigma = self.domain_fan()
        Sigma_p = self.codomain_fan()
        Sigma_0 = self.kernel_fan()
        if (Sigma.ngenerating_cones() !=
            Sigma_0.ngenerating_cones() * Sigma_p.ngenerating_cones()):
            return False    # Definitely no splitting.
        try:
            ray_index_map = self._ray_index_map()
        except ValueError:
            return False  # Rays are not mapped onto rays or the origin.
        # Figure out how Sigma_0 sits inside Sigma in terms of ray indices.
        I_0s = [Sigma.embed(sigma_0).ambient_ray_indices()
                for sigma_0 in Sigma_0]
        # We examine only generating cones, this is sufficient.
        for sigma_p in Sigma_p:
            primitive_cones = self.primitive_preimage_cones(sigma_p)
            if len(primitive_cones) != 1:   # Should be only sigma_hat.
                return False
            sigma_hat = primitive_cones[0]
            if sigma_p.dim() != sigma_hat.dim():
                return False    # sigma -> sigma_p is not a bijection
            I_p = sigma_p.ambient_ray_indices()
            I_hat = sigma_hat.ambient_ray_indices()
            if I_p != tuple(sorted(ray_index_map[i] for i in I_hat)):
                return False    # sigma -> sigma_p is not a bijection
            # Check that sigma_hat + sigma_0 is always in Sigma.
            for I_0 in I_0s:
                I = tuple(sorted(I_hat + I_0))
                if all(sigma.ambient_ray_indices() != I for sigma in Sigma):
                    return False
        return True

    @cached_method
    def is_fibration(self):
        r"""
        Check if ``self`` is a fibration.

        OUTPUT:

        - ``True`` if ``self`` is a fibration, ``False`` otherwise.

        A fan morphism `\phi: \Sigma \to \Sigma'` is a **fibration** if for any
        cone `\sigma' \in \Sigma'` and any primitive preimage cone `\sigma \in
        \Sigma` corresponding to `\sigma'` the linear map of vector spaces
        `\phi_\RR` induces a bijection between `\sigma` and `\sigma'`, and, in
        addition, `\phi_\RR: N_\RR \to N'_\RR$ is surjective.

        If a fan morphism `\phi: \Sigma \to \Sigma'` is a fibration, then the
        associated morphism between toric varieties `\tilde{\phi}: X_\Sigma \to
        X_{\Sigma'}` is a fibration in the sense that it is surjective and all
        of its fibers have the same dimension, namely `\dim X_\Sigma -
        \dim X_{\Sigma'}`. These fibers do *not* have to be isomorphic, i.e. a
        fibration is not necessarily a fiber bundle. See [HLY02]_ for more
        details.

        .. seealso:: :meth:`is_bundle`, :meth:`primitive_preimage_cones`.

        REFERENCES:

        ..  [HLY02]
            Yi Hu, Chien-Hao Liu, and Shing-Tung Yau.
            Toric morphisms and fibrations of toric Calabi-Yau hypersurfaces.
            *Adv. Theor. Math. Phys.*, 6(3):457-506, 2002.
            arXiv:math/0010082v2 [math.AG].

        EXAMPLES:

        We consider several maps between fans of a del Pezzo surface and the
        projective line::

            sage: Sigma = toric_varieties.dP8().fan()
            sage: Sigma_p = toric_varieties.P1().fan()
            sage: phi = FanMorphism(matrix([[1], [-1]]), Sigma, Sigma_p)
            sage: psi = FanMorphism(matrix([[2], [-2]]), Sigma, Sigma_p)
            sage: xi = FanMorphism(matrix([[1, 0]]), Sigma_p, Sigma)
            sage: phi.is_bundle()
            True
            sage: phi.is_fibration()
            True
            sage: phi.index()
            1
            sage: psi.is_bundle()
            False
            sage: psi.is_fibration()
            True
            sage: psi.index()
            2
            sage: xi.is_fibration()
            False
            sage: xi.index()
            +Infinity

        The first of these maps induces not only a fibration, but a fiber
        bundle structure. The second map is very similar, yet it fails to be
        a bundle, as its index is 2. The last map is not even a fibration.

        TESTS:

        We check that reviewer's example on Trac 11200 works as expected::

            sage: P1 = toric_varieties.P1()
            sage: A1 = toric_varieties.A1()
            sage: phi = FanMorphism(matrix([[1]]), A1.fan(), P1.fan())
            sage: phi.is_fibration()
            False
        """
        if not self.is_surjective():
            return False
        try:
            ray_index_map = self._ray_index_map()
        except ValueError:
            return False    # Rays are not mapped onto rays or the origin.
        Sigma_p = self.codomain_fan()
        # Rays are already checked, the origin is trivial, start with 2-cones.
        for d in range(2, Sigma_p.dim() + 1):
            for sigma_p in Sigma_p(d):
                I_p = sigma_p.ambient_ray_indices()
                for sigma in self.primitive_preimage_cones(sigma_p):
                    if sigma.dim() != d:
                        return False    # sigma -> sigma_p is not a bijection
                    I = sigma.ambient_ray_indices()
                    if I_p != tuple(sorted(ray_index_map[i] for i in I)):
                        return False    # sigma -> sigma_p is not a bijection
        return True

    @cached_method
    def is_injective(self):
        r"""
        Check if ``self`` is injective.

        OUTPUT:

        - ``True`` if ``self`` is injective, ``False`` otherwise.

        Let `\phi: \Sigma \to \Sigma'` be a fan morphism such that the
        underlying lattice morphism `\phi: N \to N'` bijectively maps `N` to a
        *saturated* sublattice of `N'`. Let `\psi: \Sigma \to \Sigma'_0` be the
        restriction of `\phi` to the image. Then `\phi` is **injective** if the
        map between cones corresponding to `\psi` (injectively) maps each cone
        of `\Sigma` to a cone of the same dimension.

        If a fan morphism `\phi: \Sigma \to \Sigma'` is injective, then the
        associated morphism between toric varieties `\tilde{\phi}: X_\Sigma \to
        X_{\Sigma'}` is injective.

        .. seealso:: :meth:`restrict_to_image`.

        EXAMPLES:

        Consider the fan of the affine plane::

            sage: A2 = toric_varieties.A(2).fan()

        We will map several fans consisting of a single ray into the interior
        of the 2-cone::

            sage: Sigma = Fan([Cone([(1,1)])])
            sage: m = identity_matrix(2)
            sage: FanMorphism(m, Sigma, A2).is_injective()
            False

        This morphism was not injective since (in the toric varieties
        interpretation) the 1-dimensional orbit corresponding to the ray was
        mapped to the 0-dimensional orbit corresponding to the 2-cone. ::

            sage: Sigma = Fan([Cone([(1,)])])
            sage: m = matrix(1, 2, [1,1])
            sage: FanMorphism(m, Sigma, A2).is_injective()
            True

        While the fans in this example are close to the previous one, here the
        ray corresponds to a 0-dimensional orbit. ::

            sage: Sigma = Fan([Cone([(1,)])])
            sage: m = matrix(1, 2, [2,2])
            sage: FanMorphism(m, Sigma, A2).is_injective()
            False

        Here the problem is that ``m`` maps the domain lattice to a
        non-saturated sublattice of the codomain. The corresponding map of the
        toric varieties is a two-sheeted cover of its image.

        We also embed the affine plane into the projective one::

            sage: P2 = toric_varieties.P(2).fan()
            sage: m = identity_matrix(2)
            sage: FanMorphism(m, A2, P2).is_injective()
            True
        """
        if self.matrix().index_in_saturation() != 1:
            return False
        if self.image().dimension() < self.codomain().dimension():
            return prod(self.factor()[1:]).is_injective()
        # Now we know that underlying lattice morphism is bijective.
        Sigma = self.domain_fan()
        return all(all(self.image_cone(sigma).dim() == d for sigma in Sigma(d))
                   for d in range(1, Sigma.dim() + 1))

    @cached_method
    def is_surjective(self):
        r"""
        Check if ``self`` is surjective.

        OUTPUT:

        - ``True`` if ``self`` is surjective, ``False`` otherwise.

        A fan morphism `\phi: \Sigma \to \Sigma'` is **surjective** if the
        corresponding map between cones is surjective, i.e. for each cone
        `\sigma' \in \Sigma'` there is at least one preimage cone `\sigma \in
        \Sigma` such that the relative interior of `\sigma` is mapped to the
        relative interior of `\sigma'` and, in addition,
        `\phi_\RR: N_\RR \to N'_\RR$ is surjective.

        If a fan morphism `\phi: \Sigma \to \Sigma'` is surjective, then the
        associated morphism between toric varieties `\tilde{\phi}: X_\Sigma \to
        X_{\Sigma'}` is surjective.

        .. seealso:: :meth:`is_bundle`, :meth:`is_fibration`,
            :meth:`preimage_cones`.

        EXAMPLES:

        We check that the blow up of the affine plane at the origin is
        surjective::

            sage: A2 = toric_varieties.A(2).fan()
            sage: Bl = A2.subdivide([(1,1)])
            sage: m = identity_matrix(2)
            sage: FanMorphism(m, Bl, A2).is_surjective()
            True

        It remains surjective if we throw away "south and north poles" of the
        exceptional divisor::

            sage: FanMorphism(m, Fan(Bl.cones(1)), A2).is_surjective()
            True

        But a single patch of the blow up does not cover the plane::

            sage: F = Fan([Bl.generating_cone(0)])
            sage: FanMorphism(m, F, A2).is_surjective()
            False

        TESTS:

        We check that reviewer's example on Trac 11200 works as expected::

            sage: P1 = toric_varieties.P1()
            sage: A1 = toric_varieties.A1()
            sage: phi = FanMorphism(matrix([[1]]), A1.fan(), P1.fan())
            sage: phi.is_surjective()
            False
        """
        if is_Infinite(self.index()):
            return False    # Not surjective between vector spaces.
        for dcones in self.codomain_fan().cones():
            for sigma_p in dcones:
                if not self.preimage_cones(sigma_p):
                    return False
        return True

    @cached_method
    def kernel_fan(self):
        r"""
        Return the subfan of the domain fan mapped into the origin.

        OUTPUT:

        - a :class:`fan <sage.geometry.fan.RationalPolyhedralFan>`.

        .. NOTE::

            The lattice of the kernel fan is the :meth:`kernel` sublattice of
            ``self``.

        .. seealso:: :meth:`preimage_fan`.

        EXAMPLES::

            sage: fan = Fan(rays=[(1,0), (1,1), (0,1)], cones=[(0,1), (1,2)])
            sage: fm = FanMorphism(matrix(2, 1, [1,-1]), fan, ToricLattice(1))
            sage: fm.kernel_fan()
            Rational polyhedral fan in Sublattice <N(1, 1)>
            sage: _.rays()
            N(1, 1)
            in Sublattice <N(1, 1)>
            sage: fm.kernel_fan().cones()
            ((0-d cone of Rational polyhedral fan in Sublattice <N(1, 1)>,),
             (1-d cone of Rational polyhedral fan in Sublattice <N(1, 1)>,))
        """
        fan = self.preimage_fan(Cone([], lattice=self.codomain()))
        return Fan((cone.ambient_ray_indices() for cone in fan), fan.rays(),
                   lattice=self.kernel(), check=False)

    def preimage_cones(self, cone):
        r"""
        Return cones of the domain fan whose :meth:`image_cone` is ``cone``.

        INPUT:

        - ``cone`` -- a :class:`cone
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` equivalent to a
          cone of the :meth:`codomain_fan` of ``self``.

        OUTPUT:

        - a :class:`tuple` of :class:`cones
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` of the
          :meth:`domain_fan` of ``self``, sorted by dimension.

        .. seealso:: :meth:`preimage_fan`.

        EXAMPLES::

            sage: quadrant = Cone([(1,0), (0,1)])
            sage: quadrant = Fan([quadrant])
            sage: quadrant_bl = quadrant.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2), quadrant_bl, quadrant)
            sage: fm.preimage_cones(Cone([(1,0)]))
            (1-d cone of Rational polyhedral fan in 2-d lattice N,)
            sage: fm.preimage_cones(Cone([(1,0), (0,1)]))
            (1-d cone of Rational polyhedral fan in 2-d lattice N,
             2-d cone of Rational polyhedral fan in 2-d lattice N,
             2-d cone of Rational polyhedral fan in 2-d lattice N)

        TESTS:

        We check that reviewer's example from Trac #9972 is handled correctly::

            sage: N1 = ToricLattice(1)
            sage: N2 = ToricLattice(2)
            sage: Hom21 = Hom(N2, N1)
            sage: pr = Hom21([N1.0,0])
            sage: P1xP1 = toric_varieties.P1xP1()
            sage: f = FanMorphism(pr, P1xP1.fan())
            sage: c = f.image_cone(Cone([(1,0), (0,1)]))
            sage: c
            1-d cone of Rational polyhedral fan in 1-d lattice N
            sage: f.preimage_cones(c)
            (1-d cone of Rational polyhedral fan in 2-d lattice N,
             2-d cone of Rational polyhedral fan in 2-d lattice N,
             2-d cone of Rational polyhedral fan in 2-d lattice N)
        """
        cone = self._codomain_fan.embed(cone)
        if cone not in self._preimage_cones:
            # "Images of preimages" must sit in all of these generating cones
            CSGI = cone.star_generator_indices()
            RISGIS = self._RISGIS()
            domain_fan = self._domain_fan
            possible_rays = frozenset(i for i in range(domain_fan.nrays())
                                        if RISGIS[i].issuperset(CSGI))
            preimage_cones = []
            for dcones in domain_fan.cones():
                for dcone in dcones:
                    if (possible_rays.issuperset(dcone.ambient_ray_indices())
                        and self.image_cone(dcone) == cone):
                        preimage_cones.append(dcone)
            self._preimage_cones[cone] = tuple(preimage_cones)
        return self._preimage_cones[cone]

    def preimage_fan(self, cone):
        r"""
        Return the subfan of the domain fan mapped into ``cone``.

        INPUT:

        - ``cone`` -- a :class:`cone
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` equivalent to a
          cone of the :meth:`codomain_fan` of ``self``.

        OUTPUT:

        - a :class:`fan <sage.geometry.fan.RationalPolyhedralFan>`.

        .. NOTE::

            The preimage fan of ``cone`` consists of all cones of the
            :meth:`domain_fan` which are mapped into ``cone``, including those
            that are mapped into its boundary. So this fan is not necessarily
            generated by :meth:`preimage_cones` of ``cone``.

        .. seealso:: :meth:`kernel_fan`, :meth:`preimage_cones`.

        EXAMPLES::

            sage: quadrant_cone = Cone([(1,0), (0,1)])
            sage: quadrant_fan = Fan([quadrant_cone])
            sage: quadrant_bl = quadrant_fan.subdivide([(1,1)])
            sage: fm = FanMorphism(identity_matrix(2),
            ...                    quadrant_bl, quadrant_fan)
            sage: fm.preimage_fan(Cone([(1,0)])).cones()
            ((0-d cone of Rational polyhedral fan in 2-d lattice N,),
             (1-d cone of Rational polyhedral fan in 2-d lattice N,))
            sage: fm.preimage_fan(quadrant_cone).ngenerating_cones()
            2
            sage: len(fm.preimage_cones(quadrant_cone))
            3
        """
        cone = self._codomain_fan.embed(cone)
        if cone not in self._preimage_fans:
            # Find all cones mapped into the given one, including its boundary.
            domain_fan = self._domain_fan
            cones = []
            for dcones in reversed(domain_fan.cones()):
                for dcone in dcones:
                    if (not any(dcone.is_face_of(other) for other in cones) and
                        self.image_cone(dcone).is_face_of(cone)):
                        cones.append(dcone)
            # Now form the fan from these cones, keeping the ray order.
            ray_indices = set(cones[0].ambient_ray_indices())
            for c in cones[1:]:
                ray_indices.update(c.ambient_ray_indices())
            self._preimage_fans[cone] = Fan(cones,
                            domain_fan.rays(sorted(ray_indices)), check=False)
        return self._preimage_fans[cone]

    def primitive_preimage_cones(self, cone):
        r"""
        Return the primitive cones of the domain fan corresponding to ``cone``.

        INPUT:

        - ``cone`` -- a :class:`cone
          <sage.geometry.cone.ConvexRationalPolyhedralCone>` equivalent to a
          cone of the :meth:`codomain_fan` of ``self``.

        OUTPUT:

        - a :class:`cone <sage.geometry.cone.ConvexRationalPolyhedralCone>`.

        Let `\phi: \Sigma \to \Sigma'` be a fan morphism, let `\sigma \in
        \Sigma`, and let `\sigma' = \phi(\sigma)`. Then `\sigma` is a
        **primitive cone corresponding to** `\sigma'` if there is no proper
        face `\tau` of `\sigma` such that `\phi(\tau) = \sigma'`.

        Primitive cones play an important role for fibration morphisms.

        .. seealso:: :meth:`is_fibration`, :meth:`preimage_cones`,
            :meth:`preimage_fan`.

        EXAMPLES:

        Consider a projection of a del Pezzo surface onto the projective line::

            sage: Sigma = toric_varieties.dP6().fan()
            sage: Sigma.rays()
            N( 0,  1),
            N(-1,  0),
            N(-1, -1),
            N( 0, -1),
            N( 1,  0),
            N( 1,  1)
            in 2-d lattice N
            sage: Sigma_p = toric_varieties.P1().fan()
            sage: phi = FanMorphism(matrix([[1], [-1]]), Sigma, Sigma_p)

        Under this map, one pair of rays is mapped to the origin, one in the
        positive direction, and one in the negative one. Also three
        2-dimensional cones are mapped in the positive direction and three in
        the negative one, so there are 5 preimage cones corresponding to either
        of the rays of the codomain fan ``Sigma_p``::

            sage: len(phi.preimage_cones(Cone([(1,)])))
            5

        Yet only rays are primitive::

            sage: phi.primitive_preimage_cones(Cone([(1,)]))
            (1-d cone of Rational polyhedral fan in 2-d lattice N,
             1-d cone of Rational polyhedral fan in 2-d lattice N)

        Since all primitive cones are mapped onto their images bijectively, we
        get a fibration::

            sage: phi.is_fibration()
            True

        But since there are several primitive cones corresponding to the same
        cone of the codomain fan, this map is not a bundle, even though its
        index is 1::

            sage: phi.is_bundle()
            False
            sage: phi.index()
            1
        """
        sigma_p = self._codomain_fan.embed(cone) # Necessary if used as a key
        if sigma_p not in self._primitive_preimage_cones:
            primitive_cones = []
            for sigma in self.preimage_cones(sigma_p):  # Sorted by dimension
                if not any(tau.is_face_of(sigma) for tau in primitive_cones):
                    primitive_cones.append(sigma)
            self._primitive_preimage_cones[sigma_p] = tuple(primitive_cones)
        return self._primitive_preimage_cones[sigma_p]

    def factor(self):
        r"""
        Factor ``self`` into injective * birational * surjective morphisms.

        OUTPUT:

        - a triple of :class:`FanMorphism` `(\phi_i, \phi_b, \phi_s)`, such that
          `\phi_s` is surjective, `\phi_b` is birational, `\phi_i` is injective,
          and ``self`` is equal to `\phi_i \circ \phi_b \circ \phi_s`.

        Intermediate fans live in the saturation of the image of ``self``
        as a map between lattices and are the image of the :meth:`domain_fan`
        and the restriction of the :meth:`codomain_fan`, i.e. if ``self`` maps
        `\Sigma \to \Sigma'`, then we have factorization into

        .. math::

            \Sigma
            \twoheadrightarrow
            \Sigma_s
            \to
            \Sigma_i
            \hookrightarrow
            \Sigma.

        .. note::

            * `\Sigma_s` is the finest fan with the smallest support that is
              compatible with ``self``: any fan morphism from `\Sigma` given by
              the same map of lattices as ``self`` factors through `\Sigma_s`.

            * `\Sigma_i` is the coarsest fan of the largest support that is
              compatible with ``self``: any fan morphism into `\Sigma'` given by
              the same map of lattices as ``self`` factors though `\Sigma_i`.

        EXAMPLES:

        We map an affine plane into a projective 3-space in such a way, that it
        becomes "a double cover of a chart of the blow up of one of the
        coordinate planes"::

            sage: A2 = toric_varieties.A2()
            sage: P3 = toric_varieties.P(3)
            sage: m = matrix([(2,0,0), (1,1,0)])
            sage: phi = A2.hom(m, P3)
            sage: phi.as_polynomial_map()
            Scheme morphism:
              From: 2-d affine toric variety
              To:   3-d CPR-Fano toric variety covered by 4 affine patches
              Defn: Defined on coordinates by sending [x : y] to
                    [x^2*y : y : 1 : 1]

        Now we will work with the underlying fan morphism::

            sage: phi = phi.fan_morphism()
            sage: phi
            Fan morphism defined by the matrix
            [2 0 0]
            [1 1 0]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in 3-d lattice N
            sage: phi.is_surjective(), phi.is_birational(), phi.is_injective()
            (False, False, False)
            sage: phi_i, phi_b, phi_s = phi.factor()
            sage: phi_s.is_surjective(), phi_b.is_birational(), phi_i.is_injective()
            (True, True, True)
            sage: prod(phi.factor()) == phi
            True

        Double cover (surjective)::

            sage: A2.fan().rays()
            N(1, 0),
            N(0, 1)
            in 2-d lattice N
            sage: phi_s
            Fan morphism defined by the matrix
            [2 0]
            [1 1]
            Domain fan: Rational polyhedral fan in 2-d lattice N
            Codomain fan: Rational polyhedral fan in Sublattice <N(1, 0, 0), N(0, 1, 0)>
            sage: phi_s.codomain_fan().rays()
            N(1, 0, 0),
            N(1, 1, 0)
            in Sublattice <N(1, 0, 0), N(0, 1, 0)>

        Blowup chart (birational)::

            sage: phi_b
            Fan morphism defined by the matrix
            [1 0]
            [0 1]
            Domain fan: Rational polyhedral fan in Sublattice <N(1, 0, 0), N(0, 1, 0)>
            Codomain fan: Rational polyhedral fan in Sublattice <N(1, 0, 0), N(0, 1, 0)>
            sage: phi_b.codomain_fan().rays()
            N( 1,  0, 0),
            N( 0,  1, 0),
            N(-1, -1, 0)
            in Sublattice <N(1, 0, 0), N(0, 1, 0)>

        Coordinate plane inclusion (injective)::

            sage: phi_i
            Fan morphism defined by the matrix
            [1 0 0]
            [0 1 0]
            Domain fan: Rational polyhedral fan in Sublattice <N(1, 0, 0), N(0, 1, 0)>
            Codomain fan: Rational polyhedral fan in 3-d lattice N
            sage: phi.codomain_fan().rays()
            N( 1,  0,  0),
            N( 0,  1,  0),
            N( 0,  0,  1),
            N(-1, -1, -1)
            in 3-d lattice N

        TESTS::

            sage: phi_s.matrix() * phi_b.matrix() * phi_i.matrix() == m
            True

            sage: phi.domain_fan() is phi_s.domain_fan()
            True
            sage: phi_s.codomain_fan() is phi_b.domain_fan()
            True
            sage: phi_b.codomain_fan() is phi_i.domain_fan()
            True
            sage: phi_i.codomain_fan() is phi.codomain_fan()
            True
        """
        L = self.image().saturation()
        d = L.dimension()
        m = self.matrix()
        m = matrix(ZZ, m.nrows(), d, (L.coordinates(c) for c in m.rows()))
        phi_s = FanMorphism(m, self.domain_fan(), L, check=False)
        Sigma_prime = self.codomain_fan()
        L_cone = Cone(sum(([g, -g] for g in L.gens()), []))
        Sigma_i = Fan(cones=(L_cone.intersection(cone) for cone in Sigma_prime),
                      lattice=L, discard_faces=True, check=False)
        phi_b = FanMorphism(identity_matrix(d), phi_s.codomain_fan(), Sigma_i,
                            check=False)
        phi_i = FanMorphism(L.basis_matrix(), Sigma_i, Sigma_prime, check=False)
        return (phi_i, phi_b, phi_s)