91a85b4208a1fd3b8903bc543b8373f4954993ad
[linpy.git] / examples / nsad2010.py
1 #!/usr/bin/env python3
2 #
3 # Copyright 2014 MINES ParisTech
4 #
5 # This file is part of LinPy.
6 #
7 # LinPy is free software: you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation, either version 3 of the License, or
10 # (at your option) any later version.
11 #
12 # LinPy is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with LinPy. If not, see <http://www.gnu.org/licenses/>.
19
20 from linpy import *
21
22
23 class Transformer:
24
25 def __new__(cls, polyhedron, range_symbols, domain_symbols):
26 self = object().__new__(cls)
27 self.polyhedron = polyhedron
28 self.range_symbols = range_symbols
29 self.domain_symbols = domain_symbols
30 return self
31
32 @property
33 def symbols(self):
34 return self.range_symbols + self.domain_symbols
35
36 def star(self):
37 delta_symbols = [symbol.asdummy() for symbol in self.range_symbols]
38 k = Dummy('k')
39 polyhedron = self.polyhedron
40 for x, xprime, dx in zip(self.range_symbols, self.domain_symbols, delta_symbols):
41 polyhedron &= Eq(dx, xprime - x)
42 polyhedron = polyhedron.project(self.symbols)
43 equalities, inequalities = [], []
44 for equality in polyhedron.equalities:
45 equality += (k-1) * equality.constant
46 equalities.append(equality)
47 for inequality in polyhedron.inequalities:
48 inequality += (k-1) * inequality.constant
49 inequalities.append(inequality)
50 polyhedron = Polyhedron(equalities, inequalities) & Ge(k, 0)
51 polyhedron = polyhedron.project([k])
52 for x, xprime, dx in zip(self.range_symbols, self.domain_symbols, delta_symbols):
53 polyhedron &= Eq(dx, xprime - x)
54 polyhedron = polyhedron.project(delta_symbols)
55 return Transformer(polyhedron, self.range_symbols, self.domain_symbols)
56
57
58 if __name__ == '__main__':
59 i, iprime, j, jprime = symbols("i i' j j'")
60 transformer = Transformer(Eq(iprime, i + 2) & Eq(jprime, j + 1),
61 [i, j], [iprime, jprime])
62 print('T =', transformer.polyhedron)
63 print('T* =', transformer.star().polyhedron)