pypol/polyhedra.py

   1 import functools
   2 import numbers
   3
   4 from . import islhelper
   5
   6 from .islhelper import mainctx, libisl
   7 from .linexprs import Expression, Rational
   8 from .domains import Domain
   9
  10
  11 __all__ = [
  12     'Polyhedron',
  13     'Lt', 'Le', 'Eq', 'Ne', 'Ge', 'Gt',
  14     'Empty', 'Universe',
  15 ]
  16
  17
  18 class Polyhedron(Domain):
  19
  20     __slots__ = (
  21         '_equalities',
  22         '_inequalities',
  23         '_constraints',
  24         '_symbols',
  25         '_dimension',
  26     )
  27
  28     def __new__(cls, equalities=None, inequalities=None):
  29         if isinstance(equalities, str):
  30             if inequalities is not None:
  31                 raise TypeError('too many arguments')
  32             return cls.fromstring(equalities)
  33         elif isinstance(equalities, Polyhedron):
  34             if inequalities is not None:
  35                 raise TypeError('too many arguments')
  36             return equalities
  37         elif isinstance(equalities, Domain):
  38             if inequalities is not None:
  39                 raise TypeError('too many arguments')
  40             return equalities.polyhedral_hull()
  41         if equalities is None:
  42             equalities = []
  43         else:
  44             for i, equality in enumerate(equalities):
  45                 if not isinstance(equality, Expression):
  46                     raise TypeError('equalities must be linear expressions')
  47                 equalities[i] = equality.scaleint()
  48         if inequalities is None:
  49             inequalities = []
  50         else:
  51             for i, inequality in enumerate(inequalities):
  52                 if not isinstance(inequality, Expression):
  53                     raise TypeError('inequalities must be linear expressions')
  54                 inequalities[i] = inequality.scaleint()
  55         symbols = cls._xsymbols(equalities + inequalities)
  56         islbset = cls._toislbasicset(equalities, inequalities, symbols)
  57         return cls._fromislbasicset(islbset, symbols)
  58
  59     @property
  60     def equalities(self):
  61         return self._equalities
  62
  63     @property
  64     def inequalities(self):
  65         return self._inequalities
  66
  67     @property
  68     def constraints(self):
  69         return self._constraints
  70
  71     @property
  72     def polyhedra(self):
  73         return self,
  74
  75     def disjoint(self):
  76         return self
  77
  78     def isuniverse(self):
  79         islbset = self._toislbasicset(self.equalities, self.inequalities,
  80             self.symbols)
  81         universe = bool(libisl.isl_basic_set_is_universe(islbset))
  82         libisl.isl_basic_set_free(islbset)
  83         return universe
  84
  85     def polyhedral_hull(self):
  86         return self
  87
  88     @classmethod
  89     def _fromislbasicset(cls, islbset, symbols):
  90         islconstraints = islhelper.isl_basic_set_constraints(islbset)
  91         equalities = []
  92         inequalities = []
  93         for islconstraint in islconstraints:
  94             constant = libisl.isl_constraint_get_constant_val(islconstraint)
  95             constant = islhelper.isl_val_to_int(constant)
  96             coefficients = {}
  97             for index, symbol in enumerate(symbols):
  98                 coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint,
  99                     libisl.isl_dim_set, index)
 100                 coefficient = islhelper.isl_val_to_int(coefficient)
 101                 if coefficient != 0:
 102                     coefficients[symbol] = coefficient
 103             expression = Expression(coefficients, constant)
 104             if libisl.isl_constraint_is_equality(islconstraint):
 105                 equalities.append(expression)
 106             else:
 107                 inequalities.append(expression)
 108         libisl.isl_basic_set_free(islbset)
 109         self = object().__new__(Polyhedron)
 110         self._equalities = tuple(equalities)
 111         self._inequalities = tuple(inequalities)
 112         self._constraints = tuple(equalities + inequalities)
 113         self._symbols = cls._xsymbols(self._constraints)
 114         self._dimension = len(self._symbols)
 115         return self
 116
 117     @classmethod
 118     def _toislbasicset(cls, equalities, inequalities, symbols):
 119         dimension = len(symbols)
 120         indices = {symbol: index for index, symbol in enumerate(symbols)}
 121         islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
 122         islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
 123         islls = libisl.isl_local_space_from_space(islsp)
 124         for equality in equalities:
 125             isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls))
 126             for symbol, coefficient in equality.coefficients():
 127                 islval = str(coefficient).encode()
 128                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 129                 index = indices[symbol]
 130                 isleq = libisl.isl_constraint_set_coefficient_val(isleq,
 131                     libisl.isl_dim_set, index, islval)
 132             if equality.constant != 0:
 133                 islval = str(equality.constant).encode()
 134                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 135                 isleq = libisl.isl_constraint_set_constant_val(isleq, islval)
 136             islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
 137         for inequality in inequalities:
 138             islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls))
 139             for symbol, coefficient in inequality.coefficients():
 140                 islval = str(coefficient).encode()
 141                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 142                 index = indices[symbol]
 143                 islin = libisl.isl_constraint_set_coefficient_val(islin,
 144                     libisl.isl_dim_set, index, islval)
 145             if inequality.constant != 0:
 146                 islval = str(inequality.constant).encode()
 147                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 148                 islin = libisl.isl_constraint_set_constant_val(islin, islval)
 149             islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
 150         return islbset
 151
 152     @classmethod
 153     def fromstring(cls, string):
 154         domain = Domain.fromstring(string)
 155         if not isinstance(domain, Polyhedron):
 156             raise ValueError('non-polyhedral expression: {!r}'.format(string))
 157         return domain
 158
 159     def __repr__(self):
 160         if self.isempty():
 161             return 'Empty'
 162         elif self.isuniverse():
 163             return 'Universe'
 164         else:
 165             strings = []
 166             for equality in self.equalities:
 167                 strings.append('Eq({}, 0)'.format(equality))
 168             for inequality in self.inequalities:
 169                 strings.append('Ge({}, 0)'.format(inequality))
 170             if len(strings) == 1:
 171                 return strings[0]
 172             else:
 173                 return 'And({})'.format(', '.join(strings))
 174
 175     @classmethod
 176     def _fromsympy(cls, expr):
 177         import sympy
 178         equalities = []
 179         inequalities = []
 180         if expr.func == sympy.And:
 181             for arg in expr.args:
 182                 arg_eqs, arg_ins = cls._fromsympy(arg)
 183                 equalities.extend(arg_eqs)
 184                 inequalities.extend(arg_ins)
 185         elif expr.func == sympy.Eq:
 186             expr = Expression.fromsympy(expr.args[0] - expr.args[1])
 187             equalities.append(expr)
 188         else:
 189             if expr.func == sympy.Lt:
 190                 expr = Expression.fromsympy(expr.args[1] - expr.args[0] - 1)
 191             elif expr.func == sympy.Le:
 192                 expr = Expression.fromsympy(expr.args[1] - expr.args[0])
 193             elif expr.func == sympy.Ge:
 194                 expr = Expression.fromsympy(expr.args[0] - expr.args[1])
 195             elif expr.func == sympy.Gt:
 196                 expr = Expression.fromsympy(expr.args[0] - expr.args[1] - 1)
 197             else:
 198                 raise ValueError('non-polyhedral expression: {!r}'.format(expr))
 199             inequalities.append(expr)
 200         return equalities, inequalities
 201
 202     @classmethod
 203     def fromsympy(cls, expr):
 204         import sympy
 205         equalities, inequalities = cls._fromsympy(expr)
 206         return cls(equalities, inequalities)
 207
 208     def tosympy(self):
 209         import sympy
 210         constraints = []
 211         for equality in self.equalities:
 212             constraints.append(sympy.Eq(equality.tosympy(), 0))
 213         for inequality in self.inequalities:
 214             constraints.append(sympy.Ge(inequality.tosympy(), 0))
 215         return sympy.And(*constraints)
 216
 217
 218 def _polymorphic(func):
 219     @functools.wraps(func)
 220     def wrapper(left, right):
 221         if isinstance(left, numbers.Rational):
 222             left = Rational(left)
 223         elif not isinstance(left, Expression):
 224             raise TypeError('left must be a a rational number '
 225                 'or a linear expression')
 226         if isinstance(right, numbers.Rational):
 227             right = Rational(right)
 228         elif not isinstance(right, Expression):
 229             raise TypeError('right must be a a rational number '
 230                 'or a linear expression')
 231         return func(left, right)
 232     return wrapper
 233
 234 @_polymorphic
 235 def Lt(left, right):
 236     return Polyhedron([], [right - left - 1])
 237
 238 @_polymorphic
 239 def Le(left, right):
 240     return Polyhedron([], [right - left])
 241
 242 @_polymorphic
 243 def Eq(left, right):
 244     return Polyhedron([left - right], [])
 245
 246 @_polymorphic
 247 def Ne(left, right):
 248     return ~Eq(left, right)
 249
 250 @_polymorphic
 251 def Gt(left, right):
 252     return Polyhedron([], [left - right - 1])
 253
 254 @_polymorphic
 255 def Ge(left, right):
 256     return Polyhedron([], [left - right])
 257
 258
 259 Empty = Eq(1, 0)
 260
 261 Universe = Polyhedron([])