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