1 # Copyright 2014 MINES ParisTech
3 # This file is part of Linpy.
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.
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.
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/>.
23 from fractions
import Fraction
25 from . import islhelper
26 from .islhelper
import mainctx
, libisl
27 from .linexprs
import Expression
, Symbol
, Rational
28 from .geometry
import GeometricObject
, Point
, Vector
37 @functools.total_ordering
38 class Domain(GeometricObject
):
46 def __new__(cls
, *polyhedra
):
47 from .polyhedra
import Polyhedron
48 if len(polyhedra
) == 1:
49 argument
= polyhedra
[0]
50 if isinstance(argument
, str):
51 return cls
.fromstring(argument
)
52 elif isinstance(argument
, GeometricObject
):
53 return argument
.aspolyhedron()
55 raise TypeError('argument must be a string '
56 'or a GeometricObject instance')
58 for polyhedron
in polyhedra
:
59 if not isinstance(polyhedron
, Polyhedron
):
60 raise TypeError('arguments must be Polyhedron instances')
61 symbols
= cls
._xsymbols
(polyhedra
)
62 islset
= cls
._toislset
(polyhedra
, symbols
)
63 return cls
._fromislset
(islset
, symbols
)
66 def _xsymbols(cls
, iterator
):
68 Return the ordered tuple of symbols present in iterator.
72 symbols
.update(item
.symbols
)
73 return tuple(sorted(symbols
, key
=Symbol
.sortkey
))
77 return self
._polyhedra
85 return self
._dimension
89 Returns this set as disjoint.
91 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
92 islset
= libisl
.isl_set_make_disjoint(mainctx
, islset
)
93 return self
._fromislset
(islset
, self
.symbols
)
97 Returns true if this set is an Empty set.
99 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
100 empty
= bool(libisl
.isl_set_is_empty(islset
))
101 libisl
.isl_set_free(islset
)
105 return not self
.isempty()
107 def isuniverse(self
):
109 Returns true if this set is the Universe set.
111 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
112 universe
= bool(libisl
.isl_set_plain_is_universe(islset
))
113 libisl
.isl_set_free(islset
)
118 Returns true if this set is bounded.
120 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
121 bounded
= bool(libisl
.isl_set_is_bounded(islset
))
122 libisl
.isl_set_free(islset
)
125 def __eq__(self
, other
):
127 Returns true if two sets are equal.
129 symbols
= self
._xsymbols
([self
, other
])
130 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
131 islset2
= other
._toislset
(other
.polyhedra
, symbols
)
132 equal
= bool(libisl
.isl_set_is_equal(islset1
, islset2
))
133 libisl
.isl_set_free(islset1
)
134 libisl
.isl_set_free(islset2
)
137 def isdisjoint(self
, other
):
139 Return True if two sets have a null intersection.
141 symbols
= self
._xsymbols
([self
, other
])
142 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
143 islset2
= self
._toislset
(other
.polyhedra
, symbols
)
144 equal
= bool(libisl
.isl_set_is_disjoint(islset1
, islset2
))
145 libisl
.isl_set_free(islset1
)
146 libisl
.isl_set_free(islset2
)
149 def issubset(self
, other
):
151 Report whether another set contains this set.
153 symbols
= self
._xsymbols
([self
, other
])
154 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
155 islset2
= self
._toislset
(other
.polyhedra
, symbols
)
156 equal
= bool(libisl
.isl_set_is_subset(islset1
, islset2
))
157 libisl
.isl_set_free(islset1
)
158 libisl
.isl_set_free(islset2
)
161 def __le__(self
, other
):
163 Returns true if this set is less than or equal to another set.
165 return self
.issubset(other
)
167 def __lt__(self
, other
):
169 Returns true if this set is less than another set.
171 symbols
= self
._xsymbols
([self
, other
])
172 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
173 islset2
= self
._toislset
(other
.polyhedra
, symbols
)
174 equal
= bool(libisl
.isl_set_is_strict_subset(islset1
, islset2
))
175 libisl
.isl_set_free(islset1
)
176 libisl
.isl_set_free(islset2
)
179 def complement(self
):
181 Returns the complement of this set.
183 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
184 islset
= libisl
.isl_set_complement(islset
)
185 return self
._fromislset
(islset
, self
.symbols
)
187 def __invert__(self
):
189 Returns the complement of this set.
191 return self
.complement()
195 Returns a set without redundant constraints.
197 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
198 islset
= libisl
.isl_set_remove_redundancies(islset
)
199 return self
._fromislset
(islset
, self
.symbols
)
201 def aspolyhedron(self
):
203 Returns polyhedral hull of set.
205 from .polyhedra
import Polyhedron
206 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
207 islbset
= libisl
.isl_set_polyhedral_hull(islset
)
208 return Polyhedron
._fromislbasicset
(islbset
, self
.symbols
)
213 def project(self
, dims
):
215 Return new set with given dimensions removed.
217 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
219 for index
, symbol
in reversed(list(enumerate(self
.symbols
))):
223 islset
= libisl
.isl_set_project_out(islset
, libisl
.isl_dim_set
, index
+ 1, n
)
226 islset
= libisl
.isl_set_project_out(islset
, libisl
.isl_dim_set
, 0, n
)
227 dims
= [symbol
for symbol
in self
.symbols
if symbol
not in dims
]
228 return Domain
._fromislset
(islset
, dims
)
232 Returns a single subset of the input.
234 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
235 islpoint
= libisl
.isl_set_sample_point(islset
)
236 if bool(libisl
.isl_point_is_void(islpoint
)):
237 libisl
.isl_point_free(islpoint
)
238 raise ValueError('domain must be non-empty')
240 for index
, symbol
in enumerate(self
.symbols
):
241 coordinate
= libisl
.isl_point_get_coordinate_val(islpoint
,
242 libisl
.isl_dim_set
, index
)
243 coordinate
= islhelper
.isl_val_to_int(coordinate
)
244 point
[symbol
] = coordinate
245 libisl
.isl_point_free(islpoint
)
248 def intersection(self
, *others
):
250 Return the intersection of two sets as a new set.
254 symbols
= self
._xsymbols
((self
,) + others
)
255 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
257 islset2
= other
._toislset
(other
.polyhedra
, symbols
)
258 islset1
= libisl
.isl_set_intersect(islset1
, islset2
)
259 return self
._fromislset
(islset1
, symbols
)
261 def __and__(self
, other
):
263 Return the intersection of two sets as a new set.
265 return self
.intersection(other
)
267 def union(self
, *others
):
269 Return the union of sets as a new set.
273 symbols
= self
._xsymbols
((self
,) + others
)
274 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
276 islset2
= other
._toislset
(other
.polyhedra
, symbols
)
277 islset1
= libisl
.isl_set_union(islset1
, islset2
)
278 return self
._fromislset
(islset1
, symbols
)
280 def __or__(self
, other
):
282 Return a new set with elements from both sets.
284 return self
.union(other
)
286 def __add__(self
, other
):
288 Return new set containing all elements in both sets.
290 return self
.union(other
)
292 def difference(self
, other
):
294 Return the difference of two sets as a new set.
296 symbols
= self
._xsymbols
([self
, other
])
297 islset1
= self
._toislset
(self
.polyhedra
, symbols
)
298 islset2
= other
._toislset
(other
.polyhedra
, symbols
)
299 islset
= libisl
.isl_set_subtract(islset1
, islset2
)
300 return self
._fromislset
(islset
, symbols
)
302 def __sub__(self
, other
):
304 Return the difference of two sets as a new set.
306 return self
.difference(other
)
310 Return a new set containing the lexicographic minimum of the elements in the set.
312 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
313 islset
= libisl
.isl_set_lexmin(islset
)
314 return self
._fromislset
(islset
, self
.symbols
)
318 Return a new set containing the lexicographic maximum of the elements in the set.
320 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
321 islset
= libisl
.isl_set_lexmax(islset
)
322 return self
._fromislset
(islset
, self
.symbols
)
325 def involves_vars(self
, vars):
327 Returns true if a set depends on given dimensions.
329 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
331 symbols
= sorted(list(self
.symbols
))
336 first
= symbols
.index(dims
[0])
342 value
= bool(libisl
.isl_set_involves_dims(islset
, libisl
.isl_dim_set
, first
, n
))
343 libisl
.isl_set_free(islset
)
346 _RE_COORDINATE
= re
.compile(r
'\((?P<num>\-?\d+)\)(/(?P<den>\d+))?')
350 Return a list of vertices for this Polygon.
352 from .polyhedra
import Polyhedron
353 if not self
.isbounded():
354 raise ValueError('domain must be bounded')
355 islbset
= self
._toislbasicset
(self
.equalities
, self
.inequalities
, self
.symbols
)
356 vertices
= libisl
.isl_basic_set_compute_vertices(islbset
);
357 vertices
= islhelper
.isl_vertices_vertices(vertices
)
359 for vertex
in vertices
:
360 expr
= libisl
.isl_vertex_get_expr(vertex
)
362 if islhelper
.isl_version
< '0.13':
363 constraints
= islhelper
.isl_basic_set_constraints(expr
)
364 for constraint
in constraints
:
365 constant
= libisl
.isl_constraint_get_constant_val(constraint
)
366 constant
= islhelper
.isl_val_to_int(constant
)
367 for index
, symbol
in enumerate(self
.symbols
):
368 coefficient
= libisl
.isl_constraint_get_coefficient_val(constraint
,
369 libisl
.isl_dim_set
, index
)
370 coefficient
= islhelper
.isl_val_to_int(coefficient
)
372 coordinate
= -Fraction(constant
, coefficient
)
373 coordinates
.append((symbol
, coordinate
))
375 string
= islhelper
.isl_multi_aff_to_str(expr
)
376 matches
= self
._RE
_COORDINATE
.finditer(string
)
377 for symbol
, match
in zip(self
.symbols
, matches
):
378 numerator
= int(match
.group('num'))
379 denominator
= match
.group('den')
380 denominator
= 1 if denominator
is None else int(denominator
)
381 coordinate
= Fraction(numerator
, denominator
)
382 coordinates
.append((symbol
, coordinate
))
383 points
.append(Point(coordinates
))
388 Returns the points contained in the set.
390 if not self
.isbounded():
391 raise ValueError('domain must be bounded')
392 from .polyhedra
import Universe
, Eq
393 islset
= self
._toislset
(self
.polyhedra
, self
.symbols
)
394 islpoints
= islhelper
.isl_set_points(islset
)
396 for islpoint
in islpoints
:
398 for index
, symbol
in enumerate(self
.symbols
):
399 coordinate
= libisl
.isl_point_get_coordinate_val(islpoint
,
400 libisl
.isl_dim_set
, index
)
401 coordinate
= islhelper
.isl_val_to_int(coordinate
)
402 coordinates
[symbol
] = coordinate
403 points
.append(Point(coordinates
))
407 def _polygon_inner_point(cls
, points
):
408 symbols
= points
[0].symbols
409 coordinates
= {symbol
: 0 for symbol
in symbols
}
411 for symbol
, coordinate
in point
.coordinates():
412 coordinates
[symbol
] += coordinate
413 for symbol
in symbols
:
414 coordinates
[symbol
] /= len(points
)
415 return Point(coordinates
)
418 def _sort_polygon_2d(cls
, points
):
421 o
= cls
._polygon
_inner
_point
(points
)
425 dx
, dy
= (coordinate
for symbol
, coordinate
in om
.coordinates())
426 angle
= math
.atan2(dy
, dx
)
428 return sorted(points
, key
=angles
.get
)
431 def _sort_polygon_3d(cls
, points
):
434 o
= cls
._polygon
_inner
_point
(points
)
445 raise ValueError('degenerate polygon')
449 normprod
= norm_oa
* om
.norm()
450 cosinus
= max(oa
.dot(om
) / normprod
, -1.)
451 sinus
= u
.dot(oa
.cross(om
)) / normprod
452 angle
= math
.acos(cosinus
)
453 angle
= math
.copysign(angle
, sinus
)
455 return sorted(points
, key
=angles
.get
)
459 Returns the vertices of the faces of a polyhedra.
462 for polyhedron
in self
.polyhedra
:
463 vertices
= polyhedron
.vertices()
464 for constraint
in polyhedron
.constraints
:
466 for vertex
in vertices
:
467 if constraint
.subs(vertex
.coordinates()) == 0:
473 def _plot_2d(self
, plot
=None, **kwargs
):
474 import matplotlib
.pyplot
as plt
475 from matplotlib
.patches
import Polygon
478 plot
= fig
.add_subplot(1, 1, 1)
479 xmin
, xmax
= plot
.get_xlim()
480 ymin
, ymax
= plot
.get_ylim()
481 for polyhedron
in self
.polyhedra
:
482 vertices
= polyhedron
._sort
_polygon
_2d
(polyhedron
.vertices())
483 xys
= [tuple(vertex
.values()) for vertex
in vertices
]
485 xmin
, xmax
= min(xmin
, float(min(xs
))), max(xmax
, float(max(xs
)))
486 ymin
, ymax
= min(ymin
, float(min(ys
))), max(ymax
, float(max(ys
)))
487 plot
.add_patch(Polygon(xys
, closed
=True, **kwargs
))
488 plot
.set_xlim(xmin
, xmax
)
489 plot
.set_ylim(ymin
, ymax
)
492 def _plot_3d(self
, plot
=None, **kwargs
):
493 import matplotlib
.pyplot
as plt
494 from mpl_toolkits
.mplot3d
import Axes3D
495 from mpl_toolkits
.mplot3d
.art3d
import Poly3DCollection
501 xmin
, xmax
= axes
.get_xlim()
502 ymin
, ymax
= axes
.get_ylim()
503 zmin
, zmax
= axes
.get_zlim()
505 for vertices
in self
.faces():
506 vertices
= self
._sort
_polygon
_3d
(vertices
)
507 vertices
.append(vertices
[0])
508 face_xyzs
= [tuple(vertex
.values()) for vertex
in vertices
]
509 xs
, ys
, zs
= zip(*face_xyzs
)
510 xmin
, xmax
= min(xmin
, float(min(xs
))), max(xmax
, float(max(xs
)))
511 ymin
, ymax
= min(ymin
, float(min(ys
))), max(ymax
, float(max(ys
)))
512 zmin
, zmax
= min(zmin
, float(min(zs
))), max(zmax
, float(max(zs
)))
513 poly_xyzs
.append(face_xyzs
)
514 collection
= Poly3DCollection(poly_xyzs
, **kwargs
)
515 axes
.add_collection3d(collection
)
516 axes
.set_xlim(xmin
, xmax
)
517 axes
.set_ylim(ymin
, ymax
)
518 axes
.set_zlim(zmin
, zmax
)
522 def plot(self
, plot
=None, **kwargs
):
524 Display plot of this set.
526 if not self
.isbounded():
527 raise ValueError('domain must be bounded')
528 elif self
.dimension
== 2:
529 return self
._plot
_2d
(plot
=plot
, **kwargs
)
530 elif self
.dimension
== 3:
531 return self
._plot
_3d
(plot
=plot
, **kwargs
)
533 raise ValueError('polyhedron must be 2 or 3-dimensional')
535 def __contains__(self
, point
):
536 for polyhedron
in self
.polyhedra
:
537 if point
in polyhedron
:
541 def subs(self
, symbol
, expression
=None):
543 Subsitute the given value into an expression and return the resulting
546 polyhedra
= [polyhedron
.subs(symbol
, expression
)
547 for polyhedron
in self
.polyhedra
]
548 return Domain(*polyhedra
)
551 def _fromislset(cls
, islset
, symbols
):
552 from .polyhedra
import Polyhedron
553 islset
= libisl
.isl_set_remove_divs(islset
)
554 islbsets
= islhelper
.isl_set_basic_sets(islset
)
555 libisl
.isl_set_free(islset
)
557 for islbset
in islbsets
:
558 polyhedron
= Polyhedron
._fromislbasicset
(islbset
, symbols
)
559 polyhedra
.append(polyhedron
)
560 if len(polyhedra
) == 0:
561 from .polyhedra
import Empty
563 elif len(polyhedra
) == 1:
566 self
= object().__new
__(Domain
)
567 self
._polyhedra
= tuple(polyhedra
)
568 self
._symbols
= cls
._xsymbols
(polyhedra
)
569 self
._dimension
= len(self
._symbols
)
573 def _toislset(cls
, polyhedra
, symbols
):
574 polyhedron
= polyhedra
[0]
575 islbset
= polyhedron
._toislbasicset
(polyhedron
.equalities
,
576 polyhedron
.inequalities
, symbols
)
577 islset1
= libisl
.isl_set_from_basic_set(islbset
)
578 for polyhedron
in polyhedra
[1:]:
579 islbset
= polyhedron
._toislbasicset
(polyhedron
.equalities
,
580 polyhedron
.inequalities
, symbols
)
581 islset2
= libisl
.isl_set_from_basic_set(islbset
)
582 islset1
= libisl
.isl_set_union(islset1
, islset2
)
586 def _fromast(cls
, node
):
587 from .polyhedra
import Polyhedron
588 if isinstance(node
, ast
.Module
) and len(node
.body
) == 1:
589 return cls
._fromast
(node
.body
[0])
590 elif isinstance(node
, ast
.Expr
):
591 return cls
._fromast
(node
.value
)
592 elif isinstance(node
, ast
.UnaryOp
):
593 domain
= cls
._fromast
(node
.operand
)
594 if isinstance(node
.operand
, ast
.invert
):
596 elif isinstance(node
, ast
.BinOp
):
597 domain1
= cls
._fromast
(node
.left
)
598 domain2
= cls
._fromast
(node
.right
)
599 if isinstance(node
.op
, ast
.BitAnd
):
600 return And(domain1
, domain2
)
601 elif isinstance(node
.op
, ast
.BitOr
):
602 return Or(domain1
, domain2
)
603 elif isinstance(node
, ast
.Compare
):
606 left
= Expression
._fromast
(node
.left
)
607 for i
in range(len(node
.ops
)):
609 right
= Expression
._fromast
(node
.comparators
[i
])
610 if isinstance(op
, ast
.Lt
):
611 inequalities
.append(right
- left
- 1)
612 elif isinstance(op
, ast
.LtE
):
613 inequalities
.append(right
- left
)
614 elif isinstance(op
, ast
.Eq
):
615 equalities
.append(left
- right
)
616 elif isinstance(op
, ast
.GtE
):
617 inequalities
.append(left
- right
)
618 elif isinstance(op
, ast
.Gt
):
619 inequalities
.append(left
- right
- 1)
624 return Polyhedron(equalities
, inequalities
)
625 raise SyntaxError('invalid syntax')
627 _RE_BRACES
= re
.compile(r
'^\{\s*|\s*\}$')
628 _RE_EQ
= re
.compile(r
'([^<=>])=([^<=>])')
629 _RE_AND
= re
.compile(r
'\band\b|,|&&|/\\|∧|∩')
630 _RE_OR
= re
.compile(r
'\bor\b|;|\|\||\\/|∨|∪')
631 _RE_NOT
= re
.compile(r
'\bnot\b|!|¬')
632 _RE_NUM_VAR
= Expression
._RE
_NUM
_VAR
633 _RE_OPERATORS
= re
.compile(r
'(&|\||~)')
636 def fromstring(cls
, string
):
637 # remove curly brackets
638 string
= cls
._RE
_BRACES
.sub(r
'', string
)
639 # replace '=' by '=='
640 string
= cls
._RE
_EQ
.sub(r
'\1==\2', string
)
641 # replace 'and', 'or', 'not'
642 string
= cls
._RE
_AND
.sub(r
' & ', string
)
643 string
= cls
._RE
_OR
.sub(r
' | ', string
)
644 string
= cls
._RE
_NOT
.sub(r
' ~', string
)
645 # add implicit multiplication operators, e.g. '5x' -> '5*x'
646 string
= cls
._RE
_NUM
_VAR
.sub(r
'\1*\2', string
)
647 # add parentheses to force precedence
648 tokens
= cls
._RE
_OPERATORS
.split(string
)
649 for i
, token
in enumerate(tokens
):
651 token
= '({})'.format(token
)
653 string
= ''.join(tokens
)
654 tree
= ast
.parse(string
, 'eval')
655 return cls
._fromast
(tree
)
658 assert len(self
.polyhedra
) >= 2
659 strings
= [repr(polyhedron
) for polyhedron
in self
.polyhedra
]
660 return 'Or({})'.format(', '.join(strings
))
662 def _repr_latex_(self
):
664 for polyhedron
in self
.polyhedra
:
665 strings
.append('({})'.format(polyhedron
._repr
_latex
_().strip('$')))
666 return '${}$'.format(' \\vee '.join(strings
))
669 def fromsympy(cls
, expr
):
671 from .polyhedra
import Lt
, Le
, Eq
, Ne
, Ge
, Gt
673 sympy
.And
: And
, sympy
.Or
: Or
, sympy
.Not
: Not
,
674 sympy
.Lt
: Lt
, sympy
.Le
: Le
,
675 sympy
.Eq
: Eq
, sympy
.Ne
: Ne
,
676 sympy
.Ge
: Ge
, sympy
.Gt
: Gt
,
678 if expr
.func
in funcmap
:
679 args
= [Domain
.fromsympy(arg
) for arg
in expr
.args
]
680 return funcmap
[expr
.func
](*args
)
681 elif isinstance(expr
, sympy
.Expr
):
682 return Expression
.fromsympy(expr
)
683 raise ValueError('non-domain expression: {!r}'.format(expr
))
687 polyhedra
= [polyhedron
.tosympy() for polyhedron
in polyhedra
]
688 return sympy
.Or(*polyhedra
)
693 Return the intersection of two sets as a new set.
695 if len(domains
) == 0:
696 from .polyhedra
import Universe
699 return domains
[0].intersection(*domains
[1:])
703 Return the union of sets as a new set.
705 if len(domains
) == 0:
706 from .polyhedra
import Empty
709 return domains
[0].union(*domains
[1:])
713 Returns the complement of this set.