Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagesmc
Path: blob/master/src/sage/geometry/triangulation/element.py
8817 views
1
""""
2
A triangulation
3
4
In Sage, the
5
:class:`~sage.geometry.triangulation.point_configuration.PointConfiguration`
6
and :class:`Triangulation` satisfy a parent/element relationship. In
7
particular, each triangulation refers back to its point
8
configuration. If you want to triangulate a point configuration, you
9
should construct a point configuration first and then use one of its
10
methods to triangulate it according to your requirements. You should
11
never have to construct a :class:`Triangulation` object directly.
12
13
EXAMPLES:
14
15
First, we select the internal implementation for enumerating
16
triangulations::
17
18
sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM
19
20
Here is a simple example of how to triangulate a point configuration::
21
22
sage: p = [[0,-1,-1],[0,0,1],[0,1,0], [1,-1,-1],[1,0,1],[1,1,0]]
23
sage: points = PointConfiguration(p)
24
sage: triang = points.triangulate(); triang
25
(<0,1,2,5>, <0,1,3,5>, <1,3,4,5>)
26
sage: triang.plot(axes=False)
27
28
See :mod:`sage.geometry.triangulation.point_configuration` for more details.
29
"""
30
31
32
########################################################################
33
# Copyright (C) 2010 Volker Braun <[email protected]>
34
#
35
# Distributed under the terms of the GNU General Public License (GPL)
36
#
37
# http://www.gnu.org/licenses/
38
########################################################################
39
40
from sage.structure.element import Element
41
from sage.rings.all import QQ, ZZ
42
from sage.modules.all import vector
43
from sage.misc.cachefunc import cached_method
44
from sage.sets.set import Set
45
from sage.graphs.graph import Graph
46
47
48
########################################################################
49
def triangulation_render_2d(triangulation, **kwds):
50
r"""
51
Return a graphical representation of a 2-d triangulation.
52
53
INPUT:
54
55
- ``triangulation`` -- a :class:`Triangulation`.
56
57
- ``**kwds`` -- keywords that are passed on to the graphics primitives.
58
59
OUTPUT:
60
61
A 2-d graphics object.
62
63
EXAMPLES::
64
65
sage: points = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
66
sage: triang = points.triangulate()
67
sage: triang.plot(axes=False, aspect_ratio=1) # indirect doctest
68
"""
69
from sage.plot.all import point2d, line2d, arrow, polygon2d
70
points = [ point.reduced_affine() for point in triangulation.point_configuration() ]
71
coord = [ [p[0], p[1]] for p in points ]
72
plot_points = sum([ point2d(p,
73
zorder=2, pointsize=10, **kwds)
74
for p in coord ])
75
76
tmp_lines = []
77
for t in triangulation:
78
if len(t)>=2:
79
tmp_lines.append([t[0], t[1]])
80
if len(t)>=3:
81
tmp_lines.append([t[0], t[2]])
82
tmp_lines.append([t[1], t[2]])
83
all_lines = []
84
interior_lines = []
85
for l in tmp_lines:
86
if l not in all_lines:
87
all_lines.append(l)
88
else:
89
interior_lines.append(l)
90
exterior_lines = [ l for l in all_lines if not l in interior_lines ]
91
92
plot_interior_lines = sum([ line2d([ coord[l[0]], coord[l[1]] ],
93
zorder=1, rgbcolor=(0,1,0), **kwds)
94
for l in interior_lines ])
95
plot_exterior_lines = sum([ line2d([ coord[l[0]], coord[l[1]] ],
96
zorder=1, rgbcolor=(0,0,1), **kwds)
97
for l in exterior_lines ])
98
99
plot_triangs = sum([ polygon2d([coord[t[0]], coord[t[1]], coord[t[2]]],
100
zorder=0, rgbcolor=(0.8, 1, 0.8), **kwds)
101
for t in triangulation if len(t)>=3 ])
102
103
return \
104
plot_points + \
105
plot_interior_lines + plot_exterior_lines + \
106
plot_triangs
107
108
109
110
111
def triangulation_render_3d(triangulation, **kwds):
112
r"""
113
Return a graphical representation of a 3-d triangulation.
114
115
INPUT:
116
117
- ``triangulation`` -- a :class:`Triangulation`.
118
119
- ``**kwds`` -- keywords that are passed on to the graphics primitives.
120
121
OUTPUT:
122
123
A 3-d graphics object.
124
125
EXAMPLES::
126
127
sage: p = [[0,-1,-1],[0,0,1],[0,1,0], [1,-1,-1],[1,0,1],[1,1,0]]
128
sage: points = PointConfiguration(p)
129
sage: triang = points.triangulate()
130
sage: triang.plot(axes=False) # indirect doctest
131
"""
132
from sage.plot.plot3d.all import point3d, line3d, arrow3d, polygon3d
133
points = [ point.reduced_affine() for point in triangulation.point_configuration() ]
134
coord = [ [p[0], p[1], p[2] ] for p in points ]
135
plot_points = sum([ point3d(p, size=15,
136
**kwds)
137
for p in coord ])
138
139
tmp_lines = []
140
for t in triangulation:
141
if len(t)>=2:
142
tmp_lines.append([t[0], t[1]])
143
if len(t)>=3:
144
tmp_lines.append([t[0], t[2]])
145
tmp_lines.append([t[1], t[2]])
146
if len(t)>=4:
147
tmp_lines.append([t[0], t[3]])
148
tmp_lines.append([t[1], t[3]])
149
tmp_lines.append([t[2], t[3]])
150
all_lines = []
151
interior_lines = []
152
for l in tmp_lines:
153
if l not in all_lines:
154
all_lines.append(l)
155
else:
156
interior_lines.append(l)
157
exterior_lines = [ l for l in all_lines if not l in interior_lines ]
158
159
from sage.plot.plot3d.texture import Texture
160
line_int = Texture(color='darkblue', ambient=1, diffuse=0)
161
line_ext = Texture(color='green', ambient=1, diffuse=0)
162
triang_int = Texture(opacity=0.3, specular=0, shininess=0, diffuse=0, ambient=1, color='yellow')
163
triang_ext = Texture(opacity=0.6, specular=0, shininess=0, diffuse=0, ambient=1, color='green')
164
165
plot_interior_lines = sum([ line3d([ coord[l[0]], coord[l[1]] ],
166
thickness=2, texture=line_int, **kwds)
167
for l in interior_lines ])
168
plot_exterior_lines = sum([ line3d([ coord[l[0]], coord[l[1]] ],
169
thickness=3, texture=line_ext, **kwds)
170
for l in exterior_lines ])
171
172
tmp_triangs = []
173
for t in triangulation:
174
if len(t)>=3:
175
tmp_triangs.append([t[0], t[1], t[2]])
176
if len(t)>=4:
177
tmp_triangs.append([t[0], t[1], t[3]])
178
tmp_triangs.append([t[0], t[2], t[3]])
179
tmp_triangs.append([t[1], t[2], t[3]])
180
all_triangs = []
181
interior_triangs = []
182
for l in tmp_triangs:
183
if l not in all_triangs:
184
all_triangs.append(l)
185
else:
186
interior_triangs.append(l)
187
exterior_triangs = [ l for l in all_triangs if not l in interior_triangs ]
188
189
plot_interior_triangs = \
190
sum([ polygon3d([coord[t[0]], coord[t[1]], coord[t[2]]],
191
texture = triang_int, **kwds)
192
for t in interior_triangs ])
193
plot_exterior_triangs = \
194
sum([ polygon3d([coord[t[0]], coord[t[1]], coord[t[2]]],
195
texture = triang_ext, **kwds)
196
for t in exterior_triangs ])
197
198
return \
199
plot_points + \
200
plot_interior_lines + plot_exterior_lines + \
201
plot_interior_triangs + plot_exterior_triangs
202
203
204
205
206
207
########################################################################
208
class Triangulation(Element):
209
"""
210
A triangulation of a
211
:class:`~sage.geometry.triangulation.point_configuration.PointConfiguration`.
212
213
.. WARNING::
214
215
You should never create :class:`Triangulation` objects
216
manually. See
217
:meth:`~sage.geometry.triangulation.point_configuration.PointConfiguration.triangulate`
218
and
219
:meth:`~sage.geometry.triangulation.point_configuration.PointConfiguration.triangulations`
220
to triangulate point configurations.
221
"""
222
223
def __init__(self, triangulation, parent, check=True):
224
"""
225
The constructor of a ``Triangulation`` object. Note that an
226
internal reference to the underlying ``PointConfiguration`` is
227
kept.
228
229
INPUT:
230
231
- ``parent`` -- a
232
:class:`~sage.geometry.triangulation.point_configuration.PointConfiguration`
233
234
- ``triangulation`` -- an iterable of integers or iterable of
235
iterables (e.g. a list of lists). In the first case, the
236
integers specify simplices via
237
:meth:`PointConfiguration.simplex_to_int`. In the second
238
case, the point indices of the maximal simplices of the
239
triangulation.
240
241
- ``check`` -- boolean. Whether to perform checks that the
242
triangulation is, indeed, a triangulation of the point
243
configuration.
244
245
NOTE:
246
247
Passing ``check=False`` allows you to create triangulations of
248
subsets of the points of the configuration, see
249
:meth:`~sage.geometry.triangulation.point_configuration.PointConfiguration.bistellar_flips`.
250
251
EXAMPLES::
252
253
sage: p = [[0,1],[0,0],[1,0]]
254
sage: points = PointConfiguration(p)
255
sage: from sage.geometry.triangulation.point_configuration import Triangulation
256
sage: Triangulation([(0,1,2)], points)
257
(<0,1,2>)
258
sage: Triangulation([1], points)
259
(<0,1,2>)
260
"""
261
Element.__init__(self, parent=parent)
262
self._point_configuration = parent
263
264
try:
265
triangulation = tuple(sorted( tuple(sorted(t)) for t in triangulation))
266
except TypeError:
267
triangulation = tuple( self.point_configuration().int_to_simplex(i)
268
for i in triangulation )
269
assert not check or all( len(t)==self.point_configuration().dim()+1
270
for t in triangulation)
271
self._triangulation = triangulation
272
273
274
def point_configuration(self):
275
"""
276
Returns the point configuration underlying the triangulation.
277
278
EXAMPLES::
279
280
sage: pconfig = PointConfiguration([[0,0],[0,1],[1,0]])
281
sage: pconfig
282
A point configuration in QQ^2 consisting of 3 points. The
283
triangulations of this point configuration are assumed to
284
be connected, not necessarily fine, not necessarily regular.
285
sage: triangulation = pconfig.triangulate()
286
sage: triangulation
287
(<0,1,2>)
288
sage: triangulation.point_configuration()
289
A point configuration in QQ^2 consisting of 3 points. The
290
triangulations of this point configuration are assumed to
291
be connected, not necessarily fine, not necessarily regular.
292
sage: pconfig == triangulation.point_configuration()
293
True
294
"""
295
return self._point_configuration
296
297
298
def __cmp__(self, right):
299
r"""
300
Compare ``self`` and ``right``.
301
302
INPUT:
303
304
- ``right`` -- anything.
305
306
OUTPUT:
307
308
- 0 if ``right`` is the same triangulation as ``self``, 1 or
309
-1 otherwise.
310
311
TESTS::
312
313
sage: pc = PointConfiguration([[0,0],[0,1],[1,0]])
314
sage: t1 = pc.triangulate()
315
sage: from sage.geometry.triangulation.point_configuration import Triangulation
316
sage: t2 = Triangulation([[2,1,0]], pc)
317
sage: t1 is t2
318
False
319
sage: cmp(t1, t2)
320
0
321
sage: t1 == t2 # indirect doctest
322
True
323
sage: abs( cmp(t1, Triangulation(((0,1),(1,2)), pc, check=False) ))
324
1
325
sage: abs( cmp(t2, "not a triangulation") )
326
1
327
"""
328
left = self
329
c = cmp(isinstance(left,Triangulation), isinstance(right,Triangulation))
330
if c: return c
331
c = cmp(left.point_configuration(), right.point_configuration())
332
if c: return c
333
return cmp(left._triangulation, right._triangulation)
334
335
336
def __iter__(self):
337
"""
338
Iterate through the simplices of the triangulation.
339
340
EXAMPLES::
341
342
sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM
343
sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
344
sage: triangulation = pc.triangulate()
345
sage: iter = triangulation.__iter__()
346
sage: iter.next()
347
(1, 3, 4)
348
sage: iter.next()
349
(2, 3, 4)
350
sage: iter.next()
351
Traceback (most recent call last):
352
...
353
StopIteration
354
"""
355
for p in self._triangulation:
356
yield p
357
358
359
def __getitem__(self, i):
360
"""
361
Access the point indices of the i-th simplex of the triangulation.
362
363
INPUT:
364
365
- ``i`` -- integer. The index of a simplex.
366
367
OUTPUT:
368
369
A tuple of integers. The vertex indices of the i-th simplex.
370
371
EXAMPLES::
372
373
sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM
374
sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
375
sage: triangulation = pc.triangulate()
376
sage: triangulation[1]
377
(2, 3, 4)
378
"""
379
return self._triangulation[i]
380
381
382
def __len__(self):
383
"""
384
Returns the length of the triangulation.
385
386
TESTS::
387
388
sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM
389
sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
390
sage: triangulation = pc.triangulations().next()
391
sage: triangulation.__len__()
392
2
393
sage: len(triangulation) # equivalent
394
2
395
"""
396
return len(self._triangulation)
397
398
399
def _repr_(self):
400
r"""
401
Return a string representation.
402
403
TESTS::
404
405
sage: PointConfiguration.set_engine('internal') # to make doctests independent of TOPCOM
406
sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1],[2,2]])
407
sage: t = pc.triangulations()
408
sage: t.next()._repr_()
409
'(<1,4,5>, <2,4,5>)'
410
"""
411
#s = 'A triangulation'
412
#s += ' in QQ^'+str(self.point_configuration().ambient_dim())
413
#s += ' consisting of '+str(len(self))+' simplices.'
414
s = '('
415
s += ', '.join([ '<'+','.join(map(str,t))+'>' for t in self._triangulation])
416
s += ')'
417
return s
418
419
420
def plot(self, **kwds):
421
r"""
422
Produce a graphical representation of the triangulation.
423
424
EXAMPLES::
425
426
sage: p = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]])
427
sage: triangulation = p.triangulate()
428
sage: triangulation
429
(<1,3,4>, <2,3,4>)
430
sage: triangulation.plot(axes=False)
431
"""
432
dim = self.point_configuration().dim()
433
434
if dim == 2:
435
return triangulation_render_2d(self, **kwds)
436
437
if dim == 3:
438
return triangulation_render_3d(self, **kwds)
439
440
raise NotImplementedError, \
441
'Plotting '+str(dim)+'-dimensional triangulations not implemented!'
442
443
444
def gkz_phi(self):
445
r"""
446
Calculate the GKZ phi vector of the triangulation.
447
448
The phi vector is a vector of length equals to the number of
449
points in the point configuration. For a fixed triangulation
450
`T`, the entry corresponding to the `i`-th point `p_i` is
451
452
.. math::
453
454
\phi_T(p_i) = \sum_{t\in T, t\owns p_i} Vol(t)
455
456
that is, the total volume of all simplices containing `p_i`.
457
See also [GKZ]_ page 220 equation 1.4.
458
459
OUTPUT:
460
461
The phi vector of self.
462
463
EXAMPLES::
464
465
sage: p = PointConfiguration([[0,0],[1,0],[2,1],[1,2],[0,1]])
466
sage: p.triangulate().gkz_phi()
467
(3, 1, 5, 2, 4)
468
sage: p.lexicographic_triangulation().gkz_phi()
469
(1, 3, 4, 2, 5)
470
"""
471
vec = vector(ZZ, self.point_configuration().n_points())
472
for simplex in self:
473
vol = self.point_configuration().volume(simplex)
474
for i in simplex:
475
vec[i] = vec[i] + vol
476
return vec
477
478
479
def enumerate_simplices(self):
480
r"""
481
Return the enumerated simplices.
482
483
OUTPUT:
484
485
A tuple of integers that uniquely specifies the triangulation.
486
487
EXAMPLES::
488
489
sage: pc = PointConfiguration(matrix([
490
... [ 0, 0, 0, 0, 0, 2, 4,-1, 1, 1, 0, 0, 1, 0],
491
... [ 0, 0, 0, 1, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0],
492
... [ 0, 2, 0, 0, 0, 0,-1, 0, 1, 0, 1, 0, 0, 1],
493
... [ 0, 1, 1, 0, 0, 1, 0,-2, 1, 0, 0,-1, 1, 1],
494
... [ 0, 0, 0, 0, 1, 0,-1, 0, 0, 0, 0, 0, 0, 0]
495
... ]).columns())
496
sage: triangulation = pc.lexicographic_triangulation()
497
sage: triangulation.enumerate_simplices()
498
(1678, 1688, 1769, 1779, 1895, 1905, 2112, 2143, 2234, 2360, 2555, 2580,
499
2610, 2626, 2650, 2652, 2654, 2661, 2663, 2667, 2685, 2755, 2757, 2759,
500
2766, 2768, 2772, 2811, 2881, 2883, 2885, 2892, 2894, 2898)
501
502
You can recreate the triangulation from this list by passing
503
it to the constructor::
504
505
sage: from sage.geometry.triangulation.point_configuration import Triangulation
506
sage: Triangulation([1678, 1688, 1769, 1779, 1895, 1905, 2112, 2143,
507
... 2234, 2360, 2555, 2580, 2610, 2626, 2650, 2652, 2654, 2661, 2663,
508
... 2667, 2685, 2755, 2757, 2759, 2766, 2768, 2772, 2811, 2881, 2883,
509
... 2885, 2892, 2894, 2898], pc)
510
(<1,3,4,7,10,13>, <1,3,4,8,10,13>, <1,3,6,7,10,13>, <1,3,6,8,10,13>,
511
<1,4,6,7,10,13>, <1,4,6,8,10,13>, <2,3,4,6,7,12>, <2,3,4,7,12,13>,
512
<2,3,6,7,12,13>, <2,4,6,7,12,13>, <3,4,5,6,9,12>, <3,4,5,8,9,12>,
513
<3,4,6,7,11,12>, <3,4,6,9,11,12>, <3,4,7,10,11,13>, <3,4,7,11,12,13>,
514
<3,4,8,9,10,12>, <3,4,8,10,12,13>, <3,4,9,10,11,12>, <3,4,10,11,12,13>,
515
<3,5,6,8,9,12>, <3,6,7,10,11,13>, <3,6,7,11,12,13>, <3,6,8,9,10,12>,
516
<3,6,8,10,12,13>, <3,6,9,10,11,12>, <3,6,10,11,12,13>, <4,5,6,8,9,12>,
517
<4,6,7,10,11,13>, <4,6,7,11,12,13>, <4,6,8,9,10,12>, <4,6,8,10,12,13>,
518
<4,6,9,10,11,12>, <4,6,10,11,12,13>)
519
"""
520
pc = self._point_configuration
521
return tuple( pc.simplex_to_int(t) for t in self )
522
523
524
def fan(self, origin=None):
525
r"""
526
Construct the fan of cones over the simplices of the triangulation.
527
528
INPUT:
529
530
- ``origin`` -- ``None`` (default) or coordinates of a
531
point. The common apex of all cones of the fan. If ``None``,
532
the triangulation must be a star triangulation and the
533
distinguished central point is used as the origin.
534
535
OUTPUT:
536
537
A :class:`~sage.geometry.fan.RationalPolyhedralFan`. The
538
coordinates of the points are shifted so that the apex of the
539
fan is the origin of the coordinate system.
540
541
.. note:: If the set of cones over the simplices is not a fan, a
542
suitable exception is raised.
543
544
EXAMPLES::
545
546
sage: pc = PointConfiguration([(0,0), (1,0), (0,1), (-1,-1)], star=0, fine=True)
547
sage: triangulation = pc.triangulate()
548
sage: fan = triangulation.fan(); fan
549
Rational polyhedral fan in 2-d lattice N
550
sage: fan.is_equivalent( toric_varieties.P2().fan() )
551
True
552
553
Toric diagrams (the `\ZZ_5` hyperconifold)::
554
555
sage: vertices=[(0, 1, 0), (0, 3, 1), (0, 2, 3), (0, 0, 2)]
556
sage: interior=[(0, 1, 1), (0, 1, 2), (0, 2, 1), (0, 2, 2)]
557
sage: points = vertices+interior
558
sage: pc = PointConfiguration(points, fine=True)
559
sage: triangulation = pc.triangulate()
560
sage: fan = triangulation.fan( (-1,0,0) )
561
sage: fan
562
Rational polyhedral fan in 3-d lattice N
563
sage: fan.rays()
564
N(1, 1, 0),
565
N(1, 3, 1),
566
N(1, 2, 3),
567
N(1, 0, 2),
568
N(1, 1, 1),
569
N(1, 1, 2),
570
N(1, 2, 1),
571
N(1, 2, 2)
572
in 3-d lattice N
573
"""
574
from sage.geometry.fan import Fan
575
if origin is None:
576
origin = self.point_configuration().star_center()
577
R = self.base_ring()
578
origin = vector(R, origin)
579
points = self.point_configuration().points()
580
return Fan(self, (vector(R, p) - origin for p in points))
581
582
583
@cached_method
584
def simplicial_complex(self):
585
r"""
586
Return a simplicial complex from a triangulation of the point
587
configuration.
588
589
OUTPUT:
590
591
A :class:`~sage.homology.simplicial_complex.SimplicialComplex`.
592
593
EXAMPLES::
594
595
sage: p = polytopes.cuboctahedron()
596
sage: sc = p.triangulate(engine='internal').simplicial_complex()
597
sage: sc
598
Simplicial complex with 12 vertices and 16 facets
599
600
Any convex set is contractable, so its reduced homology groups vanish::
601
602
sage: sc.homology()
603
{0: 0, 1: 0, 2: 0, 3: 0}
604
"""
605
from sage.homology.simplicial_complex import SimplicialComplex
606
return SimplicialComplex(self)
607
608
609
@cached_method
610
def _boundary_simplex_dictionary(self):
611
"""
612
Return facets and the simplices they bound
613
614
TESTS::
615
616
sage: triangulation = polytopes.n_cube(2).triangulate(engine='internal')
617
sage: triangulation._boundary_simplex_dictionary()
618
{(0, 1): ((0, 1, 3),),
619
(0, 3): ((0, 1, 3), (0, 2, 3)),
620
(1, 3): ((0, 1, 3),),
621
(2, 3): ((0, 2, 3),),
622
(0, 2): ((0, 2, 3),)}
623
624
sage: triangulation = polytopes.n_cube(3).triangulate(engine='internal')
625
sage: triangulation._boundary_simplex_dictionary()
626
{(1, 4, 7): ((0, 1, 4, 7), (1, 4, 5, 7)),
627
(1, 3, 7): ((1, 2, 3, 7),),
628
(0, 1, 7): ((0, 1, 2, 7), (0, 1, 4, 7)),
629
(0, 2, 7): ((0, 1, 2, 7), (0, 2, 4, 7)),
630
(0, 1, 4): ((0, 1, 4, 7),),
631
(2, 4, 6): ((2, 4, 6, 7),),
632
(0, 1, 2): ((0, 1, 2, 7),),
633
(1, 2, 7): ((0, 1, 2, 7), (1, 2, 3, 7)),
634
(2, 6, 7): ((2, 4, 6, 7),),
635
(2, 3, 7): ((1, 2, 3, 7),),
636
(1, 4, 5): ((1, 4, 5, 7),),
637
(1, 5, 7): ((1, 4, 5, 7),),
638
(4, 5, 7): ((1, 4, 5, 7),),
639
(0, 4, 7): ((0, 1, 4, 7), (0, 2, 4, 7)),
640
(2, 4, 7): ((0, 2, 4, 7), (2, 4, 6, 7)),
641
(1, 2, 3): ((1, 2, 3, 7),),
642
(4, 6, 7): ((2, 4, 6, 7),),
643
(0, 2, 4): ((0, 2, 4, 7),)}
644
"""
645
result = dict()
646
for simplex in self:
647
for i in range(len(simplex)):
648
facet = simplex[:i] + simplex[i+1:]
649
result[facet] = result.get(facet, tuple()) + (simplex,)
650
return result
651
652
653
@cached_method
654
def boundary(self):
655
"""
656
Return the boundary of the triangulation.
657
658
OUTPUT:
659
660
The outward-facing boundary simplices (of dimension `d-1`) of
661
the `d`-dimensional triangulation as a set. Each boundary is
662
returned by a tuple of point indices.
663
664
EXAMPLES::
665
666
sage: triangulation = polytopes.n_cube(3).triangulate(engine='internal')
667
sage: triangulation
668
(<0,1,2,7>, <0,1,4,7>, <0,2,4,7>, <1,2,3,7>, <1,4,5,7>, <2,4,6,7>)
669
sage: triangulation.boundary()
670
frozenset([(1, 3, 7), (4, 5, 7), (1, 2, 3), (0, 1, 2), (2, 4, 6), (2, 6, 7),
671
(2, 3, 7), (1, 5, 7), (0, 1, 4), (1, 4, 5), (4, 6, 7), (0, 2, 4)])
672
sage: triangulation.interior_facets()
673
frozenset([(1, 4, 7), (1, 2, 7), (2, 4, 7), (0, 1, 7), (0, 4, 7), (0, 2, 7)])
674
"""
675
return frozenset(facet for facet, bounded_simplices
676
in self._boundary_simplex_dictionary().iteritems()
677
if len(bounded_simplices)==1)
678
679
@cached_method
680
def interior_facets(self):
681
"""
682
Return the interior facets of the triangulation.
683
684
OUTPUT:
685
686
The inward-facing boundary simplices (of dimension `d-1`) of
687
the `d`-dimensional triangulation as a set. Each boundary is
688
returned by a tuple of point indices.
689
690
EXAMPLES::
691
692
sage: triangulation = polytopes.n_cube(3).triangulate(engine='internal')
693
sage: triangulation
694
(<0,1,2,7>, <0,1,4,7>, <0,2,4,7>, <1,2,3,7>, <1,4,5,7>, <2,4,6,7>)
695
sage: triangulation.boundary()
696
frozenset([(1, 3, 7), (4, 5, 7), (1, 2, 3), (0, 1, 2), (2, 4, 6), (2, 6, 7),
697
(2, 3, 7), (1, 5, 7), (0, 1, 4), (1, 4, 5), (4, 6, 7), (0, 2, 4)])
698
sage: triangulation.interior_facets()
699
frozenset([(1, 4, 7), (1, 2, 7), (2, 4, 7), (0, 1, 7), (0, 4, 7), (0, 2, 7)])
700
"""
701
return frozenset(facet for facet, bounded_simplices
702
in self._boundary_simplex_dictionary().iteritems()
703
if len(bounded_simplices)==2)
704
705
@cached_method
706
def normal_cone(self):
707
r"""
708
Return the (closure of the) normal cone of the triangulation.
709
710
Recall that a regular triangulation is one that equals the
711
"crease lines" of a convex piecewise-linear function. This
712
support function is not unique, for example, you can scale it
713
by a positive constant. The set of all piecewise-linear
714
functions with fixed creases forms an open cone. This cone can
715
be interpreted as the cone of normal vectors at a point of the
716
secondary polytope, which is why we call it normal cone. See
717
[GKZ]_ Section 7.1 for details.
718
719
OUTPUT:
720
721
The closure of the normal cone. The `i`-th entry equals the
722
value of the piecewise-linear function at the `i`-th point of
723
the configuration.
724
725
For an irregular triangulation, the normal cone is empty. In
726
this case, a single point (the origin) is returned.
727
728
EXAMPLES::
729
730
sage: triangulation = polytopes.n_cube(2).triangulate(engine='internal')
731
sage: triangulation
732
(<0,1,3>, <0,2,3>)
733
sage: N = triangulation.normal_cone(); N
734
4-d cone in 4-d lattice
735
sage: N.rays()
736
(-1, 0, 0, 0),
737
( 1, 0, 1, 0),
738
(-1, 0, -1, 0),
739
( 1, 0, 0, -1),
740
(-1, 0, 0, 1),
741
( 1, 1, 0, 0),
742
(-1, -1, 0, 0)
743
in Ambient free module of rank 4
744
over the principal ideal domain Integer Ring
745
sage: N.dual().rays()
746
(-1, 1, 1, -1)
747
in Ambient free module of rank 4
748
over the principal ideal domain Integer Ring
749
750
TESTS::
751
752
sage: polytopes.n_simplex(2).triangulate().normal_cone()
753
3-d cone in 3-d lattice
754
sage: _.dual().is_trivial()
755
True
756
"""
757
if not self.point_configuration().base_ring().is_subring(QQ):
758
raise NotImplementedError('Only base rings ZZ and QQ are supported')
759
from sage.libs.ppl import Variable, Constraint, Constraint_System, Linear_Expression, C_Polyhedron
760
from sage.matrix.constructor import matrix
761
from sage.misc.misc import uniq
762
from sage.rings.arith import lcm
763
pc = self.point_configuration()
764
cs = Constraint_System()
765
for facet in self.interior_facets():
766
s0, s1 = self._boundary_simplex_dictionary()[facet]
767
p = set(s0).difference(facet).pop()
768
q = set(s1).difference(facet).pop()
769
origin = pc.point(p).reduced_affine_vector()
770
base_indices = [ i for i in s0 if i!=p ]
771
base = matrix([ pc.point(i).reduced_affine_vector()-origin for i in base_indices ])
772
sol = base.solve_left( pc.point(q).reduced_affine_vector()-origin )
773
relation = [0]*pc.n_points()
774
relation[p] = sum(sol)-1
775
relation[q] = 1
776
for i, base_i in enumerate(base_indices):
777
relation[base_i] = -sol[i]
778
rel_denom = lcm([QQ(r).denominator() for r in relation])
779
relation = [ ZZ(r*rel_denom) for r in relation ]
780
ex = Linear_Expression(relation,0)
781
cs.insert(ex >= 0)
782
from sage.modules.free_module import FreeModule
783
ambient = FreeModule(ZZ, self.point_configuration().n_points())
784
if cs.empty():
785
cone = C_Polyhedron(ambient.dimension(), 'universe')
786
else:
787
cone = C_Polyhedron(cs)
788
from sage.geometry.cone import _Cone_from_PPL
789
return _Cone_from_PPL(cone, lattice=ambient)
790
791
def adjacency_graph(self):
792
"""
793
Returns a graph showing which simplices are adjacent in the
794
triangulation
795
796
OUTPUT:
797
798
A graph consisting of vertices referring to the simplices in the
799
triangulation, and edges showing which simplices are adjacent to each
800
other.
801
802
.. SEEALSO::
803
804
* To obtain the triangulation's 1-skeleton, use
805
:meth:`SimplicialComplex.graph` through
806
``MyTriangulation.simplicial_complex().graph()``.
807
808
AUTHORS:
809
810
* Stephen Farley (2013-08-10): initial version
811
812
EXAMPLES::
813
814
sage: p = PointConfiguration([[1,0,0], [0,1,0], [0,0,1], [-1,0,1],
815
....: [1,0,-1], [-1,0,0], [0,-1,0], [0,0,-1]])
816
sage: t = p.triangulate()
817
sage: t.adjacency_graph()
818
Graph on 8 vertices
819
820
"""
821
vertices = map(Set,list(self))
822
return Graph([vertices,
823
lambda x,y: len(x-y)==1])
824
825
826