Skip to content

Commit 180f61d

Browse files
mkoeppemkoeppe
committed
Merge branch 'computer_based_search_squashed'
2 parents 1c982d8 + b19c501 commit 180f61d

9 files changed

+4753
-9
lines changed

‎.gitignore

+6
Original file line numberDiff line numberDiff line change
@@ -8,3 +8,9 @@
88
desktop.ini
99
.dropbox
1010
Icon*
11+
2q_mip/*
12+
profiler/*
13+
sym_mode_2d_diagrams/*
14+
animation/*
15+
region_graphics/*
16+

‎2q_mip.sage

+298
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,298 @@
1+
# Make sure current directory is in path.
2+
# That's not true while doctesting (sage -t).
3+
if '' not in sys.path:
4+
sys.path = [''] + sys.path
5+
6+
from igp import *
7+
8+
##########################################
9+
# MIP approach to search for 2q example
10+
# use code in kslope_mip.sage
11+
##########################################
12+
13+
def print_trivial_additive_points_2q(filename, q, f, a):
14+
"""
15+
EXAMPLES::
16+
17+
sage: print_trivial_additive_points_2q(sys.stdout, 3, 2/3, 1/3)
18+
p_0_0 = 0
19+
p_0_1 = 0
20+
p_0_2 = 0
21+
p_0_3 = 0
22+
p_1_3 = 0
23+
p_2_3 = 0
24+
p_3_3 = 0
25+
p_0_2 = 0
26+
p_1_1 = 0
27+
p_2_0 = 0
28+
p_2_3 = 0
29+
p_3_2 = 0
30+
p_1_0 = 0
31+
p_1_1 = 0
32+
"""
33+
bkpt = [x/q for x in range(q+1)]
34+
# border x = 0 and border y = 0 are green
35+
for x in bkpt:
36+
print >> filename, '%s = 0' % vertex_variable(q, (0, x))
37+
for x in bkpt[1::]:
38+
print >> filename, '%s = 0' % vertex_variable(q, (x, 1))
39+
# diagonals corresponding to f
40+
for x in bkpt:
41+
if x < f:
42+
print >> filename, '%s = 0' % vertex_variable(q, (x, f - x))
43+
elif x == f:
44+
print >> filename, '%s = 0' % vertex_variable(q, (x, f - x))
45+
print >> filename, '%s = 0' % vertex_variable(q, (x, f - x + 1))
46+
elif x > f:
47+
print >> filename, '%s = 0' % vertex_variable(q, (x, f - x + 1))
48+
49+
b = f - a
50+
print >> filename, '%s = 0' % vertex_variable(q, (b - a + 1/q, a - 1/q))
51+
print >> filename, '%s = 0' % vertex_variable(q, (b - a + 1/q, a))
52+
53+
def write_lpfile_2q(q, f, a, kslopes, maxstep=None, m=0):
54+
"""
55+
EXAMPLES:
56+
57+
sage: write_lpfile_2q(37, 25/37, 11/37, 4, maxstep=2, m=4) # not tested
58+
"""
59+
if maxstep is None:
60+
maxstep = q
61+
62+
destdir = output_dir+"2q_mip/"
63+
mkdir_p(destdir)
64+
filename = open(destdir + "mip_q%s_f%s_a%s_%sslope_%smaxstep_m%s.lp" % (q, int(f*q), int(a*q), kslopes, maxstep, m), "w")
65+
faces_2d = []
66+
faces_diag = []
67+
faces_hor = []
68+
faces_ver = []
69+
faces_0d = []
70+
for xx in range(q):
71+
for yy in range(q):
72+
faces_2d.append( Face(([xx/q, (xx+1)/q], [yy/q, (yy+1)/q], [(xx+yy)/q, (xx+yy+1)/q])) )
73+
faces_2d.append( Face(([xx/q, (xx+1)/q], [yy/q, (yy+1)/q], [(xx+yy+1)/q, (xx+yy+2)/q])) )
74+
faces_diag.append( Face(([xx/q, (xx+1)/q], [yy/q, (yy+1)/q], [(xx+yy+1)/q])) )
75+
76+
for xx in range(q):
77+
for yy in range(q+1):
78+
faces_hor.append( Face(([xx/q, (xx+1)/q], [yy/q], [(xx+yy)/q, (xx+yy+1)/q])) )
79+
faces_ver.append( Face(([yy/q], [xx/q, (xx+1)/q], [(xx+yy)/q, (xx+yy+1)/q])) )
80+
81+
for xx in range(q+1):
82+
for yy in range(q+1):
83+
faces_0d.append( Face(([xx/q], [yy/q], [(xx+yy)/q])) )
84+
85+
print >> filename, '\ MIP model with q = %s, f = %s, a = %s, num of slopes = %s, maxstep of tran/refl = %s, small_m = %s' % (q, f, a, kslopes, maxstep, m)
86+
87+
print >> filename, 'Maximize'
88+
#print >> filename, 0
89+
print_obj_max_slope_slack(filename, kslopes)
90+
#print_obj_max_subadd_slack(filename, q) # is a constant!
91+
#print_obj_min_directly_covered_times(filename, q)
92+
#print_obj_min_undirectly_covered_times(filename, q)
93+
#print_obj_min_covered_times_max_subadd_slack(filename, q, maxstep=maxstep)
94+
#print_obj_5slope22(filename, q, weight=1)
95+
#print_obj_min_add_points(filename, q, weight=1)
96+
print >> filename
97+
98+
print >> filename, 'Subject to'
99+
for face in faces_2d + faces_diag + faces_hor + faces_ver:
100+
#if face.minimal_triple[0][0] <= face.minimal_triple[1][0]:
101+
print_logical_constraints(filename, q, face)
102+
103+
for face in faces_0d:
104+
if face.minimal_triple[0][0] < face.minimal_triple[1][0]:
105+
print_xy_swapped_constraints(filename, q, face)
106+
107+
# if no 1/m trick for subadditivity, set m=0 in print_fn_minimality_test.
108+
print_fn_minimality_test(filename, q, f, m)
109+
110+
print_trivial_additive_points_2q(filename, q, f, a)
111+
112+
for zz in range(q):
113+
for xx in range(q):
114+
x = xx / q
115+
z = zz / q
116+
if x != z:
117+
if maxstep > 1:
118+
print_move_constraints(filename, q, x, z)
119+
for step in range(1, maxstep):
120+
print_translation_i_constraints(filename, q, x, z, step)
121+
print_reflection_i_constraints(filename, q, x, z, step)
122+
123+
for zz in range(q):
124+
z = zz / q
125+
print_directly_covered_constraints(filename, q, z)
126+
for step in range(1, maxstep):
127+
print_undirectly_covered_i_constraints(filename, q, z, step)
128+
129+
for zz in range(q):
130+
z = zz / q
131+
if z != a - 1/q and z != f - a:
132+
print >> filename, '%s = 0' % covered_i_variable(q, z, maxstep - 1)
133+
else:
134+
print >> filename, '%s = 1' % covered_i_variable(q, z, maxstep - 1)
135+
136+
print_slope_constraints_2q(filename, q, f, a, kslopes, m)
137+
138+
print >> filename, 'Bounds'
139+
print_fn_bounds(filename, q)
140+
print_slope_bounds(filename, q, kslopes)
141+
142+
print >> filename, 'Binary'
143+
for face in faces_2d + faces_diag + faces_hor + faces_ver + faces_0d :
144+
print >> filename, face_variable(q, face),
145+
146+
for z in range(q):
147+
for step in range(maxstep):
148+
print >> filename, 'c_%s_%s' % (z, step),
149+
for z in range(q):
150+
for x in range(q):
151+
if x != z:
152+
if maxstep > 1:
153+
print >> filename, 'm_%s_%s' % (x, z),
154+
for step in range(1, maxstep):
155+
print >> filename, 't_%s_%s_%s' % (x, z, step),
156+
print >> filename, 'r_%s_%s_%s' % (x, z, step),
157+
158+
for k in range(kslopes):
159+
for j in range(q):
160+
print >> filename, '%s' % interval_slope_variable(j, k),
161+
162+
print >> filename
163+
print >> filename, 'End'
164+
filename.close()
165+
166+
def print_slope_constraints_2q(filename, q, f, a, kslopes, m=0):
167+
"""
168+
EXAMPLES::
169+
170+
sage: print_slope_constraints_2q(sys.stdout, 4, 3/4, 2/4, 3, m=0)
171+
s_0 - s_1 >= 0
172+
s_1 - s_2 >= 0
173+
s_0 - 4 fn_1 = 0
174+
i_0_s_0 = 1
175+
s_2 + 4 fn_3 = 0
176+
i_3_s_2 = 1
177+
s_0 + 4 fn_1 - 4 fn_2 + 8 i_1_s_0 <= 8
178+
s_0 + 4 fn_1 - 4 fn_2 - 8 i_1_s_0 >= -8
179+
s_1 + 4 fn_1 - 4 fn_2 + 8 i_1_s_1 <= 8
180+
s_1 + 4 fn_1 - 4 fn_2 - 8 i_1_s_1 >= -8
181+
s_2 + 4 fn_1 - 4 fn_2 + 8 i_1_s_2 <= 8
182+
s_2 + 4 fn_1 - 4 fn_2 - 8 i_1_s_2 >= -8
183+
s_0 + 4 fn_2 - 4 fn_3 + 8 i_2_s_0 <= 8
184+
s_0 + 4 fn_2 - 4 fn_3 - 8 i_2_s_0 >= -8
185+
s_1 + 4 fn_2 - 4 fn_3 + 8 i_2_s_1 <= 8
186+
s_1 + 4 fn_2 - 4 fn_3 - 8 i_2_s_1 >= -8
187+
s_2 + 4 fn_2 - 4 fn_3 + 8 i_2_s_2 <= 8
188+
s_2 + 4 fn_2 - 4 fn_3 - 8 i_2_s_2 >= -8
189+
+ i_0_s_0 + i_0_s_1 + i_0_s_2 = 1
190+
+ i_1_s_0 + i_1_s_1 + i_1_s_2 = 1
191+
+ i_2_s_0 + i_2_s_1 + i_2_s_2 = 1
192+
+ i_3_s_0 + i_3_s_1 + i_3_s_2 = 1
193+
+ i_0_s_0 + i_1_s_0 + i_2_s_0 + i_3_s_0 >= 1
194+
+ i_0_s_1 + i_1_s_1 + i_2_s_1 + i_3_s_1 >= 1
195+
+ i_0_s_2 + i_1_s_2 + i_2_s_2 + i_3_s_2 >= 1
196+
i_1_s_0 + i_0_s_0 <= 1
197+
i_1_s_0 - i_1_s_0 = 0
198+
i_1_s_0 + i_2_s_0 <= 1
199+
i_1_s_0 + i_3_s_0 <= 1
200+
i_1_s_1 + i_0_s_1 <= 1
201+
i_1_s_1 - i_1_s_1 = 0
202+
i_1_s_1 + i_2_s_1 <= 1
203+
i_1_s_1 + i_3_s_1 <= 1
204+
i_1_s_2 + i_0_s_2 <= 1
205+
i_1_s_2 - i_1_s_2 = 0
206+
i_1_s_2 + i_2_s_2 <= 1
207+
i_1_s_2 + i_3_s_2 <= 1
208+
"""
209+
# s_0 > s_1 > ... > s_kslopes-1
210+
for k in range(0, kslopes - 1):
211+
if m == 0:
212+
print >> filename, '%s - %s >= 0' % (slope_variable(k), slope_variable(k+1))
213+
else:
214+
print >> filename, '%s - %s >= %s' % (slope_variable(k), slope_variable(k+1), RR(q/m))
215+
216+
# first interval has the largest positive slope s_0
217+
print >> filename, 's_0 - %s fn_1 = 0' % q
218+
print >> filename, 'i_0_s_0 = 1'
219+
# last interval has slope s_kslopes-1
220+
print >> filename, 's_%s + %s fn_%s = 0' % (kslopes - 1, q, q - 1)
221+
print >> filename, 'i_%s_s_%s = 1' % (q - 1, kslopes - 1)
222+
# Condition: s_k + q(fn_j - fn_(j+1)) = 0 iff i_j_s_k = 1
223+
# ==> 1) s_k + q * fn_j - q * fn_(j+1) <= 2*q * (1 - i_j_s_k)
224+
# ==> 2) s_k + q * fn_j - q * fn_(j+1) >= - 2*q * (1 - i_j_s_k)
225+
# ==> 3) sum i_j_s_k over k = 1
226+
for j in range(1, q-1):
227+
for k in range(kslopes):
228+
print >> filename, 's_%s + %s fn_%s - %s fn_%s + %s %s <= %s' % (k, q, j, q, j + 1, 2*q, interval_slope_variable(j, k), 2*q)
229+
print >> filename, 's_%s + %s fn_%s - %s fn_%s - %s %s >= %s' % (k, q, j, q, j + 1, 2*q, interval_slope_variable(j, k), -2*q)
230+
for j in range(q):
231+
for k in range(kslopes):
232+
print >> filename, '+ %s' % interval_slope_variable(j, k),
233+
print >> filename, '= 1'
234+
# Condition: sum i_j_s_k over j >= 1
235+
for k in range(kslopes):
236+
for j in range(q):
237+
print >> filename, '+ %s' % interval_slope_variable(j, k),
238+
print >> filename, '>= 1'
239+
# Two special intervals have the same slope value,
240+
# which is different from the slope values of other intervals
241+
ja = int(a * q) - 1
242+
jb = int((f - a) * q)
243+
for k in range(kslopes):
244+
for j in range(q):
245+
if j == jb:
246+
print >> filename, '%s - %s = 0' % ( interval_slope_variable(ja, k), \
247+
interval_slope_variable(jb, k) )
248+
elif j != ja:
249+
print >> filename, '%s + %s <= 1' % ( interval_slope_variable(ja, k), \
250+
interval_slope_variable(j, k) )
251+
252+
def refind_function_from_lpsolution_2q(filename, q, f, a):
253+
"""
254+
EXAMPLES::
255+
256+
sage: h = refind_function_from_lpsolution_2q('solution_2q_example_m4.sol', 37, 25/37, 11/37) # not tested
257+
"""
258+
faces, fn = painted_faces_and_funciton_from_solution(filename, q)
259+
covered_intervals = generate_covered_intervals_from_faces(faces)
260+
additive_vertices = generate_additive_vertices_from_faces(q, faces)
261+
final_uncovered = [[a - 1/q, a], [f - a, f - a + 1 / q]]
262+
components = covered_intervals + [final_uncovered]
263+
fn_sym = generate_symbolic_continuous(None, components, field=QQ)
264+
ff = int(f * q)
265+
for h in generate_vertex_function(q, ff, fn_sym, additive_vertices):
266+
if not extremality_test(h): # h is not extreme
267+
logging.info("Testing the extremality of pi restricted to 1/2q.")
268+
h_2q = restrict_to_finite_group(h, f=f, oversampling=2, order=None)
269+
if extremality_test(h_2q): # but h restricted to 1/2q is extreme
270+
logging.info("Find a valid 2q_example!")
271+
return h
272+
# Comments:
273+
# def kzh_2q_example_1() is obtained by
274+
# LP model: mip_q37_f25_a11_4slope_2maxstep_m4.lp
275+
# q = 37; f = 25/37; a = 11/37;
276+
# given number of slopes = 4; maxstep = 2;
277+
# slope gap >= q/4; subadditivity slack = 1/4;
278+
# obj = max_slope_slack.
279+
#
280+
# Gurobi run on my laptop (276s, obj = 33.66393)
281+
# Optimal solution: solution_2q_example_m4.sol
282+
#
283+
# ( Or
284+
# (1) LP model: mip_q37_f25_a11_4slope_2maxstep_m12.lp
285+
# slope gap >= q/12; subadditivity slack = 1/12;
286+
# Gurobi run on my laptop (1080s, obj = 33.66393)
287+
# Optimal solution: solution_2q_example_m12.sol
288+
#
289+
# (2) LP model: mip_q37_f25_a11_4slope_2maxstep_m12slopegaponly.lp
290+
# slope gap >= q/12; no 1/m trick on subadditivity slack;
291+
# Gurobi run on point.math.ucdavis.edu (7412s, obj = 33.66393)
292+
# Optimal solution: solution_2q_example_m12slopegaponly.sol )
293+
#
294+
# The vertex function is a 4-slope non-extreme function,
295+
# whose restriction to 1/2q is extreme.
296+
#
297+
# This two_q_example is obtained by:
298+
# sage: h = refind_function_from_lpsolution_2q('solution_2q_example_m4.sol', 37, 25/37, 11/37) # not tested

‎Makefile

+7-2
Original file line numberDiff line numberDiff line change
@@ -12,9 +12,14 @@ SAGEFILES = \
1212
fast_linear.sage \
1313
functions.sage \
1414
simple_extremality_test.sage \
15-
survey_examples.sage \
15+
survey_examples.sage \
1616
extreme_functions_mlr_cpl3.sage \
17-
quasi_periodic.sage
17+
quasi_periodic.sage \
18+
kslope_ppl_mip.py \
19+
vertex_enumeration.py \
20+
kslope_pattern.sage \
21+
2q_mip.sage \
22+
kslope_mip.sage
1823

1924
all:
2025
@echo "No need to 'make' anything. Just run it in Sage; see README.rst"

‎functions.sage

+18-7
Original file line numberDiff line numberDiff line change
@@ -261,7 +261,12 @@ class Face:
261261
def is_diagonal(self):
262262
return self.is_1D() and \
263263
self.vertices[0][0] + self.vertices[0][1] == self.vertices[1][0] + self.vertices[1][1]
264-
264+
def __hash__(self):
265+
return sum([hash(x) for i in self.minimal_triple for x in i])
266+
def __cmp__(left, right):
267+
return cmp(left.minimal_triple, right.minimal_triple)
268+
269+
265270
def plot_faces(faces, **kwds):
266271
p = Graphics()
267272
for f in faces:
@@ -800,8 +805,12 @@ def generate_directly_covered_intervals(function):
800805
return function._directly_covered_intervals
801806

802807
faces = generate_maximal_additive_faces(function)
808+
covered_intervals = generate_directly_covered_intervals_from_faces(faces)
809+
function._directly_covered_intervals = covered_intervals
810+
return covered_intervals
803811

804-
covered_intervals = []
812+
def generate_directly_covered_intervals_from_faces(faces):
813+
covered_intervals = []
805814
for face in faces:
806815
if face.is_2D():
807816
component = []
@@ -822,17 +831,21 @@ def generate_directly_covered_intervals(function):
822831
covered_intervals[i] = []
823832

824833
covered_intervals = remove_empty_comp(covered_intervals)
825-
function._directly_covered_intervals = covered_intervals
826834
return covered_intervals
827835

828836
def generate_covered_intervals(function):
829837
if hasattr(function, '_covered_intervals'):
830838
return function._covered_intervals
831839

832840
logging.info("Computing covered intervals...")
833-
covered_intervals = copy(generate_directly_covered_intervals(function))
834841
faces = generate_maximal_additive_faces(function)
842+
covered_intervals = generate_covered_intervals_from_faces(faces)
843+
logging.info("Computing covered intervals... done")
844+
function._covered_intervals = covered_intervals
845+
return covered_intervals
835846

847+
def generate_covered_intervals_from_faces(faces):
848+
covered_intervals = generate_directly_covered_intervals_from_faces(faces)
836849
# debugging plot:
837850
# show(plot_covered_intervals(function, covered_intervals), \
838851
# legend_fancybox=True, \
@@ -861,9 +874,6 @@ def generate_covered_intervals(function):
861874
any_change = True
862875

863876
covered_intervals = remove_empty_comp(covered_intervals)
864-
logging.info("Computing covered intervals... done")
865-
866-
function._covered_intervals = covered_intervals
867877
return covered_intervals
868878

869879
def uncovered_intervals_from_covered_intervals(covered_intervals):
@@ -1506,6 +1516,7 @@ class FastPiecewise (PiecewisePolynomial):
15061516
15071517
sage: f1(x) = 1
15081518
sage: f2(x) = 1-x
1519+
15091520
sage: f3(x) = exp(x)
15101521
sage: f4(x) = sin(2*x)
15111522
sage: f = FastPiecewise([[(0,1),f1],

0 commit comments

Comments
 (0)