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