New example: basic implementation of ACI'10

[linpy.git] / examples / nsad2010.py

diff --git a/examples/nsad2010.py b/examples/nsad2010.py
new file mode 100755 (executable)
index 0000000..14c984c
--- /dev/null
+++ b/examples/nsad2010.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python3
+
+from pypol import *
+
+def affine_derivative_closure(T, x0s):
+
+    xs = [Symbol("{}'".format(x0.name)) for x0 in x0s]
+    dxs = [Symbol('d{}'.format(x0.name)) for x0 in x0s]
+    k = Symbol('k')
+
+    for x in T.symbols:
+        x = Symbol(x)
+        assert x in x0s + xs
+    for dx in dxs:
+        assert dx.name not in T.symbols
+    assert k.name not in T.symbols
+
+    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])
+print('T* =', Tstar)