Implement as_disjoint_union for unions of regions

blueprint/src/python/polytope.py

@@ -247,7 +247,8 @@ def intersection(polys, canonicalize=True):
     #
     # This implementation only works with finite polytopes. Extreme rays are
     # ignored by the vertex-only feasibility checks below, so unbounded inputs
-    # previously produced silently incorrect envelopes; reject them instead.
+    # would produce silently incorrect envelopes; report them as not a polytope
+    # instead, which leaves the caller with the pieces it started with.
     def try_union(polys):
         if not isinstance(polys, list) or len(polys) == 0:
             return None
@@ -259,10 +260,7 @@ def try_union(polys):
             if p.rays is None:
                 p.compute_V_rep()
             if p.rays:
-                raise ValueError(
-                    "Polytope.try_union only supports finite polytopes "
-                    "(an input has extreme rays)"
-                )
+                return None
 
         # Performance optimisation: see if intersection is nonempty first. This
         # is relative cheap compared to the full union operation, since it only

blueprint/src/python/region.py

@@ -432,7 +432,30 @@ def _as_disjoint_union_poly(self, verbose=False, track_dependencies=False) -> li
             return A
 
         if self.region_type == Region_Type.UNION:
-            raise NotImplementedError(self.region_type) # TODO: implement this
+            # The children may overlap, so each new piece is trimmed against the
+            # pieces already accepted: Polytope.set_minus splits A \\ B into
+            # pairwise disjoint parts, so the accumulated list stays disjoint and
+            # its union is unchanged.
+            result = []
+            for r in self.child:
+                for p in r._as_disjoint_union_poly(
+                    verbose=verbose,
+                    track_dependencies=track_dependencies
+                ):
+                    parts = [p]
+                    for q in result:
+                        trimmed = []
+                        for part in parts:
+                            pieces = part.set_minus(q)
+                            if track_dependencies:
+                                for piece in pieces:
+                                    piece.dependencies = part.dependencies
+                            trimmed.extend(pieces)
+                        parts = trimmed
+                        if not parts:
+                            break
+                    result.extend(parts)
+            return result
         if self.region_type == Region_Type.DISJOINT_UNION:
             result = []
             for r in self.child:

blueprint/src/python/tests/test_region.py

@@ -34,5 +34,34 @@ def test_intersect():
         x = (rd.uniform(0, 3), rd.uniform(0, 3))
         assert intersection.contains(x) == (R1.contains(x) and R2.contains(x))
 
+def test_as_disjoint_union():
+
+    # A union of overlapping rectangles, rewritten as a disjoint union: the
+    # pieces must cover exactly the same points, and no point may lie in the
+    # interior of two of them.
+    for i in range(50):
+        rects = []
+        for _ in range(3):
+            xlims = (rd.uniform(0, 3), rd.uniform(0, 3))
+            ylims = (rd.uniform(0, 3), rd.uniform(0, 3))
+            rects.append(
+                Polytope.rect((min(xlims), max(xlims)), (min(ylims), max(ylims)))
+            )
+
+        union = Region.union([Region.from_polytope(r) for r in rects])
+        pieces = union.as_disjoint_union()
+
+        for _ in range(20):
+            x = (rd.uniform(0, 3), rd.uniform(0, 3))
+            assert pieces.contains(x) == any(r.contains(x) for r in rects)
+
+        # Pairwise disjoint: any two pieces may share a face, not an interior.
+        polys = [piece.child for piece in pieces.child]
+        for j in range(len(polys)):
+            for k in range(j + 1, len(polys)):
+                overlap = polys[j].intersect(polys[k])
+                assert overlap.is_empty(include_boundary=False)
+
 test_union()
 test_intersect()
+test_as_disjoint_union()