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().ray_matrix()
    [-1  1 -1  1  0  0]
    [ 1  1 -1 -1 -1  1]    
    sage: normal.ray_matrix()
    [-1  1 -1  1]
    [ 1  1 -1 -1]
    
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 matrix, is_Matrix
from sage.misc.all import cached_method, latex, walltime
from sage.modules.free_module_morphism import (FreeModuleMorphism,
                                               is_FreeModuleMorphism)
from sage.rings.all import ZZ, 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().ray_matrix()
        [ 1 -1]
        [ 0  0]
        
    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()
        else:
            self._codomain_fan = codomain_fan
            if subdivide:
                self._subdivide_domain_fan(check, verbose)
            elif check:
                self._validate()
                
    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):
        r"""
        Construct the codomain fan as the image of the domain one.
        
        .. WARNING::
        
            This method should be called only during initialization.
            
        OUTPUT:
        
        - none, but the codomain fan of ``self`` is set to the constucted 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().ray_matrix()  # indirect doctest
            [1 0 1]
            [0 1 1]
        """
        # 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)
    
    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().ray_matrix()  # indirect doctest
            [1 0 1]
            [0 1 1]

        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().ray_matrix()
            [ 1  1  1]
            [ 1 -1  0]
            [ 0  0  0]
            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.ray_matrix()
            [ 1  1  1  1]
            [ 1  1 -1 -1]
            [ 1 -1  1 -1]
            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.ray_matrix()
            [ 1  1  1  1]
            [ 1  1  2  2]
            [ 1 -1  1 -1]
            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]
        
    @cached_method
    def index(self):
        r"""
        Return the index of ``self`` as a map between lattices.
        
        OUTPUT:
        
        - an integer or infinity.
        
        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
        """
        return self.matrix().image().index_in(self.codomain())

    @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 self.restrict_to_image().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: _.ray_matrix()
            [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.ray_matrix()
            [ 0 -1 -1  0  1  1]
            [ 1  0 -1 -1  0  1]
            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]

    @cached_method
    def restrict_to_image(self):
        r"""
        Return the fan morphism from the domain fan to the image fan.
        
        OUTPUT:
        
        - a :class:`fan morphism <FanMorphism>`.
        
        .. note::
        
            By the *image fan* of a fan morphism we mean the fan generated by
            the intersections of cones of the codomain fan with the image
            vector space of ``self``. The lattice of this fan is the saturation
            of the image of ``self``.
        
        EXAMPLES:
        
        We embed a projective line "diagonally" into the product of two lines::
        
            sage: Sigma = toric_varieties.P1().fan()
            sage: Sigmap = toric_varieties.P1xP1().fan()
            sage: m = matrix(1, 2, [1,1])
            sage: phi = FanMorphism(m, Sigma, Sigmap)
            sage: phi
            Fan morphism defined by the matrix
            [1 1]
            Domain fan: Rational polyhedral fan in 1-d lattice N
            Codomain fan: Rational polyhedral fan in 2-d lattice N
            
        Now we restrict this morphism to its image::
        
            sage: psi = phi.restrict_to_image()
            sage: psi
            Fan morphism defined by the matrix
            [1]
            Domain fan: Rational polyhedral fan in 1-d lattice N
            Codomain fan: Rational polyhedral fan in Sublattice <N(1, 1)>
            
        Note that the matrix defining a morphism to a fan in a sublattice
        operates with generators of this sublattice, so ``[1]`` in the above
        example means that `\psi` sends the only generator of the domain
        lattice to the generator of the image sublattice::
        
            sage: psi.codomain().gens()
            (N(1, 1),)
            sage: psi.codomain_fan().ray_matrix()
            [ 1 -1]
            [ 1 -1]
            
        Restriction to image returns exactly the same map if the corresponding
        map of vector spaces is surjective, e.g. in the case of double
        restriction::
        
            sage: psi.restrict_to_image() is psi
            True
        """
        L = self.image().saturation()
        d = L.dimension()
        if self.codomain().dimension() == d:
            return self
        m = self.matrix()
        m = matrix(ZZ, m.nrows(), d, (L.coordinates(c) for c in m.rows()))
        L_cone = Cone(sum(([g, -g] for g in L.gens()), []))
        Sigma = Fan(cones=discard_faces(L_cone.intersection(cone)
                                        for cone in self.codomain_fan()),
                    lattice=L, check=False)
        return FanMorphism(m, self.domain_fan(), Sigma)