# Copyright 2014 MINES ParisTech
#
# This file is part of LinPy.
#
# LinPy 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 3 of the License, or
# (at your option) any later version.
#
# LinPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with LinPy.  If not, see <http://www.gnu.org/licenses/>.

import ast
import functools
import re
import math

from fractions import Fraction

from . import islhelper
from .geometry import GeometricObject, Point, Vector
from .islhelper import libisl
from .linexprs import LinExpr, Symbol


__all__ = [
    'And',
    'Domain',
    'Not',
    'Or',
]


@functools.total_ordering
class Domain(GeometricObject):
    """
    A domain is a union of polyhedra. Unlike polyhedra, domains allow exact
    computation of union, subtraction and complementary operations.

    A domain with a unique polyhedron is automatically subclassed as a
    Polyhedron instance.
    """

    __slots__ = (
        '_polyhedra',
        '_symbols',
        '_dimension',
    )

    def __new__(cls, *polyhedra):
        """
        Return a domain from a sequence of polyhedra.

        >>> square1 = Polyhedron('0 <= x <= 2, 0 <= y <= 2')
        >>> square2 = Polyhedron('1 <= x <= 3, 1 <= y <= 3')
        >>> dom = Domain(square1, square2)
        >>> dom
        Or(And(x <= 2, 0 <= x, y <= 2, 0 <= y),
           And(x <= 3, 1 <= x, y <= 3, 1 <= y))

        It is also possible to build domains from polyhedra using arithmetic
        operators Domain.__or__(), Domain.__invert__() or functions Or() and
        Not(), using one of the following instructions:

        >>> dom = square1 | square2
        >>> dom = Or(square1, square2)

        Alternatively, a domain can be built from a string:

        >>> dom = Domain('0 <= x <= 2, 0 <= y <= 2; 1 <= x <= 3, 1 <= y <= 3')

        Finally, a domain can be built from a GeometricObject instance, calling
        the GeometricObject.asdomain() method.
        """
        from .polyhedra import Polyhedron
        if len(polyhedra) == 1:
            argument = polyhedra[0]
            if isinstance(argument, str):
                return cls.fromstring(argument)
            elif isinstance(argument, GeometricObject):
                return argument.aspolyhedron()
            else:
                raise TypeError('argument must be a string '
                                'or a GeometricObject instance')
        else:
            for polyhedron in polyhedra:
                if not isinstance(polyhedron, Polyhedron):
                    raise TypeError('arguments must be Polyhedron instances')
            symbols = cls._xsymbols(polyhedra)
            islset = cls._toislset(polyhedra, symbols)
            return cls._fromislset(islset, symbols)

    @classmethod
    def _xsymbols(cls, iterator):
        """
        Return the ordered tuple of symbols present in iterator.
        """
        symbols = set()
        for item in iterator:
            symbols.update(item.symbols)
        return tuple(sorted(symbols, key=Symbol.sortkey))

    @property
    def polyhedra(self):
        """
        The tuple of polyhedra present in the domain.
        """
        return self._polyhedra

    @property
    def symbols(self):
        """
        The tuple of symbols present in the domain equations, sorted according
        to Symbol.sortkey().
        """
        return self._symbols

    @property
    def dimension(self):
        """
        The dimension of the domain, i.e. the number of symbols present in it.
        """
        return self._dimension

    def isempty(self):
        """
        Return True if the domain is empty, that is, equal to Empty.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        empty = bool(libisl.isl_set_is_empty(islset))
        libisl.isl_set_free(islset)
        return empty

    def __bool__(self):
        """
        Return True if the domain is non-empty.
        """
        return not self.isempty()

    def isuniverse(self):
        """
        Return True if the domain is universal, that is, equal to Universe.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        universe = bool(libisl.isl_set_plain_is_universe(islset))
        libisl.isl_set_free(islset)
        return universe

    def isbounded(self):
        """
        Return True if the domain is bounded.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        bounded = bool(libisl.isl_set_is_bounded(islset))
        libisl.isl_set_free(islset)
        return bounded

    def __eq__(self, other):
        """
        Return True if two domains are equal.
        """
        if isinstance(other, Domain):
            symbols = self._xsymbols([self, other])
            islset1 = self._toislset(self.polyhedra, symbols)
            islset2 = other._toislset(other.polyhedra, symbols)
            equal = bool(libisl.isl_set_is_equal(islset1, islset2))
            libisl.isl_set_free(islset1)
            libisl.isl_set_free(islset2)
            return equal
        return NotImplemented

    def isdisjoint(self, other):
        """
        Return True if two domains have a null intersection.
        """
        if not isinstance(other, Domain):
            raise TypeError('other must be a Domain instance')
        symbols = self._xsymbols([self, other])
        islset1 = self._toislset(self.polyhedra, symbols)
        islset2 = self._toislset(other.polyhedra, symbols)
        equal = bool(libisl.isl_set_is_disjoint(islset1, islset2))
        libisl.isl_set_free(islset1)
        libisl.isl_set_free(islset2)
        return equal

    def issubset(self, other):
        """
        Report whether another domain contains the domain.
        """
        return self < other

    def __le__(self, other):
        if isinstance(other, Domain):
            symbols = self._xsymbols([self, other])
            islset1 = self._toislset(self.polyhedra, symbols)
            islset2 = self._toislset(other.polyhedra, symbols)
            equal = bool(libisl.isl_set_is_subset(islset1, islset2))
            libisl.isl_set_free(islset1)
            libisl.isl_set_free(islset2)
            return equal
        return NotImplemented
    __le__.__doc__ = issubset.__doc__

    def __lt__(self, other):
        """
        Report whether another domain is contained within the domain.
        """
        if isinstance(other, Domain):
            symbols = self._xsymbols([self, other])
            islset1 = self._toislset(self.polyhedra, symbols)
            islset2 = self._toislset(other.polyhedra, symbols)
            equal = bool(libisl.isl_set_is_strict_subset(islset1, islset2))
            libisl.isl_set_free(islset1)
            libisl.isl_set_free(islset2)
            return equal
        return NotImplemented

    def complement(self):
        """
        Return the complementary domain of the domain.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islset = libisl.isl_set_complement(islset)
        return self._fromislset(islset, self.symbols)

    def __invert__(self):
        return self.complement()
    __invert__.__doc__ = complement.__doc__

    def make_disjoint(self):
        """
        Return an equivalent domain, whose polyhedra are disjoint.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islset = libisl.isl_set_make_disjoint(islset)
        return self._fromislset(islset, self.symbols)

    def coalesce(self):
        """
        Simplify the representation of the domain by trying to combine pairs of
        polyhedra into a single polyhedron, and return the resulting domain.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islset = libisl.isl_set_coalesce(islset)
        return self._fromislset(islset, self.symbols)

    def detect_equalities(self):
        """
        Simplify the representation of the domain by detecting implicit
        equalities, and return the resulting domain.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islset = libisl.isl_set_detect_equalities(islset)
        return self._fromislset(islset, self.symbols)

    def remove_redundancies(self):
        """
        Remove redundant constraints in the domain, and return the resulting
        domain.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islset = libisl.isl_set_remove_redundancies(islset)
        return self._fromislset(islset, self.symbols)

    def aspolyhedron(self):
        from .polyhedra import Polyhedron
        islset = self._toislset(self.polyhedra, self.symbols)
        islbset = libisl.isl_set_polyhedral_hull(islset)
        return Polyhedron._fromislbasicset(islbset, self.symbols)

    def asdomain(self):
        return self

    def project(self, symbols):
        """
        Project out the sequence of symbols given in arguments, and return the
        resulting domain.
        """
        symbols = list(symbols)
        for symbol in symbols:
            if not isinstance(symbol, Symbol):
                raise TypeError('symbols must be Symbol instances')
        islset = self._toislset(self.polyhedra, self.symbols)
        n = 0
        for index, symbol in reversed(list(enumerate(self.symbols))):
            if symbol in symbols:
                n += 1
            elif n > 0:
                islset = libisl.isl_set_project_out(
                    islset, libisl.isl_dim_set, index + 1, n)
                n = 0
        if n > 0:
            islset = libisl.isl_set_project_out(
                islset, libisl.isl_dim_set, 0, n)
        symbols = [symbol for symbol in self.symbols if symbol not in symbols]
        return Domain._fromislset(islset, symbols)

    def sample(self):
        """
        Return a sample of the domain, as an integer instance of Point. If the
        domain is empty, a ValueError exception is raised.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islpoint = libisl.isl_set_sample_point(islset)
        if bool(libisl.isl_point_is_void(islpoint)):
            libisl.isl_point_free(islpoint)
            raise ValueError('domain must be non-empty')
        point = {}
        for index, symbol in enumerate(self.symbols):
            coordinate = libisl.isl_point_get_coordinate_val(
                islpoint, libisl.isl_dim_set, index)
            coordinate = islhelper.isl_val_to_int(coordinate)
            point[symbol] = coordinate
        libisl.isl_point_free(islpoint)
        return point

    def intersection(self, *others):
        """
        Return the intersection of two or more domains as a new domain. As an
        alternative, function And() can be used.
        """
        result = self
        for other in others:
            result &= other
        return result

    def __and__(self, other):
        if isinstance(other, Domain):
            symbols = self._xsymbols([self, other])
            islset1 = self._toislset(self.polyhedra, symbols)
            islset2 = other._toislset(other.polyhedra, symbols)
            islset = libisl.isl_set_intersect(islset1, islset2)
            return self._fromislset(islset, symbols)
        return NotImplemented
    __and__.__doc__ = intersection.__doc__

    def union(self, *others):
        """
        Return the union of two or more domains as a new domain. As an
        alternative, function Or() can be used.
        """
        result = self
        for other in others:
            result |= other
        return result

    def __or__(self, other):
        if isinstance(other, Domain):
            symbols = self._xsymbols([self, other])
            islset1 = self._toislset(self.polyhedra, symbols)
            islset2 = other._toislset(other.polyhedra, symbols)
            islset = libisl.isl_set_union(islset1, islset2)
            return self._fromislset(islset, symbols)
        return NotImplemented
    __or__.__doc__ = union.__doc__

    def __add__(self, other):
        return self | other
    __add__.__doc__ = union.__doc__

    def difference(self, other):
        """
        Return the difference of two domains as a new domain.
        """
        return self - other

    def __sub__(self, other):
        if isinstance(other, Domain):
            symbols = self._xsymbols([self, other])
            islset1 = self._toislset(self.polyhedra, symbols)
            islset2 = other._toislset(other.polyhedra, symbols)
            islset = libisl.isl_set_subtract(islset1, islset2)
            return self._fromislset(islset, symbols)
        return NotImplemented
    __sub__.__doc__ = difference.__doc__

    def lexmin(self):
        """
        Return the lexicographic minimum of the elements in the domain.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islset = libisl.isl_set_lexmin(islset)
        return self._fromislset(islset, self.symbols)

    def lexmax(self):
        """
        Return the lexicographic maximum of the elements in the domain.
        """
        islset = self._toislset(self.polyhedra, self.symbols)
        islset = libisl.isl_set_lexmax(islset)
        return self._fromislset(islset, self.symbols)

    if islhelper.isl_version >= '0.13':
        _RE_COORDINATE = re.compile(r'\((?P<num>\-?\d+)\)(/(?P<den>\d+))?')
    else:
        _RE_COORDINATE = None

    def vertices(self):
        """
        Return the vertices of the domain, as a list of rational instances of
        Point.
        """
        if not self.isbounded():
            raise ValueError('domain must be bounded')
        islbset = self._toislbasicset(self.equalities, self.inequalities,
                                      self.symbols)
        vertices = libisl.isl_basic_set_compute_vertices(islbset)
        vertices = islhelper.isl_vertices_vertices(vertices)
        points = []
        for vertex in vertices:
            expression = libisl.isl_vertex_get_expr(vertex)
            coordinates = []
            if self._RE_COORDINATE is None:
                constraints = islhelper.isl_basic_set_constraints(expression)
                for constraint in constraints:
                    constant = libisl.isl_constraint_get_constant_val(
                        constraint)
                    constant = islhelper.isl_val_to_int(constant)
                    for index, symbol in enumerate(self.symbols):
                        coefficient = \
                            libisl.isl_constraint_get_coefficient_val(
                                constraint, libisl.isl_dim_set, index)
                        coefficient = islhelper.isl_val_to_int(coefficient)
                        if coefficient != 0:
                            coordinate = -Fraction(constant, coefficient)
                            coordinates.append((symbol, coordinate))
            else:
                string = islhelper.isl_multi_aff_to_str(expression)
                matches = self._RE_COORDINATE.finditer(string)
                for symbol, match in zip(self.symbols, matches):
                    numerator = int(match.group('num'))
                    denominator = match.group('den')
                    denominator = \
                        1 if denominator is None else int(denominator)
                    coordinate = Fraction(numerator, denominator)
                    coordinates.append((symbol, coordinate))
            points.append(Point(coordinates))
        return points

    def points(self):
        """
        Return the integer points of a bounded domain, as a list of integer
        instances of Point. If the domain is not bounded, a ValueError
        exception is raised.
        """
        if not self.isbounded():
            raise ValueError('domain must be bounded')
        islset = self._toislset(self.polyhedra, self.symbols)
        islpoints = islhelper.isl_set_points(islset)
        points = []
        for islpoint in islpoints:
            coordinates = {}
            for index, symbol in enumerate(self.symbols):
                coordinate = libisl.isl_point_get_coordinate_val(
                    islpoint, libisl.isl_dim_set, index)
                coordinate = islhelper.isl_val_to_int(coordinate)
                coordinates[symbol] = coordinate
            points.append(Point(coordinates))
        return points

    def __contains__(self, point):
        """
        Return True if the point if contained within the domain.
        """
        for polyhedron in self.polyhedra:
            if point in polyhedron:
                return True
        return False

    @classmethod
    def _polygon_inner_point(cls, points):
        symbols = points[0].symbols
        coordinates = {symbol: 0 for symbol in symbols}
        for point in points:
            for symbol, coordinate in point.coordinates():
                coordinates[symbol] += coordinate
        for symbol in symbols:
            coordinates[symbol] /= len(points)
        return Point(coordinates)

    @classmethod
    def _sort_polygon_2d(cls, points):
        if len(points) <= 3:
            return points
        o = cls._polygon_inner_point(points)
        angles = {}
        for m in points:
            om = Vector(o, m)
            dx, dy = (coordinate for symbol, coordinate in om.coordinates())
            angle = math.atan2(dy, dx)
            angles[m] = angle
        return sorted(points, key=angles.get)

    @classmethod
    def _sort_polygon_3d(cls, points):
        if len(points) <= 3:
            return points
        o = cls._polygon_inner_point(points)
        a = points[0]
        oa = Vector(o, a)
        norm_oa = oa.norm()
        for b in points[1:]:
            ob = Vector(o, b)
            u = oa.cross(ob)
            if not u.isnull():
                u = u.asunit()
                break
        else:
            raise ValueError('degenerate polygon')
        angles = {a: 0.}
        for m in points[1:]:
            om = Vector(o, m)
            normprod = norm_oa * om.norm()
            cosinus = max(oa.dot(om) / normprod, -1.)
            sinus = u.dot(oa.cross(om)) / normprod
            angle = math.acos(cosinus)
            angle = math.copysign(angle, sinus)
            angles[m] = angle
        return sorted(points, key=angles.get)

    def faces(self):
        """
        Return the list of faces of a bounded domain. Each face is represented
        by a list of vertices, in the form of rational instances of Point. If
        the domain is not bounded, a ValueError exception is raised.
        """
        faces = []
        for polyhedron in self.polyhedra:
            vertices = polyhedron.vertices()
            for constraint in polyhedron.constraints:
                face = []
                for vertex in vertices:
                    if constraint.subs(vertex.coordinates()) == 0:
                        face.append(vertex)
                if len(face) >= 3:
                    faces.append(face)
        return faces

    def _plot_2d(self, plot=None, **kwargs):
        import matplotlib.pyplot as plt
        from matplotlib.patches import Polygon
        if plot is None:
            fig = plt.figure()
            plot = fig.add_subplot(1, 1, 1)
        xmin, xmax = plot.get_xlim()
        ymin, ymax = plot.get_ylim()
        for polyhedron in self.polyhedra:
            vertices = polyhedron._sort_polygon_2d(polyhedron.vertices())
            xys = [tuple(vertex.values()) for vertex in vertices]
            xs, ys = zip(*xys)
            xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
            ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
            plot.add_patch(Polygon(xys, closed=True, **kwargs))
        plot.set_xlim(xmin, xmax)
        plot.set_ylim(ymin, ymax)
        return plot

    def _plot_3d(self, plot=None, **kwargs):
        import matplotlib.pyplot as plt
        from mpl_toolkits.mplot3d import Axes3D
        from mpl_toolkits.mplot3d.art3d import Poly3DCollection
        if plot is None:
            fig = plt.figure()
            axes = Axes3D(fig)
        else:
            axes = plot
        xmin, xmax = axes.get_xlim()
        ymin, ymax = axes.get_ylim()
        zmin, zmax = axes.get_zlim()
        poly_xyzs = []
        for vertices in self.faces():
            vertices = self._sort_polygon_3d(vertices)
            vertices.append(vertices[0])
            face_xyzs = [tuple(vertex.values()) for vertex in vertices]
            xs, ys, zs = zip(*face_xyzs)
            xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
            ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
            zmin, zmax = min(zmin, float(min(zs))), max(zmax, float(max(zs)))
            poly_xyzs.append(face_xyzs)
        collection = Poly3DCollection(poly_xyzs, **kwargs)
        axes.add_collection3d(collection)
        axes.set_xlim(xmin, xmax)
        axes.set_ylim(ymin, ymax)
        axes.set_zlim(zmin, zmax)
        return axes

    def plot(self, plot=None, **kwargs):
        """
        Plot a 2D or 3D domain using matplotlib. Draw it to the current plot
        object if present, otherwise create a new one. options are keyword
        arguments passed to the matplotlib drawing functions, they can be used
        to set the drawing color for example. Raise ValueError is the domain is
        not 2D or 3D.
        """
        if not self.isbounded():
            raise ValueError('domain must be bounded')
        elif self.dimension == 2:
            return self._plot_2d(plot=plot, **kwargs)
        elif self.dimension == 3:
            return self._plot_3d(plot=plot, **kwargs)
        else:
            raise ValueError('domain must be two or three-dimensional')

    def subs(self, symbol, expression=None):
        """
        Substitute the given symbol by an expression in the domain constraints.
        To perform multiple substitutions at once, pass a sequence or a
        dictionary of (old, new) pairs to subs. The syntax of this function is
        similar to LinExpr.subs().
        """
        polyhedra = [polyhedron.subs(symbol, expression)
                     for polyhedron in self.polyhedra]
        return Domain(*polyhedra)

    @classmethod
    def _fromislset(cls, islset, symbols):
        from .polyhedra import Polyhedron
        islset = libisl.isl_set_remove_divs(islset)
        islbsets = islhelper.isl_set_basic_sets(islset)
        libisl.isl_set_free(islset)
        polyhedra = []
        for islbset in islbsets:
            polyhedron = Polyhedron._fromislbasicset(islbset, symbols)
            polyhedra.append(polyhedron)
        if len(polyhedra) == 0:
            from .polyhedra import Empty
            return Empty
        elif len(polyhedra) == 1:
            return polyhedra[0]
        else:
            self = object().__new__(Domain)
            self._polyhedra = tuple(polyhedra)
            self._symbols = cls._xsymbols(polyhedra)
            self._dimension = len(self._symbols)
            return self

    @classmethod
    def _toislset(cls, polyhedra, symbols):
        polyhedron = polyhedra[0]
        islbset = polyhedron._toislbasicset(
            polyhedron.equalities, polyhedron.inequalities, symbols)
        islset1 = libisl.isl_set_from_basic_set(islbset)
        for polyhedron in polyhedra[1:]:
            islbset = polyhedron._toislbasicset(
                polyhedron.equalities, polyhedron.inequalities, symbols)
            islset2 = libisl.isl_set_from_basic_set(islbset)
            islset1 = libisl.isl_set_union(islset1, islset2)
        return islset1

    @classmethod
    def _fromast(cls, node):
        from .polyhedra import Polyhedron
        if isinstance(node, ast.Module) and len(node.body) == 1:
            return cls._fromast(node.body[0])
        elif isinstance(node, ast.Expr):
            return cls._fromast(node.value)
        elif isinstance(node, ast.UnaryOp):
            domain = cls._fromast(node.operand)
            if isinstance(node.operand, ast.invert):
                return Not(domain)
        elif isinstance(node, ast.BinOp):
            domain1 = cls._fromast(node.left)
            domain2 = cls._fromast(node.right)
            if isinstance(node.op, ast.BitAnd):
                return And(domain1, domain2)
            elif isinstance(node.op, ast.BitOr):
                return Or(domain1, domain2)
        elif isinstance(node, ast.Compare):
            equalities = []
            inequalities = []
            left = LinExpr._fromast(node.left)
            for i in range(len(node.ops)):
                op = node.ops[i]
                right = LinExpr._fromast(node.comparators[i])
                if isinstance(op, ast.Lt):
                    inequalities.append(right - left - 1)
                elif isinstance(op, ast.LtE):
                    inequalities.append(right - left)
                elif isinstance(op, ast.Eq):
                    equalities.append(left - right)
                elif isinstance(op, ast.GtE):
                    inequalities.append(left - right)
                elif isinstance(op, ast.Gt):
                    inequalities.append(left - right - 1)
                else:
                    break
                left = right
            else:
                return Polyhedron(equalities, inequalities)
        raise SyntaxError('invalid syntax')

    _RE_BRACES = re.compile(r'^\{\s*|\s*\}$')
    _RE_EQ = re.compile(r'([^<=>])=([^<=>])')
    _RE_AND = re.compile(r'\band\b|,|&&|/\\|∧|∩')
    _RE_OR = re.compile(r'\bor\b|;|\|\||\\/|∨|∪')
    _RE_NOT = re.compile(r'\bnot\b|!|¬')
    _RE_NUM_VAR = LinExpr._RE_NUM_VAR
    _RE_OPERATORS = re.compile(r'(&|\||~)')

    @classmethod
    def fromstring(cls, string):
        """
        Create a domain from a string. Raise SyntaxError if the string is not
        properly formatted.
        """
        # Remove curly brackets.
        string = cls._RE_BRACES.sub(r'', string)
        # Replace '=' by '=='.
        string = cls._RE_EQ.sub(r'\1==\2', string)
        # Replace 'and', 'or', 'not'.
        string = cls._RE_AND.sub(r' & ', string)
        string = cls._RE_OR.sub(r' | ', string)
        string = cls._RE_NOT.sub(r' ~', string)
        # Add implicit multiplication operators, e.g. '5x' -> '5*x'.
        string = cls._RE_NUM_VAR.sub(r'\1*\2', string)
        # Add parentheses to force precedence.
        tokens = cls._RE_OPERATORS.split(string)
        for i, token in enumerate(tokens):
            if i % 2 == 0:
                token = '({})'.format(token)
                tokens[i] = token
        string = ''.join(tokens)
        tree = ast.parse(string, 'eval')
        return cls._fromast(tree)

    def __repr__(self):
        assert len(self.polyhedra) >= 2
        strings = [repr(polyhedron) for polyhedron in self.polyhedra]
        return 'Or({})'.format(', '.join(strings))

    @classmethod
    def fromsympy(cls, expression):
        """
        Create a domain from a SymPy expression.
        """
        import sympy
        from .polyhedra import Lt, Le, Eq, Ne, Ge, Gt
        funcmap = {
            sympy.And: And, sympy.Or: Or, sympy.Not: Not,
            sympy.Lt: Lt, sympy.Le: Le,
            sympy.Eq: Eq, sympy.Ne: Ne,
            sympy.Ge: Ge, sympy.Gt: Gt,
        }
        if expression.func in funcmap:
            args = [Domain.fromsympy(arg) for arg in expression.args]
            return funcmap[expression.func](*args)
        elif isinstance(expression, sympy.Expr):
            return LinExpr.fromsympy(expression)
        raise ValueError('non-domain expression: {!r}'.format(expression))

    def tosympy(self):
        """
        Convert the domain to a SymPy expression.
        """
        import sympy
        polyhedra = [polyhedron.tosympy() for polyhedron in self.polyhedra]
        return sympy.Or(*polyhedra)


def And(*domains):
    """
    Create the intersection domain of the domains given in arguments.
    """
    if len(domains) == 0:
        from .polyhedra import Universe
        return Universe
    else:
        return domains[0].intersection(*domains[1:])


def Or(*domains):
    """
    Create the union domain of the domains given in arguments.
    """
    if len(domains) == 0:
        from .polyhedra import Empty
        return Empty
    else:
        return domains[0].union(*domains[1:])


def Not(domain):
    """
    Create the complementary domain of the domain given in argument.
    """
    return ~domain