Cleaner implementation of Rational
[linpy.git] / pypol / polyhedra.py
1 import functools
2 import math
3 import numbers
4
5 from . import islhelper
6
7 from .islhelper import mainctx, libisl
8 from .geometry import GeometricObject, Point, Vector
9 from .linexprs import Expression, Symbol, Rational
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, GeometricObject):
36 if inequalities is not None:
37 raise TypeError('too many arguments')
38 return equalities.aspolyhedron()
39 if equalities is None:
40 equalities = []
41 else:
42 for i, equality in enumerate(equalities):
43 if not isinstance(equality, Expression):
44 raise TypeError('equalities must be linear expressions')
45 equalities[i] = equality.scaleint()
46 if inequalities is None:
47 inequalities = []
48 else:
49 for i, inequality in enumerate(inequalities):
50 if not isinstance(inequality, Expression):
51 raise TypeError('inequalities must be linear expressions')
52 inequalities[i] = inequality.scaleint()
53 symbols = cls._xsymbols(equalities + inequalities)
54 islbset = cls._toislbasicset(equalities, inequalities, symbols)
55 return cls._fromislbasicset(islbset, symbols)
56
57 @property
58 def equalities(self):
59 return self._equalities
60
61 @property
62 def inequalities(self):
63 return self._inequalities
64
65 @property
66 def constraints(self):
67 return self._constraints
68
69 @property
70 def polyhedra(self):
71 return self,
72
73 def disjoint(self):
74 return self
75
76 def isuniverse(self):
77 islbset = self._toislbasicset(self.equalities, self.inequalities,
78 self.symbols)
79 universe = bool(libisl.isl_basic_set_is_universe(islbset))
80 libisl.isl_basic_set_free(islbset)
81 return universe
82
83 def aspolyhedron(self):
84 return self
85
86 def __contains__(self, point):
87 if not isinstance(point, Point):
88 raise TypeError('point must be a Point instance')
89 if self.symbols != point.symbols:
90 raise ValueError('arguments must belong to the same space')
91 for equality in self.equalities:
92 if equality.subs(point.coordinates()) != 0:
93 return False
94 for inequality in self.inequalities:
95 if inequality.subs(point.coordinates()) < 0:
96 return False
97 return True
98
99 def subs(self, symbol, expression=None):
100 equalities = [equality.subs(symbol, expression)
101 for equality in self.equalities]
102 inequalities = [inequality.subs(symbol, expression)
103 for inequality in self.inequalities]
104 return Polyhedron(equalities, inequalities)
105
106 @classmethod
107 def _fromislbasicset(cls, islbset, symbols):
108 islconstraints = islhelper.isl_basic_set_constraints(islbset)
109 equalities = []
110 inequalities = []
111 for islconstraint in islconstraints:
112 constant = libisl.isl_constraint_get_constant_val(islconstraint)
113 constant = islhelper.isl_val_to_int(constant)
114 coefficients = {}
115 for index, symbol in enumerate(symbols):
116 coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint,
117 libisl.isl_dim_set, index)
118 coefficient = islhelper.isl_val_to_int(coefficient)
119 if coefficient != 0:
120 coefficients[symbol] = coefficient
121 expression = Expression(coefficients, constant)
122 if libisl.isl_constraint_is_equality(islconstraint):
123 equalities.append(expression)
124 else:
125 inequalities.append(expression)
126 libisl.isl_basic_set_free(islbset)
127 self = object().__new__(Polyhedron)
128 self._equalities = tuple(equalities)
129 self._inequalities = tuple(inequalities)
130 self._constraints = tuple(equalities + inequalities)
131 self._symbols = cls._xsymbols(self._constraints)
132 self._dimension = len(self._symbols)
133 return self
134
135 @classmethod
136 def _toislbasicset(cls, equalities, inequalities, symbols):
137 dimension = len(symbols)
138 indices = {symbol: index for index, symbol in enumerate(symbols)}
139 islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
140 islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
141 islls = libisl.isl_local_space_from_space(islsp)
142 for equality in equalities:
143 isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls))
144 for symbol, coefficient in equality.coefficients():
145 islval = str(coefficient).encode()
146 islval = libisl.isl_val_read_from_str(mainctx, islval)
147 index = indices[symbol]
148 isleq = libisl.isl_constraint_set_coefficient_val(isleq,
149 libisl.isl_dim_set, index, islval)
150 if equality.constant != 0:
151 islval = str(equality.constant).encode()
152 islval = libisl.isl_val_read_from_str(mainctx, islval)
153 isleq = libisl.isl_constraint_set_constant_val(isleq, islval)
154 islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
155 for inequality in inequalities:
156 islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls))
157 for symbol, coefficient in inequality.coefficients():
158 islval = str(coefficient).encode()
159 islval = libisl.isl_val_read_from_str(mainctx, islval)
160 index = indices[symbol]
161 islin = libisl.isl_constraint_set_coefficient_val(islin,
162 libisl.isl_dim_set, index, islval)
163 if inequality.constant != 0:
164 islval = str(inequality.constant).encode()
165 islval = libisl.isl_val_read_from_str(mainctx, islval)
166 islin = libisl.isl_constraint_set_constant_val(islin, islval)
167 islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
168 return islbset
169
170 @classmethod
171 def fromstring(cls, string):
172 domain = Domain.fromstring(string)
173 if not isinstance(domain, Polyhedron):
174 raise ValueError('non-polyhedral expression: {!r}'.format(string))
175 return domain
176
177 def __repr__(self):
178 if self.isempty():
179 return 'Empty'
180 elif self.isuniverse():
181 return 'Universe'
182 else:
183 strings = []
184 for equality in self.equalities:
185 strings.append('0 == {}'.format(equality))
186 for inequality in self.inequalities:
187 strings.append('0 <= {}'.format(inequality))
188 if len(strings) == 1:
189 return strings[0]
190 else:
191 return 'And({})'.format(', '.join(strings))
192
193 @classmethod
194 def fromsympy(cls, expr):
195 domain = Domain.fromsympy(expr)
196 if not isinstance(domain, Polyhedron):
197 raise ValueError('non-polyhedral expression: {!r}'.format(expr))
198 return domain
199
200 def tosympy(self):
201 import sympy
202 constraints = []
203 for equality in self.equalities:
204 constraints.append(sympy.Eq(equality.tosympy(), 0))
205 for inequality in self.inequalities:
206 constraints.append(sympy.Ge(inequality.tosympy(), 0))
207 return sympy.And(*constraints)
208
209 @classmethod
210 def _polygon_inner_point(cls, points):
211 symbols = points[0].symbols
212 coordinates = {symbol: 0 for symbol in symbols}
213 for point in points:
214 for symbol, coordinate in point.coordinates():
215 coordinates[symbol] += coordinate
216 for symbol in symbols:
217 coordinates[symbol] /= len(points)
218 return Point(coordinates)
219
220 @classmethod
221 def _sort_polygon_2d(cls, points):
222 if len(points) <= 3:
223 return points
224 o = cls._polygon_inner_point(points)
225 angles = {}
226 for m in points:
227 om = Vector(o, m)
228 dx, dy = (coordinate for symbol, coordinate in om.coordinates())
229 angle = math.atan2(dy, dx)
230 angles[m] = angle
231 return sorted(points, key=angles.get)
232
233 @classmethod
234 def _sort_polygon_3d(cls, points):
235 if len(points) <= 3:
236 return points
237 o = cls._polygon_inner_point(points)
238 a = points[0]
239 oa = Vector(o, a)
240 norm_oa = oa.norm()
241 for b in points[1:]:
242 ob = Vector(o, b)
243 u = oa.cross(ob)
244 if not u.isnull():
245 u = u.asunit()
246 break
247 else:
248 raise ValueError('degenerate polygon')
249 angles = {a: 0.}
250 for m in points[1:]:
251 om = Vector(o, m)
252 normprod = norm_oa * om.norm()
253 cosinus = oa.dot(om) / normprod
254 sinus = u.dot(oa.cross(om)) / normprod
255 angle = math.acos(cosinus)
256 angle = math.copysign(angle, sinus)
257 angles[m] = angle
258 return sorted(points, key=angles.get)
259
260 def faces(self):
261 vertices = self.vertices()
262 faces = []
263 for constraint in self.constraints:
264 face = []
265 for vertex in vertices:
266 if constraint.subs(vertex.coordinates()) == 0:
267 face.append(vertex)
268 faces.append(face)
269 return faces
270
271 def plot(self):
272 import matplotlib.pyplot as plt
273 from matplotlib.path import Path
274 import matplotlib.patches as patches
275
276 if len(self.symbols)> 3:
277 raise TypeError
278
279 elif len(self.symbols) == 2:
280 verts = self.vertices()
281 points = []
282 codes = [Path.MOVETO]
283 for vert in verts:
284 pairs = ()
285 for sym in sorted(vert, key=Symbol.sortkey):
286 num = vert.get(sym)
287 pairs = pairs + (num,)
288 points.append(pairs)
289 points.append((0.0, 0.0))
290 num = len(points)
291 while num > 2:
292 codes.append(Path.LINETO)
293 num = num - 1
294 else:
295 codes.append(Path.CLOSEPOLY)
296 path = Path(points, codes)
297 fig = plt.figure()
298 ax = fig.add_subplot(111)
299 patch = patches.PathPatch(path, facecolor='blue', lw=2)
300 ax.add_patch(patch)
301 ax.set_xlim(-5,5)
302 ax.set_ylim(-5,5)
303 plt.show()
304
305 elif len(self.symbols)==3:
306 return 0
307
308 return points
309
310
311 def _polymorphic(func):
312 @functools.wraps(func)
313 def wrapper(left, right):
314 if not isinstance(left, Expression):
315 if isinstance(left, numbers.Rational):
316 left = Rational(left)
317 else:
318 raise TypeError('left must be a a rational number '
319 'or a linear expression')
320 if not isinstance(right, Expression):
321 if isinstance(right, numbers.Rational):
322 right = Rational(right)
323 else:
324 raise TypeError('right must be a a rational number '
325 'or a linear expression')
326 return func(left, right)
327 return wrapper
328
329 @_polymorphic
330 def Lt(left, right):
331 return Polyhedron([], [right - left - 1])
332
333 @_polymorphic
334 def Le(left, right):
335 return Polyhedron([], [right - left])
336
337 @_polymorphic
338 def Eq(left, right):
339 return Polyhedron([left - right], [])
340
341 @_polymorphic
342 def Ne(left, right):
343 return ~Eq(left, right)
344
345 @_polymorphic
346 def Gt(left, right):
347 return Polyhedron([], [left - right - 1])
348
349 @_polymorphic
350 def Ge(left, right):
351 return Polyhedron([], [left - right])
352
353
354 Empty = Eq(1, 0)
355
356 Universe = Polyhedron([])