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