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([])