787e965ee6d36a506a472843fa82d68420d606ce
[linpy.git] / pypol / polyhedra.py
1 import ast
2 import functools
3 import numbers
4 import re
5
6 from . import islhelper
7
8 from .islhelper import mainctx, libisl
9 from .linexprs import Expression, Constant
10 from .domains import Domain
11
12
13 __all__ = [
14 'Polyhedron',
15 'Lt', 'Le', 'Eq', 'Ne', 'Ge', 'Gt',
16 'Empty', 'Universe',
17 ]
18
19
20 class Polyhedron(Domain):
21
22 __slots__ = (
23 '_equalities',
24 '_inequalities',
25 '_constraints',
26 '_symbols',
27 '_dimension',
28 )
29
30 def __new__(cls, equalities=None, inequalities=None):
31 if isinstance(equalities, str):
32 if inequalities is not None:
33 raise TypeError('too many arguments')
34 return cls.fromstring(equalities)
35 elif isinstance(equalities, Polyhedron):
36 if inequalities is not None:
37 raise TypeError('too many arguments')
38 return equalities
39 elif isinstance(equalities, Domain):
40 if inequalities is not None:
41 raise TypeError('too many arguments')
42 return equalities.polyhedral_hull()
43 if equalities is None:
44 equalities = []
45 else:
46 for i, equality in enumerate(equalities):
47 if not isinstance(equality, Expression):
48 raise TypeError('equalities must be linear expressions')
49 equalities[i] = equality._toint()
50 if inequalities is None:
51 inequalities = []
52 else:
53 for i, inequality in enumerate(inequalities):
54 if not isinstance(inequality, Expression):
55 raise TypeError('inequalities must be linear expressions')
56 inequalities[i] = inequality._toint()
57 symbols = cls._xsymbols(equalities + inequalities)
58 islbset = cls._toislbasicset(equalities, inequalities, symbols)
59 return cls._fromislbasicset(islbset, symbols)
60
61 @property
62 def equalities(self):
63 return self._equalities
64
65 @property
66 def inequalities(self):
67 return self._inequalities
68
69 @property
70 def constraints(self):
71 return self._constraints
72
73 @property
74 def polyhedra(self):
75 return self,
76
77 def disjoint(self):
78 return self
79
80 def isuniverse(self):
81 islbset = self._toislbasicset(self.equalities, self.inequalities,
82 self.symbols)
83 universe = bool(libisl.isl_basic_set_is_universe(islbset))
84 libisl.isl_basic_set_free(islbset)
85 return universe
86
87 def polyhedral_hull(self):
88 return self
89
90 @classmethod
91 def _fromislbasicset(cls, islbset, symbols):
92 islconstraints = islhelper.isl_basic_set_constraints(islbset)
93 equalities = []
94 inequalities = []
95 for islconstraint in islconstraints:
96 islpr = libisl.isl_printer_to_str(mainctx)
97 constant = libisl.isl_constraint_get_constant_val(islconstraint)
98 constant = islhelper.isl_val_to_int(constant)
99 coefficients = {}
100 for dim, symbol in enumerate(symbols):
101 coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint, libisl.isl_dim_set, dim)
102 coefficient = islhelper.isl_val_to_int(coefficient)
103 if coefficient != 0:
104 coefficients[symbol] = coefficient
105 expression = Expression(coefficients, constant)
106 if libisl.isl_constraint_is_equality(islconstraint):
107 equalities.append(expression)
108 else:
109 inequalities.append(expression)
110 libisl.isl_basic_set_free(islbset)
111 self = object().__new__(Polyhedron)
112 self._equalities = tuple(equalities)
113 self._inequalities = tuple(inequalities)
114 self._constraints = tuple(equalities + inequalities)
115 self._symbols = cls._xsymbols(self._constraints)
116 self._dimension = len(self._symbols)
117 return self
118
119 @classmethod
120 def _toislbasicset(cls, equalities, inequalities, symbols):
121 dimension = len(symbols)
122 islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
123 islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
124 islls = libisl.isl_local_space_from_space(islsp)
125 for equality in equalities:
126 isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls))
127 for symbol, coefficient in equality.coefficients():
128 val = str(coefficient).encode()
129 val = libisl.isl_val_read_from_str(mainctx, val)
130 sid = symbols.index(symbol)
131 isleq = libisl.isl_constraint_set_coefficient_val(isleq,
132 libisl.isl_dim_set, sid, val)
133 if equality.constant != 0:
134 val = str(equality.constant).encode()
135 val = libisl.isl_val_read_from_str(mainctx, val)
136 isleq = libisl.isl_constraint_set_constant_val(isleq, val)
137 islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
138 for inequality in inequalities:
139 islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls))
140 for symbol, coefficient in inequality.coefficients():
141 val = str(coefficient).encode()
142 val = libisl.isl_val_read_from_str(mainctx, val)
143 sid = symbols.index(symbol)
144 islin = libisl.isl_constraint_set_coefficient_val(islin,
145 libisl.isl_dim_set, sid, val)
146 if inequality.constant != 0:
147 val = str(inequality.constant).encode()
148 val = libisl.isl_val_read_from_str(mainctx, val)
149 islin = libisl.isl_constraint_set_constant_val(islin, val)
150 islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
151 return islbset
152
153 @classmethod
154 def _fromast(cls, node):
155 if isinstance(node, ast.Module) and len(node.body) == 1:
156 return cls._fromast(node.body[0])
157 elif isinstance(node, ast.Expr):
158 return cls._fromast(node.value)
159 elif isinstance(node, ast.BinOp) and isinstance(node.op, ast.BitAnd):
160 equalities1, inequalities1 = cls._fromast(node.left)
161 equalities2, inequalities2 = cls._fromast(node.right)
162 equalities = equalities1 + equalities2
163 inequalities = inequalities1 + inequalities2
164 return equalities, inequalities
165 elif isinstance(node, ast.Compare):
166 equalities = []
167 inequalities = []
168 left = Expression._fromast(node.left)
169 for i in range(len(node.ops)):
170 op = node.ops[i]
171 right = Expression._fromast(node.comparators[i])
172 if isinstance(op, ast.Lt):
173 inequalities.append(right - left - 1)
174 elif isinstance(op, ast.LtE):
175 inequalities.append(right - left)
176 elif isinstance(op, ast.Eq):
177 equalities.append(left - right)
178 elif isinstance(op, ast.GtE):
179 inequalities.append(left - right)
180 elif isinstance(op, ast.Gt):
181 inequalities.append(left - right - 1)
182 else:
183 break
184 left = right
185 else:
186 return equalities, inequalities
187 raise SyntaxError('invalid syntax')
188
189 @classmethod
190 def fromstring(cls, string):
191 string = string.strip()
192 string = re.sub(r'^\{\s*|\s*\}$', '', string)
193 string = re.sub(r'([^<=>])=([^<=>])', r'\1==\2', string)
194 string = re.sub(r'(\d+|\))\s*([^\W\d_]\w*|\()', r'\1*\2', string)
195 tokens = re.split(r',|;|and|&&|/\\|∧', string, flags=re.I)
196 tokens = ['({})'.format(token) for token in tokens]
197 string = ' & '.join(tokens)
198 tree = ast.parse(string, 'eval')
199 equalities, inequalities = cls._fromast(tree)
200 return cls(equalities, inequalities)
201
202 def __repr__(self):
203 if self.isempty():
204 return 'Empty'
205 elif self.isuniverse():
206 return 'Universe'
207 else:
208 strings = []
209 for equality in self.equalities:
210 strings.append('Eq({}, 0)'.format(equality))
211 for inequality in self.inequalities:
212 strings.append('Ge({}, 0)'.format(inequality))
213 if len(strings) == 1:
214 return strings[0]
215 else:
216 return 'And({})'.format(', '.join(strings))
217
218 @classmethod
219 def _fromsympy(cls, expr):
220 import sympy
221 equalities = []
222 inequalities = []
223 if expr.func == sympy.And:
224 for arg in expr.args:
225 arg_eqs, arg_ins = cls._fromsympy(arg)
226 equalities.extend(arg_eqs)
227 inequalities.extend(arg_ins)
228 elif expr.func == sympy.Eq:
229 expr = Expression.fromsympy(expr.args[0] - expr.args[1])
230 equalities.append(expr)
231 else:
232 if expr.func == sympy.Lt:
233 expr = Expression.fromsympy(expr.args[1] - expr.args[0] - 1)
234 elif expr.func == sympy.Le:
235 expr = Expression.fromsympy(expr.args[1] - expr.args[0])
236 elif expr.func == sympy.Ge:
237 expr = Expression.fromsympy(expr.args[0] - expr.args[1])
238 elif expr.func == sympy.Gt:
239 expr = Expression.fromsympy(expr.args[0] - expr.args[1] - 1)
240 else:
241 raise ValueError('non-polyhedral expression: {!r}'.format(expr))
242 inequalities.append(expr)
243 return equalities, inequalities
244
245 @classmethod
246 def fromsympy(cls, expr):
247 import sympy
248 equalities, inequalities = cls._fromsympy(expr)
249 return cls(equalities, inequalities)
250
251 def tosympy(self):
252 import sympy
253 constraints = []
254 for equality in self.equalities:
255 constraints.append(sympy.Eq(equality.tosympy(), 0))
256 for inequality in self.inequalities:
257 constraints.append(sympy.Ge(inequality.tosympy(), 0))
258 return sympy.And(*constraints)
259
260
261 def _polymorphic(func):
262 @functools.wraps(func)
263 def wrapper(left, right):
264 if isinstance(left, numbers.Rational):
265 left = Constant(left)
266 elif not isinstance(left, Expression):
267 raise TypeError('left must be a a rational number '
268 'or a linear expression')
269 if isinstance(right, numbers.Rational):
270 right = Constant(right)
271 elif not isinstance(right, Expression):
272 raise TypeError('right must be a a rational number '
273 'or a linear expression')
274 return func(left, right)
275 return wrapper
276
277 @_polymorphic
278 def Lt(left, right):
279 return Polyhedron([], [right - left - 1])
280
281 @_polymorphic
282 def Le(left, right):
283 return Polyhedron([], [right - left])
284
285 @_polymorphic
286 def Eq(left, right):
287 return Polyhedron([left - right], [])
288
289 @_polymorphic
290 def Ne(left, right):
291 return ~Eq(left, right)
292
293 @_polymorphic
294 def Gt(left, right):
295 return Polyhedron([], [left - right - 1])
296
297 @_polymorphic
298 def Ge(left, right):
299 return Polyhedron([], [left - right])
300
301
302 Empty = Eq(1, 0)
303
304 Universe = Polyhedron([])