1 #!/usr/bin/env python3
3 from pypol import *
5 def affine_derivative_closure(T, x0s):
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')
11 for x in T.symbols:
12 x = Symbol(x)
13 assert x in x0s + xs
14 for dx in dxs:
15 assert dx.name not in T.symbols
16 assert k.name not in T.symbols
18 T0 = T
20 T1 = T0
21 for i, x0 in enumerate(x0s):
22 x, dx = xs[i], dxs[i]
23 T1 &= Eq(dx, x - x0)
25 T2 = T1.project_out(T0.symbols)
27 T3_eqs = []
28 T3_ins = []
29 for T2_eq in T2.equalities:
30 c = T2_eq.constant
31 T3_eq = T2_eq + (k - 1) * c
32 T3_eqs.append(T3_eq)
33 for T2_in in T2.inequalities:
34 c = T2_in.constant
35 T3_in = T2_in + (k - 1) * c
36 T3_ins.append(T3_in)
37 T3 = Polyhedron(T3_eqs, T3_ins)
38 T3 &= Ge(k, 0)
40 T4 = T3.project_out([k])
41 for i, dx in enumerate(dxs):
42 x0, x = x0s[i], xs[i]
43 T4 &= Eq(dx, x - x0)
44 T4 = T4.project_out(dxs)
46 return T4
48 i0, j0, i, j = symbols(['i', 'j', "i'", "j'"])
49 T = Eq(i, i0 + 2) & Eq(j, j0 + 1)
51 print('T =', T)
52 Tstar = affine_derivative_closure(T, [i0, j0])
53 print('T* =', Tstar)