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