this repo has no description
dekker.one
1/*
2% FlatZinc built-in redefinitions for linear solvers.
3%
4% AUTHORS
5% Sebastian Brand
6% Gleb Belov (2015-)
7% cf. Belov, Stuckey, Tack, Wallace. Improved Linearization of Constraint Programming Models. CP 2016.
8*/
9
10%----------------------------- BOOL2INT --------------------------------%
11function var bool: reverse_map(var int: x) = (x==1);
12function bool: reverse_map(int: x) = (x==1);
13
14predicate mzn_reverse_map_var(var bool: b) = let {
15 var int: x = bool2int(b)
16} in true;
17
18function var int: bool2int(var bool: x) :: promise_total =
19 let { var 0..1: b2i;
20 constraint (x = reverse_map(b2i)) ::is_reverse_map ;
21 } in
22 b2i;
23
24predicate bool_eq(var bool: x, var bool: y) =
25%% trace(" bool_eq: \(x), \(y) \n") /\
26 bool2int(x)==bool2int(y);
27
28predicate bool2int(var bool: x, var int: y) = y = bool2int(x);
29
30%---------------------------- BASIC (HALF)REIFS -----------------------------%
31
32include "options.mzn";
33include "redefs_bool_reifs.mzn";
34include "redefs_bool_imp.mzn";
35
36include "domain_encodings.mzn";
37include "redefs_lin_reifs.mzn";
38include "redefs_lin_imp.mzn";
39include "redefs_lin_halfreifs.mzn";
40
41include "nosets.mzn"; %% For set_le, set_lt ... Usind std/nosets
42 %% as long as the linearization is good.
43
44%-----------------------------------------------------------------------------%
45% Strict inequality
46
47% Uncomment the following redefinition for FlatZinc MIP solver interfaces that
48% do not support strict inequality. Note that it does not preserve equivalence
49% (some solutions of the original problem may become invalid).
50
51predicate float_lt(var float: x, var float: y) =
52% (x - y) <= (-float_lt_EPS_coef__)*max(abs(ub(x - y)), abs(ub(y-x)));
53 x <= y - float_lt_EPS;
54
55predicate float_lin_lt(array[int] of float: c, array[int] of var float: x,
56 float: d) =
57 float_lt(sum(i in index_set(x))( c[i]*x[i] ), d);
58
59%-----------------------------------------------------------------------------%
60% Minimum, maximum, absolute value
61% Use unary as well? TODO
62
63predicate int_abs(var int: x, var int: z) =
64%% The simplifications seem worse on league.mzn model90-18-20.dzn:
65%% but the .lp seem to differ just by order...?? TODO
66 if lb(x)>=0 then z==x
67 elseif ub(x)<=0 then z==-x
68 else
69 let { var bool: p }
70 in
71 z >= x /\
72 z >= -x /\
73 z >= 0 /\ % This is just for preprocessor
74 z <= max([ub(x), -lb(x)]) /\ % And this
75 % z <= x \/ z <= -x %% simple
76 aux_int_le_if_1(z, x, p) /\ %% even simpler
77 aux_int_le_if_0(z, -x, p) /\
78 int_le_reif(0, x, p) % with reifs
79 %int_eq_reif(z, x, p) /\
80 %int_eq_reif(z, -x, not p)
81 endif
82 ;
83
84predicate int_min(var int: x, var int: y, var int: z) =
85 array_int_minimum(z, [x, y]);
86
87predicate int_max(var int: x, var int: y, var int: z) =
88 array_int_maximum(z, [x, y]);
89
90predicate float_abs(var float: x, var float: z) =
91 if lb(x)>=0.0 then z==x
92 elseif ub(x)<=0.0 then z==-x
93 else
94 let { var bool: p }
95 in
96 z >= x /\
97 z >= -x /\
98 z >= 0.0 /\ % This is just for preprocessor
99 z <= max([ub(x), -lb(x)]) /\ % And this
100 % z <= x \/ z <= -x
101 aux_float_le_if_1(z, x, (p)) /\
102 aux_float_le_if_0(z, -x, (p))
103 % /\
104 % float_le_reif(0.0, x, p) % with reifs - no point for floats? TODO
105 % float_eq_reif(z, x, p) /\
106 % float_eq_reif(z, -x, not p)
107 endif;
108
109predicate float_min(var float: x, var float: y, var float: z) =
110 array_float_minimum(z, [x, y]);
111
112predicate float_max(var float: x, var float: y, var float: z) =
113 array_float_maximum(z, [x, y]);
114
115predicate array_int_minimum_I(var int: m, array[int] of var int: x) =
116 let { int: n = length(x),
117 constraint assert(1 == min(index_set(x)), " array_int_minimum_I: argument indexed not from 1??"),
118 int: iMinUB = arg_min([ub(x[i]) | i in 1..n]),
119 int: MinUB = ub(x[iMinUB]),
120 set of int: sLBLess = { i | i in 1..n where lb(x[i])<MinUB },
121 set of int: sUBEqual = { i | i in 1..n where ub(x[i])==MinUB },
122 set of int: sActive = if card(sLBLess intersect sUBEqual)>0 then sLBLess
123 else sLBLess union { iMinUB } endif,
124 } in
125 if 1==card(sActive) then
126 m == x[min(sActive)]
127 elseif MZN__MinMaxGeneral then
128 fzn_array_float_minimum(m, x) %% Send to backend
129 else
130 let {
131 array[1..n] of var int: p = [
132 if i in sActive then let { var 0..1: pi } in pi else 0 endif | i in 1..n ],
133 constraint 1==sum(p),
134 constraint m >= lb_array(x),
135 constraint m <= MinUB,
136 } in
137 forall (i in index_set(x))
138 (
139 if i in sActive %% for at least 1 element
140 then
141 m<=x[i] /\ aux_int_ge_if_1(m, x[i], p[i])
142 endif ) %% -- exclude too big x[i]
143 endif;
144
145predicate array_float_minimum_I(var float: m, array[int] of var float: x) =
146 let { int: n = length(x),
147 constraint assert(1 == min(index_set(x)), " array_float_minimum_I: argument indexed not from 1??"),
148 int: iMinUB = arg_min([ub(x[i]) | i in 1..n]),
149 float: MinUB = ub(x[iMinUB]),
150 set of int: sLBLess = { i | i in 1..n where lb(x[i])<MinUB },
151 set of int: sUBEqual = { i | i in 1..n where ub(x[i])==MinUB },
152 set of int: sActive = if card(sLBLess intersect sUBEqual)>0 then sLBLess
153 else sLBLess union { iMinUB } endif,
154 } in
155 if 1==card(sActive) then
156 m == x[min(sActive)]
157 elseif MZN__MinMaxGeneral then
158 fzn_array_float_minimum(m, x) %% Send to backend
159 else
160 let {
161 array[1..n] of var int: p = [
162 if i in sActive then let { var 0..1: pi } in pi else 0 endif | i in 1..n ],
163 constraint 1==sum(p),
164 constraint m >= lb_array(x),
165 constraint m <= MinUB,
166 } in
167 forall (i in index_set(x))
168 (
169 if i in sActive %% for at least 1 element
170 then
171 m<=x[i] /\ aux_float_ge_if_1(m, x[i], p[i])
172 endif ) %% -- exclude too big x[i]
173 /\
174 if card(sActive)>1 /\ fMinimumCutsXZ then
175 let {
176 array[int] of float: AL = [ lb(x[i]) | i in 1..n],
177 array[int] of int: srt = sort_by([i | i in 1..n], AL),
178 %indices of lb in sorted order
179 array[int] of float: AL_srt = [AL[srt[i]] | i in 1..n],
180 array[int] of float: AU_srt = [ub(x[srt[i]]) | i in 1..n],
181 array[int] of float: AM_srt = AL_srt ++ [MinUB]
182 %% -- these are z-levels of extreme points
183 } in
184 forall (i in 2..n+1 where
185 AM_srt[i]<=MinUB /\ %% this is a new "start level"
186 AM_srt[i]!=AM_srt[i-1] )( %% and would produce a new cut
187 m >= AM_srt[i]
188 - sum(j in 1..i-1 where AL_srt[j]<AM_srt[i] /\ AL_srt[j]<AU_srt[j])
189 ( (AU_srt[j]-x[srt[j]]) * (AM_srt[i]-AL_srt[j]) / (AU_srt[j]-AL_srt[j]) ) :: MIP_cut
190 )
191 else true endif
192 /\
193 if card(sActive)>1 /\ fMinimumCutsXZB then
194 array_var_float_element__XBZ_lb([ -x[i] | i in sActive ], [ p[i] | i in sActive ], -m) :: MIP_cut
195 else true endif
196 endif;
197
198%-----------------------------------------------------------------------------%
199% Multiplication and division
200
201predicate int_div(var int: x, var int: y, var int: q) =
202 q == aux_int_division_modulo_fn(x,y)[1];
203
204
205predicate int_mod(var int: x, var int: y, var int: r) =
206 r == aux_int_division_modulo_fn(x,y)[2];
207
208
209function array[int] of var int: aux_int_division_modulo_fn(var int: x, var int: y) =
210 let {
211 %% Domain of q
212 set of int: dom_q =
213 if lb(y)*ub(y)>0 then
214 let {
215 set of int: EP = { ub(x) div ub(y), ub(x) div lb(y), lb(x) div ub(y), lb(x) div lb(y) },
216 } in min(EP)..max(EP)
217 else
218 let {
219 int: mm = max( abs(lb(x)), abs(ub(x)) ),
220 } in -mm..mm %% TODO case when -1 or 1 not in dom(x)
221 endif,
222 var dom_q: q;
223 int: by = max(abs(lb(y)), abs(ub(y)));
224 var -by+1..by-1: r;
225 constraint x = y * q + r,
226 constraint 0 <= x -> 0 <= r, %% which is 0 > x \/ 0 <= r
227 constraint x < 0 -> r <= 0, %% which is x >= 0 \/ r <= 0
228 % abs(r) < abs(y)
229 var 1.. max(abs(lb(y)), abs(ub(y))): w = abs(y),
230 constraint w > r /\ w > -r,
231 } in
232 [ q, r ];
233
234%% Can also have int_times(var float, var int) ......... TODO
235
236predicate int_times(int: x, var int: y, var int: z);
237predicate int_times(var int: x, int: y, var int: z) = int_times(y, x, z);
238
239predicate int_times(var int: x, var int: y, var int: z) =
240 if 0..1==dom(x) /\ 0..1==dom(y) then bool_and__INT(x,y,z)
241 elseif card(dom(x))==2 /\ card(dom(y))==2 /\ 0 in dom(x) /\ 0 in dom(y)
242 then let {
243 var 0..1: xn;
244 var 0..1: yn;
245 var 0..1: zn;
246 constraint x=xn*max(dom(x) diff {0});
247 constraint y=yn*max(dom(y) diff {0});
248 constraint z=zn*max(dom(x) diff {0})*max(dom(y) diff {0});
249 } in
250 bool_and__INT(xn,yn,zn)
251 elseif min(card(dom(x)), card(dom(y))) >= MZN__QuadrIntCard then
252 fzn_int_times(x, y, z)
253 elseif card(dom(x)) * card(dom(y)) > nMZN__UnarySizeMax_intTimes
254 \/ ( fIntTimesBool /\ (
255 %% Peter's idea for *bool. More optimal but worse values on carpet cutting.
256 (card(dom(x))==2 /\ 0 in dom(x))
257 \/ (card(dom(y))==2 /\ 0 in dom(y))
258 ) )
259 then %% PARAM
260 %% ALSO NO POINT IF <=4. TODO
261 if card(dom(x)) > card(dom(y)) \/
262 ( card(dom(x))==card(dom(y)) /\ 0 in dom(y) /\ not (0 in dom(x)) )
263 then int_times(y,x,z)
264 else
265 let {
266 set of int: s = lb(x)..ub(x),
267 set of int: r = {lb(x)*lb(y), lb(x)*ub(y), ub(x)*lb(y), ub(x)*ub(y)},
268 array[s] of var min(r)..max(r): ady = array1d(s, [
269 if d in dom(x) then d*y else min(r) endif | d in s ]) }
270 in
271 ady[x] = z %% use element()
272 endif
273 else
274 int_times_unary(x, { }, y, z)
275 endif;
276
277%% domx__ can be used to narrow domain... NOT IMPL.
278predicate int_times_unary(var int: x, set of int: domx__, var int: y, var int: z) =
279 let {
280 set of int: r = {lb(x)*lb(y), lb(x)*ub(y), ub(x)*lb(y), ub(x)*ub(y)},
281 %% set of int: domx = if card(domx__)>0 then domx__ else dom(x) endif,
282 array[int, int] of var int: pp=eq_encode(x, y)
283 } in
284 z>=min(r) /\ z<=max(r) /\
285 z==sum(i in index_set_1of2(pp), j in index_set_2of2(pp))
286 (i * j * pp[i, j]) /\
287 forall(i in index_set_1of2(pp), j in index_set_2of2(pp)
288 where not ((i*j) in dom(z))
289 )(pp[i, j]==0)
290 ;
291
292
293predicate int_times_unary__NOFN(var int: x, set of int: domx__, var int: y, var int: z) =
294 let {
295 set of int: r = {lb(x)*lb(y), lb(x)*ub(y), ub(x)*lb(y), ub(x)*ub(y)},
296 %% set of int: domx = if card(domx__)>0 then domx__ else dom(x) endif,
297 array[int] of var int: pX = eq_encode(x),
298 array[int] of var int: pY = eq_encode(y),
299 array[int] of int: valX = [ v | v in index_set(pX) ], %% NOT domx.
300 array[int] of int: valY = [ v | v in index_set(pY) ], %% -- according to eq_encode!
301 array[index_set(valX), index_set(valY)] of var 0..1: pp %% both dim 1..
302 } in
303 if is_fixed(x) \/ is_fixed(y) then
304 z==x*y
305 else
306 z>=min(r) /\ z<=max(r) /\
307 sum(pp)==1 /\
308 z==sum(i in index_set(valX), j in index_set(valY))
309 (valX[i] * valY[j] * pp[i, j]) /\
310 forall(i in index_set(valX))
311 ( pX[valX[i]] == sum(j in index_set(valY))( pp[i, j] ) ) /\
312 forall(j in index_set(valY))
313 ( pY[valY[j]] == sum(i in index_set(valX))( pp[i, j] ) )
314 endif;
315
316predicate float_times(var float: x, var float: y, var float: z) =
317 if is_fixed(x) then
318 z==fix(x)*y %%%%% Need to use fix() otherwise added to CSE & nothing happens
319 elseif is_fixed(y) then
320 z==x*fix(y)
321 elseif MZN__QuadrFloat then
322 fzn_float_times(x, y, z)
323 else
324 abort("Unable to create linear formulation for the `float_times(\(x), \(y), \(z))` constraint." ++
325 " Define QuadrFloat=true if your linear solver supports quadratic constraints," ++
326 " or use integer variables.")
327 endif;
328
329%%%Define int_pow
330predicate int_pow( var int: x, var int: y, var int: r ) =
331 let {
332 array[ int, int ] of int: x2y = array2d( lb(x)..ub(x), lb(y)..ub(y),
333 [ pow( X, Y ) | X in lb(x)..ub(x), Y in lb(y)..ub(y) ] )
334 } in
335 r == x2y[ x, y ];
336
337%%% Adding a version returning float for efficiency
338/** @group builtins.arithmetic Return \(\a x ^ {\a y}\) */
339function var float: pow_float(var int: x, var int: y) =
340 let {
341 int: yy = if is_fixed(y) then fix(y) else -1 endif;
342 } in
343 if yy = 0 then 1
344 elseif yy = 1 then x else
345 let { var float: r ::is_defined_var;
346 constraint int_pow_float(x,y,r) ::defines_var(r);
347 } in r
348 endif;
349%%%Define int_pow_float
350predicate int_pow_float( var int: x, var int: y, var float: r ) =
351 let {
352 array[ int, int ] of float: x2y = array2d( lb(x)..ub(x), lb(y)..ub(y),
353 [ pow( X, Y ) | X in lb(x)..ub(x), Y in lb(y)..ub(y) ] )
354 } in
355 r == x2y[ x, y ];
356
357
358
359%-----------------------------------------------------------------------------%
360% Array 'element' constraints
361
362predicate array_bool_element(var int: x, array[int] of bool: a, var bool: z) =
363 array_int_element(x, arrayXd(a, [bool2int(a[i]) | i in index_set(a)]), bool2int(z));
364
365predicate array_var_bool_element(var int: x, array[int] of var bool: a,
366 var bool: z) =
367 array_var_int_element(x, arrayXd(a, [bool2int(a[i]) | i in index_set(a)]), bool2int(z));
368
369predicate array_int_element(var int: i00, array[int] of int: a, var int: z) =
370 let {
371 set of int: ix = index_set(a),
372 constraint i00 in { i | i in ix where a[i] in dom(z) },
373 } in %%% Tighten domain of i00 before dMin/dMax
374 let {
375 int: dMin = min(i in dom(i00))(a[i]),
376 int: dMax = max(i in dom(i00))(a[i]),
377 } in
378 if dMin==dMax then
379 z==dMin
380 else
381 z >= dMin /\
382 z <= dMax /\
383 let {
384 int: nUBi00 = max(dom(i00)),
385 int: nLBi00 = min(dom(i00)),
386 int: nMinDist = min(i in nLBi00 .. nUBi00-1)(a[i+1]-a[i]),
387 int: nMaxDist = max(i in nLBi00 .. nUBi00-1)(a[i+1]-a[i]),
388 } in
389 if nMinDist == nMaxDist then %% The linear case
390 z == a[nLBi00] + nMinDist*(i00-nLBi00)
391 else
392 let {
393 array[int] of var int: p = eq_encode(i00) %% this needs i00 in ix
394 } %% Faster flattening than (i==i00) @2a9df1f7
395 in
396 sum(i in dom(i00))( a[i] * p[i] ) == z
397 endif
398 endif;
399
400predicate array_var_int_element(var int: i00, array[int] of var int: a,
401 var int: z) =
402 let {
403 constraint i00 in { i | i in index_set(a) where
404 0 < card(dom(a[i]) intersect dom(z)) },
405 } in %% finish domain first
406 let {
407 int: minLB=min(i in dom(i00))(lb(a[i])),
408 int: maxUB=max(i in dom(i00))(ub(a[i]))
409 } in
410 if minLB==maxUB then
411 z==minLB
412 else
413 z >= minLB /\
414 z <= maxUB /\
415 if {0,1}==dom(i00) /*ub(i00)-lb(i00)==1*/ /*2==card( dom( i00 ) )*/ then
416 aux_int_eq_if_1(z, a[lb(i00)], (ub(i00)-i00)) /\
417 aux_int_eq_if_1(z, a[ub(i00)], (i00-lb(i00)))
418 else
419 let {
420 array[int] of var int: p = eq_encode(i00) %% this needs i00 in ix
421 } %% Faster flattening than (i==i00) @2a9df1f7
422 in
423 forall (i in dom(i00))(
424 aux_int_eq_if_1(z, a[i], p[i])
425 )
426 endif
427 endif;
428
429predicate array_float_element(var int: i00, array[int] of float: a,
430 var float: z) =
431 let { set of int: ix = index_set(a),
432 constraint i00 in { i | i in ix where a[i]>=lb(z) /\ a[i]<=ub(z) },
433 } in %%% Tighten domain of i00 before dMin/dMax
434 let {
435 float: dMin = min(i in dom(i00))(a[i]),
436 float: dMax = max(i in dom(i00))(a[i]),
437 } in
438 if dMin==dMax then
439 z==dMin
440 else
441 z >= dMin /\
442 z <= dMax /\
443 let {
444 int: nUBi00 = max(dom(i00)),
445 int: nLBi00 = min(dom(i00)),
446 float: nMinDist = min(i in nLBi00 .. nUBi00-1)(a[i+1]-a[i]),
447 float: nMaxDist = max(i in nLBi00 .. nUBi00-1)(a[i+1]-a[i]),
448 } in
449 if nMinDist == nMaxDist then %% The linear case
450 z == a[nLBi00] + nMinDist*(i00-nLBi00)
451 else
452 let {
453 array[int] of var int: p = eq_encode(i00) %% this needs i00 in ix
454 } %% Faster flattening than (i==i00) @2a9df1f7
455 in
456 sum(i in dom(i00))( a[i] * p[i] ) == z
457 endif
458 endif;
459
460predicate array_var_float_element(var int: i00, array[int] of var float: a,
461 var float: z) =
462 let { set of int: ix = index_set(a),
463 constraint i00 in { i | i in ix where
464 lb(a[i])<=ub(z) /\ ub(a[i])>=lb(z)
465 },
466 } in %% finish domain first
467 let {
468 float: minLB=min(i in dom(i00))(lb(a[i])),
469 float: maxUB=max(i in dom(i00))(ub(a[i]))
470 } in
471 if minLB==maxUB then
472 z==minLB
473 else
474 z >= minLB /\
475 z <= maxUB /\
476 if {0,1}==dom(i00) /*ub(i00)-lb(i00)==1*/ /*2==card( dom( i00 ) )*/ then
477 aux_float_eq_if_1(z, a[lb(i00)], (ub(i00)-i00)) /\
478 aux_float_eq_if_1(z, a[ub(i00)], (i00-lb(i00)))
479 else
480 %%% The convexified bounds seem slow for ^2 and ^3 equations:
481 % sum(i in dom(i01))( lb(a[i]) * (i==i00) ) <= z /\ %% convexify lower bounds
482 % sum(i in dom(i01))( ub(a[i]) * (i==i00) ) >= z /\ %% convexify upper bounds
483 let {
484 array[int] of var int: p = eq_encode(i00) %% this needs i00 in ix
485 } %% Faster flattening than (i==i00) @2a9df1f7
486 in
487 forall (i in dom(i00))(
488 aux_float_eq_if_1(z, a[i], p[i])
489 )
490 %% Cuts:
491 /\
492 if fElementCutsXZ then
493 array_var_float_element__ROOF([ a[i] | i in dom(i00) ], z) :: MIP_cut %% these 2 as user cuts - too slow
494 /\ array_var_float_element__ROOF([ -a[i] | i in dom(i00) ], -z) :: MIP_cut %% or even skip them
495 else true endif
496 /\
497 if fElementCutsXZB then
498 array_var_float_element__XBZ_lb([ a[i] | i in dom(i00) ], [ p[i] | i in dom(i00) ], z) :: MIP_cut
499 /\ array_var_float_element__XBZ_lb([ -a[i] | i in dom(i00) ], [ p[i] | i in dom(i00) ], -z) :: MIP_cut
500 else true endif
501 endif
502 endif;
503
504%%% Facets on the upper surface of the z-a polytope
505%%% Possible parameter: maximal number of first cuts taken only
506predicate array_var_float_element__ROOF(array[int] of var float: a,
507 var float: z) =
508 let { set of int: ix = index_set(a),
509 int: n = length(a),
510 array[int] of float: AU = [ ub(a[i]) | i in 1..n],
511 array[int] of int: srt_ub = sort_by([i | i in 1..n], AU),
512 %indices of ub sorted up
513 array[int] of float: AU_srt_ub = [ub(a[srt_ub[i]]) | i in 1..n],
514 array[int] of float: AL_srt_ub = [lb(a[srt_ub[i]]) | i in 1..n],
515 array[int] of float: MaxLBFrom =
516 [ max(j in index_set(AL_srt_ub) where j>=i)(AL_srt_ub[j])
517 | i in 1..n ], %% direct, O(n^2)
518 array[int] of float: ULB = [
519 if 1==i then MaxLBFrom[1]
520 else max([AU_srt_ub[i-1], MaxLBFrom[i]])
521 endif | i in 1..n ]
522 } in
523 %%% "ROOF"
524 forall (i in 1..n where
525 if i==n then true else ULB[i]!=ULB[i+1] endif %% not the same base bound
526 )(
527 z <= ULB[i]
528 + sum( j in i..n where AU_srt_ub[i] != AL_srt_ub[i] ) %% not a const
529 ( (AU_srt_ub[j]-ULB[i]) * (a[srt_ub[j]]-AL_srt_ub[j]) / (AU_srt_ub[j]-AL_srt_ub[j]) ) )
530 ;
531
532predicate array_var_float_element__XBZ_lb(array[int] of var float: x, array[int] of var int: b, var float: z) =
533 if fUseXBZCutGen then
534 array_var_float_element__XBZ_lb__cutgen(x, b, z) :: MIP_cut
535 else
536 %% Adding some cuts a priori, also to make solver extract the variables
537 let {
538 int: i1 = min(index_set(x))
539 } in
540 (z <= sum(i in index_set(x))(ub(x[i]) * b[i])) %:: MIP_cut -- does not work to put them here TODO
541 /\
542 forall(i in index_set(x) intersect i1..(i1+19)) %% otherwise too many on amaze2
543 ( assert(lb(x[i]) == -ub(-x[i]) /\ ub(x[i]) == -lb(-x[i]), " negated var's bounds should swap " ) /\
544 z <= x[i] + sum(j in index_set(x) where i!=j)((ub(x[j])-lb(x[i]))*b[j])) %:: MIP_cut %% (ub_j-lb_i) * b_j
545 /\
546 forall(i in index_set(x) intersect i1..(i1+19))
547 ( z <= ub(x[i])*b[i] + sum(j in index_set(x) where i!=j)(x[j]+lb(x[j])*(b[j]-1)) ) %:: MIP_cut
548 /\
549 (z <= sum(i in index_set(x))(x[i] + lb(x[i]) * (b[i]-1))) %:: MIP_cut
550 endif;
551
552%-----------------------------------------------------------------------------%
553% Set constraints
554%% ----------------------------------------------- (NO) SETS ----------------------------------------------
555% XXX only for a fixed set here, general see below.
556% Normally not called because all plugged into the domain.
557% Might be called instead of set_in(x, var set of int s) if s gets fixed?
558predicate set_in(var int: x, set of int: s__) =
559 let {
560 set of int: s = if has_bounds(x) then s__ intersect dom(x) else s__ endif,
561 constraint min(s) <= x,
562 constraint x <= max(s),
563 } in
564 if s = min(s)..max(s) then true
565 elseif fPostprocessDomains then
566 set_in__POST(x, s)
567 else %% Update eq_encode
568 let {
569 array[int] of var int: p = eq_encode(x);
570 } in
571 forall(i in index_set(p) diff s)(p[i]==0)
572% let {
573% array[int] of int: sL = [ e | e in s where not (e - 1 in s) ];
574% array[int] of int: sR = [ e | e in s where not (e + 1 in s) ];
575% array [index_set(sR)] of var 0..1: B;
576% constraint assert(length(sR)==length(sL), "N of lb and ub of sub-intervals of a set should be equal");
577% } in
578% sum(B) = 1 %% use indicators
579% /\
580% x >= sum(i in index_set(sL))(B[i]*sL[i])
581% /\
582% x <= sum(i in index_set(sR))(B[i]*sR[i])
583 endif;
584
585%%% for a fixed set
586predicate set_in_reif(var int: x, set of int: s__, var bool: b) =
587 if is_fixed(b) then
588 if fix(b) then x in s__ else x in dom(x) diff s__ endif
589 elseif has_bounds(x) /\ not (s__ subset dom(x)) then
590 b <-> x in s__ intersect dom(x) %% Use CSE
591 else
592 let {
593 set of int: s = if has_bounds(x) then s__ intersect dom(x) else s__ endif,
594 } in
595 (
596 if dom(x) subset s then b==true
597 elseif card(dom(x) intersect s)==0 then b==false
598 elseif fPostprocessDomains then
599 set_in_reif__POST(x, s, b)
600 %% Bad. Very much so for CBC. 27.06.2019: elseif s == min(s)..max(s) then
601 %% b <-> (min(s) <= x /\ x <= max(s))
602 else
603 if card(dom(x))<=nMZN__UnaryLenMax_setInReif then %% PARAM TODO
604 let {
605 array[int] of var int: p = eq_encode(x);
606 } in
607 sum(i in s intersect dom(x))(p[i]) == bool2int(b)
608 else
609 bool2int(b) == fVarInBigSetOfInt(x, s)
610 endif
611 endif
612 )
613 endif;
614
615 % Alternative
616predicate alt_set_in_reif(var int: x, set of int: s, var bool: b) =
617 b <->
618 exists(i in 1..length([ 0 | e in s where not (e - 1 in s) ]))(
619 let { int: l = [ e | e in s where not (e - 1 in s) ][i],
620 int: r = [ e | e in s where not (e + 1 in s) ][i] }
621 in
622 l <= x /\ x <= r
623 );
624
625%%% for a fixed set
626predicate set_in_imp(var int: x, set of int: s__, var bool: b) =
627 if is_fixed(b) then
628 if fix(b) then x in s__ else true endif
629 elseif has_bounds(x) /\ not (s__ subset dom(x)) then
630 b -> x in s__ intersect dom(x) %% Use CSE
631 else
632 let {
633 set of int: s = if has_bounds(x) then s__ intersect dom(x) else s__ endif,
634 } in
635 (
636 if dom(x) subset s then true
637 elseif card(dom(x) intersect s)==0 then b==false
638 elseif s == min(s)..max(s) then
639 (b -> min(s) <= x) /\ (b -> x <= max(s))
640 else
641 if card(dom(x))<=nMZN__UnaryLenMax_setInReif then %% PARAM TODO
642 let {
643 array[int] of var int: p = eq_encode(x);
644 } in
645 sum(i in s intersect dom(x))(p[i]) >= bool2int(b)
646 else
647 bool2int(b) <= fVarInBigSetOfInt(x, s)
648 endif
649 endif
650 )
651 endif;
652
653function var 0..1: fVarInBigSetOfInt(var int: x, set of int: s)
654 = let {
655 array[int] of int: sL = [ e | e in s where not (e - 1 in s) ];
656 array[int] of int: sR = [ e | e in s where not (e + 1 in s) ];
657 constraint assert(length(sR)==length(sL), "N of lb and ub of sub-intervals of a set should be equal");
658 } in
659 sum(i in index_set(sL)) (bool2int(x>=sL[i] /\ x<=sR[i])); %% use indicators
660
661%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OTHER SET STUFF COMING FROM nosets.mzn %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
662%-----------------------------------------------------------------------------%
663%-----------------------------------------------------------------------------%
664
665annotation bool_search(array[$X] of var bool: x, ann:a1, ann:a2, ann:a3) =
666 let { array[int] of var bool: xx = array1d(x) } in
667 int_search([bool2int(xx[i]) | i in index_set(xx)],a1,a2,a3);
668
669annotation warm_start( array[int] of var bool: x, array[int] of bool: v ) =
670 warm_start( [ bool2int(x[i]) | i in index_set(x) ], [ bool2int(v[i]) | i in index_set(v) ] );
671
672
673annotation sat_goal(var bool: b) = int_max_goal(b);
674annotation int_max_goal(var int: x);
675
676
677%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DOMAIN POSTPROCESSING BUILT-INS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
678 % Single variable: x = d <-> x_eq_d[d]
679predicate equality_encoding__POST(var int: x, array[int] of var int: x_eq_d);
680
681%%%%%%% var int: b: bool2int is a reverse_map, not passed to .fzn
682predicate set_in__POST(var int: x, set of int: s__);
683predicate set_in_reif__POST(var int: x, set of int: s__, var int: b);
684
685
686
687
688%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LOGICAL CONSTRAINTS TO THE SOLVER %%%%%%%%%%%%%%%%%%%%%%%%%%
689%%%%%%% var int: b: bool2int is a reverse_map, not passed to .fzn => REPEAT TESTS. TODO
690predicate int_lin_eq_reif__IND(array[int] of int: c, array[int] of var int: x, int: d, var int: b);
691predicate int_lin_le_reif__IND(array[int] of int: c, array[int] of var int: x, int: d, var int: b);
692predicate int_lin_ne__IND(array[int] of int: c, array[int] of var int: x, int: d);
693predicate aux_int_le_zero_if_0__IND(var int: x, var int: b);
694predicate float_lin_le_reif__IND(array[int] of float: c, array[int] of var float: x, float: d, var int: b);
695predicate aux_float_eq_if_1__IND(var float: x, var float: y, var int: b);
696predicate aux_float_le_zero_if_0__IND(var float: x, var int: b);
697
698predicate array_int_minimum__IND(var int: m, array[int] of var int: x);
699predicate array_int_maximum__IND(var int: m, array[int] of var int: x);
700predicate array_float_minimum__IND(var float: m, array[int] of var float: x);
701predicate array_float_maximum__IND(var float: m, array[int] of var float: x);
702
703%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% XBZ cut generator, currently CPLEX only %%%%%%%%%%%%%%%%%%%%%%%%%%
704predicate array_var_float_element__XBZ_lb__cutgen(array[int] of var float: x, array[int] of var int: b, var float: z);
705
706%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Quadratic expressions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
707predicate fzn_float_times(var float: a, var float: b, var float: c);
708predicate fzn_int_times(var int: a, var int: b, var int: c);
709predicate fzn_array_float_minimum(var float: m, array[int] of var float: x);