diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index d50b238..2374f02 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -10,7 +10,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.7", "3.8", "3.9", "3.10"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
     steps:
       - uses: actions/checkout@v2
       - name: Set up Python
diff --git a/README.md b/README.md
index ad2aee0..fa587e7 100644
--- a/README.md
+++ b/README.md
@@ -8,7 +8,7 @@
 [![Commits](https://img.shields.io/github/last-commit/tilmann-bartsch/rportion/master)](https://github.com/tilmann-bartsch/rportion/commits/master)
 
 The `rportion` library provides data structure to represent 
-2D [rectilinear polygons](https://en.wikipedia.org/wiki/Rectilinear_polygon) (unions of 2D-intervals) in Python 3.7+.
+2D [rectilinear polygons](https://en.wikipedia.org/wiki/Rectilinear_polygon) (unions of 2D-intervals) in Python 3.9+.
 It is built upon the library [`portion`](https://github.com/AlexandreDecan/portion) and follows its concepts.
 The following features are provided:
 
@@ -55,6 +55,8 @@ pip install -e ".[test]"
 python -m unittest discover -s tests
 ```
 
+Note that `python
+
 ## Documentation & usage
 
 ### Polygon creation
diff --git a/rportion/rportion.py b/rportion/rportion.py
index 769a839..b395faf 100644
--- a/rportion/rportion.py
+++ b/rportion/rportion.py
@@ -2,11 +2,10 @@
 from typing import Iterator, List, Tuple, Callable
 
 import portion as P
-from portion.interval import Interval, open, empty, closedopen, closed, openclosed, singleton
 from sortedcontainers import SortedList
 
 
-def _sub_contained_atomics(int_1: Interval, int_2: Interval) -> Interval:
+def _sub_contained_atomics(int_1: P.Interval, int_2: P.Interval) -> P.Interval:
     """
     Remove all atomic sectors from int_1 which are contained completely in int_2, i.e.
 
@@ -15,7 +14,7 @@ def _sub_contained_atomics(int_1: Interval, int_2: Interval) -> Interval:
          int_2: ----------    ---------
         return:             ----
     """
-    ret_int = empty()
+    ret_int = P.empty()
     for y_atomic in int_1:
         if not int_2.contains(y_atomic):
             ret_int |= y_atomic
@@ -93,7 +92,7 @@ def __repr__(self):
         return f"({self.val}, {self.btype})"
 
 
-def _extend_ranges_mat(mat: List[List[Interval]],
+def _extend_ranges_mat(mat: List[List[P.Interval]],
                        x_boundaries: 'SortedList[RBoundary]',
                        rbound: RBoundary):
     n = len(mat)
@@ -105,12 +104,12 @@ def _extend_ranges_mat(mat: List[List[Interval]],
         mat[row].insert(n - index, prev_y_int)
 
 
-def _y_interval_triangle_prunable(bound_index: int, y_interval_triangle: List[List[Interval]]) -> bool:
+def _y_interval_triangle_prunable(bound_index: int, y_interval_triangle: List[List[P.Interval]]) -> bool:
     """
     Test if the given bound can be removed from y_interval_triangle.
 
     :bound_index int:
-    :y_interval_triangle list[list[Interval]]:
+    :y_interval_triangle list[list[P.Interval]]:
     :return bool:
 
     For the following triangle the boundary 4 can be pruned, because
@@ -158,19 +157,19 @@ def _y_interval_triangle_prunable(bound_index: int, y_interval_triangle: List[Li
     return True
 
 
-def _interval_triangle_prunable_indices(y_interval_triangle: List[List[Interval]]) -> List[int]:
+def _interval_triangle_prunable_indices(y_interval_triangle: List[List[P.Interval]]) -> List[int]:
     """Return a list of indices which can be removed from y_interval_triangle."""
     n = len(y_interval_triangle)
     return [b_index for b_index in range(1, n + 1)
             if _y_interval_triangle_prunable(b_index, y_interval_triangle)]
 
 
-def _remove_boundaries(bound_indices: List[int], y_interval_triangle: List[List[Interval]]):
+def _remove_boundaries(bound_indices: List[int], y_interval_triangle: List[List[P.Interval]]):
     """
     Remove rows and columnss from y_interval_triangle specified by bound_indices.
 
     :bound_indices list[int]:
-    :y_interval_triangle list[list[Interval]]: Data
+    :y_interval_triangle list[list[P.Interval]]: Data
 
     This function is meant to be performed if the boundary is prunable, see function
     `_y_interval_triangle_prunable`.
@@ -197,9 +196,9 @@ def _remove_boundaries(bound_indices: List[int], y_interval_triangle: List[List[
 
 def _update_x_boundaries_and_y_interval_triangles(
         x_boundaries: 'SortedList[RBoundary]',
-        y_interval_triangle_add: List[List[Interval]],
-        y_interval_triangle_sub: List[List[Interval]],
-        other_x_atom: Interval, other_y_interval: Interval):
+        y_interval_triangle_add: List[List[P.Interval]],
+        y_interval_triangle_sub: List[List[P.Interval]],
+        other_x_atom: P.Interval, other_y_interval: P.Interval):
     """
     Update
        - y_interval_triangle_add
@@ -208,9 +207,9 @@ def _update_x_boundaries_and_y_interval_triangles(
                 such that it represents the difference of the old polygon and the one provided as paramters
 
     :x_boundaries SortedList[RBoundary]:
-    :y_interval_triangle_add list[list[Interval]]:
-    :y_interval_triangle_sub list[list[Interval]]:
-    :new_x_atom Atomic, new_y_interval: Interval:
+    :y_interval_triangle_add list[list[P.Interval]]:
+    :y_interval_triangle_sub list[list[P.Interval]]:
+    :new_x_atom Atomic, new_y_interval: P.Interval:
     """
     assert other_x_atom.atomic
 
@@ -236,7 +235,7 @@ def _update_x_boundaries_and_y_interval_triangles(
         for j in range(n - i):
             l_bound = x_boundaries[i]
             r_bound = x_boundaries[-j - 1]
-            x_interval = Interval.from_atomic(~l_bound.btype, l_bound.val, r_bound.val, r_bound.btype)
+            x_interval = P.Interval.from_atomic(~l_bound.btype, l_bound.val, r_bound.val, r_bound.btype)
 
             if other_x_atom.contains(x_interval):
                 y_interval_triangle_add[i][j] |= other_y_interval
@@ -249,21 +248,21 @@ def _update_x_boundaries_and_y_interval_triangles(
                 left_ind = x_boundaries.bisect_left(other_x_int_l_bound)
                 col = n - left_ind
                 if l_bound < other_x_int_l_bound and (0 < col < n - i):
-                    adj_other_x_interval |= Interval.from_atomic(~l_bound.btype, l_bound.val,
+                    adj_other_x_interval |= P.Interval.from_atomic(~l_bound.btype, l_bound.val,
                                                                  adj_other_x_interval.upper,
                                                                  adj_other_x_interval.right)
                     y_int_left = y_interval_triangle_add[i][col]
                 else:
-                    y_int_left = open(-P.inf, P.inf)
+                    y_int_left = P.open(-P.inf, P.inf)
 
                 right_ind = x_boundaries.bisect_left(other_y_int_r_bound)
                 if other_y_int_r_bound < r_bound and right_ind < n - j:
-                    adj_other_x_interval |= Interval.from_atomic(adj_other_x_interval.left,
+                    adj_other_x_interval |= P.Interval.from_atomic(adj_other_x_interval.left,
                                                                  adj_other_x_interval.lower,
                                                                  r_bound.val, r_bound.btype)
                     y_int_right = y_interval_triangle_add[right_ind][j]
                 else:
-                    y_int_right = open(-P.inf, P.inf)
+                    y_int_right = P.open(-P.inf, P.inf)
 
                 if adj_other_x_interval.contains(x_interval):
                     y_interval_triangle_add[i][j] |= y_int_left & other_y_interval & y_int_right
@@ -274,7 +273,7 @@ def _update_x_boundaries_and_y_interval_triangles(
         for j in range(n - i):
             l_bound = ~x_boundaries[i]
             r_bound = x_boundaries[-j - 1]
-            x_interval = Interval.from_atomic(l_bound.btype, l_bound.val, r_bound.val, r_bound.btype)
+            x_interval = P.Interval.from_atomic(l_bound.btype, l_bound.val, r_bound.val, r_bound.btype)
             if other_x_atom.overlaps(x_interval):
                 y_interval_triangle_sub[i][j] -= other_y_interval
 
@@ -291,10 +290,10 @@ def _update_x_boundaries_and_y_interval_triangles(
 
 
 def _traverse_diagonally(boundaries: List[RBoundary],
-                         interval_triangle: List[List[Interval]],
-                         next_accumulator: Callable[[Interval, Interval, Interval], Interval],
-                         adj_y_interval: Callable[[Interval, Interval, Interval], Interval]
-                         ) -> Iterator[Tuple[Interval, Interval]]:
+                         interval_triangle: List[List[P.Interval]],
+                         next_accumulator: Callable[[P.Interval, P.Interval, P.Interval], P.Interval],
+                         adj_y_interval: Callable[[P.Interval, P.Interval, P.Interval], P.Interval]
+                         ) -> Iterator[Tuple[P.Interval, P.Interval]]:
     """
     Traverse `interval_triangle` diagonally from the top left and yield rectangles specified by the parameters
     `asdf` and `next_accumulator`.
@@ -302,10 +301,10 @@ def _traverse_diagonally(boundaries: List[RBoundary],
     The iterator *DOES NOT* return tuples where either the first or second interval is empty.
 
     :param boundaries: List[RBoundary]:
-    :param interval_triangle: List[List[Interval]]:
-    :param next_accumulator: Callable[[Interval, Interval, Interval], Interval]:
+    :param interval_triangle: List[List[P.Interval]]:
+    :param next_accumulator: Callable[[P.Interval, P.Interval, P.Interval], P.Interval]:
         Function determining how to accumulate values while traversing. See explanation below.
-    :param adj_y_interval: Callable[[Interval, Interval, Interval], Interval]:
+    :param adj_y_interval: Callable[[P.Interval, P.Interval, P.Interval], P.Interval]:
         Function determining which rectangles to return. See explanation below.
 
     `interval_triangle` represents an interval tree of the form
@@ -326,10 +325,10 @@ def _traverse_diagonally(boundaries: List[RBoundary],
     All generated rectangles which are non-empty are returned.
     """
     n = len(interval_triangle)
-    next_parent_y_intervals = [empty(), empty()]
+    next_parent_y_intervals = [P.empty(), P.empty()]
     for start_row in range(n):
         parent_y_intervals = next_parent_y_intervals
-        next_parent_y_intervals = [empty()]
+        next_parent_y_intervals = [P.empty()]
         for row, col in enumerate(range(start_row, -1, -1)):
             curr_y_int = interval_triangle[row][col]
 
@@ -337,7 +336,7 @@ def _traverse_diagonally(boundaries: List[RBoundary],
 
             l_bound = ~boundaries[row]
             r_bound = boundaries[-col - 1]
-            curr_x_int = Interval.from_atomic(l_bound.btype, l_bound.val,
+            curr_x_int = P.Interval.from_atomic(l_bound.btype, l_bound.val,
                                               r_bound.val, r_bound.btype)
             if not curr_x_int.empty and not adj_curr_y_int.empty:
                 yield curr_x_int, adj_curr_y_int
@@ -345,7 +344,7 @@ def _traverse_diagonally(boundaries: List[RBoundary],
             next_parent_y_intervals.append(next_accumulator(curr_y_int,
                                                             parent_y_intervals[row],
                                                             parent_y_intervals[row + 1]))
-        next_parent_y_intervals.append(empty())
+        next_parent_y_intervals.append(P.empty())
 
 
 class RPolygon:
@@ -374,8 +373,8 @@ def __init__(self):
         #      b2 |         ___/
         #       : |      __/
         #      bn |_____/
-        self._used_y_ranges: List[List[Interval]] = [[empty()]]
-        self._free_y_ranges: List[List[Interval]] = [[open(-P.inf, P.inf)]]
+        self._used_y_ranges: List[List[P.Interval]] = [[P.empty()]]
+        self._free_y_ranges: List[List[P.Interval]] = [[P.open(-P.inf, P.inf)]]
 
     @property
     def x_left(self):
@@ -456,36 +455,36 @@ def enclosure(self) -> 'RPolygon':
         return self.__class__.from_interval_product(self.x_enclosure_interval, self.y_enclosure_interval)
 
     @property
-    def x_enclosure_interval(self) -> Interval:
+    def x_enclosure_interval(self) -> P.Interval:
         """
         smallest y_interval enclosing the y-dimension of the polygon.
         """
         n = len(self._free_y_ranges[0])
         for l_int in range(0, n):
-            if self._free_y_ranges[0][l_int] == open(-P.inf, P.inf):
+            if self._free_y_ranges[0][l_int] == P.open(-P.inf, P.inf):
                 break
         else:
             l_int += 1
         n = len(self._free_y_ranges)
         for r_int in range(0, n):
-            if self._free_y_ranges[r_int][0] == open(-P.inf, P.inf):
+            if self._free_y_ranges[r_int][0] == P.open(-P.inf, P.inf):
                 break
         else:
             r_int += 1
         l_bound = self._x_boundaries[n - l_int]
         r_bound = self._x_boundaries[r_int]
-        return Interval.from_atomic(~l_bound.btype, l_bound.val,
+        return P.Interval.from_atomic(~l_bound.btype, l_bound.val,
                                     r_bound.val, r_bound.btype)
 
     @property
-    def y_enclosure_interval(self) -> Interval:
+    def y_enclosure_interval(self) -> P.Interval:
         """
         smallest y_interval enclosing the y-dimension of the polygon.
         """
         return ~self._free_y_ranges[0][0]
 
     @classmethod
-    def from_interval_product(cls, x_interval: Interval, y_interval: Interval) -> 'RPolygon':
+    def from_interval_product(cls, x_interval: P.Interval, y_interval: P.Interval) -> 'RPolygon':
         """
         Create a (simple) rectangular polygon as the product of two intervals.
 
@@ -562,7 +561,7 @@ def boundary(self) -> 'RPolygon':
 
         :return RPolygon: boundary of this polygon represented as another polygon
         """
-        def adj_y_interval(curr: Interval, l_parent: Interval, r_parent: Interval) -> Interval:
+        def adj_y_interval(curr: P.Interval, l_parent: P.Interval, r_parent: P.Interval) -> P.Interval:
             return curr - l_parent - r_parent
 
         iterator = _traverse_diagonally(list(self._x_boundaries), self._used_y_ranges,
@@ -570,7 +569,7 @@ def adj_y_interval(curr: Interval, l_parent: Interval, r_parent: Interval) -> In
                                         adj_y_interval)
         x_boundary_poly = rempty()
         for x_atom, y_interval in iterator:
-            x_int = singleton(x_atom.lower) | singleton(x_atom.upper)
+            x_int = P.singleton(x_atom.lower) | P.singleton(x_atom.upper)
             x_boundary_poly |= self.__class__.from_interval_product(x_int, y_interval)
 
         # Traverse the leafs of the interval tree
@@ -581,10 +580,10 @@ def adj_y_interval(curr: Interval, l_parent: Interval, r_parent: Interval) -> In
             col = n - i - 1
             l_bound = ~self._x_boundaries[i]
             r_bound = self._x_boundaries[i+1]
-            x_int = closed(l_bound.val, r_bound.val)
+            x_int = P.closed(l_bound.val, r_bound.val)
             leaf_y_int = self._used_y_ranges[row][col]
             for y_atom in leaf_y_int:
-                y_int = singleton(y_atom.lower) | singleton(y_atom.upper)
+                y_int = P.singleton(y_atom.lower) | P.singleton(y_atom.upper)
                 y_boundary_poly |= self.__class__.from_interval_product(x_int, y_int)
 
         return x_boundary_poly | y_boundary_poly
@@ -648,24 +647,24 @@ def __repr__(self):
             string.append(f"(x={repr(x_int)}, y={repr(y_int)})")
         return " | ".join(string)
 
-    def _atomic_x_rectangle_partitioning(self) -> Iterator[Tuple[Interval, Interval]]:
+    def _atomic_x_rectangle_partitioning(self) -> Iterator[Tuple[P.Interval, P.Interval]]:
         """Traverse `self._used_y_ranges` to obtain a rectangle partitioning of the polygon.
 
         This function *MUST NOT* return tuples where either the first or the second interval is empty.
         """
-        def adj_y_interval(curr: Interval, l_parent: Interval, r_parent: Interval) -> Interval:
+        def adj_y_interval(curr: P.Interval, l_parent: P.Interval, r_parent: P.Interval) -> P.Interval:
             return curr - l_parent - r_parent
 
         return _traverse_diagonally(list(self._x_boundaries), self._used_y_ranges,
                                     lambda curr, l_parent, r_parent: curr | l_parent | r_parent,
                                     adj_y_interval)
 
-    def _maximal_atomic_x_rectangles(self) -> Iterator[Tuple[Interval, Interval]]:
+    def _maximal_atomic_x_rectangles(self) -> Iterator[Tuple[P.Interval, P.Interval]]:
         """Traverse `self._used_y_ranges` to obtain the maximum contained rectangles.
 
         This function *MUST NOT* return tuples where either the first or the second interval is empty.
         """
-        def adj_y_interval(curr: Interval, l_parent: Interval, r_parent: Interval) -> Interval:
+        def adj_y_interval(curr: P.Interval, l_parent: P.Interval, r_parent: P.Interval) -> P.Interval:
             return _sub_contained_atomics(_sub_contained_atomics(curr, l_parent), r_parent)
 
         return _traverse_diagonally(list(self._x_boundaries), self._used_y_ranges,
@@ -677,17 +676,17 @@ def _invert(self):
         self._used_y_ranges = self._free_y_ranges
         self._free_y_ranges = temp
 
-    def _add_interval_product(self, x_interval: Interval, y_interval: Interval):
+    def _add_interval_product(self, x_interval: P.Interval, y_interval: P.Interval):
         if not x_interval.empty and not y_interval.empty:
             for x_atom in x_interval:
                 self._add_atomic(x_atom, y_interval)
 
-    def _sub_interval_product(self, x_interval: Interval, y_interval: Interval):
+    def _sub_interval_product(self, x_interval: P.Interval, y_interval: P.Interval):
         if not x_interval.empty and not y_interval.empty:
             for x_atom in x_interval:
                 self._sub_atomic(x_atom, y_interval)
 
-    def _add_atomic(self, add_x_atom: Interval, add_y_interval: Interval):
+    def _add_atomic(self, add_x_atom: P.Interval, add_y_interval: P.Interval):
         """
         Update self._used_y_ranges and self._free_y_ranges such that they represent the union of the old rectangular
         polygon and the provided interval product.
@@ -700,7 +699,7 @@ def _add_atomic(self, add_x_atom: Interval, add_y_interval: Interval):
             add_y_interval,
         )
 
-    def _sub_atomic(self, sub_x_atom: Interval, sub_y_interval: Interval):
+    def _sub_atomic(self, sub_x_atom: P.Interval, sub_y_interval: P.Interval):
         """
         Update self._used_y_ranges and self._free_y_ranges such that they represent the set difference of the old
         rectangular polygon and the provided interval product.
@@ -720,7 +719,7 @@ def rempty() -> RPolygon:
 
     :return RPolygon:
     """
-    return RPolygon.from_interval_product(empty(), empty())
+    return RPolygon.from_interval_product(P.empty(), P.empty())
 
 
 def rsingleton(x, y) -> RPolygon:
@@ -731,10 +730,10 @@ def rsingleton(x, y) -> RPolygon:
     :param y: y-value of the singleton
     :return RPolygon:
     """
-    return RPolygon.from_interval_product(singleton(x), singleton(y))
+    return RPolygon.from_interval_product(P.singleton(x), P.singleton(y))
 
 
-def rproduct(x_interval: Interval, y_interval: Interval):
+def rproduct(x_interval: P.Interval, y_interval: P.Interval):
     """
     Create a polygon as the (cartesian) product of two intervals.
 
@@ -755,7 +754,7 @@ def ropen(x_lower, x_upper, y_lower, y_upper) -> RPolygon:
     :param y_upper: value of the upper bound in y-dimension
     :return RPolygon:
     """
-    return RPolygon.from_interval_product(open(x_lower, x_upper), open(y_lower, y_upper))
+    return RPolygon.from_interval_product(P.open(x_lower, x_upper), P.open(y_lower, y_upper))
 
 
 def rclosed(x_lower, x_upper, y_lower, y_upper) -> RPolygon:
@@ -768,7 +767,7 @@ def rclosed(x_lower, x_upper, y_lower, y_upper) -> RPolygon:
     :param y_upper: value of the upper bound in y-dimension
     :return RPolygon:
     """
-    return RPolygon.from_interval_product(closed(x_lower, x_upper), closed(y_lower, y_upper))
+    return RPolygon.from_interval_product(P.closed(x_lower, x_upper), P.closed(y_lower, y_upper))
 
 
 def ropenclosed(x_lower, x_upper, y_lower, y_upper) -> RPolygon:
@@ -782,7 +781,7 @@ def ropenclosed(x_lower, x_upper, y_lower, y_upper) -> RPolygon:
     :param y_upper: value of the upper bound in y-dimension
     :return RPolygon:
     """
-    return RPolygon.from_interval_product(openclosed(x_lower, x_upper), openclosed(y_lower, y_upper))
+    return RPolygon.from_interval_product(P.openclosed(x_lower, x_upper), P.openclosed(y_lower, y_upper))
 
 
 def rclosedopen(x_lower, x_upper, y_lower, y_upper) -> RPolygon:
@@ -796,4 +795,4 @@ def rclosedopen(x_lower, x_upper, y_lower, y_upper) -> RPolygon:
     :param y_upper: value of the upper bound in y-dimension
     :return RPolygon:
     """
-    return RPolygon.from_interval_product(closedopen(x_lower, x_upper), closedopen(y_lower, y_upper))
+    return RPolygon.from_interval_product(P.closedopen(x_lower, x_upper), P.closedopen(y_lower, y_upper))
diff --git a/setup.py b/setup.py
index 1b85b73..84d1b05 100644
--- a/setup.py
+++ b/setup.py
@@ -9,7 +9,7 @@
 
 setup(
     name="rportion",
-    version="0.1.0",
+    version="0.2.0",
     license="MIT",
     author="Tilmann Bartsch",
     url="https://github.com/tilmann-bartsch/rportion",
@@ -27,15 +27,15 @@
     ],
     keywords="rectangle polygon interval-tree",
     packages=find_packages(include=["rportion"]),
-    python_requires="~= 3.7",
+    python_requires=">=3.9",
     install_requires=[
-        "sortedcontainers ~= 2.2", "portion ~= 2.2.0"
+        "sortedcontainers ~= 2.2", "portion ~= 2.4.2"
     ],
     extras_require={
         "test": ["pytest ~= 7.0",
                  "coverage ~= 6.0",
-                 "numpy ~= 1.21.6"],
-        "docu": ["numpy ~= 1.21.6",
+                 "numpy ~= 1.26.4"],
+        "docu": ["numpy ~= 1.26.4",
                  "matplotlib ~= 3.5.2",
                  "imageio ~= 2.19.5",
                  "pillow ~= 9.2.0",
diff --git a/tests/test_rportion.py b/tests/test_rportion.py
index c647ef3..ed11b33 100644
--- a/tests/test_rportion.py
+++ b/tests/test_rportion.py
@@ -5,7 +5,6 @@
 
 import numpy as np
 import portion as P
-from portion.interval import open, closedopen, empty, Interval, singleton, openclosed, closed
 from typing import List, Set, Tuple
 from rportion import RPolygon
 from rportion.rportion import ropen, rclosedopen, rclosed, rempty, RBoundary, rsingleton, ropenclosed, \
@@ -15,7 +14,7 @@
 
 
 def data_tree_to_string(x_boundaries: List[RBoundary],
-                        y_intervals: List[List[Interval]],
+                        y_intervals: List[List[P.Interval]],
                         spacing: int):
     n = len(y_intervals)
     msg = " " * spacing + "|"
@@ -205,14 +204,14 @@ def test_rsingleton(self):
             RBoundary(P.inf, P.OPEN),
         ])
         self.assertListEqual(p._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), singleton(y)],
-            [empty()],
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.singleton(y)],
+            [P.empty()],
         ])
         self.assertListEqual(p._free_y_ranges, [
-            [open(-P.inf, y) | open(y, P.inf), open(-P.inf, y) | open(y, P.inf), open(-P.inf, P.inf)],
-            [open(-P.inf, y) | open(y, P.inf), open(-P.inf, y) | open(y, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.open(-P.inf, y) | P.open(y, P.inf), P.open(-P.inf, y) | P.open(y, P.inf), P.open(-P.inf, P.inf)],
+            [P.open(-P.inf, y) | P.open(y, P.inf), P.open(-P.inf, y) | P.open(y, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
 
     def test_rempty(self):
@@ -221,8 +220,8 @@ def test_rempty(self):
             RBoundary(-P.inf, P.OPEN),
             RBoundary(P.inf, P.OPEN),
         ])
-        self.assertListEqual(p._used_y_ranges, [[empty()]])
-        self.assertListEqual(p._free_y_ranges, [[open(-P.inf, P.inf)]])
+        self.assertListEqual(p._used_y_ranges, [[P.empty()]])
+        self.assertListEqual(p._free_y_ranges, [[P.open(-P.inf, P.inf)]])
 
     def test_ropen(self):
         x1, x2 = 0, 3
@@ -235,16 +234,16 @@ def test_ropen(self):
             RBoundary(P.inf, P.OPEN),
         ])
         self.assertListEqual(p._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), open(y1, y2)],
-            [empty()],
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.open(y1, y2)],
+            [P.empty()],
         ])
         self.assertListEqual(p._free_y_ranges, [
-            [openclosed(-P.inf, y1) | closedopen(y2, P.inf),
-             openclosed(-P.inf, y1) | closedopen(y2, P.inf),
-             open(-P.inf, P.inf)],
-            [openclosed(-P.inf, y1) | closedopen(y2, P.inf), openclosed(-P.inf, y1) | closedopen(y2, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.closedopen(y2, P.inf),
+             P.openclosed(-P.inf, y1) | P.closedopen(y2, P.inf),
+             P.open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.closedopen(y2, P.inf), P.openclosed(-P.inf, y1) | P.closedopen(y2, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
 
     def test_rclosed(self):
@@ -258,14 +257,14 @@ def test_rclosed(self):
             RBoundary(P.inf, P.OPEN),
         ])
         self.assertListEqual(p._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), closed(y1, y2)],
-            [empty()],
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.closed(y1, y2)],
+            [P.empty()],
         ])
         self.assertListEqual(p._free_y_ranges, [
-            [open(-P.inf, y1) | open(y2, P.inf), open(-P.inf, y1) | open(y2, P.inf), open(-P.inf, P.inf)],
-            [open(-P.inf, y1) | open(y2, P.inf), open(-P.inf, y1) | open(y2, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.open(-P.inf, y1) | P.open(y2, P.inf), P.open(-P.inf, y1) | P.open(y2, P.inf), P.open(-P.inf, P.inf)],
+            [P.open(-P.inf, y1) | P.open(y2, P.inf), P.open(-P.inf, y1) | P.open(y2, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
 
     def test_ropenclosed(self):
@@ -279,14 +278,14 @@ def test_ropenclosed(self):
             RBoundary(P.inf, P.OPEN),
         ])
         self.assertListEqual(p._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), openclosed(y1, y2)],
-            [empty()],
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.openclosed(y1, y2)],
+            [P.empty()],
         ])
         self.assertListEqual(p._free_y_ranges, [
-            [openclosed(-P.inf, y1) | open(y2, P.inf), openclosed(-P.inf, y1) | open(y2, P.inf), open(-P.inf, P.inf)],
-            [openclosed(-P.inf, y1) | open(y2, P.inf), openclosed(-P.inf, y1) | open(y2, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.open(y2, P.inf), P.openclosed(-P.inf, y1) | P.open(y2, P.inf), P.open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.open(y2, P.inf), P.openclosed(-P.inf, y1) | P.open(y2, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
 
     def test_rclosedopen(self):
@@ -300,14 +299,14 @@ def test_rclosedopen(self):
             RBoundary(P.inf, P.OPEN),
         ])
         self.assertListEqual(p._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), closedopen(y1, y2)],
-            [empty()],
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.closedopen(y1, y2)],
+            [P.empty()],
         ])
         self.assertListEqual(p._free_y_ranges, [
-            [open(-P.inf, y1) | closedopen(y2, P.inf), open(-P.inf, y1) | closedopen(y2, P.inf), open(-P.inf, P.inf)],
-            [open(-P.inf, y1) | closedopen(y2, P.inf), open(-P.inf, y1) | closedopen(y2, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.open(-P.inf, y1) | P.closedopen(y2, P.inf), P.open(-P.inf, y1) | P.closedopen(y2, P.inf), P.open(-P.inf, P.inf)],
+            [P.open(-P.inf, y1) | P.closedopen(y2, P.inf), P.open(-P.inf, y1) | P.closedopen(y2, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
 
 
@@ -345,11 +344,11 @@ def test_x_enclosure_interval(self):
         y1, y2, y3 = 0, 1, 2
         self.assertEqual(
             (rclosed(x1, x3, y1, y3) | rclosed(x2, x4, y1, y2)).x_enclosure_interval,
-            closed(x1, x4)
+            P.closed(x1, x4)
         )
         self.assertEqual(
             (ropen(-P.inf, x1, y1, y2) | ropen(x2, P.inf, y1, y2)).x_enclosure_interval,
-            open(-P.inf, P.inf)
+            P.open(-P.inf, P.inf)
         )
 
     def test_y_enclosure_interval(self):
@@ -357,7 +356,7 @@ def test_y_enclosure_interval(self):
         y1, y2, y3 = 0, 1, 2
         self.assertEqual(
             (rclosed(x1, x3, y1, y3) | rclosed(x2, x4, y1, y2)).y_enclosure_interval,
-            closed(y1, y3)
+            P.closed(y1, y3)
         )
 
     def test_boundaries(self):
@@ -381,7 +380,7 @@ def __init__(self, *args, **kwargs):
     def test_from_interval_product(self):
         x1, x2 = 0, 1
         y1, y2, y3, y4 = 0, 2, 3, 4
-        p = RPolygon.from_interval_product(open(x1, x2), openclosed(y1, y2))
+        p = RPolygon.from_interval_product(P.open(x1, x2), P.openclosed(y1, y2))
         self.assertListEqual(list(p._x_boundaries), [
             RBoundary(-P.inf, P.OPEN),
             RBoundary(x1, P.CLOSED),
@@ -389,16 +388,16 @@ def test_from_interval_product(self):
             RBoundary(P.inf, P.OPEN),
         ])
         self.assertListEqual(p._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), openclosed(y1, y2)],
-            [empty()],
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.openclosed(y1, y2)],
+            [P.empty()],
         ])
         self.assertListEqual(p._free_y_ranges, [
-            [openclosed(-P.inf, y1) | open(y2, P.inf), openclosed(-P.inf, y1) | open(y2, P.inf), open(-P.inf, P.inf)],
-            [openclosed(-P.inf, y1) | open(y2, P.inf), openclosed(-P.inf, y1) | open(y2, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.open(y2, P.inf), P.openclosed(-P.inf, y1) | P.open(y2, P.inf), P.open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.open(y2, P.inf), P.openclosed(-P.inf, y1) | P.open(y2, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
-        p = RPolygon.from_interval_product(open(x1, x2), openclosed(y1, y2) | closedopen(y3, y4))
+        p = RPolygon.from_interval_product(P.open(x1, x2), P.openclosed(y1, y2) | P.closedopen(y3, y4))
         self.assertListEqual(list(p._x_boundaries), [
             RBoundary(-P.inf, P.OPEN),
             RBoundary(x1, P.CLOSED),
@@ -406,17 +405,17 @@ def test_from_interval_product(self):
             RBoundary(P.inf, P.OPEN),
         ])
         self.assertListEqual(p._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), openclosed(y1, y2) | closedopen(y3, y4)],
-            [empty()],
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.openclosed(y1, y2) | P.closedopen(y3, y4)],
+            [P.empty()],
         ])
         self.assertListEqual(p._free_y_ranges, [
-            [openclosed(-P.inf, y1) | open(y2, y3) | closedopen(y4, P.inf),
-             openclosed(-P.inf, y1) | open(y2, y3) | closedopen(y4, P.inf),
-             open(-P.inf, P.inf)],
-            [openclosed(-P.inf, y1) | open(y2, y3) | closedopen(y4, P.inf),
-             openclosed(-P.inf, y1) | open(y2, y3) | closedopen(y4, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.open(y2, y3) | P.closedopen(y4, P.inf),
+             P.openclosed(-P.inf, y1) | P.open(y2, y3) | P.closedopen(y4, P.inf),
+             P.open(-P.inf, P.inf)],
+            [P.openclosed(-P.inf, y1) | P.open(y2, y3) | P.closedopen(y4, P.inf),
+             P.openclosed(-P.inf, y1) | P.open(y2, y3) | P.closedopen(y4, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
 
 
@@ -442,178 +441,178 @@ def test_eq(self):
     def test__add_atomic(self):
         # empty x_interval
         poly = rempty()
-        x_atom = closedopen(1, 0)
-        y_interval = open(1, P.inf)
+        x_atom = P.closedopen(1, 0)
+        y_interval = P.open(1, P.inf)
         poly._add_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[empty()]])
-        self.assertListEqual(poly._free_y_ranges, [[open(-P.inf, P.inf)]])
+        self.assertListEqual(poly._used_y_ranges, [[P.empty()]])
+        self.assertListEqual(poly._free_y_ranges, [[P.open(-P.inf, P.inf)]])
 
         # empty y_interval
         poly = rempty()
-        x_atom = closedopen(-P.inf, 1)
-        y_interval = empty()
+        x_atom = P.closedopen(-P.inf, 1)
+        y_interval = P.empty()
         poly._add_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[empty()]])
-        self.assertListEqual(poly._free_y_ranges, [[open(-P.inf, P.inf)]])
+        self.assertListEqual(poly._used_y_ranges, [[P.empty()]])
+        self.assertListEqual(poly._free_y_ranges, [[P.open(-P.inf, P.inf)]])
 
         # add the whole area
         poly = rempty()
-        x_atom = closedopen(-P.inf, P.inf)
-        y_interval = open(-P.inf, P.inf)
+        x_atom = P.closedopen(-P.inf, P.inf)
+        y_interval = P.open(-P.inf, P.inf)
         poly._add_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[open(-P.inf, P.inf)]])
-        self.assertListEqual(poly._free_y_ranges, [[empty()]])
+        self.assertListEqual(poly._used_y_ranges, [[P.open(-P.inf, P.inf)]])
+        self.assertListEqual(poly._free_y_ranges, [[P.empty()]])
 
         # add half planes
         # # (a) Right half space
         poly = rempty()
-        x_atom = closedopen(0, P.inf)
-        y_interval = open(-P.inf, P.inf)
+        x_atom = P.closedopen(0, P.inf)
+        y_interval = P.open(-P.inf, P.inf)
         poly._add_atomic(x_atom, y_interval)
         self.assertListEqual(poly._used_y_ranges, [
-            [empty(), empty()],
-            [open(-P.inf, P.inf)]
+            [P.empty(), P.empty()],
+            [P.open(-P.inf, P.inf)]
         ])
         self.assertListEqual(poly._free_y_ranges, [
-            [empty(), open(-P.inf, P.inf)],
-            [empty()]
+            [P.empty(), P.open(-P.inf, P.inf)],
+            [P.empty()]
         ])
         # # (b) Left half space
         poly = rempty()
-        x_atom = closedopen(-P.inf, 0)
-        y_interval = open(-P.inf, P.inf)
+        x_atom = P.closedopen(-P.inf, 0)
+        y_interval = P.open(-P.inf, P.inf)
         poly._add_atomic(x_atom, y_interval)
         self.assertListEqual(poly._used_y_ranges, [
-            [empty(), open(-P.inf, P.inf)],
-            [empty()]
+            [P.empty(), P.open(-P.inf, P.inf)],
+            [P.empty()]
         ])
         self.assertListEqual(poly._free_y_ranges, [
-            [empty(), empty()],
-            [open(-P.inf, P.inf)]
+            [P.empty(), P.empty()],
+            [P.open(-P.inf, P.inf)]
         ])
         # # (c) Upper half space
         poly = rempty()
-        x_atom = closedopen(-P.inf, P.inf)
-        y_interval = open(-P.inf, 0)
+        x_atom = P.closedopen(-P.inf, P.inf)
+        y_interval = P.open(-P.inf, 0)
         poly._add_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[open(-P.inf, 0)]])
-        self.assertListEqual(poly._free_y_ranges, [[closedopen(0, P.inf)]])
+        self.assertListEqual(poly._used_y_ranges, [[P.open(-P.inf, 0)]])
+        self.assertListEqual(poly._free_y_ranges, [[P.closedopen(0, P.inf)]])
         # # (d) Lower half space
         poly = rempty()
-        x_atom = closedopen(-P.inf, P.inf)
-        y_interval = closedopen(0, P.inf)
+        x_atom = P.closedopen(-P.inf, P.inf)
+        y_interval = P.closedopen(0, P.inf)
         poly._add_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[closedopen(0, P.inf)]])
-        self.assertListEqual(poly._free_y_ranges, [[open(-P.inf, 0)]])
+        self.assertListEqual(poly._used_y_ranges, [[P.closedopen(0, P.inf)]])
+        self.assertListEqual(poly._free_y_ranges, [[P.open(-P.inf, 0)]])
 
         # add a single bounded polygon
         poly = rempty()
-        x_atom = closedopen(1, 3)
-        y_interval = closedopen(2, 4)
+        x_atom = P.closedopen(1, 3)
+        y_interval = P.closedopen(2, 4)
         poly._add_atomic(x_atom, y_interval)
         self.assertListEqual(poly._used_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), closedopen(2, 4)],
-            [empty()]
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.closedopen(2, 4)],
+            [P.empty()]
         ])
         self.assertListEqual(poly._free_y_ranges, [
-            [open(-P.inf, 2) | closedopen(4, P.inf),
-             open(-P.inf, 2) | closedopen(4, P.inf),
-             open(-P.inf, P.inf)],
-            [open(-P.inf, 2) | closedopen(4, P.inf),
-             open(-P.inf, 2) | closedopen(4, P.inf)],
-            [open(-P.inf, P.inf)]
+            [P.open(-P.inf, 2) | P.closedopen(4, P.inf),
+             P.open(-P.inf, 2) | P.closedopen(4, P.inf),
+             P.open(-P.inf, P.inf)],
+            [P.open(-P.inf, 2) | P.closedopen(4, P.inf),
+             P.open(-P.inf, 2) | P.closedopen(4, P.inf)],
+            [P.open(-P.inf, P.inf)]
         ])
 
     def test__sub_atomic(self):
         # Remove rectangle with empty x_interval from the whole plane
         poly = ropen(-P.inf, P.inf, -P.inf, P.inf)
-        x_atom = closedopen(1, 0)
-        y_interval = open(1, P.inf)
+        x_atom = P.closedopen(1, 0)
+        y_interval = P.open(1, P.inf)
         poly._sub_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[open(-P.inf, P.inf)]])
-        self.assertListEqual(poly._free_y_ranges, [[empty()]])
+        self.assertListEqual(poly._used_y_ranges, [[P.open(-P.inf, P.inf)]])
+        self.assertListEqual(poly._free_y_ranges, [[P.empty()]])
 
         # Remove rectangle with empty y_interval from the whole plane
         poly = ropen(-P.inf, P.inf, -P.inf, P.inf)
-        x_atom = closedopen(-P.inf, 1)
-        y_interval = empty()
+        x_atom = P.closedopen(-P.inf, 1)
+        y_interval = P.empty()
         poly._sub_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[open(-P.inf, P.inf)]])
-        self.assertListEqual(poly._free_y_ranges, [[empty()]])
+        self.assertListEqual(poly._used_y_ranges, [[P.open(-P.inf, P.inf)]])
+        self.assertListEqual(poly._free_y_ranges, [[P.empty()]])
 
         # Remove half planes from the whole plane
         # # (a) Right half space
         poly = ropen(-P.inf, P.inf, -P.inf, P.inf)
-        x_atom = closedopen(0, P.inf)
-        y_interval = open(-P.inf, P.inf)
+        x_atom = P.closedopen(0, P.inf)
+        y_interval = P.open(-P.inf, P.inf)
         poly._sub_atomic(x_atom, y_interval)
         self.assertListEqual(poly._used_y_ranges, [
-            [empty(), open(-P.inf, P.inf)],
-            [empty()]
+            [P.empty(), P.open(-P.inf, P.inf)],
+            [P.empty()]
         ])
         self.assertListEqual(poly._free_y_ranges, [
-            [empty(), empty()],
-            [open(-P.inf, P.inf)]
+            [P.empty(), P.empty()],
+            [P.open(-P.inf, P.inf)]
         ])
         # # (b) Left half space
         poly = ropen(-P.inf, P.inf, -P.inf, P.inf)
-        x_atom = closedopen(-P.inf, 0)
-        y_interval = open(-P.inf, P.inf)
+        x_atom = P.closedopen(-P.inf, 0)
+        y_interval = P.open(-P.inf, P.inf)
         poly._sub_atomic(x_atom, y_interval)
         self.assertListEqual(poly._used_y_ranges, [
-            [empty(), empty()],
-            [open(-P.inf, P.inf)]
+            [P.empty(), P.empty()],
+            [P.open(-P.inf, P.inf)]
         ])
         self.assertListEqual(poly._free_y_ranges, [
-            [empty(), open(-P.inf, P.inf)],
-            [empty()]
+            [P.empty(), P.open(-P.inf, P.inf)],
+            [P.empty()]
         ])
         # # (c) Upper half space
         poly = ropen(-P.inf, P.inf, -P.inf, P.inf)
-        x_atom = closedopen(-P.inf, P.inf)
-        y_interval = open(-P.inf, 0)
+        x_atom = P.closedopen(-P.inf, P.inf)
+        y_interval = P.open(-P.inf, 0)
         poly._sub_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[closedopen(0, P.inf)]])
-        self.assertListEqual(poly._free_y_ranges, [[open(-P.inf, 0)]])
+        self.assertListEqual(poly._used_y_ranges, [[P.closedopen(0, P.inf)]])
+        self.assertListEqual(poly._free_y_ranges, [[P.open(-P.inf, 0)]])
         # # (d) Lower half space
         poly = ropen(-P.inf, P.inf, -P.inf, P.inf)
-        x_atom = closedopen(-P.inf, P.inf)
-        y_interval = closedopen(0, P.inf)
+        x_atom = P.closedopen(-P.inf, P.inf)
+        y_interval = P.closedopen(0, P.inf)
         poly._sub_atomic(x_atom, y_interval)
-        self.assertListEqual(poly._used_y_ranges, [[open(-P.inf, 0)]])
-        self.assertListEqual(poly._free_y_ranges, [[closedopen(0, P.inf)]])
+        self.assertListEqual(poly._used_y_ranges, [[P.open(-P.inf, 0)]])
+        self.assertListEqual(poly._free_y_ranges, [[P.closedopen(0, P.inf)]])
 
         # remove a single bounded rectangle from the plane
         poly = ropen(-P.inf, P.inf, -P.inf, P.inf)
-        x_atom = closedopen(1, 3)
-        y_interval = closedopen(2, 4)
+        x_atom = P.closedopen(1, 3)
+        y_interval = P.closedopen(2, 4)
         poly._sub_atomic(x_atom, y_interval)
         self.assertListEqual(poly._used_y_ranges, [
-            [open(-P.inf, 2) | closedopen(4, P.inf),
-             open(-P.inf, 2) | closedopen(4, P.inf),
-             open(-P.inf, P.inf)],
-            [open(-P.inf, 2) | closedopen(4, P.inf),
-             open(-P.inf, 2) | closedopen(4, P.inf)],
-            [open(-P.inf, P.inf)]
+            [P.open(-P.inf, 2) | P.closedopen(4, P.inf),
+             P.open(-P.inf, 2) | P.closedopen(4, P.inf),
+             P.open(-P.inf, P.inf)],
+            [P.open(-P.inf, 2) | P.closedopen(4, P.inf),
+             P.open(-P.inf, 2) | P.closedopen(4, P.inf)],
+            [P.open(-P.inf, P.inf)]
         ])
         self.assertListEqual(poly._free_y_ranges, [
-            [empty(), empty(), empty()],
-            [empty(), closedopen(2, 4)],
-            [empty()]
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.closedopen(2, 4)],
+            [P.empty()]
         ])
 
         # Another example
         poly = rclosedopen(-P.inf, 0, -P.inf, P.inf)
-        poly._sub_atomic(open(-P.inf, P.inf), open(-P.inf, 0))
-        poly._sub_atomic(open(-P.inf, P.inf), open(1, P.inf))
+        poly._sub_atomic(P.open(-P.inf, P.inf), P.open(-P.inf, 0))
+        poly._sub_atomic(P.open(-P.inf, P.inf), P.open(1, P.inf))
         self.assertListEqual(poly._used_y_ranges, [
-            [empty(), closed(0, 1)],
-            [empty()],
+            [P.empty(), P.closed(0, 1)],
+            [P.empty()],
         ])
         self.assertListEqual(poly._free_y_ranges, [
-            [open(-P.inf, 0) | open(1, P.inf), open(-P.inf, 0) | open(1, P.inf)],
-            [open(-P.inf, P.inf)],
+            [P.open(-P.inf, 0) | P.open(1, P.inf), P.open(-P.inf, 0) | P.open(1, P.inf)],
+            [P.open(-P.inf, P.inf)],
         ])
 
     def test_intersection(self):
@@ -682,8 +681,8 @@ def test_difference(self):
         )
         self.assertListEqual(
             list((p1 - p2).maximal_rectangles()),
-            [RPolygon.from_interval_product(closed(0, 2), closedopen(0, 1)),
-             RPolygon.from_interval_product(closedopen(0, 1), closed(0, 2))]
+            [RPolygon.from_interval_product(P.closed(0, 2), P.closedopen(0, 1)),
+             RPolygon.from_interval_product(P.closedopen(0, 1), P.closed(0, 2))]
         )
 
     def test_rectangle_partition(self):
@@ -697,8 +696,8 @@ def test_rectangle_partition(self):
         poly = ropen(x0, x2, y0, y2) | ropen(x1, x3, y1, y3)
         self.assertListEqual(list(poly.rectangle_partitioning()),
                              [ropen(x0, x3, y1, y2),
-                              RPolygon.from_interval_product(open(x0, x2), openclosed(y0, y1)),
-                              RPolygon.from_interval_product(open(x1, x3), closedopen(y2, y3))])
+                              RPolygon.from_interval_product(P.open(x0, x2), P.openclosed(y0, y1)),
+                              RPolygon.from_interval_product(P.open(x1, x3), P.closedopen(y2, y3))])
 
     def test_maximal_atomic_x_rectangles(self):
         # Add two finite (in x dimension) polygons.
@@ -711,61 +710,61 @@ def test_maximal_atomic_x_rectangles(self):
         #     |  |          |   |       |      |       |   |          |  |
         #     +--+          +---+       +------+       +---+          +--+
         x1, x2, x3, x4 = (1, 2, 3, 4)
-        y_interval_1 = closedopen(-1, 1)
-        y_interval_2 = closedopen(0, 2)
+        y_interval_1 = P.closedopen(-1, 1)
+        y_interval_2 = P.closedopen(0, 2)
         # (1)
         poly = rempty()
-        poly._add_atomic(closedopen(x1, x2), y_interval_1)
-        poly._add_atomic(closedopen(x3, x4), y_interval_2)
+        poly._add_atomic(P.closedopen(x1, x2), y_interval_1)
+        poly._add_atomic(P.closedopen(x3, x4), y_interval_2)
         self.assertListEqual(list(poly._maximal_atomic_x_rectangles()), [
-            (closedopen(x1, x2), y_interval_1),
-            (closedopen(x3, x4), y_interval_2)
+            (P.closedopen(x1, x2), y_interval_1),
+            (P.closedopen(x3, x4), y_interval_2)
         ])
         self.assertListEqual(list((~poly)._maximal_atomic_x_rectangles()), [
-            (open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2)),
-            (open(-P.inf, x3),
-             Interval.from_atomic(~y_interval_1.right, y_interval_1.upper, P.inf, P.OPEN)),
-            (closedopen(x2, P.inf),
-             Interval.from_atomic(P.OPEN, -P.inf, y_interval_2.lower, ~y_interval_2.left)),
-            (open(-P.inf, x1), open(-P.inf, P.inf)),
-            (closedopen(x2, x3), open(-P.inf, P.inf)),
-            (closedopen(x4, P.inf), open(-P.inf, P.inf)),
+            (P.open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2)),
+            (P.open(-P.inf, x3),
+             P.Interval.from_atomic(~y_interval_1.right, y_interval_1.upper, P.inf, P.OPEN)),
+            (P.closedopen(x2, P.inf),
+             P.Interval.from_atomic(P.OPEN, -P.inf, y_interval_2.lower, ~y_interval_2.left)),
+            (P.open(-P.inf, x1), P.open(-P.inf, P.inf)),
+            (P.closedopen(x2, x3), P.open(-P.inf, P.inf)),
+            (P.closedopen(x4, P.inf), P.open(-P.inf, P.inf)),
         ])
         # (2)
         poly = rempty()
-        poly._add_atomic(closedopen(x1, x3), y_interval_1)
-        poly._add_atomic(closedopen(x2, x4), y_interval_2)
+        poly._add_atomic(P.closedopen(x1, x3), y_interval_1)
+        poly._add_atomic(P.closedopen(x2, x4), y_interval_2)
         self.assertListEqual(list(poly._maximal_atomic_x_rectangles()), [
-            (closedopen(x1, x4), y_interval_1 & y_interval_2),
-            (closedopen(x1, x3), y_interval_1),
-            (closedopen(x2, x4), y_interval_2),
-            (closedopen(x2, x3), y_interval_1 | y_interval_2)
+            (P.closedopen(x1, x4), y_interval_1 & y_interval_2),
+            (P.closedopen(x1, x3), y_interval_1),
+            (P.closedopen(x2, x4), y_interval_2),
+            (P.closedopen(x2, x3), y_interval_1 | y_interval_2)
         ])
         self.assertListEqual(list((~poly)._maximal_atomic_x_rectangles()), [
-            (open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2)),
-            (open(-P.inf, x2),
-             Interval.from_atomic(~y_interval_1.right, y_interval_1.upper, P.inf, P.OPEN)),
-            (closedopen(x3, P.inf),
-             Interval.from_atomic(P.OPEN, -P.inf, y_interval_2.lower, ~y_interval_2.left)),
-            (open(-P.inf, x1), open(-P.inf, P.inf)),
-            (closedopen(x4, P.inf), open(-P.inf, P.inf)),
+            (P.open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2)),
+            (P.open(-P.inf, x2),
+             P.Interval.from_atomic(~y_interval_1.right, y_interval_1.upper, P.inf, P.OPEN)),
+            (P.closedopen(x3, P.inf),
+             P.Interval.from_atomic(P.OPEN, -P.inf, y_interval_2.lower, ~y_interval_2.left)),
+            (P.open(-P.inf, x1), P.open(-P.inf, P.inf)),
+            (P.closedopen(x4, P.inf), P.open(-P.inf, P.inf)),
         ])
         # (3)
         poly = rempty()
-        poly._add_atomic(closedopen(x1, x4), y_interval_1)
-        poly._add_atomic(closedopen(x2, x3), y_interval_2)
+        poly._add_atomic(P.closedopen(x1, x4), y_interval_1)
+        poly._add_atomic(P.closedopen(x2, x3), y_interval_2)
         self.assertListEqual(list(poly._maximal_atomic_x_rectangles()), [
-            (closedopen(x1, x4), y_interval_1),
-            (closedopen(x2, x3), y_interval_1 | y_interval_2),
+            (P.closedopen(x1, x4), y_interval_1),
+            (P.closedopen(x2, x3), y_interval_1 | y_interval_2),
         ])
         self.assertListEqual(list((~poly)._maximal_atomic_x_rectangles()), [
-            (open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2)),
-            (open(-P.inf, x2),
-             Interval.from_atomic(~y_interval_1.right, y_interval_1.upper, P.inf, P.OPEN)),
-            (closedopen(x3, P.inf),
-             Interval.from_atomic(~y_interval_2.right, y_interval_1.upper, P.inf, P.OPEN)),
-            (open(-P.inf, x1), open(-P.inf, P.inf)),
-            (closedopen(x4, P.inf), open(-P.inf, P.inf)),
+            (P.open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2)),
+            (P.open(-P.inf, x2),
+             P.Interval.from_atomic(~y_interval_1.right, y_interval_1.upper, P.inf, P.OPEN)),
+            (P.closedopen(x3, P.inf),
+             P.Interval.from_atomic(~y_interval_2.right, y_interval_1.upper, P.inf, P.OPEN)),
+            (P.open(-P.inf, x1), P.open(-P.inf, P.inf)),
+            (P.closedopen(x4, P.inf), P.open(-P.inf, P.inf)),
         ])
         # skip (4) & (5)
 
@@ -779,35 +778,35 @@ def test_maximal_atomic_x_rectangles(self):
         #                              +----+            |      +-|--+
         #                              |    |            +------|-+  |
         #                              +----+                   +----+
-        y_interval_1 = closedopen(4, 6)
-        y_interval_2 = closedopen(1, 3)
-        y_interval_3 = closedopen(2, 5)
+        y_interval_1 = P.closedopen(4, 6)
+        y_interval_2 = P.closedopen(1, 3)
+        y_interval_3 = P.closedopen(2, 5)
         x1, x2, x3, x4, x5, x6 = (1, 2, 3, 4, 5, 6)
         # (1)
         poly = rempty()
-        poly._add_atomic(closedopen(x1, x3), y_interval_1)
-        poly._add_atomic(closedopen(x4, x6), y_interval_2)
-        poly._add_atomic(closedopen(x2, x5), y_interval_3)
+        poly._add_atomic(P.closedopen(x1, x3), y_interval_1)
+        poly._add_atomic(P.closedopen(x4, x6), y_interval_2)
+        poly._add_atomic(P.closedopen(x2, x5), y_interval_3)
         self.assertListEqual(list(poly._maximal_atomic_x_rectangles()), [
-            (closedopen(x1, x5), closedopen(x4, x5)),
-            (closedopen(x2, x6), closedopen(x2, x3)),
-            (closedopen(x2, x5), closedopen(x2, x5)),
-            (closedopen(x1, x3), closedopen(x4, x6)),
-            (closedopen(x4, x6), closedopen(x1, x3)),
-            (closedopen(x2, x3), closedopen(x2, x6)),
-            (closedopen(x4, x5), closedopen(x1, x5)),
+            (P.closedopen(x1, x5), P.closedopen(x4, x5)),
+            (P.closedopen(x2, x6), P.closedopen(x2, x3)),
+            (P.closedopen(x2, x5), P.closedopen(x2, x5)),
+            (P.closedopen(x1, x3), P.closedopen(x4, x6)),
+            (P.closedopen(x4, x6), P.closedopen(x1, x3)),
+            (P.closedopen(x2, x3), P.closedopen(x2, x6)),
+            (P.closedopen(x4, x5), P.closedopen(x1, x5)),
         ])
         self.assertListEqual(list((~poly)._maximal_atomic_x_rectangles()), [
-            (open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2 | y_interval_3)),
-            (open(-P.inf, x4),
-             Interval.from_atomic(P.OPEN, -P.inf, y_interval_3.lower, ~y_interval_3.left)),
-            (closedopen(x3, P.inf),
-             Interval.from_atomic(~y_interval_3.right, y_interval_3.upper, P.inf, P.OPEN)),
-            (open(-P.inf, x2), Interval.from_atomic(P.OPEN, -P.inf, y_interval_1.lower, ~y_interval_1.left)),
-            (closedopen(x5, P.inf),
-             Interval.from_atomic(~y_interval_2.right, y_interval_2.upper, P.inf, P.OPEN)),
-            (open(-P.inf, x1), open(-P.inf, P.inf)),
-            (closedopen(x6, P.inf), open(-P.inf, P.inf)),
+            (P.open(-P.inf, P.inf), ~(y_interval_1 | y_interval_2 | y_interval_3)),
+            (P.open(-P.inf, x4),
+             P.Interval.from_atomic(P.OPEN, -P.inf, y_interval_3.lower, ~y_interval_3.left)),
+            (P.closedopen(x3, P.inf),
+             P.Interval.from_atomic(~y_interval_3.right, y_interval_3.upper, P.inf, P.OPEN)),
+            (P.open(-P.inf, x2), P.Interval.from_atomic(P.OPEN, -P.inf, y_interval_1.lower, ~y_interval_1.left)),
+            (P.closedopen(x5, P.inf),
+             P.Interval.from_atomic(~y_interval_2.right, y_interval_2.upper, P.inf, P.OPEN)),
+            (P.open(-P.inf, x1), P.open(-P.inf, P.inf)),
+            (P.closedopen(x6, P.inf), P.open(-P.inf, P.inf)),
         ])
 
         # Test the following polygon addition sequence, where the third rectangle directly touches the second one.
@@ -816,12 +815,12 @@ def test_maximal_atomic_x_rectangles(self):
         #    +----+          +----+     +----+        +----+    +---------+
         #    |    |    ->    |    |     |    |    ->  |    |    |         |
         #    +----+          +----+     +----+        +----+    +---------+
-        y_interval = closedopen(0, 1)
+        y_interval = P.closedopen(0, 1)
         x1, x2, x3, x4, x5 = 1, 2, 3, 4, 5
         poly = rempty()
-        # poly._add_atomic(closedopen(x1, x2), y_interval)
-        poly._add_atomic(closedopen(x3, x4), y_interval)
-        poly._add_atomic(closedopen(x4, x5), y_interval)
+        # poly._add_atomic(P.closedopen(x1, x2), y_interval)
+        poly._add_atomic(P.closedopen(x3, x4), y_interval)
+        poly._add_atomic(P.closedopen(x4, x5), y_interval)
 
         # Add three sectors in each six possible different orders and test if the result underlying
         # data structure of RPolygon is the same for all.
@@ -831,14 +830,14 @@ def test_maximal_atomic_x_rectangles(self):
         #      |  |  |  |  |
         #      +--+--+--+--+
         x1, x2, x3, x4, x5 = (1, 2, 3, 4, 5)
-        y_interval = closedopen(0, 1)
+        y_interval = P.closedopen(0, 1)
         x_lims = [(x1, x3), (x2, x4), (x3, x5)]
         y_intervals = [y_interval, y_interval, y_interval]
         poly_list = []
         for arrangement in permutations(zip(x_lims, y_intervals)):
             poly = rempty()
             for (x_a, x_b), y_int in arrangement:
-                poly._add_atomic(closedopen(x_a, x_b), y_int)
+                poly._add_atomic(P.closedopen(x_a, x_b), y_int)
             poly_list.append(poly)
         for poly_1, poly_2 in combinations(poly_list, 2):
             self.assertListEqual(poly_1._used_y_ranges, poly_2._used_y_ranges)
@@ -956,8 +955,8 @@ def test_boundary(self):
         for constructor in [rclosed, ropen, rclosedopen, ropenclosed]:
             poly = constructor(x0, x1, y0, y1)
             self.assertEqual(poly.boundary(),
-                             rproduct(closed(x0, x1), singleton(y0) | singleton(y1))
-                             | rproduct(singleton(x0) | singleton(x1), closed(y0, y1)))
+                             rproduct(P.closed(x0, x1), P.singleton(y0) | P.singleton(y1))
+                             | rproduct(P.singleton(x0) | P.singleton(x1), P.closed(y0, y1)))
 
         # Two overlapping rectangles
         #        +--+
@@ -969,14 +968,14 @@ def test_boundary(self):
         for constructor in [rclosed, ropen, rclosedopen, ropenclosed]:
             poly = constructor(x0, x2, y0, y2) | constructor(x1, x3, y1, y3)
             self.assertEqual(poly.boundary(),
-                             rproduct(singleton(x0), closed(y0, y2))
-                             | rproduct(closed(x0, x1), singleton(y2))
-                             | rproduct(singleton(x1), closed(y2, y3))
-                             | rproduct(closed(x1, x3), singleton(y3))
-                             | rproduct(singleton(x3), closed(y1, y3))
-                             | rproduct(closed(x2, x3), singleton(y1))
-                             | rproduct(singleton(x2), closed(y0, y1))
-                             | rproduct(closed(x0, x2), singleton(y0)))
+                             rproduct(P.singleton(x0), P.closed(y0, y2))
+                             | rproduct(P.closed(x0, x1), P.singleton(y2))
+                             | rproduct(P.singleton(x1), P.closed(y2, y3))
+                             | rproduct(P.closed(x1, x3), P.singleton(y3))
+                             | rproduct(P.singleton(x3), P.closed(y1, y3))
+                             | rproduct(P.closed(x2, x3), P.singleton(y1))
+                             | rproduct(P.singleton(x2), P.closed(y0, y1))
+                             | rproduct(P.closed(x0, x2), P.singleton(y0)))
 
         # Rectangle with a hole
         #       +------+
@@ -990,10 +989,10 @@ def test_boundary(self):
             poly = constructor(x0, x3, y0, y3)
             poly -= constructor(x1, x2, y1, y2)
             self.assertEqual(poly.boundary(),
-                             rproduct(closed(x0, x3), singleton(y0) | singleton(y3))
-                             | rproduct(singleton(x0) | singleton(x3), closed(y0, y3))
-                             | rproduct(closed(x1, x2), singleton(y1) | singleton(y2))
-                             | rproduct(singleton(x1) | singleton(x2), closed(y1, y2)))
+                             rproduct(P.closed(x0, x3), P.singleton(y0) | P.singleton(y3))
+                             | rproduct(P.singleton(x0) | P.singleton(x3), P.closed(y0, y3))
+                             | rproduct(P.closed(x1, x2), P.singleton(y1) | P.singleton(y2))
+                             | rproduct(P.singleton(x1) | P.singleton(x2), P.closed(y1, y2)))
 
 
 class TestIntervalTreeFunctions(unittest.TestCase):
@@ -1005,18 +1004,18 @@ def test__traverse_diagonally(self):
         adj_y_interval = lambda curr, l_parent, r_parent: curr
 
         boundaries = [RBoundary(-P.inf, P.OPEN), RBoundary(P.inf, P.OPEN)]
-        interval_tree = [[empty()]]
+        interval_tree = [[P.empty()]]
         with self.assertRaises(StopIteration):
             next(_traverse_diagonally(boundaries, interval_tree, next_accumulator, adj_y_interval))
 
         boundaries = [RBoundary(-P.inf, P.OPEN), RBoundary(0, P.CLOSED),
                       RBoundary(1, P.OPEN), RBoundary(P.inf, P.OPEN)]
         interval_tree = [
-            [empty(), empty(), empty()],
-            [empty(), open(0, 1)],
-            [empty()]
+            [P.empty(), P.empty(), P.empty()],
+            [P.empty(), P.open(0, 1)],
+            [P.empty()]
         ]
         iterator = _traverse_diagonally(boundaries, interval_tree, next_accumulator, adj_y_interval)
-        self.assertEqual(next(iterator), (open(0, 1), open(0, 1)))
+        self.assertEqual(next(iterator), (P.open(0, 1), P.open(0, 1)))
         with self.assertRaises(StopIteration):
             next(iterator)