Implement method Polyhedron.faces()
[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
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 _sort_polygon_2d(cls, points):
211 if len(points) <= 3:
212 return points
213 o = sum((Vector(point) for point in points)) / len(points)
214 o = Point(o.coordinates())
215 angles = {}
216 for m in points:
217 om = Vector(o, m)
218 dx, dy = (coordinate for symbol, coordinates in om.coordinates())
219 angle = math.atan2(dy, dx)
220 angles[m] = angle
221 return sorted(points, key=angles.get)
222
223 @classmethod
224 def _sort_polygon_3d(cls, points):
225 if len(points) <= 3:
226 return points
227 o = sum((Vector(point) for point in points)) / len(points)
228 o = Point(o.coordinates())
229 a, b = points[:2]
230 oa = Vector(o, a)
231 ob = Vector(o, b)
232 norm_oa = oa.norm()
233 u = (oa.cross(ob)).asunit()
234 angles = {a: 0.}
235 for m in points[1:]:
236 om = Vector(o, m)
237 normprod = norm_oa * om.norm()
238 cosinus = oa.dot(om) / normprod
239 sinus = u.dot(oa.cross(om)) / normprod
240 angle = math.acos(cosinus)
241 angle = math.copysign(angle, sinus)
242 angles[m] = angle
243 return sorted(points, key=angles.get)
244
245 def faces(self):
246 vertices = self.vertices()
247 faces = []
248 for constraint in self.constraints:
249 face = []
250 for vertex in vertices:
251 if constraint.subs(vertex.coordinates()) == 0:
252 face.append(vertex)
253 faces.append(face)
254 return faces
255
256 def plot(self):
257 import matplotlib.pyplot as plt
258 from matplotlib.path import Path
259 import matplotlib.patches as patches
260
261 if len(self.symbols)> 3:
262 raise TypeError
263
264 elif len(self.symbols) == 2:
265 verts = self.vertices()
266 points = []
267 codes = [Path.MOVETO]
268 for vert in verts:
269 pairs = ()
270 for sym in sorted(vert, key=Symbol.sortkey):
271 num = vert.get(sym)
272 pairs = pairs + (num,)
273 points.append(pairs)
274 points.append((0.0, 0.0))
275 num = len(points)
276 while num > 2:
277 codes.append(Path.LINETO)
278 num = num - 1
279 else:
280 codes.append(Path.CLOSEPOLY)
281 path = Path(points, codes)
282 fig = plt.figure()
283 ax = fig.add_subplot(111)
284 patch = patches.PathPatch(path, facecolor='blue', lw=2)
285 ax.add_patch(patch)
286 ax.set_xlim(-5,5)
287 ax.set_ylim(-5,5)
288 plt.show()
289
290 elif len(self.symbols)==3:
291 return 0
292
293 return points
294
295
296 def _polymorphic(func):
297 @functools.wraps(func)
298 def wrapper(left, right):
299 if isinstance(left, numbers.Rational):
300 left = Rational(left)
301 elif not isinstance(left, Expression):
302 raise TypeError('left must be a a rational number '
303 'or a linear expression')
304 if isinstance(right, numbers.Rational):
305 right = Rational(right)
306 elif not isinstance(right, Expression):
307 raise TypeError('right must be a a rational number '
308 'or a linear expression')
309 return func(left, right)
310 return wrapper
311
312 @_polymorphic
313 def Lt(left, right):
314 return Polyhedron([], [right - left - 1])
315
316 @_polymorphic
317 def Le(left, right):
318 return Polyhedron([], [right - left])
319
320 @_polymorphic
321 def Eq(left, right):
322 return Polyhedron([left - right], [])
323
324 @_polymorphic
325 def Ne(left, right):
326 return ~Eq(left, right)
327
328 @_polymorphic
329 def Gt(left, right):
330 return Polyhedron([], [left - right - 1])
331
332 @_polymorphic
333 def Ge(left, right):
334 return Polyhedron([], [left - right])
335
336
337 Empty = Eq(1, 0)
338
339 Universe = Polyhedron([])