Update drop_dims, now works with more than 2 dims passed
[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 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
17
18 T0 = T
19
20 T1 = T0
21 for i, x0 in enumerate(x0s):
22 x, dx = xs[i], dxs[i]
23 T1 &= Eq(dx, x - x0)
24
25 T2 = T1.project_out(T0.symbols)
26
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)
39
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)
45
46 return T4
47
48 i0, j0, i, j = symbols(['i', 'j', "i'", "j'"])
49 T = Eq(i, i0 + 2) & Eq(j, j0 + 1)
50
51 print('T =', T)
52 Tstar = affine_derivative_closure(T, [i0, j0])
53 print('T* =', Tstar)