1#!/usr/bin/env python
2
3import polypy
4import polypy_test
5
6import sys
7import time
8
9polypy_test.init()
10
11[x, y, z] = [polypy.Variable(name) for name in ['x', 'y', 'z']]
12polypy.variable_order.set([z, y, x])
13
14# polypy.trace_enable("polynomial")
15# polypy.trace_enable("coefficient")
16# polypy.trace_enable("coefficient::sgn")
17# polypy.trace_enable("coefficient::roots")
18# polypy.trace_enable("value::pick");
19# polypy.trace_enable("value::cmp");
20# polypy.trace_enable("value::get_value_between")
21# polypy.trace_enable("feasibility_set")
22
23# All signs
24sgns = [polypy.SGN_LT_0,
25        polypy.SGN_LE_0,
26        polypy.SGN_EQ_0,
27        polypy.SGN_NE_0,
28        polypy.SGN_GT_0,
29        polypy.SGN_GE_0
30        ]
31
32sgn_name = {
33            polypy.SGN_LT_0 : "<  0",
34            polypy.SGN_LE_0 : "<= 0",
35            polypy.SGN_EQ_0 : "== 0",
36            polypy.SGN_NE_0 : "!= 0",
37            polypy.SGN_GT_0 : ">  0",
38            polypy.SGN_GE_0 : ">= 0"
39}
40
41polypy_test.start("Polynomial Feasibility Integer")
42
43assignment = polypy.Assignment()
44
45p1 = x**2 - 7
46S1 = p1.feasible_set(assignment, polypy.SGN_GE_0) # (-inf, -2.6], [2.6, +inf)
47p2 = x*(x - 4)
48S2 = p2.feasible_set(assignment, polypy.SGN_LT_0) # (0, 4)
49
50S = S1.intersect(S2)
51v = S.pick_value();
52polypy_test.check(v.to_double() == 3)
53
54polypy_test.start("Feasibility Intervals Intersection")
55
56assignment = polypy.Assignment()
57assignment.set_value(y, 0)
58assignment.set_value(z, 1)
59
60p1 = (x**2 - y)*(x**2 - 2*z)
61p2 = (x**2 - y)*(x**2 - 3*z)
62
63S1 = p1.feasible_set(assignment, polypy.SGN_NE_0);
64S2 = p2.feasible_set(assignment, polypy.SGN_LE_0);
65P = S1.intersect(S2);
66polypy_test.check(True)
67
68S2 = p2.feasible_set(assignment, polypy.SGN_LT_0);
69P = S1.intersect(S2);
70polypy_test.check(True)
71
72assignment = polypy.Assignment()
73assignment.set_value(y, 0)
74assignment.set_value(z, 1)
75
76p = (x - y)*(x - z)
77
78S1 = p.feasible_set(assignment, polypy.SGN_GE_0);
79S2 = p.feasible_set(assignment, polypy.SGN_EQ_0);
80
81P = S1.intersect(S2);
82polypy_test.check(True)
83
84S1 = p.feasible_set(assignment, polypy.SGN_GT_0);
85S2 = p.feasible_set(assignment, polypy.SGN_GE_0);
86
87P = S1.intersect(S2);
88polypy_test.check(True)
89
90polypy_test.start("Polynomial Feasibility Intervals")
91
92def check_feasible(p, var, assignment, expected):
93    ok = True
94    for sgn, expected in zip(sgns, expected):
95        p_feasible = p.feasible_intervals(assignment, sgn)
96        ok = len(p_feasible) == expected
97        if (not ok):
98            print("expected = {0}".format(expected))
99            print("p_feasible = {0}".format(p_feasible))
100        if ok:
101            for I in p_feasible:
102                v = I.pick_value()
103                ok = I.contains(v)
104                if (not ok):
105                    print("I = {0}".format(I))
106                    print("v = {0}".format(v))
107                    break
108        if (not ok):
109            print("p = {0}".format(p))
110            print("assignment = {0}".format(assignment))
111            print("sgn = {0}".format(sgn_name[sgn]))
112        polypy_test.check(ok)
113
114assignment = polypy.Assignment()
115assignment.set_value(y, polypy.AlgebraicNumber(x**2 - 2, 1))
116assignment.set_value(z, polypy.AlgebraicNumber(x**2 - 3, 1))
117
118# [-0.809841]
119# -1 0 1
120p_slow = ((2*z)*y**2)*x**3 + ((1*z**3)*y)*x**2 + ((1*z**2)*y**3)*x + (1*z + 4)
121p_expected = [1, 1, 1, 2, 1, 1]
122check_feasible(p_slow, x, assignment, p_expected)
123
124p_slow = (1*y**3 + (2*z**3)*y)*x**3 + ((2*z)*y**3)*x**2 + ((2*z**2)*y**2 + 3)
125
126# for random in range(1000):
127#     start = time.time()
128#     p = polypy_test.random_polynomial(3, 2, [x, y, z], 5)
129#     print(p)
130#     for sgn in sgns:
131#         p_feasible = p.feasible_intervals(assignment, sgn)
132#         # print(p_feasible)
133#     end = time.time()
134#     print(end - start)
135
136# + 0 - 0 +
137p = ((1*z)*y**2 + (1*z**2)*y)*x**2 + ((1*z**2)*y**2)*x + 1
138p_expected = [1, 1, 2, 3, 2, 2]
139check_feasible(p, x, assignment, p_expected)
140
141# [-0.686589047969039, 0.686589047969039]
142# + 0 - 0 +
143p = ((1*z**2)*y)*x**2 + (-1*y**2)
144p_expected = [1, 1, 2, 3, 2, 2]
145check_feasible(p, x, assignment, p_expected)
146
147assignment = polypy.Assignment()
148assignment.set_value(y, polypy.AlgebraicNumber(x**2 - 2, 1))
149assignment.set_value(z, polypy.AlgebraicNumber(x**2 - 2, 1))
150
151# + 0 +
152p = (y + z)*x**2 + (y - z)*x + (y - z)
153p_expected = [0, 1, 1, 2, 2, 1]
154check_feasible(p, x, assignment, p_expected)
155
156# coefficient vanish at sqrt(2)
157V = y*z -2
158
159# + 0 -
160p1 = V*x**2 - 2  *x +   (y+z)
161p1_expected = [1, 1, 1, 2, 1, 1]
162check_feasible(p1, x, assignment, p1_expected)
163
164# +
165p2 =   x**2 - 2*V*x +   (y+z)
166p2_expected = [0, 0, 0, 1, 1, 1]
167check_feasible(p2, x, assignment, p2_expected)
168
169# + 0 - 0 +
170p3 =   x**2 - 2  *x + V*(y+z)
171p3_expected = [1, 1, 2, 3, 2, 2]
172check_feasible(p3, x, assignment, p3_expected)
173
174# + 0 - 0 - 0 +
175p = (x**2 - y - z)*(x**2)
176expected = [2, 1, 3, 4, 2, 3]
177check_feasible(p, x, assignment, expected)
178
179assignment = polypy.Assignment()
180assignment.set_value(y, 1)
181assignment.set_value(z, 1)
182
183# + 0 - 0 - 0 +
184p = (x**2 - y - z)*(x**2)
185expected = [2, 1, 3, 4, 2, 3]
186check_feasible(p, x, assignment, expected)
187
188assignment = polypy.Assignment()
189assignment.set_value(y, 1)
190assignment.set_value(z, 1)
191
192# + 0 - 0 +
193p = x**2 - y - z
194expected = [1, 1, 2, 3, 2, 2]
195check_feasible(p, x, assignment, expected)
196
197assignment = polypy.Assignment()
198assignment.set_value(y, 0)
199assignment.set_value(z, 1)
200
201# - 0 + 0 -
202p = (x - y)*(x - z)*(y - z)
203expected = [2, 2, 2, 3, 1, 1]
204check_feasible(p, x, assignment, expected)
205
206assignment = polypy.Assignment()
207assignment.set_value(y, 0)
208assignment.set_value(z, 1)
209
210# +
211p = x**2 + y**2 + z**2
212expected = [0, 0, 0, 1, 1, 1]
213check_feasible(p, x, assignment, expected)
214