Ticket #9972: trac_9972_add_toric_lattice_morphisms.patch

doc/en/reference/geometry.rst

@@ -15 +15 @@
    sage/geometry/fan
    
    sage/geometry/toric_lattice
+   sage/geometry/toric_lattice_morphism

sage/geometry/toric_lattice.py

@@ -465 +465 @@
         """
         try:
             self(point)
+            return True
         except TypeError:
             return False
-        return True
+        
+    def _Hom_(self, Y, category):
+        r"""
+        Return the Hom-space from ``self`` to ``Y``.
+        
+        INPUT:
+        
+        - ``Y`` -- codomain of the Hom-space;
+        
+        - ``category`` -- category in which the Hom-space should live.
+        
+        OUTPUT:
+        
+        - :class:`toric lattce Hom-space
+          <sage.geometry.toric_lattice_morphism.ToricLatticeHomspace>`.
+          
+        TESTS::
+        
+            sage: N3 = ToricLattice(3)
+            sage: N2 = ToricLattice(2)
+            sage: N3._Hom_(N2, None)
+            Set of Morphisms from 3-d lattice N to 2-d lattice N
+            in Category of modules with basis over Integer Ring
+            sage: N3._Hom_(N2, None) == Hom(N3, N2)
+            True
+        """
+        from sage.geometry.toric_lattice_morphism import ToricLatticeHomspace
+        return ToricLatticeHomspace(self, Y, category)
         
     # We need to override this function, otherwise e.g. the sum of elements of
     # different lattices of the same dimension will live in ZZ^n.
@@ -1399 +1427 @@
             N[0, 1, 0]
         """
         return str(self.lift()).replace("(", "[", 1).replace(")", "]", 1)
-    

new file sage/geometry/toric_lattice_morphism.py

