Ticket #12159: trac_12159_normal_cone.patch

File trac_12159_normal_cone.patch, 15.4 KB (added by vbraun, 10 years ago)

Updated patch

  • sage/geometry/triangulation/element.py

    # HG changeset patch
    # User Volker Braun <vbraun@stp.dias.ie>
    # Date 1327452598 28800
    # Node ID 5b04cedf5037ddda859582d89ba2d59854d22129
    # Parent  30b758a1bed92f28e01619e8911d54d4a314d130
    Trac #12159: Placing triangulation and normal cones
    
    This patch implements the normal cone of a triangulation.
    
    diff --git a/sage/geometry/triangulation/element.py b/sage/geometry/triangulation/element.py
    a b  
    1010methods to triangulate it according to your requirements. You should
    1111never have to construct a :class:`Triangulation` object directly.
    1212
    13 EXAMPLES::
     13EXAMPLES:
     14
     15First, we select the internal implementation for enumerating
     16triangulations::
    1417
    1518    sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM
    1619
     20Here is a simple example of how to triangulate a point configuration::
     21
    1722    sage: p = [[0,-1,-1],[0,0,1],[0,1,0], [1,-1,-1],[1,0,1],[1,1,0]]
    1823    sage: points = PointConfiguration(p)
    19     sage: triang = points.triangulate()
     24    sage: triang = points.triangulate();  triang
     25    (<0,1,2,5>, <0,1,3,5>, <1,3,4,5>)
    2026    sage: triang.plot(axes=False)
    2127
    2228See :mod:`sage.geometry.triangulation.point_configuration` for more details.
     
    3440from sage.structure.element import Element
    3541from sage.rings.all import QQ, ZZ
    3642from sage.modules.all import vector
    37 
     43from sage.misc.cachefunc import cached_method
    3844
    3945
    4046########################################################################
     
    5359    A 2-d graphics object.
    5460
    5561    EXAMPLES::
    56    
     62
    5763        sage: points = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
    5864        sage: triang = points.triangulate()
    59         sage: triang.plot(axes=False)   # indirect doctest
     65        sage: triang.plot(axes=False, aspect_ratio=1)   # indirect doctest
    6066    """
    6167    from sage.plot.all import point2d, line2d, arrow, polygon2d
    6268    points = [ point.reduced_affine() for point in triangulation.point_configuration() ]
     
    115121    A 3-d graphics object.
    116122
    117123    EXAMPLES::
    118    
     124
    119125        sage: p = [[0,-1,-1],[0,0,1],[0,1,0], [1,-1,-1],[1,0,1],[1,1,0]]
    120126        sage: points = PointConfiguration(p)
    121127        sage: triang = points.triangulate()
     
    241247        :meth:`~sage.geometry.triangulation.point_configuration.PointConfiguration.bistellar_flips`.
    242248
    243249        EXAMPLES::
    244        
     250
    245251            sage: p = [[0,1],[0,0],[1,0]]
    246252            sage: points = PointConfiguration(p)
    247253            sage: from sage.geometry.triangulation.point_configuration import Triangulation
     
    268274        Returns the point configuration underlying the triangulation.
    269275
    270276        EXAMPLES::
    271        
     277
    272278            sage: pconfig = PointConfiguration([[0,0],[0,1],[1,0]])
    273279            sage: pconfig
    274280            A point configuration in QQ^2 consisting of 3 points. The
     
    331337
    332338        EXAMPLES::
    333339
     340            sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM
    334341            sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
    335342            sage: triangulation = pc.triangulate()
    336343            sage: iter = triangulation.__iter__()
     
    360367        A tuple of integers. The vertex indices of the i-th simplex.
    361368
    362369        EXAMPLES::
    363        
     370
     371            sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM
    364372            sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
    365373            sage: triangulation = pc.triangulate()
    366374            sage: triangulation[1]
     
    375383
    376384        TESTS::
    377385
     386            sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM
    378387            sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
    379388            sage: triangulation = pc.triangulations().next()
    380389            sage: triangulation.__len__()
     
    391400       
    392401        TESTS::
    393402
     403            sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM
    394404            sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1],[2,2]])
    395405            sage: t = pc.triangulations()
    396406            sage: t.next()._repr_()
     
    410420        Produce a graphical representation of the triangulation.
    411421
    412422        EXAMPLES::
    413        
     423
    414424            sage: p = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
    415425            sage: triangulation = p.triangulate()
    416426            sage: triangulation
     
    433443        r"""
    434444        Calculate the GKZ phi vector of the triangulation. 
    435445       
     446        The phi vector is a vector of length equals to the number of
     447        points in the point configuration. For a fixed triangulation
     448        `T`, the entry corresponding to the `i`-th point `p_i` is
     449
     450        .. math::
     451
     452            \phi_T(p_i) = \sum_{t\in T, t\owns p_i} Vol(t)
     453
     454        that is, the total volume of all simplices containing `p_i`.
     455        See also [GKZ]_ page 220 equation 1.4.
     456
    436457        OUTPUT:
    437458
    438         Vector -- the phi vector of self.
     459        The phi vector of self.
    439460
    440461        EXAMPLES::
    441462
     
    444465            (3, 1, 5, 2, 4)
    445466            sage: p.lexicographic_triangulation().gkz_phi()
    446467            (1, 3, 4, 2, 5)
    447 
    448         NOTE:
    449 
    450         For a definition of the phi vector, see [GKZ]_ page 220 equation 1.4.
    451468        """
    452469        vec = vector(ZZ, self.point_configuration().n_points())
    453470        for simplex in self:
     
    525542        EXAMPLES::
    526543       
    527544            sage: pc = PointConfiguration([(0,0), (1,0), (0,1), (-1,-1)], star=0, fine=True)
    528             sage: pc.set_engine('internal')   # to make doctests independent of TOPCOM
    529545            sage: triangulation = pc.triangulate()
    530546            sage: fan = triangulation.fan(); fan
    531547            Rational polyhedral fan in 2-d lattice N
     
    538554            sage: interior=[(0, 1, 1), (0, 1, 2), (0, 2, 1), (0, 2, 2)]
    539555            sage: points = vertices+interior
    540556            sage: pc = PointConfiguration(points, fine=True)
    541             sage: pc.set_engine('internal')   # to make doctests independent of TOPCOM
    542557            sage: triangulation = pc.triangulate()
    543558            sage: fan = triangulation.fan( (-1,0,0) )
    544559            sage: fan
     
    556571        return Fan(self, (vector(R, p) - origin for p in points))
    557572
    558573
     574    @cached_method
    559575    def simplicial_complex(self):
    560576        r"""
    561577        Return a simplicial complex from a triangulation of the point
     
    583599        return SimplicialComplex(vertex_set = vertex_set,
    584600                                 maximal_faces = self)
    585601
     602
     603    @cached_method
     604    def _boundary_simplex_dictionary(self):
     605        """
     606        Return facets and the simplices they bound
     607
     608        TESTS::
     609
     610            sage: triangulation = polytopes.n_cube(2).triangulate(engine='internal')
     611            sage: triangulation._boundary_simplex_dictionary()
     612            {(0, 1): ((0, 1, 3),),
     613             (0, 3): ((0, 1, 3), (0, 2, 3)),
     614             (1, 3): ((0, 1, 3),),
     615             (2, 3): ((0, 2, 3),),
     616             (0, 2): ((0, 2, 3),)}
     617
     618            sage: triangulation = polytopes.n_cube(3).triangulate(engine='internal')
     619            sage: triangulation._boundary_simplex_dictionary()
     620            {(1, 4, 7): ((0, 1, 4, 7), (1, 4, 5, 7)),
     621             (1, 3, 7): ((1, 2, 3, 7),),
     622             (0, 1, 7): ((0, 1, 2, 7), (0, 1, 4, 7)),
     623             (0, 2, 7): ((0, 1, 2, 7), (0, 2, 4, 7)),
     624             (0, 1, 4): ((0, 1, 4, 7),),
     625             (2, 4, 6): ((2, 4, 6, 7),),
     626             (0, 1, 2): ((0, 1, 2, 7),),
     627             (1, 2, 7): ((0, 1, 2, 7), (1, 2, 3, 7)),
     628             (2, 6, 7): ((2, 4, 6, 7),),
     629             (2, 3, 7): ((1, 2, 3, 7),),
     630             (1, 4, 5): ((1, 4, 5, 7),),
     631             (1, 5, 7): ((1, 4, 5, 7),),
     632             (4, 5, 7): ((1, 4, 5, 7),),
     633             (0, 4, 7): ((0, 1, 4, 7), (0, 2, 4, 7)),
     634             (2, 4, 7): ((0, 2, 4, 7), (2, 4, 6, 7)),
     635             (1, 2, 3): ((1, 2, 3, 7),),
     636             (4, 6, 7): ((2, 4, 6, 7),),
     637             (0, 2, 4): ((0, 2, 4, 7),)}
     638        """
     639        result = dict()
     640        for simplex in self:
     641            for i in range(len(simplex)):
     642                facet = simplex[:i] + simplex[i+1:]
     643                result[facet] = result.get(facet, tuple()) + (simplex,)
     644        return result
     645
     646
     647    @cached_method
     648    def boundary(self):
     649        """
     650        Return the boundary of the triangulation.
     651       
     652        OUTPUT:
     653
     654        The outward-facing boundary simplices (of dimension `d-1`) of
     655        the `d`-dimensional triangulation as a set. Each boundary is
     656        returned by a tuple of point indices.
     657
     658        EXAMPLES::
     659
     660            sage: triangulation = polytopes.n_cube(3).triangulate(engine='internal')
     661            sage: triangulation
     662            (<0,1,2,7>, <0,1,4,7>, <0,2,4,7>, <1,2,3,7>, <1,4,5,7>, <2,4,6,7>)
     663            sage: triangulation.boundary()
     664            frozenset([(1, 3, 7), (4, 5, 7), (1, 2, 3), (0, 1, 2), (2, 4, 6), (2, 6, 7),
     665                       (2, 3, 7), (1, 5, 7), (0, 1, 4), (1, 4, 5), (4, 6, 7), (0, 2, 4)])
     666            sage: triangulation.interior_facets()
     667            frozenset([(1, 4, 7), (1, 2, 7), (2, 4, 7), (0, 1, 7), (0, 4, 7), (0, 2, 7)])
     668        """
     669        return frozenset(facet for facet, bounded_simplices
     670                         in self._boundary_simplex_dictionary().iteritems()
     671                         if len(bounded_simplices)==1)
     672
     673    @cached_method
     674    def interior_facets(self):
     675        """
     676        Return the interior facets of the triangulation.
     677       
     678        OUTPUT:
     679
     680        The inward-facing boundary simplices (of dimension `d-1`) of
     681        the `d`-dimensional triangulation as a set. Each boundary is
     682        returned by a tuple of point indices.
     683
     684        EXAMPLES::
     685
     686            sage: triangulation = polytopes.n_cube(3).triangulate(engine='internal')
     687            sage: triangulation
     688            (<0,1,2,7>, <0,1,4,7>, <0,2,4,7>, <1,2,3,7>, <1,4,5,7>, <2,4,6,7>)
     689            sage: triangulation.boundary()
     690            frozenset([(1, 3, 7), (4, 5, 7), (1, 2, 3), (0, 1, 2), (2, 4, 6), (2, 6, 7),
     691                       (2, 3, 7), (1, 5, 7), (0, 1, 4), (1, 4, 5), (4, 6, 7), (0, 2, 4)])
     692            sage: triangulation.interior_facets()
     693            frozenset([(1, 4, 7), (1, 2, 7), (2, 4, 7), (0, 1, 7), (0, 4, 7), (0, 2, 7)])
     694        """
     695        return frozenset(facet for facet, bounded_simplices
     696                         in self._boundary_simplex_dictionary().iteritems()
     697                         if len(bounded_simplices)==2)
     698       
     699    @cached_method
     700    def normal_cone(self):
     701        r"""
     702        Return the (closure of the) normal cone of the triangulation.
     703
     704        Recall that a regular triangulation is one that equals the
     705        "crease lines" of a convex piecewise-linear function. This
     706        support function is not unique, for example, you can scale it
     707        by a positive constant. The set of all piecewise-linear
     708        functions with fixed creases forms an open cone. This cone can
     709        be interpreted as the cone of normal vectors at a point of the
     710        secondary polytope, which is why we call it normal cone. See
     711        [GKZ]_ Section 7.1 for details.
     712
     713        OUTPUT:
     714
     715        The closure of the normal cone. The `i`-th entry equals the
     716        value of the piecewise-linear function at the `i`-th point of
     717        the configuration.
     718
     719        For an irregular triangulation, the normal cone is empty. In
     720        this case, a single point (the origin) is returned.
     721
     722        EXAMPLES::
     723
     724            sage: triangulation = polytopes.n_cube(2).triangulate(engine='internal')
     725            sage: triangulation
     726            (<0,1,3>, <0,2,3>)
     727            sage: N = triangulation.normal_cone();  N
     728            4-d cone in 4-d lattice
     729            sage: N.rays()
     730            ((-1, 0, 0, 0), (1, 0, 1, 0), (-1, 0, -1, 0), (1, 0, 0, -1),
     731             (-1, 0, 0, 1), (1, 1, 0, 0), (-1, -1, 0, 0))
     732            sage: N.dual().rays()
     733            ((-1, 1, 1, -1),)
     734
     735        TESTS::
     736       
     737            sage: polytopes.n_simplex(2).triangulate().normal_cone()
     738            3-d cone in 3-d lattice
     739            sage: _.dual().is_trivial()
     740            True
     741        """
     742        if not self.point_configuration().base_ring().is_subring(QQ):
     743            raise NotImplementedError('Only base rings ZZ and QQ are supported')
     744        from sage.libs.ppl import Variable, Constraint, Constraint_System, Linear_Expression, C_Polyhedron
     745        from sage.matrix.constructor import matrix
     746        from sage.misc.misc import uniq
     747        from sage.rings.arith import lcm
     748        pc = self.point_configuration()
     749        cs = Constraint_System()
     750        for facet in self.interior_facets():
     751            s0, s1 = self._boundary_simplex_dictionary()[facet]
     752            p = set(s0).difference(facet).pop()
     753            q = set(s1).difference(facet).pop()
     754            origin = pc.point(p).reduced_affine_vector()
     755            base_indices = [ i for i in s0 if i!=p ]
     756            base = matrix([ pc.point(i).reduced_affine_vector()-origin for i in base_indices ])
     757            sol = base.solve_left( pc.point(q).reduced_affine_vector()-origin )
     758            relation = [0]*pc.n_points()
     759            relation[p] = sum(sol)-1
     760            relation[q] = 1
     761            for i, base_i in enumerate(base_indices):
     762                relation[base_i] = -sol[i]
     763            rel_denom = lcm([QQ(r).denominator() for r in relation])
     764            relation = [ ZZ(r*rel_denom) for r in relation ]
     765            ex = Linear_Expression(relation,0)
     766            cs.insert(ex >= 0)
     767        from sage.modules.free_module import FreeModule
     768        ambient = FreeModule(ZZ, self.point_configuration().n_points())
     769        if cs.empty():
     770            cone = C_Polyhedron(ambient.dimension(), 'universe')
     771        else:
     772            cone = C_Polyhedron(cs)
     773        from sage.geometry.cone import _Cone_from_PPL
     774        return _Cone_from_PPL(cone, lattice=ambient)
     775           
     776           
     777       
  • sage/geometry/triangulation/point_configuration.py

    diff --git a/sage/geometry/triangulation/point_configuration.py b/sage/geometry/triangulation/point_configuration.py
    a b  
    19911991           
    19921992    pushing_triangulation = placing_triangulation
    19931993
     1994    @cached_method
     1995    def Gale_transform(self, points=None):
     1996        r"""
     1997        Return the Gale transform of ``self``.
    19941998
     1999        INPUT:
     2000
     2001        - ``points`` -- a tuple of points or point indices or ``None``
     2002          (default). A subset of points for which to compute the Gale
     2003          transform. By default, all points are used.
     2004
     2005        OUTPUT:
     2006
     2007        A matrix over :meth:`base_ring`.
     2008
     2009        EXAMPLES::
     2010       
     2011            sage: pc = PointConfiguration([(0,0),(1,0),(2,1),(1,1),(0,1)])
     2012            sage: pc.Gale_transform()
     2013            [ 1 -1  0  1 -1]
     2014            [ 0  0  1 -2  1]
     2015
     2016            sage: pc.Gale_transform((0,1,3,4))
     2017            [ 1 -1  1 -1]
     2018
     2019            sage: points = (pc.point(0), pc.point(1), pc.point(3), pc.point(4))
     2020            sage: pc.Gale_transform(points)
     2021            [ 1 -1  1 -1]
     2022        """
     2023        self._assert_is_affine()
     2024        if points is None:
     2025            points = self.points()
     2026        else:
     2027            try:
     2028                points = [ self.point(ZZ(i)) for i in points ]
     2029            except TypeError:
     2030                pass
     2031        m = matrix([ (1,) + p.affine() for p in points])
     2032        return m.left_kernel().matrix()