238248b63440ce4110721b1e74222ecdce672c89
[linpy.git] / linpy / domains.py
1 # Copyright 2014 MINES ParisTech
2 #
3 # This file is part of LinPy.
4 #
5 # LinPy is free software: you can redistribute it and/or modify
6 # it under the terms of the GNU General Public License as published by
7 # the Free Software Foundation, either version 3 of the License, or
8 # (at your option) any later version.
9 #
10 # LinPy is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 # GNU General Public License for more details.
14 #
15 # You should have received a copy of the GNU General Public License
16 # along with LinPy. If not, see <http://www.gnu.org/licenses/>.
17
18 """
19 Polyhedral domains
20
21 This module provides classes and functions to deal with polyhedral
22 domains, i.e. unions of polyhedra.
23 """
24
25 import ast
26 import functools
27 import re
28 import math
29
30 from fractions import Fraction
31
32 from . import islhelper
33 from .islhelper import mainctx, libisl
34 from .linexprs import LinExpr, Symbol, Rational
35 from .geometry import GeometricObject, Point, Vector
36
37
38 __all__ = [
39 'Domain',
40 'And', 'Or', 'Not',
41 ]
42
43
44 @functools.total_ordering
45 class Domain(GeometricObject):
46 """
47 This class represents polyhedral domains, i.e. unions of polyhedra.
48 """
49
50 __slots__ = (
51 '_polyhedra',
52 '_symbols',
53 '_dimension',
54 )
55
56 def __new__(cls, *polyhedra):
57 """
58 Create and return a new domain from a string or a list of polyhedra.
59 """
60 from .polyhedra import Polyhedron
61 if len(polyhedra) == 1:
62 argument = polyhedra[0]
63 if isinstance(argument, str):
64 return cls.fromstring(argument)
65 elif isinstance(argument, GeometricObject):
66 return argument.aspolyhedron()
67 else:
68 raise TypeError('argument must be a string '
69 'or a GeometricObject instance')
70 else:
71 for polyhedron in polyhedra:
72 if not isinstance(polyhedron, Polyhedron):
73 raise TypeError('arguments must be Polyhedron instances')
74 symbols = cls._xsymbols(polyhedra)
75 islset = cls._toislset(polyhedra, symbols)
76 return cls._fromislset(islset, symbols)
77
78 @classmethod
79 def _xsymbols(cls, iterator):
80 """
81 Return the ordered tuple of symbols present in iterator.
82 """
83 symbols = set()
84 for item in iterator:
85 symbols.update(item.symbols)
86 return tuple(sorted(symbols, key=Symbol.sortkey))
87
88 @property
89 def polyhedra(self):
90 """
91 The tuple of polyhedra which constitute the domain.
92 """
93 return self._polyhedra
94
95 @property
96 def symbols(self):
97 """
98 The tuple of symbols present in the domain equations.
99 """
100 return self._symbols
101
102 @property
103 def dimension(self):
104 """
105 The dimension of the domain, i.e. the number of symbols.
106 """
107 return self._dimension
108
109 def make_disjoint(self):
110 """
111 Return an equivalent domain, whose polyhedra are disjoint.
112 """
113 islset = self._toislset(self.polyhedra, self.symbols)
114 islset = libisl.isl_set_make_disjoint(mainctx, islset)
115 return self._fromislset(islset, self.symbols)
116
117 def isempty(self):
118 """
119 Return True if the domain is empty.
120 """
121 islset = self._toislset(self.polyhedra, self.symbols)
122 empty = bool(libisl.isl_set_is_empty(islset))
123 libisl.isl_set_free(islset)
124 return empty
125
126 def __bool__(self):
127 """
128 Return True if the domain is non-empty.
129 """
130 return not self.isempty()
131
132 def isuniverse(self):
133 """
134 Return True if the domain is universal, i.e. with no constraint.
135 """
136 islset = self._toislset(self.polyhedra, self.symbols)
137 universe = bool(libisl.isl_set_plain_is_universe(islset))
138 libisl.isl_set_free(islset)
139 return universe
140
141 def isbounded(self):
142 """
143 Return True if the domain is bounded.
144 """
145 islset = self._toislset(self.polyhedra, self.symbols)
146 bounded = bool(libisl.isl_set_is_bounded(islset))
147 libisl.isl_set_free(islset)
148 return bounded
149
150 def __eq__(self, other):
151 """
152 Return True if the two domains are equal.
153 """
154 symbols = self._xsymbols([self, other])
155 islset1 = self._toislset(self.polyhedra, symbols)
156 islset2 = other._toislset(other.polyhedra, symbols)
157 equal = bool(libisl.isl_set_is_equal(islset1, islset2))
158 libisl.isl_set_free(islset1)
159 libisl.isl_set_free(islset2)
160 return equal
161
162 def isdisjoint(self, other):
163 """
164 Return True if two domains have a null intersection.
165 """
166 symbols = self._xsymbols([self, other])
167 islset1 = self._toislset(self.polyhedra, symbols)
168 islset2 = self._toislset(other.polyhedra, symbols)
169 equal = bool(libisl.isl_set_is_disjoint(islset1, islset2))
170 libisl.isl_set_free(islset1)
171 libisl.isl_set_free(islset2)
172 return equal
173
174 def issubset(self, other):
175 """
176 Report whether another domain contains the domain.
177 """
178 symbols = self._xsymbols([self, other])
179 islset1 = self._toislset(self.polyhedra, symbols)
180 islset2 = self._toislset(other.polyhedra, symbols)
181 equal = bool(libisl.isl_set_is_subset(islset1, islset2))
182 libisl.isl_set_free(islset1)
183 libisl.isl_set_free(islset2)
184 return equal
185
186 def __le__(self, other):
187 return self.issubset(other)
188 __le__.__doc__ = issubset.__doc__
189
190 def __lt__(self, other):
191 """
192 Report whether another domain is contained within the domain.
193 """
194 symbols = self._xsymbols([self, other])
195 islset1 = self._toislset(self.polyhedra, symbols)
196 islset2 = self._toislset(other.polyhedra, symbols)
197 equal = bool(libisl.isl_set_is_strict_subset(islset1, islset2))
198 libisl.isl_set_free(islset1)
199 libisl.isl_set_free(islset2)
200 return equal
201
202 def complement(self):
203 """
204 Return the complementary domain of the domain.
205 """
206 islset = self._toislset(self.polyhedra, self.symbols)
207 islset = libisl.isl_set_complement(islset)
208 return self._fromislset(islset, self.symbols)
209
210 def __invert__(self):
211 return self.complement()
212 __invert__.__doc__ = complement.__doc__
213
214 def coalesce(self):
215 """
216 Simplify the representation of the domain by trying to combine pairs of
217 polyhedra into a single polyhedron.
218 """
219 islset = self._toislset(self.polyhedra, self.symbols)
220 islset = libisl.isl_set_coalesce(islset)
221 return self._fromislset(islset, self.symbols)
222
223 def detect_equalities(self):
224 """
225 Simplify the representation of the domain by detecting implicit
226 equalities.
227 """
228 islset = self._toislset(self.polyhedra, self.symbols)
229 islset = libisl.isl_set_detect_equalities(islset)
230 return self._fromislset(islset, self.symbols)
231
232 def remove_redundancies(self):
233 """
234 Remove redundant constraints in the domain.
235 """
236 islset = self._toislset(self.polyhedra, self.symbols)
237 islset = libisl.isl_set_remove_redundancies(islset)
238 return self._fromislset(islset, self.symbols)
239
240 def aspolyhedron(self):
241 """
242 Return the polyhedral hull of the domain.
243 """
244 from .polyhedra import Polyhedron
245 islset = self._toislset(self.polyhedra, self.symbols)
246 islbset = libisl.isl_set_polyhedral_hull(islset)
247 return Polyhedron._fromislbasicset(islbset, self.symbols)
248
249 def asdomain(self):
250 return self
251
252 def project(self, symbols):
253 """
254 Project out the symbols given in arguments.
255 """
256 islset = self._toislset(self.polyhedra, self.symbols)
257 n = 0
258 for index, symbol in reversed(list(enumerate(self.symbols))):
259 if symbol in symbols:
260 n += 1
261 elif n > 0:
262 islset = libisl.isl_set_project_out(islset,
263 libisl.isl_dim_set, index + 1, n)
264 n = 0
265 if n > 0:
266 islset = libisl.isl_set_project_out(islset, libisl.isl_dim_set, 0, n)
267 symbols = [symbol for symbol in self.symbols if symbol not in symbols]
268 return Domain._fromislset(islset, symbols)
269
270 def sample(self):
271 """
272 Return a sample of the domain.
273 """
274 islset = self._toislset(self.polyhedra, self.symbols)
275 islpoint = libisl.isl_set_sample_point(islset)
276 if bool(libisl.isl_point_is_void(islpoint)):
277 libisl.isl_point_free(islpoint)
278 raise ValueError('domain must be non-empty')
279 point = {}
280 for index, symbol in enumerate(self.symbols):
281 coordinate = libisl.isl_point_get_coordinate_val(islpoint,
282 libisl.isl_dim_set, index)
283 coordinate = islhelper.isl_val_to_int(coordinate)
284 point[symbol] = coordinate
285 libisl.isl_point_free(islpoint)
286 return point
287
288 def intersection(self, *others):
289 """
290 Return the intersection of two domains as a new domain.
291 """
292 if len(others) == 0:
293 return self
294 symbols = self._xsymbols((self,) + others)
295 islset1 = self._toislset(self.polyhedra, symbols)
296 for other in others:
297 islset2 = other._toislset(other.polyhedra, symbols)
298 islset1 = libisl.isl_set_intersect(islset1, islset2)
299 return self._fromislset(islset1, symbols)
300
301 def __and__(self, other):
302 return self.intersection(other)
303 __and__.__doc__ = intersection.__doc__
304
305 def union(self, *others):
306 """
307 Return the union of two domains as a new domain.
308 """
309 if len(others) == 0:
310 return self
311 symbols = self._xsymbols((self,) + others)
312 islset1 = self._toislset(self.polyhedra, symbols)
313 for other in others:
314 islset2 = other._toislset(other.polyhedra, symbols)
315 islset1 = libisl.isl_set_union(islset1, islset2)
316 return self._fromislset(islset1, symbols)
317
318 def __or__(self, other):
319 return self.union(other)
320 __or__.__doc__ = union.__doc__
321
322 def __add__(self, other):
323 return self.union(other)
324 __add__.__doc__ = union.__doc__
325
326 def difference(self, other):
327 """
328 Return the difference of two domains as a new domain.
329 """
330 symbols = self._xsymbols([self, other])
331 islset1 = self._toislset(self.polyhedra, symbols)
332 islset2 = other._toislset(other.polyhedra, symbols)
333 islset = libisl.isl_set_subtract(islset1, islset2)
334 return self._fromislset(islset, symbols)
335
336 def __sub__(self, other):
337 return self.difference(other)
338 __sub__.__doc__ = difference.__doc__
339
340 def lexmin(self):
341 """
342 Return the lexicographic minimum of the elements in the domain.
343 """
344 islset = self._toislset(self.polyhedra, self.symbols)
345 islset = libisl.isl_set_lexmin(islset)
346 return self._fromislset(islset, self.symbols)
347
348 def lexmax(self):
349 """
350 Return the lexicographic maximum of the elements in the domain.
351 """
352 islset = self._toislset(self.polyhedra, self.symbols)
353 islset = libisl.isl_set_lexmax(islset)
354 return self._fromislset(islset, self.symbols)
355
356 _RE_COORDINATE = re.compile(r'\((?P<num>\-?\d+)\)(/(?P<den>\d+))?')
357
358 def vertices(self):
359 """
360 Return the vertices of the domain.
361 """
362 from .polyhedra import Polyhedron
363 if not self.isbounded():
364 raise ValueError('domain must be bounded')
365 islbset = self._toislbasicset(self.equalities, self.inequalities,
366 self.symbols)
367 vertices = libisl.isl_basic_set_compute_vertices(islbset);
368 vertices = islhelper.isl_vertices_vertices(vertices)
369 points = []
370 for vertex in vertices:
371 expr = libisl.isl_vertex_get_expr(vertex)
372 coordinates = []
373 if islhelper.isl_version < '0.13':
374 constraints = islhelper.isl_basic_set_constraints(expr)
375 for constraint in constraints:
376 constant = libisl.isl_constraint_get_constant_val(constraint)
377 constant = islhelper.isl_val_to_int(constant)
378 for index, symbol in enumerate(self.symbols):
379 coefficient = libisl.isl_constraint_get_coefficient_val(constraint,
380 libisl.isl_dim_set, index)
381 coefficient = islhelper.isl_val_to_int(coefficient)
382 if coefficient != 0:
383 coordinate = -Fraction(constant, coefficient)
384 coordinates.append((symbol, coordinate))
385 else:
386 string = islhelper.isl_multi_aff_to_str(expr)
387 matches = self._RE_COORDINATE.finditer(string)
388 for symbol, match in zip(self.symbols, matches):
389 numerator = int(match.group('num'))
390 denominator = match.group('den')
391 denominator = 1 if denominator is None else int(denominator)
392 coordinate = Fraction(numerator, denominator)
393 coordinates.append((symbol, coordinate))
394 points.append(Point(coordinates))
395 return points
396
397 def points(self):
398 """
399 Return the points with integer coordinates contained in the domain.
400 """
401 if not self.isbounded():
402 raise ValueError('domain must be bounded')
403 from .polyhedra import Universe, Eq
404 islset = self._toislset(self.polyhedra, self.symbols)
405 islpoints = islhelper.isl_set_points(islset)
406 points = []
407 for islpoint in islpoints:
408 coordinates = {}
409 for index, symbol in enumerate(self.symbols):
410 coordinate = libisl.isl_point_get_coordinate_val(islpoint,
411 libisl.isl_dim_set, index)
412 coordinate = islhelper.isl_val_to_int(coordinate)
413 coordinates[symbol] = coordinate
414 points.append(Point(coordinates))
415 return points
416
417 @classmethod
418 def _polygon_inner_point(cls, points):
419 symbols = points[0].symbols
420 coordinates = {symbol: 0 for symbol in symbols}
421 for point in points:
422 for symbol, coordinate in point.coordinates():
423 coordinates[symbol] += coordinate
424 for symbol in symbols:
425 coordinates[symbol] /= len(points)
426 return Point(coordinates)
427
428 @classmethod
429 def _sort_polygon_2d(cls, points):
430 if len(points) <= 3:
431 return points
432 o = cls._polygon_inner_point(points)
433 angles = {}
434 for m in points:
435 om = Vector(o, m)
436 dx, dy = (coordinate for symbol, coordinate in om.coordinates())
437 angle = math.atan2(dy, dx)
438 angles[m] = angle
439 return sorted(points, key=angles.get)
440
441 @classmethod
442 def _sort_polygon_3d(cls, points):
443 if len(points) <= 3:
444 return points
445 o = cls._polygon_inner_point(points)
446 a = points[0]
447 oa = Vector(o, a)
448 norm_oa = oa.norm()
449 for b in points[1:]:
450 ob = Vector(o, b)
451 u = oa.cross(ob)
452 if not u.isnull():
453 u = u.asunit()
454 break
455 else:
456 raise ValueError('degenerate polygon')
457 angles = {a: 0.}
458 for m in points[1:]:
459 om = Vector(o, m)
460 normprod = norm_oa * om.norm()
461 cosinus = max(oa.dot(om) / normprod, -1.)
462 sinus = u.dot(oa.cross(om)) / normprod
463 angle = math.acos(cosinus)
464 angle = math.copysign(angle, sinus)
465 angles[m] = angle
466 return sorted(points, key=angles.get)
467
468 def faces(self):
469 """
470 Return the vertices of the domain, grouped by face.
471 """
472 faces = []
473 for polyhedron in self.polyhedra:
474 vertices = polyhedron.vertices()
475 for constraint in polyhedron.constraints:
476 face = []
477 for vertex in vertices:
478 if constraint.subs(vertex.coordinates()) == 0:
479 face.append(vertex)
480 if len(face) >= 3:
481 faces.append(face)
482 return faces
483
484 def _plot_2d(self, plot=None, **kwargs):
485 import matplotlib.pyplot as plt
486 from matplotlib.patches import Polygon
487 if plot is None:
488 fig = plt.figure()
489 plot = fig.add_subplot(1, 1, 1)
490 xmin, xmax = plot.get_xlim()
491 ymin, ymax = plot.get_ylim()
492 for polyhedron in self.polyhedra:
493 vertices = polyhedron._sort_polygon_2d(polyhedron.vertices())
494 xys = [tuple(vertex.values()) for vertex in vertices]
495 xs, ys = zip(*xys)
496 xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
497 ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
498 plot.add_patch(Polygon(xys, closed=True, **kwargs))
499 plot.set_xlim(xmin, xmax)
500 plot.set_ylim(ymin, ymax)
501 return plot
502
503 def _plot_3d(self, plot=None, **kwargs):
504 import matplotlib.pyplot as plt
505 from mpl_toolkits.mplot3d import Axes3D
506 from mpl_toolkits.mplot3d.art3d import Poly3DCollection
507 if plot is None:
508 fig = plt.figure()
509 axes = Axes3D(fig)
510 else:
511 axes = plot
512 xmin, xmax = axes.get_xlim()
513 ymin, ymax = axes.get_ylim()
514 zmin, zmax = axes.get_zlim()
515 poly_xyzs = []
516 for vertices in self.faces():
517 vertices = self._sort_polygon_3d(vertices)
518 vertices.append(vertices[0])
519 face_xyzs = [tuple(vertex.values()) for vertex in vertices]
520 xs, ys, zs = zip(*face_xyzs)
521 xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
522 ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
523 zmin, zmax = min(zmin, float(min(zs))), max(zmax, float(max(zs)))
524 poly_xyzs.append(face_xyzs)
525 collection = Poly3DCollection(poly_xyzs, **kwargs)
526 axes.add_collection3d(collection)
527 axes.set_xlim(xmin, xmax)
528 axes.set_ylim(ymin, ymax)
529 axes.set_zlim(zmin, zmax)
530 return axes
531
532 def plot(self, plot=None, **kwargs):
533 """
534 Plot the domain using matplotlib.
535 """
536 if not self.isbounded():
537 raise ValueError('domain must be bounded')
538 elif self.dimension == 2:
539 return self._plot_2d(plot=plot, **kwargs)
540 elif self.dimension == 3:
541 return self._plot_3d(plot=plot, **kwargs)
542 else:
543 raise ValueError('polyhedron must be 2 or 3-dimensional')
544
545 def __contains__(self, point):
546 """
547 Return True if point if contained within the domain.
548 """
549 for polyhedron in self.polyhedra:
550 if point in polyhedron:
551 return True
552 return False
553
554 def subs(self, symbol, expression=None):
555 """
556 Subsitute symbol by expression in equations and return the resulting
557 domain.
558 """
559 polyhedra = [polyhedron.subs(symbol, expression)
560 for polyhedron in self.polyhedra]
561 return Domain(*polyhedra)
562
563 @classmethod
564 def _fromislset(cls, islset, symbols):
565 from .polyhedra import Polyhedron
566 islset = libisl.isl_set_remove_divs(islset)
567 islbsets = islhelper.isl_set_basic_sets(islset)
568 libisl.isl_set_free(islset)
569 polyhedra = []
570 for islbset in islbsets:
571 polyhedron = Polyhedron._fromislbasicset(islbset, symbols)
572 polyhedra.append(polyhedron)
573 if len(polyhedra) == 0:
574 from .polyhedra import Empty
575 return Empty
576 elif len(polyhedra) == 1:
577 return polyhedra[0]
578 else:
579 self = object().__new__(Domain)
580 self._polyhedra = tuple(polyhedra)
581 self._symbols = cls._xsymbols(polyhedra)
582 self._dimension = len(self._symbols)
583 return self
584
585 @classmethod
586 def _toislset(cls, polyhedra, symbols):
587 polyhedron = polyhedra[0]
588 islbset = polyhedron._toislbasicset(polyhedron.equalities,
589 polyhedron.inequalities, symbols)
590 islset1 = libisl.isl_set_from_basic_set(islbset)
591 for polyhedron in polyhedra[1:]:
592 islbset = polyhedron._toislbasicset(polyhedron.equalities,
593 polyhedron.inequalities, symbols)
594 islset2 = libisl.isl_set_from_basic_set(islbset)
595 islset1 = libisl.isl_set_union(islset1, islset2)
596 return islset1
597
598 @classmethod
599 def _fromast(cls, node):
600 from .polyhedra import Polyhedron
601 if isinstance(node, ast.Module) and len(node.body) == 1:
602 return cls._fromast(node.body[0])
603 elif isinstance(node, ast.Expr):
604 return cls._fromast(node.value)
605 elif isinstance(node, ast.UnaryOp):
606 domain = cls._fromast(node.operand)
607 if isinstance(node.operand, ast.invert):
608 return Not(domain)
609 elif isinstance(node, ast.BinOp):
610 domain1 = cls._fromast(node.left)
611 domain2 = cls._fromast(node.right)
612 if isinstance(node.op, ast.BitAnd):
613 return And(domain1, domain2)
614 elif isinstance(node.op, ast.BitOr):
615 return Or(domain1, domain2)
616 elif isinstance(node, ast.Compare):
617 equalities = []
618 inequalities = []
619 left = LinExpr._fromast(node.left)
620 for i in range(len(node.ops)):
621 op = node.ops[i]
622 right = LinExpr._fromast(node.comparators[i])
623 if isinstance(op, ast.Lt):
624 inequalities.append(right - left - 1)
625 elif isinstance(op, ast.LtE):
626 inequalities.append(right - left)
627 elif isinstance(op, ast.Eq):
628 equalities.append(left - right)
629 elif isinstance(op, ast.GtE):
630 inequalities.append(left - right)
631 elif isinstance(op, ast.Gt):
632 inequalities.append(left - right - 1)
633 else:
634 break
635 left = right
636 else:
637 return Polyhedron(equalities, inequalities)
638 raise SyntaxError('invalid syntax')
639
640 _RE_BRACES = re.compile(r'^\{\s*|\s*\}$')
641 _RE_EQ = re.compile(r'([^<=>])=([^<=>])')
642 _RE_AND = re.compile(r'\band\b|,|&&|/\\|∧|∩')
643 _RE_OR = re.compile(r'\bor\b|;|\|\||\\/|∨|∪')
644 _RE_NOT = re.compile(r'\bnot\b|!|¬')
645 _RE_NUM_VAR = LinExpr._RE_NUM_VAR
646 _RE_OPERATORS = re.compile(r'(&|\||~)')
647
648 @classmethod
649 def fromstring(cls, string):
650 """
651 Convert a string into a domain.
652 """
653 # remove curly brackets
654 string = cls._RE_BRACES.sub(r'', string)
655 # replace '=' by '=='
656 string = cls._RE_EQ.sub(r'\1==\2', string)
657 # replace 'and', 'or', 'not'
658 string = cls._RE_AND.sub(r' & ', string)
659 string = cls._RE_OR.sub(r' | ', string)
660 string = cls._RE_NOT.sub(r' ~', string)
661 # add implicit multiplication operators, e.g. '5x' -> '5*x'
662 string = cls._RE_NUM_VAR.sub(r'\1*\2', string)
663 # add parentheses to force precedence
664 tokens = cls._RE_OPERATORS.split(string)
665 for i, token in enumerate(tokens):
666 if i % 2 == 0:
667 token = '({})'.format(token)
668 tokens[i] = token
669 string = ''.join(tokens)
670 tree = ast.parse(string, 'eval')
671 return cls._fromast(tree)
672
673 def __repr__(self):
674 """
675 Return repr(self).
676 """
677 assert len(self.polyhedra) >= 2
678 strings = [repr(polyhedron) for polyhedron in self.polyhedra]
679 return 'Or({})'.format(', '.join(strings))
680
681 def _repr_latex_(self):
682 strings = []
683 for polyhedron in self.polyhedra:
684 strings.append('({})'.format(polyhedron._repr_latex_().strip('$')))
685 return '${}$'.format(' \\vee '.join(strings))
686
687 @classmethod
688 def fromsympy(cls, expr):
689 """
690 Convert a SymPy expression into a domain.
691 """
692 import sympy
693 from .polyhedra import Lt, Le, Eq, Ne, Ge, Gt
694 funcmap = {
695 sympy.And: And, sympy.Or: Or, sympy.Not: Not,
696 sympy.Lt: Lt, sympy.Le: Le,
697 sympy.Eq: Eq, sympy.Ne: Ne,
698 sympy.Ge: Ge, sympy.Gt: Gt,
699 }
700 if expr.func in funcmap:
701 args = [Domain.fromsympy(arg) for arg in expr.args]
702 return funcmap[expr.func](*args)
703 elif isinstance(expr, sympy.Expr):
704 return LinExpr.fromsympy(expr)
705 raise ValueError('non-domain expression: {!r}'.format(expr))
706
707 def tosympy(self):
708 """
709 Convert a domain into a SymPy expression.
710 """
711 import sympy
712 polyhedra = [polyhedron.tosympy() for polyhedron in polyhedra]
713 return sympy.Or(*polyhedra)
714
715
716 def And(*domains):
717 if len(domains) == 0:
718 from .polyhedra import Universe
719 return Universe
720 else:
721 return domains[0].intersection(*domains[1:])
722 And.__doc__ = Domain.intersection.__doc__
723
724 def Or(*domains):
725 if len(domains) == 0:
726 from .polyhedra import Empty
727 return Empty
728 else:
729 return domains[0].union(*domains[1:])
730 Or.__doc__ = Domain.union.__doc__
731
732 def Not(domain):
733 return ~domain
734 Not.__doc__ = Domain.complement.__doc__