@@ +1 @@
+r"""
+Morphisms between toric lattices
+
+This module was designed as a part of the framework for toric varieties
+(:mod:`~sage.schemes.generic.toric_variety`,
+:mod:`~sage.schemes.generic.fano_toric_variety`). Its main purpose is to
+provide support for working with morphisms compatible with fans, see
+:meth:`~ToricLatticeMorphism.is_compatible_with` and
+:meth:`~ToricLatticeMorphism.make_compatible_with`.
+
+AUTHORS:
+
+- Andrey Novoseltsev (2010-09-22): initial version.
+
+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: phi.is_compatible_with(normal, face)
+    False
+    
+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: subdivision = phi.make_compatible_with(normal, face)
+    sage: subdivision.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 <novoselt@gmail.com>
+#       Copyright (C) 2010 William Stein <wstein@gmail.com>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+
+import operator
+
+from sage.geometry.cone import Cone, is_Cone
+from sage.geometry.fan import Fan, is_Fan
+from sage.matrix.all import is_Matrix, matrix
+from sage.misc.all import walltime
+from sage.modules.free_module_homspace import FreeModuleHomspace
+from sage.modules.free_module_morphism import FreeModuleMorphism
+from sage.rings.all import ZZ
+
+
+class ToricLatticeHomspace(FreeModuleHomspace):
+    r"""
+    Create a space of homomorphisms between toric lattices.
+    
+    INPUT:
+    
+    - same as for
+      :class:`~sage.modules.free_module_homspace.FreeModuleHomspace`.
+      
+    OUTPUT:
+    
+    - a space of homomorphisms between toric lattices.
+    
+    EXAMPLES::
+    
+        sage: N = ToricLattice(3)
+        sage: M = N.dual()
+        sage: H = Hom(N, M)  # indirect doctest
+        sage: H
+        Set of Morphisms from 3-d lattice N to 3-d lattice M
+        in Category of modules with basis over Integer Ring
+    """
+    # sage: TestSuite(H).run() is not tested since the base class fails it
+    
+    # We could be happy with the base class call, but it return objects of the
+    # wrong class
+    def __call__(self, A, check=True):
+        """
+        Construct a morphism corresponding to ``A``.
+        
+        INPUT:
+        
+        - ``A`` -- a matrix or a list of images of generators of the domain of
+          ``self``;
+          
+        - ``check`` -- boolean (default: ``True``).
+
+        If ``A`` is a matrix, then it is the matrix of this linear
+        transformation, with respect to the bases for the domain and codomain
+        of ``self``.
+
+        EXAMPLES::
+        
+            sage: N3 = ToricLattice(3, "N3")
+            sage: N2 = ToricLattice(2, "N2")
+            sage: H = Hom(N3, N2)
+            sage: H
+            Set of Morphisms from 3-d lattice N3 to 2-d lattice N2
+            in Category of modules with basis over Integer Ring
+            sage: phi = H([N2.0, N2.1, N2.0])   # indirect doctest
+            sage: phi
+            Free module morphism defined by the matrix
+            [1 0]
+            [0 1]
+            [1 0]
+            Domain: 3-d lattice N3
+            Codomain: 2-d lattice N2
+            sage: phi(N3(1,2,3))
+            N2(4, 2)
+        """
+        if not is_Matrix(A):
+            # Compute the matrix of the morphism that sends the
+            # generators of the domain to the elements of A.
+            codomain = self.codomain()
+            try:
+                A = matrix(ZZ, [codomain.coordinates(codomain(a)) for a in A])
+            except TypeError:
+                pass
+        return ToricLatticeMorphism(self, A)
+
+
+class ToricLatticeMorphism(FreeModuleMorphism):
+    r"""
+    Create a morphism between toric lattices.
+    
+    INPUT:
+    
+    - ``parent`` -- a :class:`space of homomorphisms between toric lattices
+      <ToricLatticeHomspace>`;
+    
+    - ``A`` -- an integral matrix defining the morphism.
+
+    OUTPUT:
+    
+    - a morphism between toric lattices.
+    
+    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])  # indirect doctest
+        sage: phi
+        Free module morphism defined by the matrix
+        [1 0]
+        [0 1]
+        [1 0]
+        Domain: 3-d lattice N3
+        Codomain: 2-d lattice N2
+    """
+    
+    def _chambers(self, fan):
+        r"""
+        Compute chambers in the domain of ``self`` corresponding to ``fan``.
+        
+        This function is useful for automatic refinement of fans to make them
+        compatible with ``self``, see :meth:`make_compatible_with`.
+        
+        INPUT:
+        
+        - ``fan`` -- a :class:`fan <sage.geometry.fan.RationalPolyhedralFan>`
+          in the codomain of ``self``.
+        
+        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 ``fan``.
+                    
+        TESTS::
+        
+            sage: F = NormalFan(lattice_polytope.octahedron(2))
+            sage: N = F.lattice()
+            sage: H = End(N)
+            sage: phi = H([N.0, 0])
+            sage: phi._chambers(F)
+            ([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])            
+        """
+        # It probably does not make much sense to cache chambers, since their
+        # computation for fan subdivision is likely to be much shorter than
+        # the subdivision itself.
+        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 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 is_compatible_with(self, x, y):
+        r"""
+        Check if ``self`` is compatible with given cones or fans.
+        
+        INPUT:
+        
+        - ``x`` -- a :class:`cone
+          <sage.geometry.cone.ConvexRationalPolyhedralCone` or a :class:`fan
+          <sage.geometry.fan.RationalPolyhedralFan` in the domain of ``self``;
+                  
+        - ``y`` -- a :class:`cone
+          <sage.geometry.cone.ConvexRationalPolyhedralCone` or a :class:`fan
+          <sage.geometry.fan.RationalPolyhedralFan` in the codomain of
+          ``self``.
+          
+        OUTPUT:
+        
+        - ``True`` if the image of every cone of ``x`` is completely contained
+          in a single cone of ``y``, ``False`` otherwise.
+          
+        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: C1 = Cone([(1,0), (0,1), (0,-1)], lattice=N2)
+            sage: phi.is_compatible_with(F1, C1)
+            True
+            sage: phi.is_compatible_with(F1.generating_cone(0), C1)
+            True
+            sage: C2 = Cone([(1,0), (0,1)], lattice=N2)
+            sage: phi.is_compatible_with(F1, C2)
+            False
+            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: phi.is_compatible_with(F2, C2)
+            True
+            sage: F3 = Fan(cones=[(0,1), (1,2)],
+            ...            rays=[(1,0), (2,1), (0,1)],
+            ...            lattice=N2)
+            sage: phi.is_compatible_with(F2, F3)
+            True
+            sage: phi.is_compatible_with(F1, F3)
+            False
+        """
+        if is_Cone(y):
+            return all(self(ray) in y for ray in x.rays())
+        if is_Cone(x):
+            if x.is_strictly_convex():
+                x = Fan(cones=[x], check=False)
+            else:
+                return False
+        assert is_Fan(x) and is_Fan(y)        
+        try:
+            ray_images = [y.cone_containing(self(ray)) for ray in x.rays()]
+        except ValueError:
+            return False
+        ray_images = [set(image.star_generator_indices())
+                      for image in ray_images]
+        for c in x:
+            if not reduce(operator.and_, (ray_images[i]
+                                          for i in c.ambient_ray_indices())):
+                return False
+        return True
+
+    def make_compatible_with(self, domain_fan, codomain_fan, check=True,
+                             verbose=False):
+        r"""
+        Subdivide ``domain_fan`` to make it compatible with ``codomain_fan``.
+        
+        These fans are compatible with ``self`` if the image of every cone of
+        ``domain_fan`` is completely contained within a single cone of
+        ``codomain_fan``. If this is not the case, one can either subdivide 
+        ``domain_fan`` or combine some of the cones in ``codomain_fan``. The
+        latter is not always possible (e.g. if some of the cones of
+        ``domain_fan`` are mapped to non-strictly convex cones) and is not as 
+        nice from the point of view of toric varieties as the subdivision of 
+        ``domain_fan``.
+       
+        This function will intersect the generating cones of ``domain_fan``
+        with preimages of the generating cones of ``codomain_fan`` and return
+        the resulting fan in the domain. ``ValueError`` exception is raised if
+        the image of ``domain_fan`` does not lie in the support of
+        ``codomain_fan`` (which means that it is impossible to fix the
+        incompatibility between these fans using only subdivision).
+        
+        INPUT:
+        
+        - ``domain_fan`` -- a :class:`fan
+          <sage.geometry.fan.RationalPolyhedralFan` in the domain of ``self``;
+                  
+        - ``codomain_fan`` -- a :class:`fan
+          <sage.geometry.fan.RationalPolyhedralFan` in the codomain of
+          ``self``;
+          
+        - ``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:
+        
+        - a :class:`fan
+          <sage.geometry.fan.RationalPolyhedralFan` in the domain of ``self``.
+          
+        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: phi.is_compatible_with(face, normal)
+            True
+            sage: phi.make_compatible_with(face, normal) 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: phi.is_compatible_with(normal, face)
+            False
+            sage: subdivision = phi.make_compatible_with(normal, face)
+            sage: subdivision is normal
+            False
+            sage: subdivision.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.
+        
+        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: phi.is_compatible_with(F1, F2)
+            False
+            sage: phi.make_compatible_with(F1, F2)
+            Traceback (most recent call last):
+            ...
+            ValueError: Free module morphism defined by the matrix
+            [1 0]
+            [0 1]
+            Domain: 2-d lattice N
+            Codomain: 2-d lattice N
+            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 ``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
+        ``codomain_fan``, but is **not** a subdivision of ``domain_fan``. You 
+        can test it on the example above::
+        
+            sage: F3 = phi.make_compatible_with(
+            ...            F1, F2, check=False, verbose=True)
+            Placing ray images...
+            Computing chambers...
+            Subdividing cone 1 of 1...            
+            sage: F3.is_equivalent(F2)
+            True
+        """
+        assert is_Fan(domain_fan) and is_Fan(codomain_fan)
+        if verbose:
+            start = walltime()
+            print "Placing ray images...",
+        # Figure out where 1-dimensional cones (i.e. rays) are mapped.
+        try:
+            ray_images = [codomain_fan.cone_containing(self(ray))
+                          for ray in domain_fan.rays()]
+        except ValueError:
+            raise ValueError("%s\ndoes not map\n%s\ninto the support of\n%s!"
+                             % (self, domain_fan, codomain_fan))
+        ray_images = [set(image.star_generator_indices()) 
+                      for image in ray_images]
+        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_, (ray_images[i]
+                                for i in domain_cone.ambient_ray_indices())):
+                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(codomain_fan)
+                if verbose:
+                    print "%.3f ms" % walltime(start)
+            # Subdivide domain_cone.
+            if verbose:
+                start = walltime()
+                print ("Subdividing cone %d of %d..."
+                       % (cone_index + 1, domain_fan.ngenerating_cones())),
+            # Only these chambers intersect domain_cone.
+            containing_chambers = (cone_to_chamber[j]
+                for j in reduce(operator.or_, (ray_images[i]
+                                for i in domain_cone.ambient_ray_indices())))
+            # We don't care about chambers of small dimension.
+            containing_chambers = (chambers[i]
+                                   for i in set(containing_chambers)
+                                   if chambers[i].dim() >= dim)
+            parts = (domain_cone.intersection(chamber)
+                     for chamber in containing_chambers)
+            # We cannot leave parts as a generator since we use them twice.
+            parts = [part for part in parts if part.dim() == dim]
+            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())):
+                            raise ValueError("%s\ndoes not map\n%s\ninto the "
+                                             "support of\n%s!" % (self,
+                                             domain_fan, codomain_fan))
+            new_cones.extend(parts)
+            if verbose:
+                print "%.3f ms" % walltime(start)
+        if len(new_cones) == domain_fan.ngenerating_cones():
+            return domain_fan
+        # 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)
+        return Fan(new_cones, new_rays, domain_cone.lattice(), check=False)