#!/usr/bin/env python3
from pypol import *
def affine_derivative_closure(T, x0s, xs):
dxs = [x0.asdummy() for x0 in x0s]
k = Dummy('k')
for x in T.symbols:
assert x in x0s + xs
T0 = T
T1 = T0
for i, x0 in enumerate(x0s):
x, dx = xs[i], dxs[i]
T1 &= Eq(dx, x - x0)
T2 = T1.project_out(T0.symbols)
T3_eqs = []
T3_ins = []
for T2_eq in T2.equalities:
c = T2_eq.constant
T3_eq = T2_eq + (k - 1) * c
T3_eqs.append(T3_eq)
for T2_in in T2.inequalities:
c = T2_in.constant
T3_in = T2_in + (k - 1) * c
T3_ins.append(T3_in)
T3 = Polyhedron(T3_eqs, T3_ins)
T3 &= Ge(k, 0)
T4 = T3.project_out([k])
for i, dx in enumerate(dxs):
x0, x = x0s[i], xs[i]
T4 &= Eq(dx, x - x0)
T4 = T4.project_out(dxs)
return T4
i0, j0, i, j = symbols(['i', 'j', "i'", "j'"])
T = Eq(i, i0 + 2) & Eq(j, j0 + 1)
print('T =', T)
Tstar = affine_derivative_closure(T, [i0, j0], [i, j])
print('T* =', Tstar)