Working test suite
[linpy.git] / pypol / linear.py
1 import ctypes, ctypes.util
2 import functools
3 import numbers
4
5 from fractions import Fraction, gcd
6
7 from . import isl
8 from .isl import libisl
9
10
11 __all__ = [
12 'Expression', 'Constant', 'Symbol', 'symbols',
13 'eq', 'le', 'lt', 'ge', 'gt',
14 'Polyhedron',
15 'empty', 'universe'
16 ]
17
18
19 def _polymorphic_method(func):
20 @functools.wraps(func)
21 def wrapper(a, b):
22 if isinstance(b, Expression):
23 return func(a, b)
24 if isinstance(b, numbers.Rational):
25 b = Constant(b)
26 return func(a, b)
27 return NotImplemented
28 return wrapper
29
30 def _polymorphic_operator(func):
31 # A polymorphic operator should call a polymorphic method, hence we just
32 # have to test the left operand.
33 @functools.wraps(func)
34 def wrapper(a, b):
35 if isinstance(a, numbers.Rational):
36 a = Constant(a)
37 return func(a, b)
38 elif isinstance(a, Expression):
39 return func(a, b)
40 raise TypeError('arguments must be linear expressions')
41 return wrapper
42
43
44 _main_ctx = isl.Context()
45
46
47 class Expression:
48 """
49 This class implements linear expressions.
50 """
51
52 def __new__(cls, coefficients=None, constant=0):
53 if isinstance(coefficients, str):
54 if constant:
55 raise TypeError('too many arguments')
56 return cls.fromstring(coefficients)
57 if isinstance(coefficients, dict):
58 coefficients = coefficients.items()
59 if coefficients is None:
60 return Constant(constant)
61 coefficients = [(symbol, coefficient)
62 for symbol, coefficient in coefficients if coefficient != 0]
63 if len(coefficients) == 0:
64 return Constant(constant)
65 elif len(coefficients) == 1 and constant == 0:
66 symbol, coefficient = coefficients[0]
67 if coefficient == 1:
68 return Symbol(symbol)
69 self = object().__new__(cls)
70 self._coefficients = {}
71 for symbol, coefficient in coefficients:
72 if isinstance(symbol, Symbol):
73 symbol = str(symbol)
74 elif not isinstance(symbol, str):
75 raise TypeError('symbols must be strings or Symbol instances')
76 if isinstance(coefficient, Constant):
77 coefficient = coefficient.constant
78 if not isinstance(coefficient, numbers.Rational):
79 raise TypeError('coefficients must be rational numbers or Constant instances')
80 self._coefficients[symbol] = coefficient
81 if isinstance(constant, Constant):
82 constant = constant.constant
83 if not isinstance(constant, numbers.Rational):
84 raise TypeError('constant must be a rational number or a Constant instance')
85 self._constant = constant
86 self._symbols = tuple(sorted(self._coefficients))
87 self._dimension = len(self._symbols)
88 return self
89
90 @classmethod
91 def fromstring(cls, string):
92 raise NotImplementedError
93
94 @property
95 def symbols(self):
96 return self._symbols
97
98 @property
99 def dimension(self):
100 return self._dimension
101
102 def coefficient(self, symbol):
103 if isinstance(symbol, Symbol):
104 symbol = str(symbol)
105 elif not isinstance(symbol, str):
106 raise TypeError('symbol must be a string or a Symbol instance')
107 try:
108 return self._coefficients[symbol]
109 except KeyError:
110 return 0
111
112 __getitem__ = coefficient
113
114 def coefficients(self):
115 for symbol in self.symbols:
116 yield symbol, self.coefficient(symbol)
117
118 @property
119 def constant(self):
120 return self._constant
121
122 def isconstant(self):
123 return False
124
125 def values(self):
126 for symbol in self.symbols:
127 yield self.coefficient(symbol)
128 yield self.constant
129
130 @property
131 def symbol(self):
132 raise ValueError('not a symbol: {}'.format(self))
133
134 def issymbol(self):
135 return False
136
137 def __bool__(self):
138 return True
139
140 def __pos__(self):
141 return self
142
143 def __neg__(self):
144 return self * -1
145
146 @_polymorphic_method
147 def __add__(self, other):
148 coefficients = dict(self.coefficients())
149 for symbol, coefficient in other.coefficients():
150 if symbol in coefficients:
151 coefficients[symbol] += coefficient
152 else:
153 coefficients[symbol] = coefficient
154 constant = self.constant + other.constant
155 return Expression(coefficients, constant)
156
157 __radd__ = __add__
158
159 @_polymorphic_method
160 def __sub__(self, other):
161 coefficients = dict(self.coefficients())
162 for symbol, coefficient in other.coefficients():
163 if symbol in coefficients:
164 coefficients[symbol] -= coefficient
165 else:
166 coefficients[symbol] = -coefficient
167 constant = self.constant - other.constant
168 return Expression(coefficients, constant)
169
170 def __rsub__(self, other):
171 return -(self - other)
172
173 @_polymorphic_method
174 def __mul__(self, other):
175 if other.isconstant():
176 coefficients = dict(self.coefficients())
177 for symbol in coefficients:
178 coefficients[symbol] *= other.constant
179 constant = self.constant * other.constant
180 return Expression(coefficients, constant)
181 if isinstance(other, Expression) and not self.isconstant():
182 raise ValueError('non-linear expression: '
183 '{} * {}'.format(self._parenstr(), other._parenstr()))
184 return NotImplemented
185
186 __rmul__ = __mul__
187
188 @_polymorphic_method
189 def __truediv__(self, other):
190 if other.isconstant():
191 coefficients = dict(self.coefficients())
192 for symbol in coefficients:
193 coefficients[symbol] = \
194 Fraction(coefficients[symbol], other.constant)
195 constant = Fraction(self.constant, other.constant)
196 return Expression(coefficients, constant)
197 if isinstance(other, Expression):
198 raise ValueError('non-linear expression: '
199 '{} / {}'.format(self._parenstr(), other._parenstr()))
200 return NotImplemented
201
202 def __rtruediv__(self, other):
203 if isinstance(other, self):
204 if self.isconstant():
205 constant = Fraction(other, self.constant)
206 return Expression(constant=constant)
207 else:
208 raise ValueError('non-linear expression: '
209 '{} / {}'.format(other._parenstr(), self._parenstr()))
210 return NotImplemented
211
212 def __str__(self):
213 string = ''
214 i = 0
215 for symbol in self.symbols:
216 coefficient = self.coefficient(symbol)
217 if coefficient == 1:
218 if i == 0:
219 string += symbol
220 else:
221 string += ' + {}'.format(symbol)
222 elif coefficient == -1:
223 if i == 0:
224 string += '-{}'.format(symbol)
225 else:
226 string += ' - {}'.format(symbol)
227 else:
228 if i == 0:
229 string += '{}*{}'.format(coefficient, symbol)
230 elif coefficient > 0:
231 string += ' + {}*{}'.format(coefficient, symbol)
232 else:
233 assert coefficient < 0
234 coefficient *= -1
235 string += ' - {}*{}'.format(coefficient, symbol)
236 i += 1
237 constant = self.constant
238 if constant != 0 and i == 0:
239 string += '{}'.format(constant)
240 elif constant > 0:
241 string += ' + {}'.format(constant)
242 elif constant < 0:
243 constant *= -1
244 string += ' - {}'.format(constant)
245 if string == '':
246 string = '0'
247 return string
248
249 def _parenstr(self, always=False):
250 string = str(self)
251 if not always and (self.isconstant() or self.issymbol()):
252 return string
253 else:
254 return '({})'.format(string)
255
256 def __repr__(self):
257 string = '{}({{'.format(self.__class__.__name__)
258 for i, (symbol, coefficient) in enumerate(self.coefficients()):
259 if i != 0:
260 string += ', '
261 string += '{!r}: {!r}'.format(symbol, coefficient)
262 string += '}}, {!r})'.format(self.constant)
263 return string
264
265 @_polymorphic_method
266 def __eq__(self, other):
267 # "normal" equality
268 # see http://docs.sympy.org/dev/tutorial/gotchas.html#equals-signs
269 return isinstance(other, Expression) and \
270 self._coefficients == other._coefficients and \
271 self.constant == other.constant
272
273 def __hash__(self):
274 return hash((tuple(sorted(self._coefficients.items())), self._constant))
275
276 def _toint(self):
277 lcm = functools.reduce(lambda a, b: a*b // gcd(a, b),
278 [value.denominator for value in self.values()])
279 return self * lcm
280
281 @_polymorphic_method
282 def _eq(self, other):
283 return Polyhedron(equalities=[(self - other)._toint()])
284
285 @_polymorphic_method
286 def __le__(self, other):
287 return Polyhedron(inequalities=[(other - self)._toint()])
288
289 @_polymorphic_method
290 def __lt__(self, other):
291 return Polyhedron(inequalities=[(other - self)._toint() - 1])
292
293 @_polymorphic_method
294 def __ge__(self, other):
295 return Polyhedron(inequalities=[(self - other)._toint()])
296
297 @_polymorphic_method
298 def __gt__(self, other):
299 return Polyhedron(inequalities=[(self - other)._toint() - 1])
300
301
302 class Constant(Expression):
303
304 def __new__(cls, numerator=0, denominator=None):
305 self = object().__new__(cls)
306 if denominator is None:
307 if isinstance(numerator, numbers.Rational):
308 self._constant = numerator
309 elif isinstance(numerator, Constant):
310 self._constant = numerator.constant
311 else:
312 raise TypeError('constant must be a rational number or a Constant instance')
313 else:
314 self._constant = Fraction(numerator, denominator)
315 self._coefficients = {}
316 self._symbols = ()
317 self._dimension = 0
318 return self
319
320 def isconstant(self):
321 return True
322
323 def __bool__(self):
324 return bool(self.constant)
325
326 def __repr__(self):
327 return '{}({!r})'.format(self.__class__.__name__, self._constant)
328
329
330 class Symbol(Expression):
331
332 def __new__(cls, name):
333 if isinstance(name, Symbol):
334 name = name.symbol
335 elif not isinstance(name, str):
336 raise TypeError('name must be a string or a Symbol instance')
337 self = object().__new__(cls)
338 self._coefficients = {name: 1}
339 self._constant = 0
340 self._symbols = tuple(name)
341 self._symbol = name
342 self._dimension = 1
343 return self
344
345 @property
346 def symbol(self):
347 return self._symbol
348
349 def issymbol(self):
350 return True
351
352 def __repr__(self):
353 return '{}({!r})'.format(self.__class__.__name__, self._symbol)
354
355 def symbols(names):
356 if isinstance(names, str):
357 names = names.replace(',', ' ').split()
358 return (Symbol(name) for name in names)
359
360
361 @_polymorphic_operator
362 def eq(a, b):
363 return a._eq(b)
364
365 @_polymorphic_operator
366 def le(a, b):
367 return a <= b
368
369 @_polymorphic_operator
370 def lt(a, b):
371 return a < b
372
373 @_polymorphic_operator
374 def ge(a, b):
375 return a >= b
376
377 @_polymorphic_operator
378 def gt(a, b):
379 return a > b
380
381
382 class Polyhedron:
383 """
384 This class implements polyhedrons.
385 """
386
387 def __new__(cls, equalities=None, inequalities=None):
388 if isinstance(equalities, str):
389 if inequalities is not None:
390 raise TypeError('too many arguments')
391 return cls.fromstring(equalities)
392 self = super().__new__(cls)
393 self._equalities = []
394 if equalities is not None:
395 for constraint in equalities:
396 for value in constraint.values():
397 if value.denominator != 1:
398 raise TypeError('non-integer constraint: '
399 '{} == 0'.format(constraint))
400 self._equalities.append(constraint)
401 self._equalities = tuple(self._equalities)
402 self._inequalities = []
403 if inequalities is not None:
404 for constraint in inequalities:
405 for value in constraint.values():
406 if value.denominator != 1:
407 raise TypeError('non-integer constraint: '
408 '{} <= 0'.format(constraint))
409 self._inequalities.append(constraint)
410 self._inequalities = tuple(self._inequalities)
411 self._constraints = self._equalities + self._inequalities
412 self._symbols = set()
413 for constraint in self._constraints:
414 self.symbols.update(constraint.symbols)
415 self._symbols = tuple(sorted(self._symbols))
416 return self
417
418 @classmethod
419 def fromstring(cls, string):
420 raise NotImplementedError
421
422 @property
423 def equalities(self):
424 return self._equalities
425
426 @property
427 def inequalities(self):
428 return self._inequalities
429
430 @property
431 def constraints(self):
432 return self._constraints
433
434 @property
435 def symbols(self):
436 return self._symbols
437
438 @property
439 def dimension(self):
440 return len(self.symbols)
441
442 def __bool__(self):
443 return not self.is_empty()
444
445 def __contains__(self, value):
446 # is the value in the polyhedron?
447 raise NotImplementedError
448
449 def __eq__(self, other):
450 raise NotImplementedError
451
452 def isempty(self):
453 bset = self._toisl()
454 return bool(libisl.isl_basic_set_is_empty(bset))
455
456 def isuniverse(self):
457 raise NotImplementedError
458
459 def isdisjoint(self, other):
460 # return true if the polyhedron has no elements in common with other
461 raise NotImplementedError
462
463 def issubset(self, other):
464 raise NotImplementedError
465
466 def __le__(self, other):
467 return self.issubset(other)
468
469 def __lt__(self, other):
470 raise NotImplementedError
471
472 def issuperset(self, other):
473 # test whether every element in other is in the polyhedron
474 raise NotImplementedError
475
476 def __ge__(self, other):
477 return self.issuperset(other)
478
479 def __gt__(self, other):
480 raise NotImplementedError
481
482 def union(self, *others):
483 # return a new polyhedron with elements from the polyhedron and all
484 # others (convex union)
485 raise NotImplementedError
486
487 def __or__(self, other):
488 return self.union(other)
489
490 def intersection(self, *others):
491 # return a new polyhedron with elements common to the polyhedron and all
492 # others
493 # a poor man's implementation could be:
494 # equalities = list(self.equalities)
495 # inequalities = list(self.inequalities)
496 # for other in others:
497 # equalities.extend(other.equalities)
498 # inequalities.extend(other.inequalities)
499 # return self.__class__(equalities, inequalities)
500 raise NotImplementedError
501
502 def __and__(self, other):
503 return self.intersection(other)
504
505 def difference(self, *others):
506 # return a new polyhedron with elements in the polyhedron that are not
507 # in the others
508 raise NotImplementedError
509
510 def __sub__(self, other):
511 return self.difference(other)
512
513 def __str__(self):
514 constraints = []
515 for constraint in self.equalities:
516 constraints.append('{} == 0'.format(constraint))
517 for constraint in self.inequalities:
518 constraints.append('{} >= 0'.format(constraint))
519 return '{{{}}}'.format(', '.join(constraints))
520
521 def __repr__(self):
522 equalities = list(self.equalities)
523 inequalities = list(self.inequalities)
524 return '{}(equalities={!r}, inequalities={!r})' \
525 ''.format(self.__class__.__name__, equalities, inequalities)
526
527 def _symbolunion(self, *others):
528 symbols = set(self.symbols)
529 for other in others:
530 symbols.update(other.symbols)
531 return sorted(symbols)
532
533 def _toisl(self, symbols=None):
534 if symbols is None:
535 symbols = self.symbols
536 num_coefficients = len(symbols)
537 space = libisl.isl_space_set_alloc(_main_ctx, 0, num_coefficients)
538 bset = libisl.isl_basic_set_universe(libisl.isl_space_copy(space))
539 ls = libisl.isl_local_space_from_space(space)
540 ceq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(ls))
541 cin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(ls))
542 '''if there are equalities/inequalities, take each constant and coefficient and add as a constraint to the basic set'''
543 if list(self.equalities): #check if any equalities exist
544 for eq in self.equalities:
545 coeff_eq = dict(eq.coefficients())
546 if eq.constant:
547 value = eq.constant
548 ceq = libisl.isl_constraint_set_constant_si(ceq, value)
549 for eq in coeff_eq:
550 num = coeff_eq.get(eq)
551 iden = symbols.index(eq)
552 ceq = libisl.isl_constraint_set_coefficient_si(ceq, libisl.isl_dim_set, iden, num) #use 3 for type isl_dim_set
553 bset = libisl.isl_basic_set_add_constraint(bset, ceq)
554 if list(self.inequalities): #check if any inequalities exist
555 for ineq in self.inequalities:
556 coeff_in = dict(ineq.coefficients())
557 if ineq.constant:
558 value = ineq.constant
559 cin = libisl.isl_constraint_set_constant_si(cin, value)
560 for ineq in coeff_in:
561 num = coeff_in.get(ineq)
562 iden = symbols.index(ineq)
563 cin = libisl.isl_constraint_set_coefficient_si(cin, libisl.isl_dim_set, iden, num) #use 3 for type isl_dim_set
564 bset = libisl.isl_basic_set_add_constraint(bset, cin)
565 bset = isl.BasicSet(bset)
566 return bset
567
568 @classmethod
569 def _fromisl(cls, bset):
570 raise NotImplementedError
571 equalities = ...
572 inequalities = ...
573 return cls(equalities, inequalities)
574 '''takes basic set in isl form and puts back into python version of polyhedron
575 isl example code gives isl form as:
576 "{[i] : exists (a : i = 2a and i >= 10 and i <= 42)}")
577 our printer is giving form as:
578 b'{ [i0] : 1 = 0 }' '''
579 #bset = self
580 # if self._equalities:
581 # constraints = libisl.isl_basic_set_equalities_matrix(bset, 3)
582 # elif self._inequalities:
583 # constraints = libisl.isl_basic_set_inequalities_matrix(bset, 3)
584 # print(constraints)
585 # return constraints
586
587 empty = None #eq(0,1)
588 universe = None #Polyhedron()
589
590
591 if __name__ == '__main__':
592 ex1 = Expression(coefficients={'a': 1, 'x': 2}, constant=2)
593 ex2 = Expression(coefficients={'a': 3 , 'b': 2}, constant=3)
594 p = Polyhedron(inequalities=[ex1, ex2])
595 bs = p._toisl()
596 print(bs)
597 print('empty ?', p.isempty())
598 print('empty ?', eq(0, 1).isempty())