e2cacc7aa07461a9b2612ad1cb62481bb4502d1e
[linpy.git] / examples / nsad2010.py
1 #!/usr/bin/env python3
2
3 from pypol import *
4
5 def affine_derivative_closure(T, x0s):
6
7 xs = [Symbol("{}'".format(x0.name)) for x0 in x0s]
8 dxs = [Symbol('d{}'.format(x0.name)) for x0 in x0s]
9 k = Symbol('k')
10
11 for x in T.symbols:
12 assert x in x0s + xs
13 for dx in dxs:
14 assert dx.name not in T.symbols
15 assert k.name not in T.symbols
16
17 T0 = T
18
19 T1 = T0
20 for i, x0 in enumerate(x0s):
21 x, dx = xs[i], dxs[i]
22 T1 &= Eq(dx, x - x0)
23
24 T2 = T1.project_out(T0.symbols)
25
26 T3_eqs = []
27 T3_ins = []
28 for T2_eq in T2.equalities:
29 c = T2_eq.constant
30 T3_eq = T2_eq + (k - 1) * c
31 T3_eqs.append(T3_eq)
32 for T2_in in T2.inequalities:
33 c = T2_in.constant
34 T3_in = T2_in + (k - 1) * c
35 T3_ins.append(T3_in)
36 T3 = Polyhedron(T3_eqs, T3_ins)
37 T3 &= Ge(k, 0)
38
39 T4 = T3.project_out([k])
40 for i, dx in enumerate(dxs):
41 x0, x = x0s[i], xs[i]
42 T4 &= Eq(dx, x - x0)
43 T4 = T4.project_out(dxs)
44
45 return T4
46
47 i0, j0, i, j = symbols(['i', 'j', "i'", "j'"])
48 T = Eq(i, i0 + 2) & Eq(j, j0 + 1)
49
50 print('T =', T)
51 Tstar = affine_derivative_closure(T, [i0, j0])
52 print('T* =', Tstar)