Ticket #9972: trac_9972_add_toric_lattice_morphisms.patch

File trac_9972_add_toric_lattice_morphisms.patch, 23.0 KB (added by novoselt, 9 years ago)
  • doc/en/reference/geometry.rst

    # HG changeset patch
    # User Andrey Novoseltsev <novoselt@gmail.com>
    # Date 1285207580 21600
    # Node ID ae67d4d01eca3f3f6da37e31f4cb38fcb02d55ea
    # Parent  91fbc1c694603c16d3776ae0dc753588a7d35d10
    Trac 9972: Add toric lattice morphisms.
    
    diff -r 91fbc1c69460 -r ae67d4d01eca doc/en/reference/geometry.rst
    a b  
    1515   sage/geometry/fan
    1616   
    1717   sage/geometry/toric_lattice
     18   sage/geometry/toric_lattice_morphism
  • sage/geometry/toric_lattice.py

    diff -r 91fbc1c69460 -r ae67d4d01eca sage/geometry/toric_lattice.py
    a b  
    465465        """
    466466        try:
    467467            self(point)
     468            return True
    468469        except TypeError:
    469470            return False
    470         return True
     471       
     472    def _Hom_(self, Y, category):
     473        r"""
     474        Return the Hom-space from ``self`` to ``Y``.
     475       
     476        INPUT:
     477       
     478        - ``Y`` -- codomain of the Hom-space;
     479       
     480        - ``category`` -- category in which the Hom-space should live.
     481       
     482        OUTPUT:
     483       
     484        - :class:`toric lattce Hom-space
     485          <sage.geometry.toric_lattice_morphism.ToricLatticeHomspace>`.
     486         
     487        TESTS::
     488       
     489            sage: N3 = ToricLattice(3)
     490            sage: N2 = ToricLattice(2)
     491            sage: N3._Hom_(N2, None)
     492            Set of Morphisms from 3-d lattice N to 2-d lattice N
     493            in Category of modules with basis over Integer Ring
     494            sage: N3._Hom_(N2, None) == Hom(N3, N2)
     495            True
     496        """
     497        from sage.geometry.toric_lattice_morphism import ToricLatticeHomspace
     498        return ToricLatticeHomspace(self, Y, category)
    471499       
    472500    # We need to override this function, otherwise e.g. the sum of elements of
    473501    # different lattices of the same dimension will live in ZZ^n.
     
    13991427            N[0, 1, 0]
    14001428        """
    14011429        return str(self.lift()).replace("(", "[", 1).replace(")", "]", 1)
    1402    
  • new file sage/geometry/toric_lattice_morphism.py

    diff -r 91fbc1c69460 -r ae67d4d01eca sage/geometry/toric_lattice_morphism.py
    - +  
     1r"""
     2Morphisms between toric lattices
     3
     4This module was designed as a part of the framework for toric varieties
     5(:mod:`~sage.schemes.generic.toric_variety`,
     6:mod:`~sage.schemes.generic.fano_toric_variety`). Its main purpose is to
     7provide support for working with morphisms compatible with fans, see
     8:meth:`~ToricLatticeMorphism.is_compatible_with` and
     9:meth:`~ToricLatticeMorphism.make_compatible_with`.
     10
     11AUTHORS:
     12
     13- Andrey Novoseltsev (2010-09-22): initial version.
     14
     15EXAMPLES:
     16
     17Let's consider the face and normal fans of the "diamond" and the projection
     18to the `x`-axis::
     19
     20    sage: diamond = lattice_polytope.octahedron(2)
     21    sage: face = FaceFan(diamond)
     22    sage: normal = NormalFan(diamond)
     23    sage: N = face.lattice()
     24    sage: H = End(N)
     25    sage: phi = H([N.0, 0])
     26    sage: phi     
     27    Free module morphism defined by the matrix
     28    [1 0]
     29    [0 0]
     30    Domain: 2-d lattice N
     31    Codomain: 2-d lattice N
     32    sage: phi.is_compatible_with(normal, face)
     33    False
     34   
     35Some of the cones of the normal fan fail to be mapped to a single cone of the
     36face fan. We can rectify the situation in the following way::
     37
     38    sage: subdivision = phi.make_compatible_with(normal, face)
     39    sage: subdivision.ray_matrix()
     40    [-1  1 -1  1  0  0]
     41    [ 1  1 -1 -1 -1  1]   
     42    sage: normal.ray_matrix()
     43    [-1  1 -1  1]
     44    [ 1  1 -1 -1]
     45   
     46As you see, it was necessary to insert two new rays (to prevent "upper" and
     47"lower" cones of the normal fan from being mapped to the whole `x`-axis).
     48"""
     49
     50
     51#*****************************************************************************
     52#       Copyright (C) 2010 Andrey Novoseltsev <novoselt@gmail.com>
     53#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     54#
     55#  Distributed under the terms of the GNU General Public License (GPL)
     56#
     57#                  http://www.gnu.org/licenses/
     58#*****************************************************************************
     59
     60
     61import operator
     62
     63from sage.geometry.cone import Cone, is_Cone
     64from sage.geometry.fan import Fan, is_Fan
     65from sage.matrix.all import is_Matrix, matrix
     66from sage.misc.all import walltime
     67from sage.modules.free_module_homspace import FreeModuleHomspace
     68from sage.modules.free_module_morphism import FreeModuleMorphism
     69from sage.rings.all import ZZ
     70
     71
     72class ToricLatticeHomspace(FreeModuleHomspace):
     73    r"""
     74    Create a space of homomorphisms between toric lattices.
     75   
     76    INPUT:
     77   
     78    - same as for
     79      :class:`~sage.modules.free_module_homspace.FreeModuleHomspace`.
     80     
     81    OUTPUT:
     82   
     83    - a space of homomorphisms between toric lattices.
     84   
     85    EXAMPLES::
     86   
     87        sage: N = ToricLattice(3)
     88        sage: M = N.dual()
     89        sage: H = Hom(N, M)  # indirect doctest
     90        sage: H
     91        Set of Morphisms from 3-d lattice N to 3-d lattice M
     92        in Category of modules with basis over Integer Ring
     93    """
     94    # sage: TestSuite(H).run() is not tested since the base class fails it
     95   
     96    # We could be happy with the base class call, but it return objects of the
     97    # wrong class
     98    def __call__(self, A, check=True):
     99        """
     100        Construct a morphism corresponding to ``A``.
     101       
     102        INPUT:
     103       
     104        - ``A`` -- a matrix or a list of images of generators of the domain of
     105          ``self``;
     106         
     107        - ``check`` -- boolean (default: ``True``).
     108
     109        If ``A`` is a matrix, then it is the matrix of this linear
     110        transformation, with respect to the bases for the domain and codomain
     111        of ``self``.
     112
     113        EXAMPLES::
     114       
     115            sage: N3 = ToricLattice(3, "N3")
     116            sage: N2 = ToricLattice(2, "N2")
     117            sage: H = Hom(N3, N2)
     118            sage: H
     119            Set of Morphisms from 3-d lattice N3 to 2-d lattice N2
     120            in Category of modules with basis over Integer Ring
     121            sage: phi = H([N2.0, N2.1, N2.0])   # indirect doctest
     122            sage: phi
     123            Free module morphism defined by the matrix
     124            [1 0]
     125            [0 1]
     126            [1 0]
     127            Domain: 3-d lattice N3
     128            Codomain: 2-d lattice N2
     129            sage: phi(N3(1,2,3))
     130            N2(4, 2)
     131        """
     132        if not is_Matrix(A):
     133            # Compute the matrix of the morphism that sends the
     134            # generators of the domain to the elements of A.
     135            codomain = self.codomain()
     136            try:
     137                A = matrix(ZZ, [codomain.coordinates(codomain(a)) for a in A])
     138            except TypeError:
     139                pass
     140        return ToricLatticeMorphism(self, A)
     141
     142
     143class ToricLatticeMorphism(FreeModuleMorphism):
     144    r"""
     145    Create a morphism between toric lattices.
     146   
     147    INPUT:
     148   
     149    - ``parent`` -- a :class:`space of homomorphisms between toric lattices
     150      <ToricLatticeHomspace>`;
     151   
     152    - ``A`` -- an integral matrix defining the morphism.
     153
     154    OUTPUT:
     155   
     156    - a morphism between toric lattices.
     157   
     158    EXAMPLES::
     159   
     160        sage: N3 = ToricLattice(3, "N3")
     161        sage: N2 = ToricLattice(2, "N2")
     162        sage: H = Hom(N3, N2)
     163        sage: phi = H([N2.0, N2.1, N2.0])  # indirect doctest
     164        sage: phi
     165        Free module morphism defined by the matrix
     166        [1 0]
     167        [0 1]
     168        [1 0]
     169        Domain: 3-d lattice N3
     170        Codomain: 2-d lattice N2
     171    """
     172   
     173    def _chambers(self, fan):
     174        r"""
     175        Compute chambers in the domain of ``self`` corresponding to ``fan``.
     176       
     177        This function is useful for automatic refinement of fans to make them
     178        compatible with ``self``, see :meth:`make_compatible_with`.
     179       
     180        INPUT:
     181       
     182        - ``fan`` -- a :class:`fan <sage.geometry.fan.RationalPolyhedralFan>`
     183          in the codomain of ``self``.
     184       
     185        OUTPUT:
     186       
     187        - a :class:`tuple` ``(chambers, cone_to_chamber)``, where
     188       
     189        - ``chambers`` is a :class:`list` of :class:`cones
     190        <sage.geometry.cone.ConvexRationalPolyhedralCone>` in the domain of
     191        ``self``;
     192       
     193        - ``cone_to_chamber`` is a :class:`list` of integers, if its `i`-th
     194          element is `j`, then the `j`-th element of ``chambers`` is the
     195          inverse image of the `i`-th generating cone of ``fan``.
     196                   
     197        TESTS::
     198       
     199            sage: F = NormalFan(lattice_polytope.octahedron(2))
     200            sage: N = F.lattice()
     201            sage: H = End(N)
     202            sage: phi = H([N.0, 0])
     203            sage: phi._chambers(F)
     204            ([2-d cone in 2-d lattice N,
     205              1-d cone in 2-d lattice N,
     206              2-d cone in 2-d lattice N],
     207             [0, 1, 2, 1])           
     208        """
     209        # It probably does not make much sense to cache chambers, since their
     210        # computation for fan subdivision is likely to be much shorter than
     211        # the subdivision itself.
     212        kernel_rays = []
     213        for ray in self.kernel().basis():
     214            kernel_rays.append(ray)
     215            kernel_rays.append(-ray)
     216        image_rays = []
     217        for ray in self.image().basis():
     218            image_rays.append(ray)
     219            image_rays.append(-ray)
     220        image = Cone(image_rays)
     221        chambers = []
     222        cone_to_chamber = []
     223        for cone in fan:
     224            chamber = Cone([self.lift(ray) for ray in cone.intersection(image)]
     225                           + kernel_rays, lattice=self.domain())
     226            cone_to_chamber.append(len(chambers))
     227            for i, old_chamber in enumerate(chambers):
     228                if old_chamber.is_equivalent(chamber):
     229                    cone_to_chamber[-1] = i
     230                    break
     231            if cone_to_chamber[-1] == len(chambers):
     232                chambers.append(chamber)
     233        return (chambers, cone_to_chamber)
     234   
     235    def is_compatible_with(self, x, y):
     236        r"""
     237        Check if ``self`` is compatible with given cones or fans.
     238       
     239        INPUT:
     240       
     241        - ``x`` -- a :class:`cone
     242          <sage.geometry.cone.ConvexRationalPolyhedralCone` or a :class:`fan
     243          <sage.geometry.fan.RationalPolyhedralFan` in the domain of ``self``;
     244                 
     245        - ``y`` -- a :class:`cone
     246          <sage.geometry.cone.ConvexRationalPolyhedralCone` or a :class:`fan
     247          <sage.geometry.fan.RationalPolyhedralFan` in the codomain of
     248          ``self``.
     249         
     250        OUTPUT:
     251       
     252        - ``True`` if the image of every cone of ``x`` is completely contained
     253          in a single cone of ``y``, ``False`` otherwise.
     254         
     255        EXAMPLES::
     256       
     257            sage: N3 = ToricLattice(3, "N3")
     258            sage: N2 = ToricLattice(2, "N2")
     259            sage: H = Hom(N3, N2)
     260            sage: phi = H([N2.0, N2.1, N2.0])
     261            sage: F1 = Fan(cones=[(0,1,2), (1,2,3)],
     262            ...            rays=[(1,1,1), (1,1,-1), (1,-1,1), (1,-1,-1)],
     263            ...            lattice=N3)
     264            sage: F1.ray_matrix()
     265            [ 1  1  1  1]
     266            [ 1  1 -1 -1]
     267            [ 1 -1  1 -1]
     268            sage: [phi(ray) for ray in F1.rays()]
     269            [N2(2, 1), N2(0, 1), N2(2, -1), N2(0, -1)]
     270            sage: C1 = Cone([(1,0), (0,1), (0,-1)], lattice=N2)
     271            sage: phi.is_compatible_with(F1, C1)
     272            True
     273            sage: phi.is_compatible_with(F1.generating_cone(0), C1)
     274            True
     275            sage: C2 = Cone([(1,0), (0,1)], lattice=N2)
     276            sage: phi.is_compatible_with(F1, C2)
     277            False
     278            sage: F2 = Fan(cones=[(0,1,2), (1,2,3)],
     279            ...            rays=[(1,1,1), (1,1,-1), (1,2,1), (1,2,-1)],
     280            ...            lattice=N3)
     281            sage: F2.ray_matrix()
     282            [ 1  1  1  1]
     283            [ 1  1  2  2]
     284            [ 1 -1  1 -1]
     285            sage: [phi(ray) for ray in F2.rays()]
     286            [N2(2, 1), N2(0, 1), N2(2, 2), N2(0, 2)]
     287            sage: phi.is_compatible_with(F2, C2)
     288            True
     289            sage: F3 = Fan(cones=[(0,1), (1,2)],
     290            ...            rays=[(1,0), (2,1), (0,1)],
     291            ...            lattice=N2)
     292            sage: phi.is_compatible_with(F2, F3)
     293            True
     294            sage: phi.is_compatible_with(F1, F3)
     295            False
     296        """
     297        if is_Cone(y):
     298            return all(self(ray) in y for ray in x.rays())
     299        if is_Cone(x):
     300            if x.is_strictly_convex():
     301                x = Fan(cones=[x], check=False)
     302            else:
     303                return False
     304        assert is_Fan(x) and is_Fan(y)       
     305        try:
     306            ray_images = [y.cone_containing(self(ray)) for ray in x.rays()]
     307        except ValueError:
     308            return False
     309        ray_images = [set(image.star_generator_indices())
     310                      for image in ray_images]
     311        for c in x:
     312            if not reduce(operator.and_, (ray_images[i]
     313                                          for i in c.ambient_ray_indices())):
     314                return False
     315        return True
     316
     317    def make_compatible_with(self, domain_fan, codomain_fan, check=True,
     318                             verbose=False):
     319        r"""
     320        Subdivide ``domain_fan`` to make it compatible with ``codomain_fan``.
     321       
     322        These fans are compatible with ``self`` if the image of every cone of
     323        ``domain_fan`` is completely contained within a single cone of
     324        ``codomain_fan``. If this is not the case, one can either subdivide
     325        ``domain_fan`` or combine some of the cones in ``codomain_fan``. The
     326        latter is not always possible (e.g. if some of the cones of
     327        ``domain_fan`` are mapped to non-strictly convex cones) and is not as
     328        nice from the point of view of toric varieties as the subdivision of
     329        ``domain_fan``.
     330       
     331        This function will intersect the generating cones of ``domain_fan``
     332        with preimages of the generating cones of ``codomain_fan`` and return
     333        the resulting fan in the domain. ``ValueError`` exception is raised if
     334        the image of ``domain_fan`` does not lie in the support of
     335        ``codomain_fan`` (which means that it is impossible to fix the
     336        incompatibility between these fans using only subdivision).
     337       
     338        INPUT:
     339       
     340        - ``domain_fan`` -- a :class:`fan
     341          <sage.geometry.fan.RationalPolyhedralFan` in the domain of ``self``;
     342                 
     343        - ``codomain_fan`` -- a :class:`fan
     344          <sage.geometry.fan.RationalPolyhedralFan` in the codomain of
     345          ``self``;
     346         
     347        - ``check`` -- (default: ``True``) if ``False``, some of the
     348          consistency checks will be omitted, which saves time but can
     349          potentially lead to wrong results. Currently, with
     350          ``check=False`` option there will be no check that ``domain_fan``
     351          maps to ``codomain_fan``;
     352         
     353        - ``verbose`` -- (default: ``False``) if ``True``, some timing
     354          information will be printed in the process.
     355         
     356        OUTPUT:
     357       
     358        - a :class:`fan
     359          <sage.geometry.fan.RationalPolyhedralFan` in the domain of ``self``.
     360         
     361        EXAMPLES:
     362       
     363        Here we consider the face and normal fans of the "diamond" and the
     364        projection to the `x`-axis::
     365       
     366            sage: diamond = lattice_polytope.octahedron(2)
     367            sage: face = FaceFan(diamond)
     368            sage: normal = NormalFan(diamond)
     369            sage: N = face.lattice()
     370            sage: H = End(N)
     371            sage: phi = H([N.0, 0])
     372            sage: phi     
     373            Free module morphism defined by the matrix
     374            [1 0]
     375            [0 0]
     376            Domain: 2-d lattice N
     377            Codomain: 2-d lattice N
     378            sage: phi.is_compatible_with(face, normal)
     379            True
     380            sage: phi.make_compatible_with(face, normal) is face
     381            True
     382           
     383        Note, that since ``phi`` is compatible with these fans, the returned
     384        fan is exactly the same object as the initial ``domain_fan``. ::
     385       
     386            sage: phi.is_compatible_with(normal, face)
     387            False
     388            sage: subdivision = phi.make_compatible_with(normal, face)
     389            sage: subdivision is normal
     390            False
     391            sage: subdivision.ngenerating_cones()
     392            6
     393           
     394        We had to subdivide two of the four cones of the normal fan, since
     395        they were mapped by ``phi`` into non-strictly convex cones.
     396       
     397        Now we demonstrate a more subtle example. We take the first quadrant
     398        as our ``domain_fan``. Then we divide the first quadrant into three
     399        cones, throw away the middle one and take the other two as our
     400        ``codomain_fan``. These fans are incompatible with the identity
     401        lattice morphism since the image of ``domain_fan`` is out of the
     402        support of ``codomain_fan``::
     403       
     404            sage: N = ToricLattice(2)
     405            sage: phi = End(N).identity()
     406            sage: F1 = Fan(cones=[(0,1)], rays=[(1,0), (0,1)])
     407            sage: F2 = Fan(cones=[(0,1), (2,3)],
     408            ...            rays=[(1,0), (2,1), (1,2), (0,1)])
     409            sage: phi.is_compatible_with(F1, F2)
     410            False
     411            sage: phi.make_compatible_with(F1, F2)
     412            Traceback (most recent call last):
     413            ...
     414            ValueError: Free module morphism defined by the matrix
     415            [1 0]
     416            [0 1]
     417            Domain: 2-d lattice N
     418            Codomain: 2-d lattice N
     419            does not map
     420            Rational polyhedral fan in 2-d lattice N
     421            into the support of
     422            Rational polyhedral fan in 2-d lattice N!
     423           
     424        The problem was detected and handled correctly (i.e. an exception was
     425        raised). However, the used algorithm requires extra checks for this
     426        situation after constructing a potential subdivision and this can take
     427        significant time. You can save about half the time using
     428        ``check=False`` option, if you know in advance that it is possible to
     429        make fans compatible with the morphism by subdividing ``domain_fan``.
     430        Of course, if your assumption was incorrect, the result will be wrong
     431        and you will get a fan which *does* map into the support of
     432        ``codomain_fan``, but is **not** a subdivision of ``domain_fan``. You
     433        can test it on the example above::
     434       
     435            sage: F3 = phi.make_compatible_with(
     436            ...            F1, F2, check=False, verbose=True)
     437            Placing ray images...
     438            Computing chambers...
     439            Subdividing cone 1 of 1...           
     440            sage: F3.is_equivalent(F2)
     441            True
     442        """
     443        assert is_Fan(domain_fan) and is_Fan(codomain_fan)
     444        if verbose:
     445            start = walltime()
     446            print "Placing ray images...",
     447        # Figure out where 1-dimensional cones (i.e. rays) are mapped.
     448        try:
     449            ray_images = [codomain_fan.cone_containing(self(ray))
     450                          for ray in domain_fan.rays()]
     451        except ValueError:
     452            raise ValueError("%s\ndoes not map\n%s\ninto the support of\n%s!"
     453                             % (self, domain_fan, codomain_fan))
     454        ray_images = [set(image.star_generator_indices())
     455                      for image in ray_images]
     456        if verbose:
     457            print "%.3f ms" % walltime(start)
     458        # Subdivide cones that require it.
     459        chambers = None # preimages of codomain cones, computed if necessary
     460        new_cones = []
     461        for cone_index, domain_cone in enumerate(domain_fan):
     462            if reduce(operator.and_, (ray_images[i]
     463                                for i in domain_cone.ambient_ray_indices())):
     464                new_cones.append(domain_cone)
     465                continue
     466            dim = domain_cone.dim()
     467            if chambers is None:
     468                if verbose:
     469                    start = walltime()
     470                    print "Computing chambers...",
     471                chambers, cone_to_chamber = self._chambers(codomain_fan)
     472                if verbose:
     473                    print "%.3f ms" % walltime(start)
     474            # Subdivide domain_cone.
     475            if verbose:
     476                start = walltime()
     477                print ("Subdividing cone %d of %d..."
     478                       % (cone_index + 1, domain_fan.ngenerating_cones())),
     479            # Only these chambers intersect domain_cone.
     480            containing_chambers = (cone_to_chamber[j]
     481                for j in reduce(operator.or_, (ray_images[i]
     482                                for i in domain_cone.ambient_ray_indices())))
     483            # We don't care about chambers of small dimension.
     484            containing_chambers = (chambers[i]
     485                                   for i in set(containing_chambers)
     486                                   if chambers[i].dim() >= dim)
     487            parts = (domain_cone.intersection(chamber)
     488                     for chamber in containing_chambers)
     489            # We cannot leave parts as a generator since we use them twice.
     490            parts = [part for part in parts if part.dim() == dim]
     491            if check:
     492                # Check if the subdivision is complete, i.e. there are no
     493                # missing pieces of domain_cone. To do this, we construct a
     494                # fan from the obtained parts and check that interior points
     495                # of boundary cones of this fan are in the interior of the
     496                # original cone. In any case we know that we are constructing
     497                # a valid fan, so passing check=False to Fan(...) is OK.
     498                if verbose:
     499                    print "%.3f ms" % walltime(start)
     500                    start = walltime()
     501                    print "Checking for missing pieces... ",
     502                cone_subdivision = Fan(parts, check=False)
     503                for cone in cone_subdivision(dim - 1):
     504                    if len(cone.star_generators()) == 1:
     505                        if domain_cone.relative_interior_contains(
     506                                                            sum(cone.rays())):
     507                            raise ValueError("%s\ndoes not map\n%s\ninto the "
     508                                             "support of\n%s!" % (self,
     509                                             domain_fan, codomain_fan))
     510            new_cones.extend(parts)
     511            if verbose:
     512                print "%.3f ms" % walltime(start)
     513        if len(new_cones) == domain_fan.ngenerating_cones():
     514            return domain_fan
     515        # Construct a new fan keeping old rays in the same order
     516        new_rays = list(domain_fan.rays())
     517        for cone in new_cones:
     518            for ray in cone:
     519                if ray not in new_rays:
     520                    new_rays.append(ray)
     521        return Fan(new_cones, new_rays, domain_cone.lattice(), check=False)