Trac #32013: Initialize a Set from a ConvexSet_base instance

A `Polyhedron` is not a `Parent`, so `Set` refuses to construct the set
of its elements.

We change this by creating a new abstract base class for non-
necessarily-parent sets with methods `union`, `intersection`, etc.

`ConvexSet_base` (from #31919) and `RealSet` now both inherit from
`Set_base`.  To complete the implementation of the `Set_base` protocol,
we add an implementation of `RealSet.symmetric_difference`.

So we can now do the following things:
{{{
sage: Set(polytopes.cube())
Set of elements of A 3-dimensional polyhedron in ZZ^3 defined as the
convex hull of 8 vertices
sage: polytopes.cube().union(polytopes.tetrahedron())
Set-theoretic union of Set of elements of A 3-dimensional polyhedron in
ZZ^3 defined as the convex hull of 8 vertices and Set of elements of A
3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4
vertices
}}}

URL: https://trac.sagemath.org/32013
Reported by: mkoeppe
Ticket author(s): Matthias Koeppe
Reviewer(s): Travis Scrimshaw

src/sage/geometry/cone.py

@@ -2006,6 +2006,27 @@ def _repr_(self):
             result += " face of %s" % self.ambient()
         return result
 
+    def _some_elements_(self):
+        r"""
+        Generate some points of ``self``.
+
+        EXAMPLE::
+
+            sage: K = cones.nonnegative_orthant(3)
+            sage: K.some_elements()  # indirect doctest
+            [(0, 0, 0), (1/2, 0, 0), (1/4, 1/2, 0), (1/8, 1/4, 1/2)]
+        """
+        V = self.ambient_vector_space()
+        r_iter = iter(self._rays)
+        p = V(0)
+        yield p
+        for i in range(5):
+            try:
+                p = (p + next(r_iter)) / 2
+            except StopIteration:
+                return
+            yield p
+
     def _sort_faces(self,  faces):
         r"""
         Return sorted (if necessary) ``faces`` as a tuple.

src/sage/geometry/convex_set.py

@@ -13,10 +13,13 @@
 # ****************************************************************************
 
 from sage.structure.sage_object import SageObject
+from sage.sets.set import Set_base
 from sage.categories.sets_cat import EmptySetError
 from sage.misc.abstract_method import abstract_method
+from sage.rings.infinity import infinity
+from sage.rings.integer_ring import ZZ
 
-class ConvexSet_base(SageObject):
+class ConvexSet_base(SageObject, Set_base):
     """
     Abstract base class for convex sets.
     """
@@ -38,6 +41,60 @@ def is_empty(self):
         """
         return self.dim() < 0
 
+    def is_finite(self):
+        r"""
+        Test whether ``self`` is a finite set.
+
+        OUTPUT:
+
+        Boolean.
+
+        EXAMPLES::
+
+            sage: p = LatticePolytope([], lattice=ToricLattice(3).dual()); p
+            -1-d lattice polytope in 3-d lattice M
+            sage: p.is_finite()
+            True
+            sage: q = Polyhedron(ambient_dim=2); q
+            The empty polyhedron in ZZ^2
+            sage: q.is_finite()
+            True
+            sage: r = Polyhedron(rays=[(1, 0)]); r
+            A 1-dimensional polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 ray
+            sage: r.is_finite()
+            False
+        """
+        return self.dim() < 1
+
+    def cardinality(self):
+        """
+        Return the cardinality of this set.
+
+        OUTPUT:
+
+        Either an integer or ``Infinity``.
+
+        EXAMPLES::
+
+            sage: p = LatticePolytope([], lattice=ToricLattice(3).dual()); p
+            -1-d lattice polytope in 3-d lattice M
+            sage: p.cardinality()
+            0
+            sage: q = Polyhedron(ambient_dim=2); q
+            The empty polyhedron in ZZ^2
+            sage: q.cardinality()
+            0
+            sage: r = Polyhedron(rays=[(1, 0)]); r
+            A 1-dimensional polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 ray
+            sage: r.cardinality()
+            +Infinity
+        """
+        if self.dim() < 0:
+            return ZZ(0)
+        if self.dim() == 0:
+            return ZZ(1)
+        return infinity
+
     def is_universe(self):
         r"""
         Test whether ``self`` is the whole ambient space.
@@ -370,7 +427,7 @@ def _test_convex_set(self, tester=None, **options):
             ....:         return 42
             ....:     def ambient_dim(self):
             ....:         return 91
-            sage: TestSuite(FaultyConvexSet()).run(skip=('_test_pickling', '_test_contains'))
+            sage: TestSuite(FaultyConvexSet()).run(skip=('_test_pickling', '_test_contains', '_test_as_set_object'))
             Failure in _test_convex_set:
             ...
             The following tests failed: _test_convex_set
@@ -384,7 +441,7 @@ def _test_convex_set(self, tester=None, **options):
             ....:         return QQ^3
             ....:     def ambient_dim(self):
             ....:         return 3
-            sage: TestSuite(BiggerOnTheInside()).run(skip=('_test_pickling', '_test_contains'))
+            sage: TestSuite(BiggerOnTheInside()).run(skip=('_test_pickling', '_test_contains', '_test_as_set_object'))
             Failure in _test_convex_set:
             ...
             The following tests failed: _test_convex_set
@@ -619,29 +676,6 @@ def _test_contains(self, tester=None, **options):
                     tester.assertTrue(self.contains(point))
                     tester.assertTrue(point in self)
 
-    @abstract_method(optional=True)
-    def intersection(self, other):
-        r"""
-        Return the intersection of ``self`` and ``other``.
-
-        INPUT:
-
-        - ``other`` -- another convex set
-
-        OUTPUT:
-
-        The intersection.
-
-        TESTS::
-
-            sage: from sage.geometry.convex_set import ConvexSet_base
-            sage: C = ConvexSet_base()
-            sage: C.intersection(C)
-            Traceback (most recent call last):
-            ...
-            TypeError: 'NotImplementedType' object is not callable
-        """
-
 
 class ConvexSet_closed(ConvexSet_base):
     r"""

src/sage/geometry/lattice_polytope.py

@@ -3794,6 +3794,36 @@ def points(self, *args, **kwds):
         else:
             return self._points
 
+    def _some_elements_(self):
+        r"""
+        Generate some points of ``self`` as a convex polytope.
+
+        In contrast to :meth:`points`, these are not necessarily lattice points.
+
+        EXAMPLE::
+
+            sage: o = lattice_polytope.cross_polytope(3)
+            sage: o.some_elements()  # indirect doctest
+            [(1, 0, 0),
+            (1/2, 1/2, 0),
+            (1/4, 1/4, 1/2),
+            (-3/8, 1/8, 1/4),
+            (-3/16, -7/16, 1/8),
+            (-3/32, -7/32, -7/16)]
+        """
+        if not self._vertices:
+            return
+        V = self.ambient_vector_space()
+        v_iter = iter(self._vertices)
+        p = V(next(v_iter))
+        yield p
+        for i in range(5):
+            try:
+                p = (p + next(v_iter)) / 2
+            except StopIteration:
+                return
+            yield p
+
     def polar(self):
         r"""
         Return the polar polytope, if this polytope is reflexive.

src/sage/geometry/polyhedron/base.py

@@ -8909,11 +8909,22 @@ def contains(self, point):
             True
             sage: full.contains([0])
             False
+
+        TESTS:
+
+        Passing non-iterable objects does not cause an exception, see :trac:`32013`::
+
+            sage: None in Polyhedron(vertices=[(0,0)], rays=[(1,0)], base_ring=QQ)
+            False
         """
         try:
             p = vector(point)
         except TypeError:  # point not iterable or no common ring for elements
-            if len(point) > 0:
+            try:
+                l = len(point)
+            except TypeError:
+                return False
+            if l > 0:
                 return False
             else:
                 p = vector(self.base_ring(), [])
@@ -9012,7 +9023,11 @@ def interior_contains(self, point):
         try:
             p = vector(point)
         except TypeError:  # point not iterable or no common ring for elements
-            if len(point) > 0:
+            try:
+                l = len(point)
+            except TypeError:
+                return False
+            if l > 0:
                 return False
             else:
                 p = vector(self.base_ring(), [])
@@ -9129,7 +9144,11 @@ def relative_interior_contains(self, point):
         try:
             p = vector(point)
         except TypeError:  # point not iterable or no common ring for elements
-            if len(point) > 0:
+            try:
+                l = len(point)
+            except TypeError:
+                return False
+            if l > 0:
                 return False
             else:
                 p = vector(self.base_ring(), [])

src/sage/sets/real_set.py

@@ -89,6 +89,8 @@ class RealSet.
 from sage.structure.parent import Parent
 from sage.structure.unique_representation import UniqueRepresentation
 from sage.categories.topological_spaces import TopologicalSpaces
+from sage.categories.sets_cat import Sets
+from sage.sets.set import Set_base, Set_boolean_operators, Set_add_sub_operators
 from sage.rings.all import ZZ
 from sage.rings.real_lazy import LazyFieldElement, RLF
 from sage.rings.infinity import infinity, minus_infinity
@@ -791,7 +793,8 @@ def __rmul__(self, other):
         return self * other
 
 @richcmp_method
-class RealSet(UniqueRepresentation, Parent):
+class RealSet(UniqueRepresentation, Parent, Set_base,
+              Set_boolean_operators, Set_add_sub_operators):
 
     @staticmethod
     def __classcall__(cls, *args):
@@ -1621,9 +1624,6 @@ def union(self, *other):
         intervals = self._intervals + other._intervals
         return RealSet(*intervals)
 
-    __or__ = union
-    __add__ = union
-
     def intersection(self, *other):
         """
         Return the intersection of the two sets
@@ -1668,8 +1668,6 @@ def intersection(self, *other):
                 intervals.append(i1.intersection(i2))
         return RealSet(*intervals)
 
-    __and__ = intersection
-
     def inf(self):
         """
         Return the infimum
@@ -1794,7 +1792,30 @@ def difference(self, *other):
         other = RealSet(*other)
         return self.intersection(other.complement())
 
-    __sub__ = difference
+    def symmetric_difference(self, *other):
+        r"""
+        Returns the symmetric difference of ``self`` and ``other``.
+
+        INPUT:
+
+        - ``other`` -- a :class:`RealSet` or data that defines one.
+
+        OUTPUT:
+
+        The set-theoretic symmetric difference of ``self`` and ``other``
+        as a new :class:`RealSet`.
+
+        EXAMPLES::
+
+            sage: s1 = RealSet(0,2); s1
+            (0, 2)
+            sage: s2 = RealSet.unbounded_above_open(1); s2
+            (1, +oo)
+            sage: s1.symmetric_difference(s2)
+            (0, 1] ∪ [2, +oo)
+        """
+        other = RealSet(*other)
+        return self.difference(other).union(other.difference(self))
 
     def contains(self, x):
         """