This repository has been archived by the owner on Jan 30, 2023. It is now read-only.
/
base.py
10700 lines (8504 loc) · 400 KB
/
base.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
r"""
Base class for polyhedra
"""
# ****************************************************************************
# Copyright (C) 2008-2012 Marshall Hampton <hamptonio@gmail.com>
# Copyright (C) 2011-2015 Volker Braun <vbraun.name@gmail.com>
# Copyright (C) 2012-2018 Frederic Chapoton
# Copyright (C) 2013 Andrey Novoseltsev
# Copyright (C) 2014-2017 Moritz Firsching
# Copyright (C) 2014-2019 Thierry Monteil
# Copyright (C) 2015 Nathann Cohen
# Copyright (C) 2015-2017 Jeroen Demeyer
# Copyright (C) 2015-2017 Vincent Delecroix
# Copyright (C) 2015-2018 Dima Pasechnik
# Copyright (C) 2015-2020 Jean-Philippe Labbe <labbe at math.huji.ac.il>
# Copyright (C) 2015-2021 Matthias Koeppe
# Copyright (C) 2016-2019 Daniel Krenn
# Copyright (C) 2017 Marcelo Forets
# Copyright (C) 2017-2018 Mark Bell
# Copyright (C) 2019 Julian Ritter
# Copyright (C) 2019-2020 Laith Rastanawi
# Copyright (C) 2019-2020 Sophia Elia
# Copyright (C) 2019-2021 Jonathan Kliem <jonathan.kliem@fu-berlin.de>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
# https://www.gnu.org/licenses/
# ****************************************************************************
from dataclasses import dataclass
from typing import Any
import itertools
from sage.structure.element import Element, coerce_binop, is_Vector, is_Matrix
from sage.structure.richcmp import rich_to_bool, op_NE
from sage.cpython.string import bytes_to_str
from sage.misc.all import cached_method, prod
from sage.misc.randstate import current_randstate
from sage.misc.superseded import deprecated_function_alias
from sage.rings.all import QQ, ZZ, AA
from sage.rings.real_double import RDF
from sage.modules.free_module_element import vector
from sage.modules.vector_space_morphism import linear_transformation
from sage.matrix.constructor import matrix
from sage.functions.other import sqrt, floor, ceil
from sage.groups.matrix_gps.finitely_generated import MatrixGroup
from sage.graphs.graph import Graph
from sage.geometry.convex_set import ConvexSet_closed
from .constructor import Polyhedron
from sage.geometry.relative_interior import RelativeInterior
from sage.categories.sets_cat import EmptySetError
#########################################################################
# Notes if you want to implement your own backend:
#
# * derive from Polyhedron_base
#
# * you must implement _init_from_Vrepresentation and
# _init_from_Hrepresentation
#
# * You might want to override _init_empty_polyhedron
#
# * You may implement _init_from_Vrepresentation_and_Hrepresentation
#
# * You can of course also override any other method for which you
# have a faster implementation.
#########################################################################
#########################################################################
def is_Polyhedron(X):
"""
Test whether ``X`` is a Polyhedron.
INPUT:
- ``X`` -- anything.
OUTPUT:
Boolean.
EXAMPLES::
sage: p = polytopes.hypercube(2)
sage: from sage.geometry.polyhedron.base import is_Polyhedron
sage: is_Polyhedron(p)
True
sage: is_Polyhedron(123456)
False
"""
return isinstance(X, Polyhedron_base)
#########################################################################
class Polyhedron_base(Element, ConvexSet_closed):
"""
Base class for Polyhedron objects
INPUT:
- ``parent`` -- the parent, an instance of
:class:`~sage.geometry.polyhedron.parent.Polyhedra`.
- ``Vrep`` -- a list ``[vertices, rays, lines]`` or ``None``. The
V-representation of the polyhedron. If ``None``, the polyhedron
is determined by the H-representation.
- ``Hrep`` -- a list ``[ieqs, eqns]`` or ``None``. The
H-representation of the polyhedron. If ``None``, the polyhedron
is determined by the V-representation.
- ``Vrep_minimal`` (optional) -- see below
- ``Hrep_minimal`` (optional) -- see below
- ``pref_rep`` -- string (default: ``None``);
one of``Vrep`` or ``Hrep`` to pick this in case the backend
cannot initialize from complete double description
If both ``Vrep`` and ``Hrep`` are provided, then
``Vrep_minimal`` and ``Hrep_minimal`` must be set to ``True``.
TESTS::
sage: p = Polyhedron()
sage: TestSuite(p).run()
::
sage: p = Polyhedron(vertices=[(1,0), (0,1)], rays=[(1,1)], base_ring=ZZ)
sage: TestSuite(p).run()
::
sage: p=polytopes.flow_polytope(digraphs.DeBruijn(3,2))
sage: TestSuite(p).run()
::
sage: TestSuite(Polyhedron([[]])).run()
sage: TestSuite(Polyhedron([[0]])).run()
sage: TestSuite(Polyhedron([[1]])).run()
::
sage: P = polytopes.permutahedron(3) * Polyhedron(rays=[[0,0,1],[0,1,1],[1,2,3]])
sage: TestSuite(P).run()
::
sage: P = polytopes.permutahedron(3)*Polyhedron(rays=[[0,0,1],[0,1,1]], lines=[[1,0,0]])
sage: TestSuite(P).run()
::
sage: M = random_matrix(ZZ, 5, 5, distribution='uniform')
sage: while True:
....: M = random_matrix(ZZ, 5, 5, distribution='uniform')
....: if M.rank() != 5:
....: break
....:
sage: P = Polyhedron(M)
sage: TestSuite(P).run()
"""
def __init__(self, parent, Vrep, Hrep, Vrep_minimal=None, Hrep_minimal=None, pref_rep=None, **kwds):
"""
Initializes the polyhedron.
See :class:`Polyhedron_base` for a description of the input
data.
TESTS::
sage: p = Polyhedron() # indirect doctests
sage: from sage.geometry.polyhedron.backend_field import Polyhedron_field
sage: from sage.geometry.polyhedron.parent import Polyhedra_field
sage: parent = Polyhedra_field(AA, 1, 'field')
sage: Vrep = [[[0], [1/2], [1]], [], []]
sage: Hrep = [[[0, 1], [1, -1]], []]
sage: p = Polyhedron_field(parent, Vrep, Hrep,
....: Vrep_minimal=False, Hrep_minimal=True)
Traceback (most recent call last):
...
ValueError: if both Vrep and Hrep are provided, they must be minimal...
Illustration of ``pref_rep``.
Note that ``ppl`` doesn't support precomputed data::
sage: from sage.geometry.polyhedron.backend_ppl import Polyhedron_QQ_ppl
sage: from sage.geometry.polyhedron.parent import Polyhedra_QQ_ppl
sage: parent = Polyhedra_QQ_ppl(QQ, 1, 'ppl')
sage: p = Polyhedron_QQ_ppl(parent, Vrep, 'nonsense',
....: Vrep_minimal=True, Hrep_minimal=True, pref_rep='Vrep')
sage: p = Polyhedron_QQ_ppl(parent, 'nonsense', Hrep,
....: Vrep_minimal=True, Hrep_minimal=True, pref_rep='Hrep')
sage: p = Polyhedron_QQ_ppl(parent, 'nonsense', Hrep,
....: Vrep_minimal=True, Hrep_minimal=True, pref_rep='Vrepresentation')
Traceback (most recent call last):
...
ValueError: ``pref_rep`` must be one of ``(None, 'Vrep', 'Hrep')``
If the backend supports precomputed data, ``pref_rep`` is ignored::
sage: p = Polyhedron_field(parent, Vrep, 'nonsense', # py3
....: Vrep_minimal=True, Hrep_minimal=True, pref_rep='Vrep')
Traceback (most recent call last):
...
TypeError: _init_Hrepresentation() takes 3 positional arguments but 9 were given
sage: p = Polyhedron_field(parent, Vrep, 'nonsense', # py2
....: Vrep_minimal=True, Hrep_minimal=True, pref_rep='Vrep')
Traceback (most recent call last):
...
TypeError: _init_Hrepresentation() takes exactly 3 arguments (9 given)
The empty polyhedron is detected when the Vrepresentation is given with generator;
see :trac:`29899`::
sage: from sage.geometry.polyhedron.backend_cdd import Polyhedron_QQ_cdd
sage: from sage.geometry.polyhedron.parent import Polyhedra_QQ_cdd
sage: parent = Polyhedra_QQ_cdd(QQ, 0, 'cdd')
sage: p = Polyhedron_QQ_cdd(parent, [iter([]), iter([]), iter([])], None)
"""
Element.__init__(self, parent=parent)
if Vrep is not None and Hrep is not None:
if not (Vrep_minimal is True and Hrep_minimal is True):
raise ValueError("if both Vrep and Hrep are provided, they must be minimal"
" and Vrep_minimal and Hrep_minimal must both be True")
if hasattr(self, "_init_from_Vrepresentation_and_Hrepresentation"):
self._init_from_Vrepresentation_and_Hrepresentation(Vrep, Hrep)
return
else:
if pref_rep is None:
# Initialize from Hrepresentation if this seems simpler.
Vrep = [tuple(Vrep[0]), tuple(Vrep[1]), Vrep[2]]
Hrep = [tuple(Hrep[0]), Hrep[1]]
if len(Hrep[0]) < len(Vrep[0]) + len(Vrep[1]):
pref_rep = 'Hrep'
else:
pref_rep = 'Vrep'
if pref_rep == 'Vrep':
Hrep = None
elif pref_rep == 'Hrep':
Vrep = None
else:
raise ValueError("``pref_rep`` must be one of ``(None, 'Vrep', 'Hrep')``")
if Vrep is not None:
vertices, rays, lines = Vrep
# We build tuples out of generators now to detect the empty polyhedron.
# The damage is limited:
# The backend will have to obtain all elements from the generator anyway.
# The generators are mainly for saving time with initializing from
# Vrepresentation and Hrepresentation.
# If we dispose of one of them (see above), it is wasteful to have generated it.
# E.g. the dilate will be set up with new Vrepresentation and Hrepresentation
# regardless of the backend along with the argument ``pref_rep``.
# As we only use generators, there is no penalty to this approach
# (and the method ``dilation`` does not have to distinguish by backend).
if not isinstance(vertices, (tuple, list)):
vertices = tuple(vertices)
if not isinstance(rays, (tuple, list)):
rays = tuple(rays)
if not isinstance(lines, (tuple, list)):
lines = tuple(lines)
if vertices or rays or lines:
self._init_from_Vrepresentation(vertices, rays, lines, **kwds)
else:
self._init_empty_polyhedron()
elif Hrep is not None:
ieqs, eqns = Hrep
self._init_from_Hrepresentation(ieqs, eqns, **kwds)
else:
self._init_empty_polyhedron()
def __hash__(self):
r"""
TESTS::
sage: K.<a> = QuadraticField(2)
sage: p = Polyhedron(vertices=[(0,1,a),(3,a,5)],
....: rays=[(a,2,3), (0,0,1)],
....: base_ring=K)
sage: q = Polyhedron(vertices=[(3,a,5),(0,1,a)],
....: rays=[(0,0,1), (a,2,3)],
....: base_ring=K)
sage: hash(p) == hash(q)
True
"""
# TODO: find something better *but* fast
return hash((self.dim(),
self.ambient_dim(),
self.n_Hrepresentation(),
self.n_Vrepresentation(),
self.n_equations(),
self.n_facets(),
self.n_inequalities(),
self.n_lines(),
self.n_rays(),
self.n_vertices()))
def _sage_input_(self, sib, coerced):
"""
Return Sage command to reconstruct ``self``.
See :mod:`sage.misc.sage_input` for details.
.. TODO::
Add the option ``preparse`` to the method.
EXAMPLES::
sage: P = Polyhedron(vertices = [[1, 0], [0, 1]], rays = [[1, 1]], backend='ppl')
sage: sage_input(P)
Polyhedron(backend='ppl', base_ring=QQ, rays=[(QQ(1), QQ(1))], vertices=[(QQ(0), QQ(1)), (QQ(1), QQ(0))])
sage: P = Polyhedron(vertices = [[1, 0], [0, 1]], rays = [[1, 1]], backend='normaliz') # optional - pynormaliz
sage: sage_input(P) # optional - pynormaliz
Polyhedron(backend='normaliz', base_ring=QQ, rays=[(QQ(1), QQ(1))], vertices=[(QQ(0), QQ(1)), (QQ(1), QQ(0))])
sage: P = Polyhedron(vertices = [[1, 0], [0, 1]], rays = [[1, 1]], backend='polymake') # optional - polymake
sage: sage_input(P) # optional - polymake
Polyhedron(backend='polymake', base_ring=QQ, rays=[(QQ(1), QQ(1))], vertices=[(QQ(1), QQ(0)), (QQ(0), QQ(1))])
"""
kwds = dict()
kwds['base_ring'] = sib(self.base_ring())
kwds['backend'] = sib(self.backend())
if self.n_vertices() > 0:
kwds['vertices'] = [sib(tuple(v)) for v in self.vertices()]
if self.n_rays() > 0:
kwds['rays'] = [sib(tuple(r)) for r in self.rays()]
if self.n_lines() > 0:
kwds['lines'] = [sib(tuple(l)) for l in self.lines()]
return sib.name('Polyhedron')(**kwds)
def _init_from_Vrepresentation(self, vertices, rays, lines, **kwds):
"""
Construct polyhedron from V-representation data.
INPUT:
- ``vertices`` -- list of point. Each point can be specified
as any iterable container of
:meth:`~sage.geometry.polyhedron.base.base_ring` elements.
- ``rays`` -- list of rays. Each ray can be specified as any
iterable container of
:meth:`~sage.geometry.polyhedron.base.base_ring` elements.
- ``lines`` -- list of lines. Each line can be specified as
any iterable container of
:meth:`~sage.geometry.polyhedron.base.base_ring` elements.
EXAMPLES::
sage: p = Polyhedron()
sage: from sage.geometry.polyhedron.base import Polyhedron_base
sage: Polyhedron_base._init_from_Vrepresentation(p, [], [], [])
Traceback (most recent call last):
...
NotImplementedError: a derived class must implement this method
"""
raise NotImplementedError('a derived class must implement this method')
def _init_from_Hrepresentation(self, ieqs, eqns, **kwds):
"""
Construct polyhedron from H-representation data.
INPUT:
- ``ieqs`` -- list of inequalities. Each line can be specified
as any iterable container of
:meth:`~sage.geometry.polyhedron.base.base_ring` elements.
- ``eqns`` -- list of equalities. Each line can be specified
as any iterable container of
:meth:`~sage.geometry.polyhedron.base.base_ring` elements.
EXAMPLES::
sage: p = Polyhedron()
sage: from sage.geometry.polyhedron.base import Polyhedron_base
sage: Polyhedron_base._init_from_Hrepresentation(p, [], [])
Traceback (most recent call last):
...
NotImplementedError: a derived class must implement this method
"""
raise NotImplementedError('a derived class must implement this method')
def _init_empty_polyhedron(self):
"""
Initializes an empty polyhedron.
TESTS::
sage: empty = Polyhedron(); empty
The empty polyhedron in ZZ^0
sage: empty.Vrepresentation()
()
sage: empty.Hrepresentation()
(An equation -1 == 0,)
sage: Polyhedron(vertices = [])
The empty polyhedron in ZZ^0
sage: Polyhedron(vertices = [])._init_empty_polyhedron()
sage: from sage.geometry.polyhedron.parent import Polyhedra
sage: Polyhedra(QQ,7)()
A 0-dimensional polyhedron in QQ^7 defined as the convex hull of 1 vertex
"""
self._Vrepresentation = []
self._Hrepresentation = []
self.parent()._make_Equation(self, [-1] + [0]*self.ambient_dim())
self._Vrepresentation = tuple(self._Vrepresentation)
self._Hrepresentation = tuple(self._Hrepresentation)
V_matrix = matrix(ZZ, 0, 0, 0)
V_matrix.set_immutable()
self.vertex_adjacency_matrix.set_cache(V_matrix)
H_matrix = matrix(ZZ, 1, 1, 0)
H_matrix.set_immutable()
self.facet_adjacency_matrix.set_cache(H_matrix)
def _facet_adjacency_matrix(self):
"""
Compute the facet adjacency matrix in case it has not been
computed during initialization.
EXAMPLES::
sage: p = Polyhedron(vertices=[(0,0),(1,0),(0,1)])
sage: p._facet_adjacency_matrix()
[0 1 1]
[1 0 1]
[1 1 0]
Checks that :trac:`22455` is fixed::
sage: s = polytopes.simplex(2)
sage: s._facet_adjacency_matrix()
[0 1 1]
[1 0 1]
[1 1 0]
"""
# TODO: This implementation computes the whole face lattice,
# which is much more information than necessary.
M = matrix(ZZ, self.n_facets(), self.n_facets(), 0)
codim = self.ambient_dim()-self.dim()
def set_adjacent(h1, h2):
if h1 is h2:
return
i = h1.index() - codim
j = h2.index() - codim
M[i, j] = 1
M[j, i] = 1
for face in self.faces(self.dim()-2):
Hrep = face.ambient_Hrepresentation()
assert(len(Hrep) == codim+2)
set_adjacent(Hrep[-2], Hrep[-1])
M.set_immutable()
return M
def _vertex_adjacency_matrix(self):
"""
Compute the vertex adjacency matrix in case it has not been
computed during initialization.
EXAMPLES::
sage: p = Polyhedron(vertices=[(0,0),(1,0),(0,1)])
sage: p._vertex_adjacency_matrix()
[0 1 1]
[1 0 1]
[1 1 0]
"""
# TODO: This implementation computes the whole face lattice,
# which is much more information than necessary.
M = matrix(ZZ, self.n_Vrepresentation(), self.n_Vrepresentation(), 0)
def set_adjacent(v1, v2):
if v1 is v2:
return
i = v1.index()
j = v2.index()
M[i, j] = 1
M[j, i] = 1
face_lattice = self.face_lattice()
for face in face_lattice:
Vrep = face.ambient_Vrepresentation()
if len(Vrep) == 2:
set_adjacent(Vrep[0], Vrep[1])
M.set_immutable()
return M
def _delete(self):
"""
Delete this polyhedron.
This speeds up creation of new polyhedra by reusing
objects. After recycling a polyhedron object, it is not in a
consistent state any more and neither the polyhedron nor its
H/V-representation objects may be used any more.
.. SEEALSO::
:meth:`~sage.geometry.polyhedron.parent.Polyhedra_base.recycle`
EXAMPLES::
sage: p = Polyhedron([(0,0),(1,0),(0,1)])
sage: p._delete()
sage: vertices = [(0,0,0,0),(1,0,0,0),(0,1,0,0),(1,1,0,0),(0,0,1,0),(0,0,0,1)]
sage: def loop_polyhedra():
....: for i in range(100):
....: p = Polyhedron(vertices)
sage: timeit('loop_polyhedra()') # not tested - random
5 loops, best of 3: 79.5 ms per loop
sage: def loop_polyhedra_with_recycling():
....: for i in range(100):
....: p = Polyhedron(vertices)
....: p._delete()
sage: timeit('loop_polyhedra_with_recycling()') # not tested - random
5 loops, best of 3: 57.3 ms per loop
"""
self.parent().recycle(self)
def _test_basic_properties(self, tester=None, **options):
"""
Run some basic tests to see, that some general assertion on polyhedra hold.
TESTS::
sage: polytopes.cross_polytope(3)._test_basic_properties()
"""
if tester is None:
tester = self._tester(**options)
tester.assertEqual(self.n_vertices() + self.n_rays() + self.n_lines(), self.n_Vrepresentation())
tester.assertEqual(self.n_inequalities() + self.n_equations(), self.n_Hrepresentation())
if self.n_vertices():
# Depending on the backend, this does not hold for the empty polyhedron.
tester.assertEqual(self.dim() + self.n_equations(), self.ambient_dim())
tester.assertTrue(all(len(v[::]) == self.ambient_dim() for v in self.Vrep_generator()))
tester.assertTrue(all(len(h[::]) == self.ambient_dim() + 1 for h in self.Hrep_generator()))
if self.n_vertices() + self.n_rays() < 40:
tester.assertEqual(self, Polyhedron(vertices=self.vertices(), rays=self.rays(), lines=self.lines(), ambient_dim=self.ambient_dim()))
if self.n_inequalities() < 40:
tester.assertEqual(self, Polyhedron(ieqs=self.inequalities(), eqns=self.equations(), ambient_dim=self.ambient_dim()))
def base_extend(self, base_ring, backend=None):
"""
Return a new polyhedron over a larger base ring.
This method can also be used to change the backend.
INPUT:
- ``base_ring`` -- the new base ring
- ``backend`` -- the new backend, see
:func:`~sage.geometry.polyhedron.constructor.Polyhedron`.
If ``None`` (the default), attempt to keep the same backend.
Otherwise, use the same defaulting behavior
as described there.
OUTPUT:
The same polyhedron, but over a larger base ring and possibly with a changed backend.
EXAMPLES::
sage: P = Polyhedron(vertices=[(1,0), (0,1)], rays=[(1,1)], base_ring=ZZ); P
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 2 vertices and 1 ray
sage: P.base_extend(QQ)
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 2 vertices and 1 ray
sage: P.base_extend(QQ) == P
True
TESTS:
Test that :trac:`22575` is fixed::
sage: Q = P.base_extend(ZZ, backend='field')
sage: Q.backend()
'field'
"""
new_parent = self.parent().base_extend(base_ring, backend)
return new_parent(self)
def change_ring(self, base_ring, backend=None):
"""
Return the polyhedron obtained by coercing the entries of the
vertices/lines/rays of this polyhedron into the given ring.
This method can also be used to change the backend.
INPUT:
- ``base_ring`` -- the new base ring
- ``backend`` -- the new backend or ``None`` (default), see
:func:`~sage.geometry.polyhedron.constructor.Polyhedron`.
If ``None`` (the default), attempt to keep the same backend.
Otherwise, use the same defaulting behavior
as described there.
EXAMPLES::
sage: P = Polyhedron(vertices=[(1,0), (0,1)], rays=[(1,1)], base_ring=QQ); P
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 2 vertices and 1 ray
sage: P.change_ring(ZZ)
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 2 vertices and 1 ray
sage: P.change_ring(ZZ) == P
True
sage: P = Polyhedron(vertices=[(-1.3,0), (0,2.3)], base_ring=RDF); P.vertices()
(A vertex at (-1.3, 0.0), A vertex at (0.0, 2.3))
sage: P.change_ring(QQ).vertices()
(A vertex at (-13/10, 0), A vertex at (0, 23/10))
sage: P == P.change_ring(QQ)
True
sage: P.change_ring(ZZ)
Traceback (most recent call last):
...
TypeError: cannot change the base ring to the Integer Ring
sage: P = polytopes.regular_polygon(3); P
A 2-dimensional polyhedron in AA^2 defined as the convex hull of 3 vertices
sage: P.vertices()
(A vertex at (0.?e-16, 1.000000000000000?),
A vertex at (0.866025403784439?, -0.500000000000000?),
A vertex at (-0.866025403784439?, -0.500000000000000?))
sage: P.change_ring(QQ)
Traceback (most recent call last):
...
TypeError: cannot change the base ring to the Rational Field
.. WARNING::
The base ring ``RDF`` should be used with care. As it is
not an exact ring, certain computations may break or
silently produce wrong results, for example changing the
base ring from an exact ring into ``RDF`` may cause a
loss of data::
sage: P = Polyhedron([[2/3,0],[6666666666666667/10^16,0]], base_ring=AA); P
A 1-dimensional polyhedron in AA^2 defined as the convex hull of 2 vertices
sage: Q = P.change_ring(RDF); Q
A 0-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex
sage: P.n_vertices() == Q.n_vertices()
False
"""
from sage.categories.all import Rings
if base_ring not in Rings:
raise ValueError("invalid base ring")
try:
vertices = [[base_ring(x) for x in vertex] for vertex in self.vertices_list()]
rays = [[base_ring(x) for x in ray] for ray in self.rays_list()]
lines = [[base_ring(x) for x in line] for line in self.lines_list()]
except (TypeError, ValueError):
raise TypeError("cannot change the base ring to the {0}".format(base_ring))
new_parent = self.parent().change_ring(base_ring, backend)
return new_parent([vertices, rays, lines], None)
def _richcmp_(self, other, op):
"""
Compare ``self`` and ``other``.
INPUT:
- ``other`` -- a polyhedron
OUTPUT:
If ``other`` is a polyhedron, then the comparison
operator "less or equal than" means "is contained in", and
"less than" means "is strictly contained in".
EXAMPLES::
sage: P = Polyhedron(vertices=[(1,0), (0,1)], rays=[(1,1)])
sage: Q = Polyhedron(vertices=[(1,0), (0,1)])
sage: P >= Q
True
sage: Q <= P
True
sage: P == P
True
The polytope ``Q`` is strictly contained in ``P``::
sage: P > Q
True
sage: P < Q
False
sage: P == Q
False
"""
if self._Vrepresentation is None or other._Vrepresentation is None:
raise RuntimeError('some V representation is missing')
# make sure deleted polyhedra are not used in cache
if self.ambient_dim() != other.ambient_dim():
return op == op_NE
c0 = self._is_subpolyhedron(other)
c1 = other._is_subpolyhedron(self)
if c0 and c1:
return rich_to_bool(op, 0)
if c0:
return rich_to_bool(op, -1)
else:
return rich_to_bool(op, 1)
@coerce_binop
def _is_subpolyhedron(self, other):
"""
Test whether ``self`` is a (not necessarily strict)
sub-polyhedron of ``other``.
INPUT:
- ``other`` -- a :class:`Polyhedron`
OUTPUT:
Boolean
EXAMPLES::
sage: P = Polyhedron(vertices=[(1,0), (0,1)], rays=[(1,1)])
sage: Q = Polyhedron(vertices=[(1,0), (0,1)])
sage: P._is_subpolyhedron(Q)
False
sage: Q._is_subpolyhedron(P)
True
"""
return all(other_H.contains(self_V)
for other_H in other.Hrepresentation()
for self_V in self.Vrepresentation())
@cached_method
def vertex_facet_graph(self, labels=True):
r"""
Return the vertex-facet graph.
This function constructs a directed bipartite graph.
The nodes of the graph correspond to the vertices of the polyhedron
and the facets of the polyhedron. There is an directed edge
from a vertex to a face if and only if the vertex is incident to the face.
INPUT:
- ``labels`` -- boolean (default: ``True``); decide how the nodes
of the graph are labelled. Either with the original vertices/facets
of the Polyhedron or with integers.
OUTPUT:
- a bipartite DiGraph. If ``labels`` is ``True``, then the nodes
of the graph will actually be the vertices and facets of ``self``,
otherwise they will be integers.
.. SEEALSO::
:meth:`combinatorial_automorphism_group`,
:meth:`is_combinatorially_isomorphic`.
EXAMPLES::
sage: P = polytopes.cube()
sage: G = P.vertex_facet_graph(); G
Digraph on 14 vertices
sage: G.vertices(key = lambda v: str(v))
[A vertex at (-1, -1, -1),
A vertex at (-1, -1, 1),
A vertex at (-1, 1, -1),
A vertex at (-1, 1, 1),
A vertex at (1, -1, -1),
A vertex at (1, -1, 1),
A vertex at (1, 1, -1),
A vertex at (1, 1, 1),
An inequality (-1, 0, 0) x + 1 >= 0,
An inequality (0, -1, 0) x + 1 >= 0,
An inequality (0, 0, -1) x + 1 >= 0,
An inequality (0, 0, 1) x + 1 >= 0,
An inequality (0, 1, 0) x + 1 >= 0,
An inequality (1, 0, 0) x + 1 >= 0]
sage: G.automorphism_group().is_isomorphic(P.hasse_diagram().automorphism_group())
True
sage: O = polytopes.octahedron(); O
A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 6 vertices
sage: O.vertex_facet_graph()
Digraph on 14 vertices
sage: H = O.vertex_facet_graph()
sage: G.is_isomorphic(H)
False
sage: G2 = copy(G)
sage: G2.reverse_edges(G2.edges())
sage: G2.is_isomorphic(H)
True
TESTS:
Check that :trac:`28828` is fixed::
sage: G._immutable
True
Check that :trac:`29188` is fixed::
sage: P = polytopes.cube()
sage: P.vertex_facet_graph().is_isomorphic(P.vertex_facet_graph(False))
True
"""
return self.combinatorial_polyhedron().vertex_facet_graph(names=labels)
def plot(self,
point=None, line=None, polygon=None, # None means unspecified by the user
wireframe='blue', fill='green',
position=None,
orthonormal=True, # whether to use orthonormal projections
**kwds):
"""
Return a graphical representation.
INPUT:
- ``point``, ``line``, ``polygon`` -- Parameters to pass to
point (0d), line (1d), and polygon (2d) plot commands.
Allowed values are:
* A Python dictionary to be passed as keywords to the plot
commands.
* A string or triple of numbers: The color. This is
equivalent to passing the dictionary ``{'color':...}``.
* ``False``: Switches off the drawing of the corresponding
graphics object
- ``wireframe``, ``fill`` -- Similar to ``point``, ``line``,
and ``polygon``, but ``fill`` is used for the graphics
objects in the dimension of the polytope (or of dimension 2
for higher dimensional polytopes) and ``wireframe`` is used
for all lower-dimensional graphics objects
(default: 'green' for ``fill`` and 'blue' for ``wireframe``)
- ``position`` -- positive number; the position to take the projection
point in Schlegel diagrams.
- ``orthonormal`` -- Boolean (default: True); whether to use
orthonormal projections.
- ``**kwds`` -- optional keyword parameters that are passed to
all graphics objects.
OUTPUT:
A (multipart) graphics object.
EXAMPLES::
sage: square = polytopes.hypercube(2)
sage: point = Polyhedron([[1,1]])
sage: line = Polyhedron([[1,1],[2,1]])
sage: cube = polytopes.hypercube(3)
sage: hypercube = polytopes.hypercube(4)
By default, the wireframe is rendered in blue and the fill in green::
sage: square.plot()
Graphics object consisting of 6 graphics primitives
sage: point.plot()
Graphics object consisting of 1 graphics primitive
sage: line.plot()
Graphics object consisting of 2 graphics primitives
sage: cube.plot()
Graphics3d Object
sage: hypercube.plot()
Graphics3d Object
Draw the lines in red and nothing else::
sage: square.plot(point=False, line='red', polygon=False)
Graphics object consisting of 4 graphics primitives
sage: point.plot(point=False, line='red', polygon=False)
Graphics object consisting of 0 graphics primitives
sage: line.plot(point=False, line='red', polygon=False)
Graphics object consisting of 1 graphics primitive
sage: cube.plot(point=False, line='red', polygon=False)
Graphics3d Object
sage: hypercube.plot(point=False, line='red', polygon=False)
Graphics3d Object
Draw points in red, no lines, and a blue polygon::
sage: square.plot(point={'color':'red'}, line=False, polygon=(0,0,1))
Graphics object consisting of 2 graphics primitives
sage: point.plot(point={'color':'red'}, line=False, polygon=(0,0,1))
Graphics object consisting of 1 graphics primitive
sage: line.plot(point={'color':'red'}, line=False, polygon=(0,0,1))
Graphics object consisting of 1 graphics primitive
sage: cube.plot(point={'color':'red'}, line=False, polygon=(0,0,1))
Graphics3d Object
sage: hypercube.plot(point={'color':'red'}, line=False, polygon=(0,0,1))
Graphics3d Object
If we instead use the ``fill`` and ``wireframe`` options, the
coloring depends on the dimension of the object::
sage: square.plot(fill='green', wireframe='red')
Graphics object consisting of 6 graphics primitives
sage: point.plot(fill='green', wireframe='red')
Graphics object consisting of 1 graphics primitive
sage: line.plot(fill='green', wireframe='red')
Graphics object consisting of 2 graphics primitives
sage: cube.plot(fill='green', wireframe='red')
Graphics3d Object
sage: hypercube.plot(fill='green', wireframe='red')
Graphics3d Object
It is possible to draw polyhedra up to dimension 4, no matter what the
ambient dimension is::
sage: hcube = polytopes.hypercube(5)
sage: facet = hcube.facets()[0].as_polyhedron();facet
A 4-dimensional polyhedron in ZZ^5 defined as the convex hull of 16 vertices
sage: facet.plot()
Graphics3d Object
TESTS::
sage: for p in square.plot():
....: print("{} {}".format(p.options()['rgbcolor'], p))
blue Point set defined by 4 point(s)
blue Line defined by 2 points
blue Line defined by 2 points
blue Line defined by 2 points
blue Line defined by 2 points
green Polygon defined by 4 points
sage: for p in line.plot():
....: print("{} {}".format(p.options()['rgbcolor'], p))
blue Point set defined by 2 point(s)
green Line defined by 2 points
sage: for p in point.plot():
....: print("{} {}".format(p.options()['rgbcolor'], p))
green Point set defined by 1 point(s)
Draw the lines in red and nothing else::
sage: for p in square.plot(point=False, line='red', polygon=False):
....: print("{} {}".format(p.options()['rgbcolor'], p))
red Line defined by 2 points
red Line defined by 2 points
red Line defined by 2 points
red Line defined by 2 points
Draw vertices in red, no lines, and a blue polygon::
sage: for p in square.plot(point={'color':'red'}, line=False, polygon=(0,0,1)):
....: print("{} {}".format(p.options()['rgbcolor'], p))
red Point set defined by 4 point(s)
(0, 0, 1) Polygon defined by 4 points
sage: for p in line.plot(point={'color':'red'}, line=False, polygon=(0,0,1)):
....: print("{} {}".format(p.options()['rgbcolor'], p))
red Point set defined by 2 point(s)
sage: for p in point.plot(point={'color':'red'}, line=False, polygon=(0,0,1)):
....: print("{} {}".format(p.options()['rgbcolor'], p))