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