122d428f7e63dd1056a357aaae8e9d85e97939f4
6 from fractions
import Fraction
8 from . import islhelper
9 from .islhelper
import mainctx
, libisl
10 from .linexprs
import Expression
, Symbol
, Rational
11 from .geometry
import GeometricObject
, Point
, Vector
20 @functools.total_ordering
21 class Domain(GeometricObject
):
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()
38 raise TypeError('argument must be a string '
39 'or a GeometricObject instance')
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
)
49 def _xsymbols(cls
, iterator
):
51 Return the ordered tuple of symbols present in iterator.
55 symbols
.update(item
.symbols
)
56 return tuple(sorted(symbols
, key
=Symbol
.sortkey
))
60 return self
._polyhedra
68 return self
._dimension
72 Returns this set as disjoint.
74 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
75 islset
= libisl
.isl_set_make_disjoint(mainctx
, islset
)
76 return self
._fromislset
(islset
, self
.symbols
)
80 Returns true if this set is an Empty set.
82 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
83 empty
= bool(libisl
.isl_set_is_empty(islset
))
84 libisl
.isl_set_free(islset
)
88 return not self
.isempty()
92 Returns true if this set is the Universe set.
94 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
95 universe
= bool(libisl
.isl_set_plain_is_universe(islset
))
96 libisl
.isl_set_free(islset
)
101 Returns true if this set is bounded.
103 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
104 bounded
= bool(libisl
.isl_set_is_bounded(islset
))
105 libisl
.isl_set_free(islset
)
108 def __eq__(self
, other
):
110 Returns true if two sets are equal.
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
)
120 def isdisjoint(self
, other
):
122 Return True if two sets have a null intersection.
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
)
132 def issubset(self
, other
):
134 Report whether another set contains this set.
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
)
144 def __le__(self
, other
):
146 Returns true if this set is less than or equal to another set.
148 return self
.issubset(other
)
150 def __lt__(self
, other
):
152 Returns true if this set is less than another set.
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
)
162 def complement(self
):
164 Returns the complement of this set.
166 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
167 islset
= libisl
.isl_set_complement(islset
)
168 return self
._fromislset
(islset
, self
.symbols
)
170 def __invert__(self
):
172 Returns the complement of this set.
174 return self
.complement()
178 Returns a set without redundant constraints.
180 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
181 islset
= libisl
.isl_set_remove_redundancies(islset
)
182 return self
._fromislset
(islset
, self
.symbols
)
184 def aspolyhedron(self
):
186 Returns polyhedral hull of set.
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
)
196 def project(self
, dims
):
198 Return new set with given dimensions removed.
200 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
202 for index
, symbol
in reversed(list(enumerate(self
.symbols
))):
206 islset
= libisl
.isl_set_project_out(islset
, libisl
.isl_dim_set
, index
+ 1, n
)
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
)
215 Returns a single subset of the input.
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')
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
)
231 def intersection(self
, *others
):
233 Return the intersection of two sets as a new set.
237 symbols
= self
._xsymbols
((self
,) + others
)
238 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
240 islset2
= other
._toislset
(other
.polyhedra
, symbols
)
241 islset1
= libisl
.isl_set_intersect(islset1
, islset2
)
242 return self
._fromislset
(islset1
, symbols
)
244 def __and__(self
, other
):
246 Return the intersection of two sets as a new set.
248 return self
.intersection(other
)
250 def union(self
, *others
):
252 Return the union of sets as a new set.
256 symbols
= self
._xsymbols
((self
,) + others
)
257 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
259 islset2
= other
._toislset
(other
.polyhedra
, symbols
)
260 islset1
= libisl
.isl_set_union(islset1
, islset2
)
261 return self
._fromislset
(islset1
, symbols
)
263 def __or__(self
, other
):
265 Return a new set with elements from both sets.
267 return self
.union(other
)
269 def __add__(self
, other
):
271 Return new set containing all elements in both sets.
273 return self
.union(other
)
275 def difference(self
, other
):
277 Return the difference of two sets as a new set.
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
)
285 def __sub__(self
, other
):
287 Return the difference of two sets as a new set.
289 return self
.difference(other
)
293 Return a new set containing the lexicographic minimum of the elements in the set.
295 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
296 islset
= libisl
.isl_set_lexmin(islset
)
297 return self
._fromislset
(islset
, self
.symbols
)
301 Return a new set containing the lexicographic maximum of the elements in the set.
303 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
304 islset
= libisl
.isl_set_lexmax(islset
)
305 return self
._fromislset
(islset
, self
.symbols
)
307 def num_parameters(self
):
309 Return the total number of parameters, input, output or set dimensions.
311 islbset
= self
._toislbasicset
(self
.equalities
, self
.inequalities
, self
.symbols
)
312 num
= libisl
.isl_basic_set_dim(islbset
, libisl
.isl_dim_set
)
315 def involves_dims(self
, dims
):
317 Returns true if set depends on given dimensions.
319 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
321 symbols
= sorted(list(self
.symbols
))
326 first
= symbols
.index(dims
[0])
332 value
= bool(libisl
.isl_set_involves_dims(islset
, libisl
.isl_dim_set
, first
, n
))
333 libisl
.isl_set_free(islset
)
336 _RE_COORDINATE
= re
.compile(r
'\((?P<num>\-?\d+)\)(/(?P<den>\d+))?')
340 Return a list of vertices for this Polygon.
342 from .polyhedra
import Polyhedron
343 if not self
.isbounded():
344 raise ValueError('domain must be bounded')
345 islbset
= self
._toislbasicset
(self
.equalities
, self
.inequalities
, self
.symbols
)
346 vertices
= libisl
.isl_basic_set_compute_vertices(islbset
);
347 vertices
= islhelper
.isl_vertices_vertices(vertices
)
349 for vertex
in vertices
:
350 expr
= libisl
.isl_vertex_get_expr(vertex
)
352 if islhelper
.isl_version
< '0.13':
353 constraints
= islhelper
.isl_basic_set_constraints(expr
)
354 for constraint
in constraints
:
355 constant
= libisl
.isl_constraint_get_constant_val(constraint
)
356 constant
= islhelper
.isl_val_to_int(constant
)
357 for index
, symbol
in enumerate(self
.symbols
):
358 coefficient
= libisl
.isl_constraint_get_coefficient_val(constraint
,
359 libisl
.isl_dim_set
, index
)
360 coefficient
= islhelper
.isl_val_to_int(coefficient
)
362 coordinate
= -Fraction(constant
, coefficient
)
363 coordinates
.append((symbol
, coordinate
))
365 string
= islhelper
.isl_multi_aff_to_str(expr
)
366 matches
= self
._RE
_COORDINATE
.finditer(string
)
367 for symbol
, match
in zip(self
.symbols
, matches
):
368 numerator
= int(match
.group('num'))
369 denominator
= match
.group('den')
370 denominator
= 1 if denominator
is None else int(denominator
)
371 coordinate
= Fraction(numerator
, denominator
)
372 coordinates
.append((symbol
, coordinate
))
373 points
.append(Point(coordinates
))
378 Returns the points contained in the set.
380 if not self
.isbounded():
381 raise ValueError('domain must be bounded')
382 from .polyhedra
import Universe
, Eq
383 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
384 islpoints
= islhelper
.isl_set_points(islset
)
386 for islpoint
in islpoints
:
388 for index
, symbol
in enumerate(self
.symbols
):
389 coordinate
= libisl
.isl_point_get_coordinate_val(islpoint
,
390 libisl
.isl_dim_set
, index
)
391 coordinate
= islhelper
.isl_val_to_int(coordinate
)
392 coordinates
[symbol
] = coordinate
393 points
.append(Point(coordinates
))
397 def _polygon_inner_point(cls
, points
):
398 symbols
= points
[0].symbols
399 coordinates
= {symbol
: 0 for symbol
in symbols
}
401 for symbol
, coordinate
in point
.coordinates():
402 coordinates
[symbol
] += coordinate
403 for symbol
in symbols
:
404 coordinates
[symbol
] /= len(points
)
405 return Point(coordinates
)
408 def _sort_polygon_2d(cls
, points
):
411 o
= cls
._polygon
_inner
_point
(points
)
415 dx
, dy
= (coordinate
for symbol
, coordinate
in om
.coordinates())
416 angle
= math
.atan2(dy
, dx
)
418 return sorted(points
, key
=angles
.get
)
421 def _sort_polygon_3d(cls
, points
):
424 o
= cls
._polygon
_inner
_point
(points
)
435 raise ValueError('degenerate polygon')
439 normprod
= norm_oa
* om
.norm()
440 cosinus
= max(oa
.dot(om
) / normprod
, -1.)
441 sinus
= u
.dot(oa
.cross(om
)) / normprod
442 angle
= math
.acos(cosinus
)
443 angle
= math
.copysign(angle
, sinus
)
445 return sorted(points
, key
=angles
.get
)
449 Returns the vertices of the faces of a polyhedra.
452 for polyhedron
in self
.polyhedra
:
453 vertices
= polyhedron
.vertices()
454 for constraint
in polyhedron
.constraints
:
456 for vertex
in vertices
:
457 if constraint
.subs(vertex
.coordinates()) == 0:
463 def _plot_2d(self
, plot
=None, **kwargs
):
464 import matplotlib
.pyplot
as plt
465 from matplotlib
.patches
import Polygon
468 plot
= fig
.add_subplot(1, 1, 1)
469 xmin
, xmax
= plot
.get_xlim()
470 ymin
, ymax
= plot
.get_ylim()
471 for polyhedron
in self
.polyhedra
:
472 vertices
= polyhedron
._sort
_polygon
_2d
(polyhedron
.vertices())
473 xys
= [tuple(vertex
.values()) for vertex
in vertices
]
475 xmin
, xmax
= min(xmin
, float(min(xs
))), max(xmax
, float(max(xs
)))
476 ymin
, ymax
= min(ymin
, float(min(ys
))), max(ymax
, float(max(ys
)))
477 plot
.add_patch(Polygon(xys
, closed
=True, **kwargs
))
478 plot
.set_xlim(xmin
, xmax
)
479 plot
.set_ylim(ymin
, ymax
)
482 def _plot_3d(self
, plot
=None, **kwargs
):
483 import matplotlib
.pyplot
as plt
484 from mpl_toolkits
.mplot3d
import Axes3D
485 from mpl_toolkits
.mplot3d
.art3d
import Poly3DCollection
491 xmin
, xmax
= axes
.get_xlim()
492 ymin
, ymax
= axes
.get_ylim()
493 zmin
, zmax
= axes
.get_zlim()
495 for vertices
in self
.faces():
496 vertices
= self
._sort
_polygon
_3d
(vertices
)
497 vertices
.append(vertices
[0])
498 face_xyzs
= [tuple(vertex
.values()) for vertex
in vertices
]
499 xs
, ys
, zs
= zip(*face_xyzs
)
500 xmin
, xmax
= min(xmin
, float(min(xs
))), max(xmax
, float(max(xs
)))
501 ymin
, ymax
= min(ymin
, float(min(ys
))), max(ymax
, float(max(ys
)))
502 zmin
, zmax
= min(zmin
, float(min(zs
))), max(zmax
, float(max(zs
)))
503 poly_xyzs
.append(face_xyzs
)
504 collection
= Poly3DCollection(poly_xyzs
, **kwargs
)
505 axes
.add_collection3d(collection
)
506 axes
.set_xlim(xmin
, xmax
)
507 axes
.set_ylim(ymin
, ymax
)
508 axes
.set_zlim(zmin
, zmax
)
512 def plot(self
, plot
=None, **kwargs
):
514 Display plot of this set.
516 if not self
.isbounded():
517 raise ValueError('domain must be bounded')
518 elif self
.dimension
== 2:
519 return self
._plot
_2d
(plot
=plot
, **kwargs
)
520 elif self
.dimension
== 3:
521 return self
._plot
_3d
(plot
=plot
, **kwargs
)
523 raise ValueError('polyhedron must be 2 or 3-dimensional')
525 def __contains__(self
, point
):
526 for polyhedron
in self
.polyhedra
:
527 if point
in polyhedron
:
531 def subs(self
, symbol
, expression
=None):
533 Subsitute the given value into an expression and return the resulting
536 polyhedra
= [polyhedron
.subs(symbol
, expression
)
537 for polyhedron
in self
.polyhedra
]
538 return Domain(*polyhedra
)
541 def _fromislset(cls
, islset
, symbols
):
542 from .polyhedra
import Polyhedron
543 islset
= libisl
.isl_set_remove_divs(islset
)
544 islbsets
= islhelper
.isl_set_basic_sets(islset
)
545 libisl
.isl_set_free(islset
)
547 for islbset
in islbsets
:
548 polyhedron
= Polyhedron
._fromislbasicset
(islbset
, symbols
)
549 polyhedra
.append(polyhedron
)
550 if len(polyhedra
) == 0:
551 from .polyhedra
import Empty
553 elif len(polyhedra
) == 1:
556 self
= object().__new
__(Domain
)
557 self
._polyhedra
= tuple(polyhedra
)
558 self
._symbols
= cls
._xsymbols
(polyhedra
)
559 self
._dimension
= len(self
._symbols
)
563 def _toislset(cls
, polyhedra
, symbols
):
564 polyhedron
= polyhedra
[0]
565 islbset
= polyhedron
._toislbasicset
(polyhedron
.equalities
,
566 polyhedron
.inequalities
, symbols
)
567 islset1
= libisl
.isl_set_from_basic_set(islbset
)
568 for polyhedron
in polyhedra
[1:]:
569 islbset
= polyhedron
._toislbasicset
(polyhedron
.equalities
,
570 polyhedron
.inequalities
, symbols
)
571 islset2
= libisl
.isl_set_from_basic_set(islbset
)
572 islset1
= libisl
.isl_set_union(islset1
, islset2
)
576 def _fromast(cls
, node
):
577 from .polyhedra
import Polyhedron
578 if isinstance(node
, ast
.Module
) and len(node
.body
) == 1:
579 return cls
._fromast
(node
.body
[0])
580 elif isinstance(node
, ast
.Expr
):
581 return cls
._fromast
(node
.value
)
582 elif isinstance(node
, ast
.UnaryOp
):
583 domain
= cls
._fromast
(node
.operand
)
584 if isinstance(node
.operand
, ast
.invert
):
586 elif isinstance(node
, ast
.BinOp
):
587 domain1
= cls
._fromast
(node
.left
)
588 domain2
= cls
._fromast
(node
.right
)
589 if isinstance(node
.op
, ast
.BitAnd
):
590 return And(domain1
, domain2
)
591 elif isinstance(node
.op
, ast
.BitOr
):
592 return Or(domain1
, domain2
)
593 elif isinstance(node
, ast
.Compare
):
596 left
= Expression
._fromast
(node
.left
)
597 for i
in range(len(node
.ops
)):
599 right
= Expression
._fromast
(node
.comparators
[i
])
600 if isinstance(op
, ast
.Lt
):
601 inequalities
.append(right
- left
- 1)
602 elif isinstance(op
, ast
.LtE
):
603 inequalities
.append(right
- left
)
604 elif isinstance(op
, ast
.Eq
):
605 equalities
.append(left
- right
)
606 elif isinstance(op
, ast
.GtE
):
607 inequalities
.append(left
- right
)
608 elif isinstance(op
, ast
.Gt
):
609 inequalities
.append(left
- right
- 1)
614 return Polyhedron(equalities
, inequalities
)
615 raise SyntaxError('invalid syntax')
617 _RE_BRACES
= re
.compile(r
'^\{\s*|\s*\}$')
618 _RE_EQ
= re
.compile(r
'([^<=>])=([^<=>])')
619 _RE_AND
= re
.compile(r
'\band\b|,|&&|/\\|∧|∩')
620 _RE_OR
= re
.compile(r
'\bor\b|;|\|\||\\/|∨|∪')
621 _RE_NOT
= re
.compile(r
'\bnot\b|!|¬')
622 _RE_NUM_VAR
= Expression
._RE
_NUM
_VAR
623 _RE_OPERATORS
= re
.compile(r
'(&|\||~)')
626 def fromstring(cls
, string
):
627 # remove curly brackets
628 string
= cls
._RE
_BRACES
.sub(r
'', string
)
629 # replace '=' by '=='
630 string
= cls
._RE
_EQ
.sub(r
'\1==\2', string
)
631 # replace 'and', 'or', 'not'
632 string
= cls
._RE
_AND
.sub(r
' & ', string
)
633 string
= cls
._RE
_OR
.sub(r
' | ', string
)
634 string
= cls
._RE
_NOT
.sub(r
' ~', string
)
635 # add implicit multiplication operators, e.g. '5x' -> '5*x'
636 string
= cls
._RE
_NUM
_VAR
.sub(r
'\1*\2', string
)
637 # add parentheses to force precedence
638 tokens
= cls
._RE
_OPERATORS
.split(string
)
639 for i
, token
in enumerate(tokens
):
641 token
= '({})'.format(token
)
643 string
= ''.join(tokens
)
644 tree
= ast
.parse(string
, 'eval')
645 return cls
._fromast
(tree
)
648 assert len(self
.polyhedra
) >= 2
649 strings
= [repr(polyhedron
) for polyhedron
in self
.polyhedra
]
650 return 'Or({})'.format(', '.join(strings
))
652 def _repr_latex_(self
):
654 for polyhedron
in self
.polyhedra
:
655 strings
.append('({})'.format(polyhedron
._repr
_latex
_().strip('$')))
656 return '${}$'.format(' \\vee '.join(strings
))
659 def fromsympy(cls
, expr
):
661 from .polyhedra
import Lt
, Le
, Eq
, Ne
, Ge
, Gt
663 sympy
.And
: And
, sympy
.Or
: Or
, sympy
.Not
: Not
,
664 sympy
.Lt
: Lt
, sympy
.Le
: Le
,
665 sympy
.Eq
: Eq
, sympy
.Ne
: Ne
,
666 sympy
.Ge
: Ge
, sympy
.Gt
: Gt
,
668 if expr
.func
in funcmap
:
669 args
= [Domain
.fromsympy(arg
) for arg
in expr
.args
]
670 return funcmap
[expr
.func
](*args
)
671 elif isinstance(expr
, sympy
.Expr
):
672 return Expression
.fromsympy(expr
)
673 raise ValueError('non-domain expression: {!r}'.format(expr
))
677 polyhedra
= [polyhedron
.tosympy() for polyhedron
in polyhedra
]
678 return sympy
.Or(*polyhedra
)
683 Return the intersection of two sets as a new set.
685 if len(domains
) == 0:
686 from .polyhedra
import Universe
689 return domains
[0].intersection(*domains
[1:])
693 Return the union of sets as a new set.
695 if len(domains
) == 0:
696 from .polyhedra
import Empty
699 return domains
[0].union(*domains
[1:])
703 Returns the complement of this set.