Fix error message in Domain.plot()
[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 return self - other
367
368 def __sub__(self, other):
369 if isinstance(other, Domain):
370 symbols = self._xsymbols([self, other])
371 islset1 = self._toislset(self.polyhedra, symbols)
372 islset2 = other._toislset(other.polyhedra, symbols)
373 islset = libisl.isl_set_subtract(islset1, islset2)
374 return self._fromislset(islset, symbols)
375 return NotImplemented
376 __sub__.__doc__ = difference.__doc__
377
378 def lexmin(self):
379 """
380 Return the lexicographic minimum of the elements in the domain.
381 """
382 islset = self._toislset(self.polyhedra, self.symbols)
383 islset = libisl.isl_set_lexmin(islset)
384 return self._fromislset(islset, self.symbols)
385
386 def lexmax(self):
387 """
388 Return the lexicographic maximum of the elements in the domain.
389 """
390 islset = self._toislset(self.polyhedra, self.symbols)
391 islset = libisl.isl_set_lexmax(islset)
392 return self._fromislset(islset, self.symbols)
393
394 if islhelper.isl_version >= '0.13':
395 _RE_COORDINATE = re.compile(r'\((?P<num>\-?\d+)\)(/(?P<den>\d+))?')
396 else:
397 _RE_COORDINATE = None
398
399 def vertices(self):
400 """
401 Return the vertices of the domain, as a list of rational instances of
402 Point.
403 """
404 from .polyhedra import Polyhedron
405 if not self.isbounded():
406 raise ValueError('domain must be bounded')
407 islbset = self._toislbasicset(self.equalities, self.inequalities,
408 self.symbols)
409 vertices = libisl.isl_basic_set_compute_vertices(islbset);
410 vertices = islhelper.isl_vertices_vertices(vertices)
411 points = []
412 for vertex in vertices:
413 expr = libisl.isl_vertex_get_expr(vertex)
414 coordinates = []
415 if self._RE_COORDINATE is None:
416 constraints = islhelper.isl_basic_set_constraints(expr)
417 for constraint in constraints:
418 constant = libisl.isl_constraint_get_constant_val(constraint)
419 constant = islhelper.isl_val_to_int(constant)
420 for index, symbol in enumerate(self.symbols):
421 coefficient = libisl.isl_constraint_get_coefficient_val(constraint,
422 libisl.isl_dim_set, index)
423 coefficient = islhelper.isl_val_to_int(coefficient)
424 if coefficient != 0:
425 coordinate = -Fraction(constant, coefficient)
426 coordinates.append((symbol, coordinate))
427 else:
428 string = islhelper.isl_multi_aff_to_str(expr)
429 matches = self._RE_COORDINATE.finditer(string)
430 for symbol, match in zip(self.symbols, matches):
431 numerator = int(match.group('num'))
432 denominator = match.group('den')
433 denominator = 1 if denominator is None else int(denominator)
434 coordinate = Fraction(numerator, denominator)
435 coordinates.append((symbol, coordinate))
436 points.append(Point(coordinates))
437 return points
438
439 def points(self):
440 """
441 Return the integer points of a bounded domain, as a list of integer
442 instances of Point. If the domain is not bounded, a ValueError exception
443 is raised.
444 """
445 if not self.isbounded():
446 raise ValueError('domain must be bounded')
447 from .polyhedra import Universe, Eq
448 islset = self._toislset(self.polyhedra, self.symbols)
449 islpoints = islhelper.isl_set_points(islset)
450 points = []
451 for islpoint in islpoints:
452 coordinates = {}
453 for index, symbol in enumerate(self.symbols):
454 coordinate = libisl.isl_point_get_coordinate_val(islpoint,
455 libisl.isl_dim_set, index)
456 coordinate = islhelper.isl_val_to_int(coordinate)
457 coordinates[symbol] = coordinate
458 points.append(Point(coordinates))
459 return points
460
461 def __contains__(self, point):
462 """
463 Return True if the point if contained within the domain.
464 """
465 for polyhedron in self.polyhedra:
466 if point in polyhedron:
467 return True
468 return False
469
470 @classmethod
471 def _polygon_inner_point(cls, points):
472 symbols = points[0].symbols
473 coordinates = {symbol: 0 for symbol in symbols}
474 for point in points:
475 for symbol, coordinate in point.coordinates():
476 coordinates[symbol] += coordinate
477 for symbol in symbols:
478 coordinates[symbol] /= len(points)
479 return Point(coordinates)
480
481 @classmethod
482 def _sort_polygon_2d(cls, points):
483 if len(points) <= 3:
484 return points
485 o = cls._polygon_inner_point(points)
486 angles = {}
487 for m in points:
488 om = Vector(o, m)
489 dx, dy = (coordinate for symbol, coordinate in om.coordinates())
490 angle = math.atan2(dy, dx)
491 angles[m] = angle
492 return sorted(points, key=angles.get)
493
494 @classmethod
495 def _sort_polygon_3d(cls, points):
496 if len(points) <= 3:
497 return points
498 o = cls._polygon_inner_point(points)
499 a = points[0]
500 oa = Vector(o, a)
501 norm_oa = oa.norm()
502 for b in points[1:]:
503 ob = Vector(o, b)
504 u = oa.cross(ob)
505 if not u.isnull():
506 u = u.asunit()
507 break
508 else:
509 raise ValueError('degenerate polygon')
510 angles = {a: 0.}
511 for m in points[1:]:
512 om = Vector(o, m)
513 normprod = norm_oa * om.norm()
514 cosinus = max(oa.dot(om) / normprod, -1.)
515 sinus = u.dot(oa.cross(om)) / normprod
516 angle = math.acos(cosinus)
517 angle = math.copysign(angle, sinus)
518 angles[m] = angle
519 return sorted(points, key=angles.get)
520
521 def faces(self):
522 """
523 Return the list of faces of a bounded domain. Each face is represented
524 by a list of vertices, in the form of rational instances of Point. If
525 the domain is not bounded, a ValueError exception is raised.
526 """
527 faces = []
528 for polyhedron in self.polyhedra:
529 vertices = polyhedron.vertices()
530 for constraint in polyhedron.constraints:
531 face = []
532 for vertex in vertices:
533 if constraint.subs(vertex.coordinates()) == 0:
534 face.append(vertex)
535 if len(face) >= 3:
536 faces.append(face)
537 return faces
538
539 def _plot_2d(self, plot=None, **kwargs):
540 import matplotlib.pyplot as plt
541 from matplotlib.patches import Polygon
542 if plot is None:
543 fig = plt.figure()
544 plot = fig.add_subplot(1, 1, 1)
545 xmin, xmax = plot.get_xlim()
546 ymin, ymax = plot.get_ylim()
547 for polyhedron in self.polyhedra:
548 vertices = polyhedron._sort_polygon_2d(polyhedron.vertices())
549 xys = [tuple(vertex.values()) for vertex in vertices]
550 xs, ys = zip(*xys)
551 xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
552 ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
553 plot.add_patch(Polygon(xys, closed=True, **kwargs))
554 plot.set_xlim(xmin, xmax)
555 plot.set_ylim(ymin, ymax)
556 return plot
557
558 def _plot_3d(self, plot=None, **kwargs):
559 import matplotlib.pyplot as plt
560 from mpl_toolkits.mplot3d import Axes3D
561 from mpl_toolkits.mplot3d.art3d import Poly3DCollection
562 if plot is None:
563 fig = plt.figure()
564 axes = Axes3D(fig)
565 else:
566 axes = plot
567 xmin, xmax = axes.get_xlim()
568 ymin, ymax = axes.get_ylim()
569 zmin, zmax = axes.get_zlim()
570 poly_xyzs = []
571 for vertices in self.faces():
572 vertices = self._sort_polygon_3d(vertices)
573 vertices.append(vertices[0])
574 face_xyzs = [tuple(vertex.values()) for vertex in vertices]
575 xs, ys, zs = zip(*face_xyzs)
576 xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
577 ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
578 zmin, zmax = min(zmin, float(min(zs))), max(zmax, float(max(zs)))
579 poly_xyzs.append(face_xyzs)
580 collection = Poly3DCollection(poly_xyzs, **kwargs)
581 axes.add_collection3d(collection)
582 axes.set_xlim(xmin, xmax)
583 axes.set_ylim(ymin, ymax)
584 axes.set_zlim(zmin, zmax)
585 return axes
586
587 def plot(self, plot=None, **kwargs):
588 """
589 Plot a 2D or 3D domain using matplotlib. Draw it to the current plot
590 object if present, otherwise create a new one. options are keyword
591 arguments passed to the matplotlib drawing functions, they can be used
592 to set the drawing color for example. Raise ValueError is the domain is
593 not 2D or 3D.
594 """
595 if not self.isbounded():
596 raise ValueError('domain must be bounded')
597 elif self.dimension == 2:
598 return self._plot_2d(plot=plot, **kwargs)
599 elif self.dimension == 3:
600 return self._plot_3d(plot=plot, **kwargs)
601 else:
602 raise ValueError('domain must be 2 or 3-dimensional')
603
604 def subs(self, symbol, expression=None):
605 """
606 Substitute the given symbol by an expression in the domain constraints.
607 To perform multiple substitutions at once, pass a sequence or a
608 dictionary of (old, new) pairs to subs. The syntax of this function is
609 similar to LinExpr.subs().
610 """
611 polyhedra = [polyhedron.subs(symbol, expression)
612 for polyhedron in self.polyhedra]
613 return Domain(*polyhedra)
614
615 @classmethod
616 def _fromislset(cls, islset, symbols):
617 from .polyhedra import Polyhedron
618 islset = libisl.isl_set_remove_divs(islset)
619 islbsets = islhelper.isl_set_basic_sets(islset)
620 libisl.isl_set_free(islset)
621 polyhedra = []
622 for islbset in islbsets:
623 polyhedron = Polyhedron._fromislbasicset(islbset, symbols)
624 polyhedra.append(polyhedron)
625 if len(polyhedra) == 0:
626 from .polyhedra import Empty
627 return Empty
628 elif len(polyhedra) == 1:
629 return polyhedra[0]
630 else:
631 self = object().__new__(Domain)
632 self._polyhedra = tuple(polyhedra)
633 self._symbols = cls._xsymbols(polyhedra)
634 self._dimension = len(self._symbols)
635 return self
636
637 @classmethod
638 def _toislset(cls, polyhedra, symbols):
639 polyhedron = polyhedra[0]
640 islbset = polyhedron._toislbasicset(polyhedron.equalities,
641 polyhedron.inequalities, symbols)
642 islset1 = libisl.isl_set_from_basic_set(islbset)
643 for polyhedron in polyhedra[1:]:
644 islbset = polyhedron._toislbasicset(polyhedron.equalities,
645 polyhedron.inequalities, symbols)
646 islset2 = libisl.isl_set_from_basic_set(islbset)
647 islset1 = libisl.isl_set_union(islset1, islset2)
648 return islset1
649
650 @classmethod
651 def _fromast(cls, node):
652 from .polyhedra import Polyhedron
653 if isinstance(node, ast.Module) and len(node.body) == 1:
654 return cls._fromast(node.body[0])
655 elif isinstance(node, ast.Expr):
656 return cls._fromast(node.value)
657 elif isinstance(node, ast.UnaryOp):
658 domain = cls._fromast(node.operand)
659 if isinstance(node.operand, ast.invert):
660 return Not(domain)
661 elif isinstance(node, ast.BinOp):
662 domain1 = cls._fromast(node.left)
663 domain2 = cls._fromast(node.right)
664 if isinstance(node.op, ast.BitAnd):
665 return And(domain1, domain2)
666 elif isinstance(node.op, ast.BitOr):
667 return Or(domain1, domain2)
668 elif isinstance(node, ast.Compare):
669 equalities = []
670 inequalities = []
671 left = LinExpr._fromast(node.left)
672 for i in range(len(node.ops)):
673 op = node.ops[i]
674 right = LinExpr._fromast(node.comparators[i])
675 if isinstance(op, ast.Lt):
676 inequalities.append(right - left - 1)
677 elif isinstance(op, ast.LtE):
678 inequalities.append(right - left)
679 elif isinstance(op, ast.Eq):
680 equalities.append(left - right)
681 elif isinstance(op, ast.GtE):
682 inequalities.append(left - right)
683 elif isinstance(op, ast.Gt):
684 inequalities.append(left - right - 1)
685 else:
686 break
687 left = right
688 else:
689 return Polyhedron(equalities, inequalities)
690 raise SyntaxError('invalid syntax')
691
692 _RE_BRACES = re.compile(r'^\{\s*|\s*\}$')
693 _RE_EQ = re.compile(r'([^<=>])=([^<=>])')
694 _RE_AND = re.compile(r'\band\b|,|&&|/\\|∧|∩')
695 _RE_OR = re.compile(r'\bor\b|;|\|\||\\/|∨|∪')
696 _RE_NOT = re.compile(r'\bnot\b|!|¬')
697 _RE_NUM_VAR = LinExpr._RE_NUM_VAR
698 _RE_OPERATORS = re.compile(r'(&|\||~)')
699
700 @classmethod
701 def fromstring(cls, string):
702 """
703 Create a domain from a string. Raise SyntaxError if the string is not
704 properly formatted.
705 """
706 # remove curly brackets
707 string = cls._RE_BRACES.sub(r'', string)
708 # replace '=' by '=='
709 string = cls._RE_EQ.sub(r'\1==\2', string)
710 # replace 'and', 'or', 'not'
711 string = cls._RE_AND.sub(r' & ', string)
712 string = cls._RE_OR.sub(r' | ', string)
713 string = cls._RE_NOT.sub(r' ~', string)
714 # add implicit multiplication operators, e.g. '5x' -> '5*x'
715 string = cls._RE_NUM_VAR.sub(r'\1*\2', string)
716 # add parentheses to force precedence
717 tokens = cls._RE_OPERATORS.split(string)
718 for i, token in enumerate(tokens):
719 if i % 2 == 0:
720 token = '({})'.format(token)
721 tokens[i] = token
722 string = ''.join(tokens)
723 tree = ast.parse(string, 'eval')
724 return cls._fromast(tree)
725
726 def __repr__(self):
727 assert len(self.polyhedra) >= 2
728 strings = [repr(polyhedron) for polyhedron in self.polyhedra]
729 return 'Or({})'.format(', '.join(strings))
730
731 def _repr_latex_(self):
732 strings = []
733 for polyhedron in self.polyhedra:
734 strings.append('({})'.format(polyhedron._repr_latex_().strip('$')))
735 return '${}$'.format(' \\vee '.join(strings))
736
737 @classmethod
738 def fromsympy(cls, expr):
739 """
740 Create a domain from a sympy expression.
741 """
742 import sympy
743 from .polyhedra import Lt, Le, Eq, Ne, Ge, Gt
744 funcmap = {
745 sympy.And: And, sympy.Or: Or, sympy.Not: Not,
746 sympy.Lt: Lt, sympy.Le: Le,
747 sympy.Eq: Eq, sympy.Ne: Ne,
748 sympy.Ge: Ge, sympy.Gt: Gt,
749 }
750 if expr.func in funcmap:
751 args = [Domain.fromsympy(arg) for arg in expr.args]
752 return funcmap[expr.func](*args)
753 elif isinstance(expr, sympy.Expr):
754 return LinExpr.fromsympy(expr)
755 raise ValueError('non-domain expression: {!r}'.format(expr))
756
757 def tosympy(self):
758 """
759 Convert the domain to a sympy expression.
760 """
761 import sympy
762 polyhedra = [polyhedron.tosympy() for polyhedron in polyhedra]
763 return sympy.Or(*polyhedra)
764
765
766 def And(*domains):
767 """
768 Create the intersection domain of the domains given in arguments.
769 """
770 if len(domains) == 0:
771 from .polyhedra import Universe
772 return Universe
773 else:
774 return domains[0].intersection(*domains[1:])
775 And.__doc__ = Domain.intersection.__doc__
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 Or.__doc__ = Domain.union.__doc__
787
788 def Not(domain):
789 """
790 Create the complementary domain of the domain given in argument.
791 """
792 return ~domain
793 Not.__doc__ = Domain.complement.__doc__