Ticket #12159: trac_12159_normal_cone.patch
File trac_12159_normal_cone.patch, 15.4 KB (added by , 10 years ago) 


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 10 10 methods to triangulate it according to your requirements. You should 11 11 never have to construct a :class:`Triangulation` object directly. 12 12 13 EXAMPLES:: 13 EXAMPLES: 14 15 First, we select the internal implementation for enumerating 16 triangulations:: 14 17 15 18 sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM 16 19 20 Here is a simple example of how to triangulate a point configuration:: 21 17 22 sage: p = [[0,1,1],[0,0,1],[0,1,0], [1,1,1],[1,0,1],[1,1,0]] 18 23 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>) 20 26 sage: triang.plot(axes=False) 21 27 22 28 See :mod:`sage.geometry.triangulation.point_configuration` for more details. … … 34 40 from sage.structure.element import Element 35 41 from sage.rings.all import QQ, ZZ 36 42 from sage.modules.all import vector 37 43 from sage.misc.cachefunc import cached_method 38 44 39 45 40 46 ######################################################################## … … 53 59 A 2d graphics object. 54 60 55 61 EXAMPLES:: 56 62 57 63 sage: points = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[1,1]]) 58 64 sage: triang = points.triangulate() 59 sage: triang.plot(axes=False ) # indirect doctest65 sage: triang.plot(axes=False, aspect_ratio=1) # indirect doctest 60 66 """ 61 67 from sage.plot.all import point2d, line2d, arrow, polygon2d 62 68 points = [ point.reduced_affine() for point in triangulation.point_configuration() ] … … 115 121 A 3d graphics object. 116 122 117 123 EXAMPLES:: 118 124 119 125 sage: p = [[0,1,1],[0,0,1],[0,1,0], [1,1,1],[1,0,1],[1,1,0]] 120 126 sage: points = PointConfiguration(p) 121 127 sage: triang = points.triangulate() … … 241 247 :meth:`~sage.geometry.triangulation.point_configuration.PointConfiguration.bistellar_flips`. 242 248 243 249 EXAMPLES:: 244 250 245 251 sage: p = [[0,1],[0,0],[1,0]] 246 252 sage: points = PointConfiguration(p) 247 253 sage: from sage.geometry.triangulation.point_configuration import Triangulation … … 268 274 Returns the point configuration underlying the triangulation. 269 275 270 276 EXAMPLES:: 271 277 272 278 sage: pconfig = PointConfiguration([[0,0],[0,1],[1,0]]) 273 279 sage: pconfig 274 280 A point configuration in QQ^2 consisting of 3 points. The … … 331 337 332 338 EXAMPLES:: 333 339 340 sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM 334 341 sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[1,1]]) 335 342 sage: triangulation = pc.triangulate() 336 343 sage: iter = triangulation.__iter__() … … 360 367 A tuple of integers. The vertex indices of the ith simplex. 361 368 362 369 EXAMPLES:: 363 370 371 sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM 364 372 sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[1,1]]) 365 373 sage: triangulation = pc.triangulate() 366 374 sage: triangulation[1] … … 375 383 376 384 TESTS:: 377 385 386 sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM 378 387 sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[1,1]]) 379 388 sage: triangulation = pc.triangulations().next() 380 389 sage: triangulation.__len__() … … 391 400 392 401 TESTS:: 393 402 403 sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM 394 404 sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[1,1],[2,2]]) 395 405 sage: t = pc.triangulations() 396 406 sage: t.next()._repr_() … … 410 420 Produce a graphical representation of the triangulation. 411 421 412 422 EXAMPLES:: 413 423 414 424 sage: p = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[1,1]]) 415 425 sage: triangulation = p.triangulate() 416 426 sage: triangulation … … 433 443 r""" 434 444 Calculate the GKZ phi vector of the triangulation. 435 445 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 436 457 OUTPUT: 437 458 438 Vector  the phi vector of self.459 The phi vector of self. 439 460 440 461 EXAMPLES:: 441 462 … … 444 465 (3, 1, 5, 2, 4) 445 466 sage: p.lexicographic_triangulation().gkz_phi() 446 467 (1, 3, 4, 2, 5) 447 448 NOTE:449 450 For a definition of the phi vector, see [GKZ]_ page 220 equation 1.4.451 468 """ 452 469 vec = vector(ZZ, self.point_configuration().n_points()) 453 470 for simplex in self: … … 525 542 EXAMPLES:: 526 543 527 544 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 TOPCOM529 545 sage: triangulation = pc.triangulate() 530 546 sage: fan = triangulation.fan(); fan 531 547 Rational polyhedral fan in 2d lattice N … … 538 554 sage: interior=[(0, 1, 1), (0, 1, 2), (0, 2, 1), (0, 2, 2)] 539 555 sage: points = vertices+interior 540 556 sage: pc = PointConfiguration(points, fine=True) 541 sage: pc.set_engine('internal') # to make doctests independent of TOPCOM542 557 sage: triangulation = pc.triangulate() 543 558 sage: fan = triangulation.fan( (1,0,0) ) 544 559 sage: fan … … 556 571 return Fan(self, (vector(R, p)  origin for p in points)) 557 572 558 573 574 @cached_method 559 575 def simplicial_complex(self): 560 576 r""" 561 577 Return a simplicial complex from a triangulation of the point … … 583 599 return SimplicialComplex(vertex_set = vertex_set, 584 600 maximal_faces = self) 585 601 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 outwardfacing boundary simplices (of dimension `d1`) 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 inwardfacing boundary simplices (of dimension `d1`) 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 piecewiselinear function. This 706 support function is not unique, for example, you can scale it 707 by a positive constant. The set of all piecewiselinear 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 piecewiselinear 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 4d cone in 4d 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 3d cone in 3d 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 1991 1991 1992 1992 pushing_triangulation = placing_triangulation 1993 1993 1994 @cached_method 1995 def Gale_transform(self, points=None): 1996 r""" 1997 Return the Gale transform of ``self``. 1994 1998 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()