dekker.one
The models, scripts, and results of the benchmarks performed for a Half Reification Journal paper
1/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */
2
3/*
4 * Main authors:
5 * Gleb Belov <gleb.belov@monash.edu>
6 */
7
8/* This Source Code Form is subject to the terms of the Mozilla Public
9 * License, v. 2.0. If a copy of the MPL was ! distributed with this
10 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
11
12#ifdef _MSC_VER
13#define _CRT_SECURE_NO_WARNINGS
14#endif
15
16#include <minizinc/MIPdomains.hh>
17#include <minizinc/astexception.hh>
18#include <minizinc/astiterator.hh>
19#include <minizinc/copy.hh>
20#include <minizinc/eval_par.hh>
21#include <minizinc/flatten.hh>
22#include <minizinc/flatten_internal.hh>
23#include <minizinc/hash.hh>
24
25// temporary
26#include <minizinc/prettyprinter.hh>
27//#include <ostream>
28
29#include <map>
30#include <unordered_map>
31#include <unordered_set>
32
33/// TODOs
34/// TODO Not going to work for float vars because of round-offs in the domain interval sorting...
35/// set_in etc. are ! propagated between views
36/// CLEANUP after work: ~destructor
37/// Also check initexpr of all vars? DONE
38/// In case of only_range_domains we'd need to register inequalities
39/// - so better turn that off TODO
40/// CSE for lineq coefs TODO
41
42/// TODO use integer division instead of INT_EPS
43#define INT_EPS 1e-5 // the absolute epsilon for integrality of integer vars.
44
45#define MZN_MIPDOMAINS_PRINTMORESTATS
46#define MZN_DBG_CHECK_ITER_CUTOUT
47
48// #define MZN_DBGOUT_MIPDOMAINS
49#ifdef MZN_DBGOUT_MIPDOMAINS
50#define DBGOUT_MIPD(s) std::cerr << s << std::endl
51#define DBGOUT_MIPD_FLUSH(s) std::cerr << s << std::flush
52#define DBGOUT_MIPD_SELF(op) op
53#else
54#define DBGOUT_MIPD(s) \
55 do { \
56 } while (false)
57#define DBGOUT_MIPD_FLUSH(s) \
58 do { \
59 } while (false)
60#define DBGOUT_MIPD_SELF(op) \
61 do { \
62 } while (false)
63#endif
64
65namespace MiniZinc {
66
67enum EnumStatIdx_MIPD {
68 N_POSTs_all, // N all POSTs in the model
69 N_POSTs_intCmpReif,
70 N_POSTs_floatCmpReif, // in detail
71 N_POSTs_intNE,
72 N_POSTs_floatNE,
73 N_POSTs_setIn,
74 N_POSTs_domain,
75 N_POSTs_setInReif,
76 N_POSTs_eq_encode,
77 N_POSTs_intAux,
78 N_POSTs_floatAux,
79 // Kind of equality connections between involved variables
80 N_POSTs_eq2intlineq,
81 N_POSTs_eq2floatlineq,
82 N_POSTs_int2float,
83 N_POSTs_internalvarredef,
84 N_POSTs_initexpr1id,
85 N_POSTs_initexpr1linexp,
86 N_POSTs_initexprN,
87 N_POSTs_eqNlineq,
88 N_POSTs_eqNmapsize,
89 // other
90 N_POSTs_varsDirect,
91 N_POSTs_varsInvolved,
92 N_POSTs_NSubintvMin,
93 N_POSTs_NSubintvSum,
94 N_POSTs_NSubintvMax, // as N subintervals
95 N_POSTs_SubSizeMin,
96 N_POSTs_SubSizeSum,
97 N_POSTs_SubSizeMax, // subintv. size
98 N_POSTs_linCoefMin,
99 N_POSTs_linCoefMax,
100 N_POSTs_cliquesWithEqEncode,
101 N_POSTs_clEEEnforced,
102 N_POSTs_clEEFound,
103 N_POSTs_size
104};
105extern std::vector<double> MIPD_stats;
106
107enum EnumReifType { RIT_None, RIT_Static, RIT_Reif, RIT_Halfreif };
108enum EnumConstrType { CT_None, CT_Comparison, CT_SetIn, CT_Encode };
109enum EnumCmpType {
110 CMPT_None = 0,
111 CMPT_LE = -4,
112 CMPT_GE = 4,
113 CMPT_EQ = 1,
114 CMPT_NE = 3,
115 CMPT_LT = -5,
116 CMPT_GT = 5,
117 CMPT_LE_0 = -6,
118 CMPT_GE_0 = 6,
119 CMPT_EQ_0 = 2,
120 CMPT_LT_0 = -7,
121 CMPT_GT_0 = 7
122};
123enum EnumVarType { VT_None, VT_Int, VT_Float };
124
125/// struct DomainCallType describes & characterizes a possible domain constr call
126struct DCT {
127 const char* sFuncName = nullptr;
128 const std::vector<Type>& aParams;
129 // unsigned iItem; // call's item number in the flat
130 EnumReifType nReifType = RIT_None; // 0/static/halfreif/reif
131 EnumConstrType nConstrType = CT_None; //
132 EnumCmpType nCmpType = CMPT_None;
133 EnumVarType nVarType = VT_None;
134 FunctionI*& pfi;
135 // double dEps = -1.0;
136 DCT(const char* fn, const std::vector<Type>& prm, EnumReifType er, EnumConstrType ec,
137 EnumCmpType ecmp, EnumVarType ev, FunctionI*& pfi_)
138 : sFuncName(fn),
139 aParams(prm),
140 nReifType(er),
141 nConstrType(ec),
142 nCmpType(ecmp),
143 nVarType(ev),
144 pfi(pfi_) {}
145};
146
147template <class N>
148struct Interval {
149 N left = infMinus(), right = infPlus();
150 mutable VarDecl* varFlag = nullptr;
151 /*constexpr*/ static N infMinus() {
152 return (std::numeric_limits<N>::has_infinity) ? -std::numeric_limits<N>::infinity()
153 : std::numeric_limits<N>::lowest();
154 }
155 /*constexpr*/ static N infPlus() {
156 return (std::numeric_limits<N>::has_infinity) ? std::numeric_limits<N>::infinity()
157 : std::numeric_limits<N>::max();
158 }
159 Interval(N a = infMinus(), N b = infPlus()) : left(a), right(b) {}
160 bool operator<(const Interval& intv) const { return left < intv.left; }
161};
162typedef Interval<double> IntvReal;
163
164template <class N>
165std::ostream& operator<<(std::ostream& os, const Interval<N>& ii) {
166 os << "[ " << ii.left << ", " << ii.right << " ] ";
167 return os;
168}
169
170template <class N>
171class SetOfIntervals : public std::multiset<Interval<N> > {
172public:
173 using Intv = Interval<N>;
174 typedef std::multiset<Interval<N> > Base;
175 typedef typename Base::iterator iterator;
176 SetOfIntervals() : Base() {}
177 SetOfIntervals(std::initializer_list<Interval<N> > il) : Base(il) {}
178 template <class Iter>
179 SetOfIntervals(Iter i1, Iter i2) : Base(i1, i2) {}
180 /// Number of integer values in all the intervals
181 /// Assumes the interval bounds are ints
182 int cardInt() const;
183 /// Max interval length
184 N maxInterval() const;
185 /// Special insert function: check if interval is ok
186 iterator insert(const Interval<N>& iv) {
187 if (iv.left > iv.right) {
188 DBGOUT_MIPD("Interval " << iv.left << ".." << iv.right
189 << " is empty, difference: " << (iv.right - iv.left) << ". Skipping");
190 return Base::end();
191 }
192 return Base::insert(iv);
193 }
194 template <class N1>
195 void intersect(const SetOfIntervals<N1>& s2);
196 /// Assumes open intervals to cut out from closed
197 template <class N1>
198 void cutDeltas(const SetOfIntervals<N1>& s2, N1 delta);
199 template <class N1>
200 void cutDeltas(N1 left, N1 right, N1 delta) {
201 SetOfIntervals<N1> soi;
202 soi.insert(Interval<N1>(left, right));
203 cutDeltas(soi, delta);
204 }
205 /// Cut out an open interval from a set of closed ones (except for infinities)
206 void cutOut(const Interval<N>& intv);
207 typedef std::pair<iterator, iterator> SplitResult;
208 SplitResult split(iterator& it, N pos);
209 bool checkFiniteBounds();
210 /// Check there are no useless interval splittings
211 bool checkDisjunctStrict();
212 Interval<N> getBounds() const;
213 /// Split domain into the integer values
214 /// May assume integer bounds
215 void split2Bits();
216}; // class SetOfIntervals
217typedef SetOfIntervals<double> SetOfIntvReal;
218
219template <class N>
220std::ostream& operator<<(std::ostream& os, const SetOfIntervals<N>& soi) {
221 os << "[[ ";
222 for (auto& ii : soi) {
223 os << "[ " << ii.left << ", " << ii.right;
224 if (ii.varFlag) {
225 os << " @" << ii.varFlag;
226 }
227 os << " ] ";
228 }
229 os << "]]";
230 return os;
231}
232
233template <class Coefs, class Vars>
234class LinEqHelper {
235public:
236 Coefs coefs;
237 Vars vd;
238 double rhs;
239};
240
241template <class Coefs, class Vars>
242static std::ostream& operator<<(std::ostream& os, LinEqHelper<Coefs, Vars>& led) {
243 os << "( [";
244 for (auto c : led.coefs) {
245 os << c << ' ';
246 }
247 os << " ] * [ ";
248 for (auto v : led.vd) {
249 os << v->id()->str() << ' ';
250 }
251 os << " ] ) == " << led.rhs;
252 return os;
253}
254
255typedef LinEqHelper<std::array<double, 2>, std::array<VarDecl*, 2> > LinEq2Vars;
256typedef LinEqHelper<std::vector<double>, std::vector<VarDecl*> > LinEq;
257// struct LinEq2Vars {
258// std::array<double, 2> coefs;
259// std::array<PVarDecl, 2> vd = { { 0, 0 } };
260// double rhs;
261// };
262//
263// struct LinEq {
264// std::vector<double> coefs;
265// std::vector<VarDecl*> vd;
266// double rhs;
267// };
268
269std::vector<double> MIPD_stats(N_POSTs_size);
270
271template <class T>
272static std::vector<T> make_vec(T t1, T t2) {
273 T c_array[] = {t1, t2};
274 std::vector<T> result(c_array, c_array + sizeof(c_array) / sizeof(c_array[0]));
275 return result;
276}
277template <class T>
278static std::vector<T> make_vec(T t1, T t2, T t3) {
279 T c_array[] = {t1, t2, t3};
280 std::vector<T> result(c_array, c_array + sizeof(c_array) / sizeof(c_array[0]));
281 return result;
282}
283template <class T>
284static std::vector<T> make_vec(T t1, T t2, T t3, T t4) {
285 T c_array[] = {t1, t2, t3, t4};
286 std::vector<T> result(c_array, c_array + sizeof(c_array) / sizeof(c_array[0]));
287 return result;
288}
289
290class MIPD {
291public:
292 MIPD(Env* env, bool fV, int nmi, double dmd)
293 : nMaxIntv2Bits(nmi), dMaxNValueDensity(dmd), _env(env) {
294 getEnv();
295 fVerbose = fV;
296 }
297 static bool fVerbose;
298 const int nMaxIntv2Bits = 0; // Maximal interval length to enforce equality encoding
299 const double dMaxNValueDensity = 3.0; // Maximal ratio cardInt() / size() of a domain
300 // to enforce ee
301 bool doMIPdomains() {
302 MIPD_stats[N_POSTs_NSubintvMin] = 1e100;
303 MIPD_stats[N_POSTs_SubSizeMin] = 1e100;
304
305 if (!registerLinearConstraintDecls()) {
306 return true;
307 }
308 if (!registerPOSTConstraintDecls()) { // not declared => no conversions
309 return true;
310 }
311 registerPOSTVariables();
312 if (_vVarDescr.empty()) {
313 return true;
314 }
315 constructVarViewCliques();
316 if (!decomposeDomains()) {
317 return false;
318 }
319 if (fVerbose) {
320 printStats(std::cerr);
321 }
322 return true;
323 }
324
325private:
326 Env* _env = nullptr;
327 Env* getEnv() {
328 MZN_MIPD_assert_hard(_env);
329 return _env;
330 }
331
332 typedef VarDecl* PVarDecl;
333
334 // NOLINTNEXTLINE(readability-identifier-naming)
335 FunctionI* int_lin_eq;
336 // NOLINTNEXTLINE(readability-identifier-naming)
337 FunctionI* int_lin_le;
338 // NOLINTNEXTLINE(readability-identifier-naming)
339 FunctionI* float_lin_eq;
340 // NOLINTNEXTLINE(readability-identifier-naming)
341 FunctionI* float_lin_le;
342 // NOLINTNEXTLINE(readability-identifier-naming)
343 FunctionI* int2float;
344 // NOLINTNEXTLINE(readability-identifier-naming)
345 FunctionI* lin_exp_int;
346 // NOLINTNEXTLINE(readability-identifier-naming)
347 FunctionI* lin_exp_float;
348
349 // NOLINTNEXTLINE(readability-identifier-naming)
350 std::vector<Type> int_lin_eq_t = make_vec(Type::parint(1), Type::varint(1), Type::parint());
351 // NOLINTNEXTLINE(readability-identifier-naming)
352 std::vector<Type> float_lin_eq_t =
353 make_vec(Type::parfloat(1), Type::varfloat(1), Type::parfloat());
354 // NOLINTNEXTLINE(readability-identifier-naming)
355 std::vector<Type> t_VIVF = make_vec(Type::varint(), Type::varfloat());
356
357 // double float_lt_EPS_coef_ = 1e-5;
358
359 bool registerLinearConstraintDecls() {
360 EnvI& env = getEnv()->envi();
361 GCLock lock;
362
363 int_lin_eq = env.model->matchFn(env, constants().ids.int_.lin_eq, int_lin_eq_t, false);
364 DBGOUT_MIPD(" int_lin_eq = " << int_lin_eq);
365 // MZN_MIPD_assert_hard(fi);
366 // int_lin_eq = (fi && fi->e()) ? fi : NULL;
367 int_lin_le = env.model->matchFn(env, constants().ids.int_.lin_le, int_lin_eq_t, false);
368 float_lin_eq = env.model->matchFn(env, constants().ids.float_.lin_eq, float_lin_eq_t, false);
369 float_lin_le = env.model->matchFn(env, constants().ids.float_.lin_le, float_lin_eq_t, false);
370 int2float = env.model->matchFn(env, constants().ids.int2float, t_VIVF, false);
371
372 lin_exp_int = env.model->matchFn(env, constants().ids.lin_exp, int_lin_eq_t, false);
373 lin_exp_float = env.model->matchFn(env, constants().ids.lin_exp, float_lin_eq_t, false);
374
375 return (int_lin_eq != nullptr) && (int_lin_le != nullptr) && (float_lin_eq != nullptr) &&
376 (float_lin_le != nullptr);
377 // say something...
378
379 // std::cerr << " lin_exp_int=" << lin_exp_int << std::endl;
380 // std::cerr << " lin_exp_float=" << lin_exp_float << std::endl;
381 // For this to work, need to define a function, see mzn_only_range_domains()
382 // {
383 // GCLock lock;
384 // Call* call_EPS_for_LT =
385 // new Call(Location(),"mzn_float_lt_EPS_coef__", std::vector<Expression*>());
386 // call_EPS_for_LT->type(Type::parfloat());
387 // call_EPS_for_LT->decl(env.model->matchFn(getEnv()->envi(), call_EPS_for_LT));
388 // float_lt_EPS_coef_ = eval_float(getEnv()->envi(), call_EPS_for_LT);
389 // }
390 }
391 // bool matchAndMarkFunction();
392 // std::set<FunctionI*> funcs;
393
394 /// Possible function param sets
395 // NOLINTNEXTLINE(readability-identifier-naming)
396 std::vector<Type> t_VII = make_vec(Type::varint(), Type::parint());
397 // NOLINTNEXTLINE(readability-identifier-naming)
398 std::vector<Type> t_VIVI = make_vec(Type::varint(), Type::varint());
399 // NOLINTNEXTLINE(readability-identifier-naming)
400 std::vector<Type> t_VIIVI = make_vec(Type::varint(), Type::parint(), Type::varint());
401 // NOLINTNEXTLINE(readability-identifier-naming)
402 std::vector<Type> t_VFVI = make_vec(Type::varfloat(), Type::varint());
403 // NOLINTNEXTLINE(readability-identifier-naming)
404 std::vector<Type> t_VFVF = make_vec(Type::varfloat(), Type::varfloat());
405 // NOLINTNEXTLINE(readability-identifier-naming)
406 std::vector<Type> t_VFFVI = make_vec(Type::varfloat(), Type::parfloat(), Type::varint());
407 // NOLINTNEXTLINE(readability-identifier-naming)
408 std::vector<Type> t_VFFVIF =
409 make_vec(Type::varfloat(), Type::parfloat(), Type::varint(), Type::parfloat());
410 // NOLINTNEXTLINE(readability-identifier-naming)
411 std::vector<Type> t_VFVIF = make_vec(Type::varfloat(), Type::varint(), Type::parfloat());
412 // NOLINTNEXTLINE(readability-identifier-naming)
413 std::vector<Type> t_VFVIVF = make_vec(Type::varfloat(), Type::varint(), Type::varfloat());
414 // NOLINTNEXTLINE(readability-identifier-naming)
415 std::vector<Type> t_VFVIVFF =
416 make_vec(Type::varfloat(), Type::varint(), Type::varfloat(), Type::parfloat());
417 // NOLINTNEXTLINE(readability-identifier-naming)
418 std::vector<Type> t_VFVFF = make_vec(Type::varfloat(), Type::varfloat(), Type::parfloat());
419 // NOLINTNEXTLINE(readability-identifier-naming)
420 std::vector<Type> t_VFFF = make_vec(Type::varfloat(), Type::parfloat(), Type::parfloat());
421 // std::vector<Type> t_VFVFVIF({ Type::varfloat(), Type::varfloat(), Type::varint(),
422 // Type::parfloat() });
423
424 // NOLINTNEXTLINE(readability-identifier-naming)
425 std::vector<Type> t_VIAVI = make_vec(Type::varint(), Type::varint(1));
426 // NOLINTNEXTLINE(readability-identifier-naming)
427 std::vector<Type> t_VISI = make_vec(Type::varint(), Type::parsetint());
428 // NOLINTNEXTLINE(readability-identifier-naming)
429 std::vector<Type> t_VISIVI = make_vec(Type::varint(), Type::parsetint(), Type::varint());
430
431 // std::vector<Type> t_intarray(1);
432 // t_intarray[0] = Type::parint(-1);
433
434 typedef std::unordered_map<FunctionI*, DCT*> M_POSTCallTypes;
435 M_POSTCallTypes _mCallTypes; // actually declared in the input
436 std::vector<DCT> _aCT; // all possible
437
438 // Fails:
439 // DomainCallType a = { NULL, t_VII, RIT_Halfreif, CT_Comparison, CMPT_EQ, VT_Float };
440
441 /// struct VarDescr stores some info about variables involved in domain constr
442 struct VarDescr {
443 typedef unsigned char boolShort;
444 VarDescr(VarDecl* vd_, boolShort fi, double l_ = 0.0, double u_ = 0.0)
445 : lb(l_), ub(u_), vd(vd_), fInt(fi) {}
446 double lb, ub;
447 VarDecl* vd = nullptr;
448 int nClique = -1; // clique number
449 // std::vector<Call*> aCalls;
450 std::vector<ConstraintI*> aCalls;
451 boolShort fInt = 0;
452 ConstraintI* pEqEncoding = nullptr;
453 boolShort fDomainConstrProcessed = 0;
454 // boolShort fPropagatedViews=0;
455 // boolShort fPropagatedLargerEqns=0;
456 };
457
458 std::vector<VarDescr> _vVarDescr;
459
460 // NOLINTNEXTLINE(readability-identifier-naming)
461 FunctionI* int_le_reif_POST = nullptr;
462 // NOLINTNEXTLINE(readability-identifier-naming)
463 FunctionI* int_ge_reif_POST = nullptr;
464 // NOLINTNEXTLINE(readability-identifier-naming)
465 FunctionI* int_eq_reif_POST = nullptr;
466 // NOLINTNEXTLINE(readability-identifier-naming)
467 FunctionI* int_ne_POST = nullptr;
468 // NOLINTNEXTLINE(readability-identifier-naming)
469 FunctionI* float_le_reif_POST = nullptr;
470 // NOLINTNEXTLINE(readability-identifier-naming)
471 FunctionI* float_ge_reif_POST = nullptr;
472 // NOLINTNEXTLINE(readability-identifier-naming)
473 FunctionI* aux_float_lt_zero_iff_1_POST = nullptr;
474 // NOLINTNEXTLINE(readability-identifier-naming)
475 FunctionI* float_eq_reif_POST = nullptr;
476 // NOLINTNEXTLINE(readability-identifier-naming)
477 FunctionI* float_ne_POST = nullptr;
478 // NOLINTNEXTLINE(readability-identifier-naming)
479 FunctionI* aux_float_eq_zero_if_1_POST = nullptr;
480 // NOLINTNEXTLINE(readability-identifier-naming)
481 FunctionI* aux_int_le_zero_if_1_POST = nullptr;
482 // NOLINTNEXTLINE(readability-identifier-naming)
483 FunctionI* aux_float_le_zero_if_1_POST = nullptr;
484 // NOLINTNEXTLINE(readability-identifier-naming)
485 FunctionI* aux_float_lt_zero_if_1_POST = nullptr;
486 // NOLINTNEXTLINE(readability-identifier-naming)
487 FunctionI* equality_encoding_POST = nullptr;
488 // NOLINTNEXTLINE(readability-identifier-naming)
489 FunctionI* set_in_POST = nullptr;
490 // NOLINTNEXTLINE(readability-identifier-naming)
491 FunctionI* set_in_reif_POST = nullptr;
492
493 bool registerPOSTConstraintDecls() {
494 EnvI& env = getEnv()->envi();
495 GCLock lock;
496
497 _aCT.clear();
498 _aCT.emplace_back("int_le_reif__POST", t_VIIVI, RIT_Reif, CT_Comparison, CMPT_LE, VT_Int,
499 int_le_reif_POST);
500 _aCT.emplace_back("int_ge_reif__POST", t_VIIVI, RIT_Reif, CT_Comparison, CMPT_GE, VT_Int,
501 int_ge_reif_POST);
502 _aCT.emplace_back("int_eq_reif__POST", t_VIIVI, RIT_Reif, CT_Comparison, CMPT_EQ, VT_Int,
503 int_eq_reif_POST);
504 _aCT.emplace_back("int_ne__POST", t_VII, RIT_Static, CT_Comparison, CMPT_NE, VT_Int,
505 int_ne_POST);
506
507 _aCT.emplace_back("float_le_reif__POST", t_VFFVIF, RIT_Reif, CT_Comparison, CMPT_LE, VT_Float,
508 float_le_reif_POST);
509 _aCT.emplace_back("float_ge_reif__POST", t_VFFVIF, RIT_Reif, CT_Comparison, CMPT_GE, VT_Float,
510 float_ge_reif_POST);
511 _aCT.emplace_back("aux_float_lt_zero_iff_1__POST", t_VFVIF, RIT_Reif, CT_Comparison, CMPT_LT,
512 VT_Float, aux_float_lt_zero_iff_1_POST);
513 _aCT.emplace_back("float_eq_reif__POST", t_VFFVIF, RIT_Reif, CT_Comparison, CMPT_EQ, VT_Float,
514 float_eq_reif_POST);
515 _aCT.emplace_back("float_ne__POST", t_VFFF, RIT_Static, CT_Comparison, CMPT_NE, VT_Float,
516 float_ne_POST);
517
518 _aCT.emplace_back("aux_float_eq_zero_if_1__POST", t_VFVIVF, RIT_Halfreif, CT_Comparison,
519 CMPT_EQ_0, VT_Float, aux_float_eq_zero_if_1_POST);
520 _aCT.emplace_back("aux_int_le_zero_if_1__POST", t_VIVI, RIT_Halfreif, CT_Comparison, CMPT_LE_0,
521 VT_Int, aux_int_le_zero_if_1_POST);
522 _aCT.emplace_back("aux_float_le_zero_if_1__POST", t_VFVIVF, RIT_Halfreif, CT_Comparison,
523 CMPT_LE_0, VT_Float, aux_float_le_zero_if_1_POST);
524 _aCT.emplace_back("aux_float_lt_zero_if_1__POST", t_VFVIVFF, RIT_Halfreif, CT_Comparison,
525 CMPT_LT_0, VT_Float, aux_float_lt_zero_if_1_POST);
526
527 _aCT.emplace_back("equality_encoding__POST", t_VIAVI, RIT_Static, CT_Encode, CMPT_None, VT_Int,
528 equality_encoding_POST);
529 _aCT.emplace_back("set_in__POST", t_VISI, RIT_Static, CT_SetIn, CMPT_None, VT_Int, set_in_POST);
530 _aCT.emplace_back("set_in_reif__POST", t_VISIVI, RIT_Reif, CT_SetIn, CMPT_None, VT_Int,
531 set_in_reif_POST);
532 /// Registering all declared & compatible _POST constraints
533 /// (First, cleanup FunctionIs' payload: -- ! doing now)
534 for (int i = 0; i < _aCT.size(); ++i) {
535 FunctionI* fi = env.model->matchFn(env, ASTString(_aCT[i].sFuncName), _aCT[i].aParams, false);
536 if (fi != nullptr) {
537 _mCallTypes[fi] = _aCT.data() + i;
538 _aCT[i].pfi = fi;
539 // fi->pPayload = (void*)this;
540 // std::cerr << " FOund declaration: " << _aCT[i].sFuncName << std::endl;
541 } else {
542 _aCT[i].pfi = nullptr;
543 DBGOUT_MIPD(" MIssing declaration: " << _aCT[i].sFuncName);
544 return false;
545 }
546 }
547 return true;
548 }
549
550 /// Registering all _POST calls' domain-constrained variables
551 void registerPOSTVariables() {
552 EnvI& env = getEnv()->envi();
553 GCLock lock;
554 Model& mFlat = *getEnv()->flat();
555 // First, cleanup VarDecls' payload which stores index in _vVarDescr
556 for (VarDeclIterator ivd = mFlat.vardecls().begin(); ivd != mFlat.vardecls().end(); ++ivd) {
557 ivd->e()->payload(-1);
558 }
559 // Now add variables with non-contiguous domain
560 for (VarDeclIterator ivd = mFlat.vardecls().begin(); ivd != mFlat.vardecls().end(); ++ivd) {
561 VarDecl* vd0 = ivd->e();
562 bool fNonCtg = false;
563 if (vd0->type().isint()) { // currently only for int vars TODO
564 if (Expression* eDom = vd0->ti()->domain()) {
565 IntSetVal* dom = eval_intset(env, eDom);
566 fNonCtg = (dom->size() > 1);
567 }
568 }
569 if (fNonCtg) {
570 DBGOUT_MIPD(" Variable " << vd0->id()->str() << ": non-contiguous domain "
571 << (*(vd0->ti()->domain())));
572 if (vd0->payload() == -1) { // ! yet visited
573 vd0->payload(static_cast<int>(_vVarDescr.size()));
574 _vVarDescr.emplace_back(vd0, vd0->type().isint()); // can use /prmTypes/ as well
575 if (vd0->e() != nullptr) {
576 checkInitExpr(vd0);
577 }
578 } else {
579 DBGOUT_MIPD_FLUSH(" (already touched)");
580 }
581 ++MIPD_stats[N_POSTs_domain];
582 ++MIPD_stats[N_POSTs_all];
583 }
584 }
585 // Iterate thru original _POST constraints to mark constrained vars:
586 for (ConstraintIterator ic = mFlat.constraints().begin(); ic != mFlat.constraints().end();
587 ++ic) {
588 if (ic->removed()) {
589 continue;
590 }
591 if (Call* c = ic->e()->dynamicCast<Call>()) {
592 auto ipct = _mCallTypes.find(c->decl());
593 if (ipct != _mCallTypes.end()) {
594 // No ! here because might be deleted immediately in later versions.
595 // ic->remove(); // mark removed at once
596 MZN_MIPD_assert_hard(c->argCount() > 1);
597 ++MIPD_stats[N_POSTs_all];
598 VarDecl* vd0 = expr2VarDecl(c->arg(0));
599 if (nullptr == vd0) {
600 DBGOUT_MIPD_FLUSH(" Call " << *c
601 << ": 1st arg not a VarDecl, removing if eq_encoding...");
602 /// Only allow literals as main argument for equality_encoding
603 if (equality_encoding_POST ==
604 ipct->first) { // was MZN_MIPD_assert_hard before MZN 2017
605 ic->remove();
606 }
607 continue; // ignore this call
608 }
609 DBGOUT_MIPD_FLUSH(" Call " << c->id().str() << " uses variable " << vd0->id()->str());
610 if (vd0->payload() == -1) { // ! yet visited
611 vd0->payload(static_cast<int>(_vVarDescr.size()));
612 _vVarDescr.emplace_back(vd0, vd0->type().isint()); // can use /prmTypes/ as well
613 // bounds/domains later for each involved var TODO
614 if (vd0->e() != nullptr) {
615 checkInitExpr(vd0);
616 }
617 } else {
618 DBGOUT_MIPD_FLUSH(" (already touched)");
619 }
620 DBGOUT_MIPD("");
621 if (equality_encoding_POST == c->decl()) {
622 MZN_MIPD_assert_hard(!_vVarDescr[vd0->payload()].pEqEncoding);
623 _vVarDescr[vd0->payload()].pEqEncoding = &*ic;
624 DBGOUT_MIPD(" Variable " << vd0->id()->str() << " has eq_encode.");
625 } // + if has aux_ constraints?
626 else {
627 _vVarDescr[vd0->payload()].aCalls.push_back(&*ic);
628 }
629 }
630 }
631 }
632 MIPD_stats[N_POSTs_varsDirect] = static_cast<double>(_vVarDescr.size());
633 }
634
635 // Should only be called on a newly added variable
636 // OR when looking thru all non-touched vars
637 /// Checks init expr of a variable
638 /// Return true IFF new connection
639 /// The bool param requires RHS to be POST-touched
640 // Guido: can! be recursive in FZN
641 bool checkInitExpr(VarDecl* vd, bool fCheckArg = false) {
642 MZN_MIPD_assert_hard(vd->e());
643 if (!vd->type().isint() && !vd->type().isfloat()) {
644 return false;
645 }
646 if (!fCheckArg) {
647 MZN_MIPD_assert_hard(vd->payload() >= 0);
648 }
649 if (Id* id = vd->e()->dynamicCast<Id>()) {
650 // const int f1 = ( vd->payload()>=0 );
651 // const int f2 = ( id->decl()->payload()>=0 );
652 if (!fCheckArg || (id->decl()->payload() >= 0)) {
653 DBGOUT_MIPD_FLUSH(" Checking init expr ");
654 DBGOUT_MIPD_SELF(debugprint(vd));
655 LinEq2Vars led;
656 // FAILS:
657 // led.vd = { vd, expr2VarDecl(id->decl()->e()) };
658 led.vd = {{vd, expr2VarDecl(vd->e())}};
659 led.coefs = {{1.0, -1.0}};
660 led.rhs = 0.0;
661 put2VarsConnection(led, false);
662 ++MIPD_stats[N_POSTs_initexpr1id];
663 if (id->decl()->e() != nullptr) { // no initexpr for initexpr FAILS on cc-base.mzn
664 checkInitExpr(id->decl());
665 }
666 return true; // in any case
667 }
668 } else if (Call* c = vd->e()->dynamicCast<Call>()) {
669 if (lin_exp_int == c->decl() || lin_exp_float == c->decl()) {
670 // std::cerr << " !E call " << std::flush;
671 // debugprint(c);
672 MZN_MIPD_assert_hard(c->argCount() == 3);
673 // ArrayLit* al = c->args()[1]->dynamicCast<ArrayLit>();
674 auto* al = follow_id(c->arg(1))->cast<ArrayLit>();
675 MZN_MIPD_assert_hard(al);
676 MZN_MIPD_assert_hard(al->size() >= 1);
677 if (al->size() == 1) { // 1-term scalar product in the rhs
678 LinEq2Vars led;
679 led.vd = {{vd, expr2VarDecl((*al)[0])}};
680 // const int f1 = ( vd->payload()>=0 );
681 // const int f2 = ( led.vd[1]->payload()>=0 );
682 if (!fCheckArg || (led.vd[1]->payload() >= 0)) {
683 // Can use a!her map here:
684 // if ( _sCallLinEq2.end() != _sCallLinEq2.find(c) )
685 // continue;
686 // _sCallLinEq2.insert(c); // memorize this call
687 DBGOUT_MIPD_FLUSH(" REG 1-LINEXP ");
688 DBGOUT_MIPD_SELF(debugprint(vd));
689 std::array<double, 1> coef0;
690 expr2Array(c->arg(0), coef0);
691 led.coefs = {{-1.0, coef0[0]}};
692 led.rhs = -expr2Const(c->arg(2)); // MINUS
693 put2VarsConnection(led, false);
694 ++MIPD_stats[N_POSTs_initexpr1linexp];
695 if (led.vd[1]->e() != nullptr) { // no initexpr for initexpr FAILS TODO
696 checkInitExpr(led.vd[1]);
697 }
698 return true; // in any case
699 }
700 } else if (true) { // NOLINT: check larger views always. OK? TODO
701 // if ( vd->payload()>=0 ) { // larger views
702 // TODO should be here?
703 // std::cerr << " LE_" << al->v().size() << ' ' << std::flush;
704 DBGOUT_MIPD(" REG N-LINEXP ");
705 DBGOUT_MIPD_SELF(debugprint(vd));
706 // Checking all but adding only touched defined vars?
707 return findOrAddDefining(vd->id(), c);
708 }
709 }
710 }
711 return false;
712 }
713
714 /// Build var cliques (i.e. of var pairs viewing each other)
715 void constructVarViewCliques() {
716 // std::cerr << " Model: " << std::endl;
717 // debugprint(getEnv()->flat());
718
719 // TAgenda agenda(_vVarDescr.size()), agendaNext;
720 // for ( int i=0; i<agenda.size(); ++i )
721 // agenda[i] = i;
722 bool fChanges;
723 do {
724 fChanges = false;
725 propagateViews(fChanges);
726 propagateImplViews(fChanges);
727 } while (fChanges);
728
729 MIPD_stats[N_POSTs_varsInvolved] = static_cast<double>(_vVarDescr.size());
730 }
731
732 void propagateViews(bool& fChanges) {
733 // EnvI& env = getEnv()->envi();
734 GCLock lock;
735
736 // Iterate thru original 2-variable equalities to mark views:
737 Model& mFlat = *getEnv()->flat();
738
739 DBGOUT_MIPD(" CHECK ALL INITEXPR if they access a touched variable:");
740 for (VarDeclIterator ivd = mFlat.vardecls().begin(); ivd != mFlat.vardecls().end(); ++ivd) {
741 if (ivd->removed()) {
742 continue;
743 }
744 if ((ivd->e()->e() != nullptr) && ivd->e()->payload() < 0 // untouched
745 && (ivd->e()->type().isint() || ivd->e()->type().isfloat())) { // scalars
746 if (checkInitExpr(ivd->e(), true)) {
747 fChanges = true;
748 }
749 }
750 }
751
752 DBGOUT_MIPD(" CHECK ALL CONSTRAINTS for 2-var equations:");
753 for (ConstraintIterator ic = mFlat.constraints().begin(); ic != mFlat.constraints().end();
754 ++ic) {
755 // std::cerr << " SEE constraint: " << " ";
756 // debugprint(&*ic);
757 // debugprint(c->decl());
758 if (ic->removed()) {
759 continue;
760 }
761 if (Call* c = ic->e()->dynamicCast<Call>()) {
762 const bool fIntLinEq = int_lin_eq == c->decl();
763 const bool fFloatLinEq = float_lin_eq == c->decl();
764 if (fIntLinEq || fFloatLinEq) {
765 // std::cerr << " !E call " << std::flush;
766 // debugprint(c);
767 MZN_MIPD_assert_hard(c->argCount() == 3);
768 auto* al = follow_id(c->arg(1))->cast<ArrayLit>();
769 MZN_MIPD_assert_hard(al);
770 if (al->size() == 2) { // 2-term eqn
771 LinEq2Vars led;
772 expr2DeclArray(c->arg(1), led.vd);
773 // At least 1 touched var:
774 if (nullptr != led.vd[0] && nullptr != led.vd[1]) {
775 if (led.vd[0]->payload() >= 0 || led.vd[1]->payload() >= 0) {
776 if (_sCallLinEq2.end() != _sCallLinEq2.find(c)) {
777 continue;
778 }
779 _sCallLinEq2.insert(c); // memorize this call
780 DBGOUT_MIPD(" REG 2-call ");
781 DBGOUT_MIPD_SELF(debugprint(c));
782 led.rhs = expr2Const(c->arg(2));
783 expr2Array(c->arg(0), led.coefs);
784 MZN_MIPD_assert_hard(2 == led.coefs.size());
785 fChanges = true;
786 put2VarsConnection(led);
787 ++MIPD_stats[fIntLinEq ? N_POSTs_eq2intlineq : N_POSTs_eq2floatlineq];
788 }
789 }
790 } else if (al->size() == 1) {
791 static int nn = 0;
792 if (++nn <= 7) {
793 std::cerr << " MIPD: LIN_EQ with 1 variable::: " << std::flush;
794 std::cerr << (*c) << std::endl;
795 }
796 } else { // larger eqns
797 // TODO should be here?
798 auto* eVD = get_annotation(c->ann(), constants().ann.defines_var);
799 if (eVD != nullptr) {
800 if (_sCallLinEqN.end() != _sCallLinEqN.find(c)) {
801 continue;
802 }
803 _sCallLinEqN.insert(c); // memorize this call
804 DBGOUT_MIPD(" REG N-call ");
805 DBGOUT_MIPD_SELF(debugprint(c));
806 Call* pC = eVD->dynamicCast<Call>();
807 MZN_MIPD_assert_hard(pC);
808 MZN_MIPD_assert_hard(pC->argCount());
809 // Checking all but adding only touched defined vars? Seems too long.
810 VarDecl* vd = expr2VarDecl(pC->arg(0));
811 if ((vd != nullptr) && vd->payload() >= 0) { // only if touched
812 if (findOrAddDefining(pC->arg(0), c)) {
813 fChanges = true;
814 }
815 }
816 }
817 }
818 } else
819 // const bool fI2F = (int2float==c->decl());
820 // const bool fIVR = (constants().varRedef==c->decl());
821 // if ( fI2F || fIVR ) {
822 if (int2float == c->decl() || constants().varRedef == c->decl()) {
823 // std::cerr << " !E call " << std::flush;
824 // debugprint(c);
825 MZN_MIPD_assert_hard(c->argCount() == 2);
826 LinEq2Vars led;
827 // led.vd.resize(2);
828 led.vd[0] = expr2VarDecl(c->arg(0));
829 led.vd[1] = expr2VarDecl(c->arg(1));
830 // At least 1 touched var:
831 if (led.vd[0]->payload() >= 0 || led.vd[1]->payload() >= 0) {
832 if (_sCallInt2Float.end() != _sCallInt2Float.find(c)) {
833 continue;
834 }
835 _sCallInt2Float.insert(c); // memorize this call
836 DBGOUT_MIPD(" REG call ");
837 DBGOUT_MIPD_SELF(debugprint(c));
838 led.rhs = 0.0;
839 led.coefs = {{1.0, -1.0}};
840 fChanges = true;
841 put2VarsConnection(led);
842 ++MIPD_stats[int2float == c->decl() ? N_POSTs_int2float : N_POSTs_internalvarredef];
843 }
844 }
845 }
846 }
847 }
848
849 /// This vector stores the linear part of a general view
850 /// x = <linear part> + rhs
851 typedef std::vector<std::pair<VarDecl*, float> > TLinExpLin;
852 /// This struct has data describing the rest of a general view
853 struct NViewData {
854 VarDecl* pVarDefined = nullptr;
855 double coef0 = 1.0;
856 double rhs;
857 };
858 typedef std::map<TLinExpLin, NViewData> NViewMap;
859 NViewMap _mNViews;
860
861 /// compare to an existing defining linexp, || just add it to the map
862 /// adds only touched defined vars
863 /// return true iff new linear connection
864 // linexp: z = a^T x+b
865 // _lin_eq: a^T x == b
866 bool findOrAddDefining(Expression* exp, Call* pC) {
867 Id* pId = exp->dynamicCast<Id>();
868 MZN_MIPD_assert_hard(pId);
869 VarDecl* vd = pId->decl();
870 MZN_MIPD_assert_hard(vd);
871 MZN_MIPD_assert_hard(pC->argCount() == 3);
872
873 TLinExpLin rhsLin;
874 NViewData nVRest;
875 nVRest.pVarDefined = vd;
876 nVRest.rhs = expr2Const(pC->arg(2));
877
878 std::vector<VarDecl*> vars;
879 expr2DeclArray(pC->arg(1), vars);
880 std::vector<double> coefs;
881 expr2Array(pC->arg(0), coefs);
882 MZN_MIPD_assert_hard(vars.size() == coefs.size());
883
884 int nVD = 0;
885 for (int i = 0; i < vars.size(); ++i) {
886 // MZN_MIPD_assert_hard( 0.0!=std::fabs
887 if (vd ==
888 vars[i]) { // when int/float_lin_eq :: defines_var(vd) "Recursive definition of " << *vd
889 nVRest.coef0 = -coefs[i];
890 nVRest.rhs = -nVRest.rhs;
891 ++nVD;
892 } else {
893 rhsLin.push_back(std::make_pair(vars[i], coefs[i]));
894 }
895 }
896 MZN_MIPD_assert_hard(1 >= nVD);
897 std::sort(rhsLin.begin(), rhsLin.end());
898
899 // Divide the equation by the 1st coef
900 const double coef1 = rhsLin.begin()->second;
901 MZN_MIPD_assert_hard(0.0 != std::fabs(coef1));
902 nVRest.coef0 /= coef1;
903 nVRest.rhs /= coef1;
904 for (auto& rhsL : rhsLin) {
905 rhsL.second /= coef1;
906 }
907
908 auto it = _mNViews.find(rhsLin);
909 if (_mNViews.end() != it &&
910 nVRest.pVarDefined != it->second.pVarDefined) { // don't connect to itself
911 LinEq2Vars leq;
912 leq.vd = {{nVRest.pVarDefined, it->second.pVarDefined}};
913 leq.coefs = {{nVRest.coef0, -it->second.coef0}}; // +, -
914 leq.rhs = nVRest.rhs - it->second.rhs;
915 put2VarsConnection(leq, false);
916 ++MIPD_stats[nVD != 0 ? N_POSTs_eqNlineq : N_POSTs_initexprN];
917 return true;
918 }
919 if (vd->payload() >= 0) { // only touched
920 _mNViews[rhsLin] = nVRest;
921 return true; // can lead to a new connection
922 }
923
924 return false;
925 }
926
927 static void propagateImplViews(bool& fChanges) {
928 // EnvI& env = getEnv()->envi();
929 GCLock lock;
930
931 // TODO
932 }
933
934 /// Could be better to mark the calls instead:
935 std::unordered_set<Call*> _sCallLinEq2, _sCallInt2Float, _sCallLinEqN;
936
937 class TClique : public std::vector<LinEq2Vars> { // need more info?
938 public:
939 /// This function takes the 1st variable && relates all to it
940 /// Return false if contrad / disconnected graph
941 // bool findRelations0() {
942 // return true;
943 // }
944 };
945 typedef std::vector<TClique> TCLiqueList;
946 TCLiqueList _aCliques;
947
948 /// register a 2-variable lin eq
949 /// add it to the var clique, joining the participants' cliques if needed
950 void put2VarsConnection(LinEq2Vars& led, bool fCheckinitExpr = true) {
951 MZN_MIPD_assert_hard(led.coefs.size() == led.vd.size());
952 MZN_MIPD_assert_hard(led.vd.size() == 2);
953 DBGOUT_MIPD_FLUSH(" Register 2-var connection: " << led);
954 /// Check it's not same 2 vars
955 if (led.vd[0] == led.vd[1]) {
956 MZN_MIPD_assert_soft(
957 0, "MIPD: STRANGE: registering var connection to itself: " << led << ", skipping");
958 MZN_MIPD_ASSERT_FOR_SAT(fabs(led.coefs[0] + led.coefs[1]) < 1e-6, // TODO param
959 getEnv()->envi(), led.vd[0]->loc(),
960 "Var connection to itself seems to indicate UNSAT: " << led);
961 return;
962 }
963 // register if new variables
964 // std::vector<bool> fHaveClq(led.vd.size(), false);
965 int nCliqueAvailable = -1;
966 for (auto* vd : led.vd) {
967 if (vd->payload() < 0) { // ! yet visited
968 vd->payload(static_cast<int>(_vVarDescr.size()));
969 _vVarDescr.emplace_back(vd, vd->type().isint()); // can use /prmTypes/ as well
970 if (fCheckinitExpr && (vd->e() != nullptr)) {
971 checkInitExpr(vd);
972 }
973 } else {
974 int nMaybeClq = _vVarDescr[vd->payload()].nClique;
975 if (nMaybeClq >= 0) {
976 nCliqueAvailable = nMaybeClq;
977 }
978 // MZN_MIPD_assert_hard( nCliqueAvailable>=0 );
979 // fHaveClq[i] = true;
980 }
981 }
982 if (nCliqueAvailable < 0) { // no clique found
983 nCliqueAvailable = static_cast<int>(_aCliques.size());
984 _aCliques.resize(_aCliques.size() + 1);
985 }
986 DBGOUT_MIPD(" ...adding to clique " << nCliqueAvailable << " of size "
987 << _aCliques[nCliqueAvailable].size());
988 TClique& clqNew = _aCliques[nCliqueAvailable];
989 clqNew.push_back(led);
990 for (auto* vd : led.vd) { // merging cliques
991 int& nMaybeClq = _vVarDescr[vd->payload()].nClique;
992 if (nMaybeClq >= 0 && nMaybeClq != nCliqueAvailable) {
993 TClique& clqOld = _aCliques[nMaybeClq];
994 MZN_MIPD_assert_hard(clqOld.size());
995 for (auto& eq2 : clqOld) {
996 for (auto* vd : eq2.vd) { // point all the variables to the new clique
997 _vVarDescr[vd->payload()].nClique = nCliqueAvailable;
998 }
999 }
1000 clqNew.insert(clqNew.end(), clqOld.begin(), clqOld.end());
1001 clqOld.clear(); // Can use C++11 move TODO
1002 DBGOUT_MIPD(" +++ Joining cliques");
1003 }
1004 nMaybeClq = nCliqueAvailable; // Could mark as 'unused' TODO
1005 }
1006 }
1007
1008 /// Finds a clique variable to which all domain constr are related
1009 class TCliqueSorter {
1010 MIPD& _mipd;
1011 const int _iVarStart; // this is the first var to which all others are related
1012 public:
1013 // VarDecl* varRef0=0; // this is the first var to which all others are related
1014 VarDecl* varRef1 = nullptr; // this is the 2nd main reference.
1015 // it is a var with eq_encode, ||
1016 // an (integer if any) variable with the least rel. factor
1017 bool fRef1HasEqEncode = false;
1018 /// This map stores the relations y = ax+b of all the clique's vars to y
1019 typedef std::unordered_map<VarDecl*, std::pair<double, double> > TMapVars;
1020 TMapVars mRef0, mRef1; // to the main var 0, 1
1021
1022 class TMatrixVars : public std::unordered_map<VarDecl*, TMapVars> {
1023 public:
1024 /// Check existing connection
1025 template <class IVarDecl>
1026 bool checkExistingArc(IVarDecl begV, double A, double B, bool fReportRepeat = true) {
1027 auto it1 = this->find(*begV);
1028 if (this->end() != it1) {
1029 auto it2 = it1->second.find(*(begV + 1));
1030 if (it1->second.end() != it2) {
1031 MZN_MIPD_assert_hard(std::fabs(it2->second.first - A) <
1032 1e-6 * std::max(std::fabs(it2->second.first), std::fabs(A)));
1033 MZN_MIPD_assert_hard(std::fabs(it2->second.second - B) <
1034 1e-6 * std::max(std::fabs(it2->second.second), std::fabs(B)) +
1035 1e-6);
1036 MZN_MIPD_assert_hard(std::fabs(A) != 0.0);
1037 MZN_MIPD_assert_soft(!fVerbose || std::fabs(A) > 1e-12,
1038 " Very small coef: " << (*begV)->id()->str() << " = " << A << " * "
1039 << (*(begV + 1))->id()->str() << " + " << B);
1040 if (fReportRepeat) {
1041 MZN_MIPD_assert_soft(!fVerbose, "LinEqGraph: eqn between "
1042 << (*begV)->id()->str() << " && "
1043 << (*(begV + 1))->id()->str()
1044 << " is repeated. ");
1045 }
1046 return true;
1047 }
1048 }
1049 return false;
1050 }
1051 };
1052
1053 class LinEqGraph : public TMatrixVars {
1054 public:
1055 static double dCoefMin, dCoefMax;
1056 /// Stores the arc (x1, x2) as x1 = a*x2 + b
1057 /// so that a constraint on x2, say x2<=c <-> f,
1058 /// is equivalent to one for x1: x1 <=/>= a*c+b <-> f
1059 //// ( the other way involves division:
1060 //// so that a constraint on x1, say x1<=c <-> f,
1061 //// can easily be converted into one for x2 as a*x2 <= c-b <-> f
1062 //// <=> x2 (care for sign) (c-b)/a <-> f )
1063
1064 template <class ICoef, class IVarDecl>
1065 void addArc(ICoef begC, IVarDecl begV, double rhs) {
1066 MZN_MIPD_assert_soft(!fVerbose || std::fabs(*begC) >= 1e-10,
1067 " Vars " << (*begV)->id()->str() << " to "
1068 << (*(begV + 1))->id()->str() << ": coef=" << (*begC));
1069 // Transform Ax+By=C into x = -B/Ay+C/A
1070 const double negBA = -(*(begC + 1)) / (*begC);
1071 const double CA = rhs / (*begC);
1072 checkExistingArc(begV, negBA, CA);
1073 (*this)[*begV][*(begV + 1)] = std::make_pair(negBA, CA);
1074 const double dCoefAbs = std::fabs(negBA);
1075 if (dCoefAbs < dCoefMin) {
1076 dCoefMin = dCoefAbs;
1077 }
1078 if (dCoefAbs > dCoefMax) {
1079 dCoefMax = dCoefAbs;
1080 }
1081 }
1082 void addEdge(const LinEq2Vars& led) {
1083 addArc(led.coefs.begin(), led.vd.begin(), led.rhs);
1084 addArc(led.coefs.rbegin(), led.vd.rbegin(), led.rhs);
1085 }
1086 /// Propagate linear relations from the given variable
1087 void propagate(iterator itStart, TMapVars& mWhereStore) {
1088 MZN_MIPD_assert_hard(this->end() != itStart);
1089 TMatrixVars mTemp;
1090 mTemp[itStart->first] = itStart->second; // init with existing
1091 DBGOUT_MIPD("Propagation started from " << itStart->first->id()->str() << " having "
1092 << itStart->second.size() << " connections");
1093 propagate2(itStart, itStart, std::make_pair(1.0, 0.0), mTemp);
1094 mWhereStore = mTemp.begin()->second;
1095 MZN_MIPD_assert_hard_msg(
1096 mWhereStore.size() == this->size() - 1,
1097 "Variable " << (*(mTemp.begin()->first))
1098 << " should be connected to all others in the clique, but "
1099 << "|edges|==" << mWhereStore.size() << ", |all nodes|==" << this->size());
1100 }
1101 /// Propagate linear relations from it1 via it2
1102 void propagate2(iterator itSrc, iterator itVia, std::pair<double, double> rel,
1103 TMatrixVars& mWhereStore) {
1104 for (auto itDst = itVia->second.begin(); itDst != itVia->second.end(); ++itDst) {
1105 // Transform x1=A1x2+B1, x2=A2x3+B2 into x1=A1A2x3+A1B2+B1
1106 if (itDst->first == itSrc->first) {
1107 continue;
1108 }
1109 const double A1A2 = rel.first * itDst->second.first;
1110 const double A1B2plusB1 = rel.first * itDst->second.second + rel.second;
1111 bool fDive = true;
1112 if (itSrc != itVia) {
1113 PVarDecl vd[2] = {itSrc->first, itDst->first};
1114 if (!mWhereStore.checkExistingArc(vd, A1A2, A1B2plusB1, false)) {
1115 mWhereStore[vd[0]][vd[1]] = std::make_pair(A1A2, A1B2plusB1);
1116 DBGOUT_MIPD(" PROPAGATING: " << vd[0]->id()->str() << " = " << A1A2 << " * "
1117 << vd[1]->id()->str() << " + " << A1B2plusB1);
1118 } else {
1119 fDive = false;
1120 }
1121 }
1122 if (fDive) {
1123 auto itDST = this->find(itDst->first);
1124 MZN_MIPD_assert_hard(this->end() != itDST);
1125 propagate2(itSrc, itDST, std::make_pair(A1A2, A1B2plusB1), mWhereStore);
1126 }
1127 }
1128 }
1129 };
1130 LinEqGraph leg;
1131
1132 TCliqueSorter(MIPD* pm, int iv) : _mipd(*pm), _iVarStart(iv) {}
1133 void doRelate() {
1134 MZN_MIPD_assert_hard(_mipd._vVarDescr[_iVarStart].nClique >= 0);
1135 const TClique& clq = _mipd._aCliques[_mipd._vVarDescr[_iVarStart].nClique];
1136 for (const auto& eq2 : clq) {
1137 leg.addEdge(eq2);
1138 }
1139 DBGOUT_MIPD(" Clique " << _mipd._vVarDescr[_iVarStart].nClique << ": " << leg.size()
1140 << " variables, " << clq.size() << " connections.");
1141 for (auto& it1 : leg) {
1142 _mipd._vVarDescr[it1.first->payload()].fDomainConstrProcessed = 1U;
1143 }
1144
1145 // Propagate the 1st var's relations:
1146 leg.propagate(leg.begin(), mRef0);
1147
1148 // Find a best main variable according to:
1149 // 1. isInt 2. hasEqEncode 3. abs linFactor to ref0
1150 varRef1 = leg.begin()->first;
1151 std::array<double, 3> aCrit = {
1152 {(double)_mipd._vVarDescr[varRef1->payload()].fInt,
1153 static_cast<double>(_mipd._vVarDescr[varRef1->payload()].pEqEncoding != nullptr), 1.0}};
1154 for (auto& it2 : mRef0) {
1155 VarDescr& vard = _mipd._vVarDescr[it2.first->payload()];
1156 std::array<double, 3> aCrit1 = {{(double)vard.fInt,
1157 static_cast<double>(vard.pEqEncoding != nullptr),
1158 std::fabs(it2.second.first)}};
1159 if (aCrit1 > aCrit) {
1160 varRef1 = it2.first;
1161 aCrit = aCrit1;
1162 }
1163 }
1164 leg.propagate(leg.find(varRef1), mRef1);
1165 }
1166 }; // class TCliqueSorter
1167
1168 /// Build a domain decomposition for a clique
1169 /// a clique can consist of just 1 var without a clique object
1170 class DomainDecomp {
1171 public:
1172 MIPD& mipd;
1173 const int iVarStart;
1174 TCliqueSorter cls;
1175 SetOfIntvReal sDomain;
1176
1177 DomainDecomp(MIPD* pm, int iv) : mipd(*pm), iVarStart(iv), cls(pm, iv) {
1178 sDomain.insert(IntvReal()); // the decomposed domain. Init to +-inf
1179 }
1180 void doProcess() {
1181 // Choose the main variable && relate all others to it
1182 const int nClique = mipd._vVarDescr[iVarStart].nClique;
1183 if (nClique >= 0) {
1184 cls.doRelate();
1185 } else {
1186 cls.varRef1 = mipd._vVarDescr[iVarStart].vd;
1187 }
1188 // Adding itself:
1189 cls.mRef1[cls.varRef1] = std::make_pair(1.0, 0.0);
1190
1191 int iVarRef1 = cls.varRef1->payload();
1192 MZN_MIPD_assert_hard(nClique == mipd._vVarDescr[iVarRef1].nClique);
1193 cls.fRef1HasEqEncode = (mipd._vVarDescr[iVarRef1].pEqEncoding != nullptr);
1194
1195 // First, construct the domain decomposition in any case
1196 // projectVariableConstr( cls.varRef1, std::make_pair(1.0, 0.0) );
1197 // if ( nClique >= 0 ) {
1198 for (auto& iRef1 : cls.mRef1) {
1199 projectVariableConstr(iRef1.first, iRef1.second);
1200 }
1201
1202 DBGOUT_MIPD("Clique " << nClique << ": main ref var " << cls.varRef1->id()->str()
1203 << ", domain dec: " << sDomain);
1204
1205 MZN_MIPD_ASSERT_FOR_SAT(!sDomain.empty(), mipd.getEnv()->envi(), cls.varRef1->loc(),
1206 "clique " << nClique << ": main ref var " << *cls.varRef1->id()
1207 << ", domain decomposition seems empty: " << sDomain);
1208
1209 MZN_MIPD_FLATTENING_ERROR_IF_NOT(sDomain.checkFiniteBounds(), mipd.getEnv()->envi(),
1210 cls.varRef1->loc(),
1211 "variable " << *cls.varRef1->id()
1212 << " needs finite bounds for linearisation."
1213 " Or, use indicator constraints. "
1214 << "Current domain is " << sDomain);
1215
1216 MZN_MIPD_assert_hard(sDomain.checkDisjunctStrict());
1217
1218 makeRangeDomains();
1219
1220 // Then, use equality_encoding if available
1221 if (cls.fRef1HasEqEncode) {
1222 syncWithEqEncoding();
1223 syncOtherEqEncodings();
1224 } else { // ! cls.fRef1HasEqEncode
1225 if (sDomain.size() >= 2) { // need to simplify stuff otherwise
1226 considerDenseEncoding();
1227 createDomainFlags();
1228 }
1229 }
1230 implementPOSTs();
1231
1232 // Statistics
1233 if (sDomain.size() < MIPD_stats[N_POSTs_NSubintvMin]) {
1234 MIPD_stats[N_POSTs_NSubintvMin] = static_cast<double>(sDomain.size());
1235 }
1236 MIPD_stats[N_POSTs_NSubintvSum] += sDomain.size();
1237 if (sDomain.size() > MIPD_stats[N_POSTs_NSubintvMax]) {
1238 MIPD_stats[N_POSTs_NSubintvMax] = static_cast<double>(sDomain.size());
1239 }
1240 for (const auto& intv : sDomain) {
1241 const auto nSubSize = intv.right - intv.left;
1242 if (nSubSize < MIPD_stats[N_POSTs_SubSizeMin]) {
1243 MIPD_stats[N_POSTs_SubSizeMin] = nSubSize;
1244 }
1245 MIPD_stats[N_POSTs_SubSizeSum] += nSubSize;
1246 if (nSubSize > MIPD_stats[N_POSTs_SubSizeMax]) {
1247 MIPD_stats[N_POSTs_SubSizeMax] = nSubSize;
1248 }
1249 }
1250 if (cls.fRef1HasEqEncode) {
1251 ++MIPD_stats[N_POSTs_cliquesWithEqEncode];
1252 }
1253 }
1254
1255 /// Project the domain-related constraints of a variable into the clique
1256 /// Deltas should be scaled but to a minimum of the target's discr
1257 /// COmparison sense changes on negated vars
1258 void projectVariableConstr(VarDecl* vd, std::pair<double, double> eq1) {
1259 DBGOUT_MIPD_FLUSH(" MIPD: projecting variable ");
1260 DBGOUT_MIPD_SELF(debugprint(vd));
1261 // Always check if domain becomes empty? TODO
1262 const double A = eq1.first; // vd = A*arg + B. conversion
1263 const double B = eq1.second;
1264 // process domain info
1265 double lb = B;
1266 double ub = A + B; // projected bounds for bool
1267 if (vd->ti()->domain() != nullptr) {
1268 if (vd->type().isint() || vd->type().isfloat()) { // INT VAR OR FLOAT VAR
1269 SetOfIntvReal sD1;
1270 convertIntSet(vd->ti()->domain(), sD1, cls.varRef1, A, B);
1271 sDomain.intersect(sD1);
1272 DBGOUT_MIPD(" Clique domain after proj of the init. domain "
1273 << sD1 << " of " << (vd->type().isint() ? "varint" : "varfloat") << A << " * "
1274 << vd->id()->str() << " + " << B << ": " << sDomain);
1275 auto bnds = sD1.getBounds();
1276 lb = bnds.left;
1277 ub = bnds.right;
1278 } else {
1279 MZN_MIPD_FLATTENING_ERROR_IF_NOT(0, mipd.getEnv()->envi(), cls.varRef1->loc(),
1280 "Variable " << vd->id()->str() << " of type "
1281 << vd->type().toString(mipd._env->envi())
1282 << " has a domain.");
1283 }
1284 // /// Deleting var domain:
1285 // vd->ti()->domain( NULL );
1286 } else {
1287 if (nullptr == vd->ti()->domain() && !vd->type().isbool()) {
1288 lb = IntvReal::infMinus();
1289 ub = IntvReal::infPlus();
1290 }
1291 }
1292 auto bnds = sDomain.getBounds(); // can change TODO
1293 // process calls. Can use the constr type info.
1294 auto& aCalls = mipd._vVarDescr[vd->payload()].aCalls;
1295 for (Item* pItem : aCalls) {
1296 auto* pCI = pItem->dynamicCast<ConstraintI>();
1297 MZN_MIPD_assert_hard(pCI != nullptr);
1298 Call* pCall = pCI->e()->dynamicCast<Call>();
1299 MZN_MIPD_assert_hard(pCall != nullptr);
1300 DBGOUT_MIPD_FLUSH("PROPAG CALL ");
1301 DBGOUT_MIPD_SELF(debugprint(pCall));
1302 // check the bounds for bool in reifs? TODO
1303 auto ipct = mipd._mCallTypes.find(pCall->decl());
1304 MZN_MIPD_assert_hard(mipd._mCallTypes.end() != ipct);
1305 const DCT& dct = *ipct->second;
1306 int nCmpType_ADAPTED = dct.nCmpType;
1307 if (A < 0.0) { // negative factor
1308 if (std::abs(nCmpType_ADAPTED) >= 4) { // inequality
1309 nCmpType_ADAPTED = -nCmpType_ADAPTED;
1310 }
1311 }
1312 switch (dct.nConstrType) {
1313 case CT_SetIn: {
1314 SetOfIntvReal SS;
1315 convertIntSet(pCall->arg(1), SS, cls.varRef1, A, B);
1316 if (RIT_Static == dct.nReifType) {
1317 sDomain.intersect(SS);
1318 ++MIPD_stats[N_POSTs_setIn];
1319 } else {
1320 sDomain.cutDeltas(SS, std::max(1.0, std::fabs(A))); // deltas to scale
1321 ++MIPD_stats[N_POSTs_setInReif];
1322 }
1323 } break;
1324 case CT_Comparison:
1325 if (RIT_Reif == dct.nReifType) {
1326 const double rhs = (mipd.aux_float_lt_zero_iff_1_POST == pCall->decl())
1327 ? B /* + A*0.0, relating to 0 */
1328 // The 2nd argument is constant:
1329 : A * MIPD::expr2Const(pCall->arg(1)) + B;
1330 const double rhsUp = rndUpIfInt(cls.varRef1, rhs);
1331 const double rhsDown = rndDownIfInt(cls.varRef1, rhs);
1332 const double rhsRnd = rndIfInt(cls.varRef1, rhs);
1333 /// Strictly, for delta we should finish domain reductions first... TODO?
1334 const double delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 3);
1335 switch (nCmpType_ADAPTED) {
1336 case CMPT_LE:
1337 sDomain.cutDeltas(IntvReal::infMinus(), rhsDown, delta);
1338 break;
1339 case CMPT_GE:
1340 sDomain.cutDeltas(rhsUp, IntvReal::infPlus(), delta);
1341 break;
1342 case CMPT_LT_0:
1343 sDomain.cutDeltas(IntvReal::infMinus(), rhsDown - delta, delta);
1344 break;
1345 case CMPT_GT_0:
1346 sDomain.cutDeltas(rhsUp + delta, IntvReal::infPlus(), delta);
1347 break;
1348 case CMPT_EQ:
1349 if (!(cls.varRef1->type().isint() && // skip if int target var
1350 std::fabs(rhs - rhsRnd) > INT_EPS)) { // && fract value
1351 sDomain.cutDeltas(rhsRnd, rhsRnd, delta);
1352 }
1353 break;
1354 default:
1355 MZN_MIPD_assert_hard_msg(0, " No other reified cmp type ");
1356 }
1357 ++MIPD_stats[(vd->ti()->type().isint()) ? N_POSTs_intCmpReif : N_POSTs_floatCmpReif];
1358 } else if (RIT_Static == dct.nReifType) {
1359 // _ne, later maybe static ineq TODO
1360 MZN_MIPD_assert_hard(CMPT_NE == dct.nCmpType);
1361 const double rhs = A * MIPD::expr2Const(pCall->arg(1)) + B;
1362 const double rhsRnd = rndIfInt(cls.varRef1, rhs);
1363 bool fSkipNE = (cls.varRef1->type().isint() && std::fabs(rhs - rhsRnd) > INT_EPS);
1364 if (!fSkipNE) {
1365 const double delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 2);
1366 sDomain.cutOut({rhsRnd - delta, rhsRnd + delta});
1367 }
1368 ++MIPD_stats[(vd->ti()->type().isint()) ? N_POSTs_intNE : N_POSTs_floatNE];
1369 } else { // aux_ relate to 0.0
1370 // But we don't modify domain splitting for them currently
1371 ++MIPD_stats[(vd->ti()->type().isint()) ? N_POSTs_intAux : N_POSTs_floatAux];
1372 MZN_MIPD_assert_hard(RIT_Halfreif == dct.nReifType);
1373 // const double rhs = B; // + A*0
1374 // const double delta = vd->type().isint() ? 1.0 : 1e-5; //
1375 // TODO : eps
1376 }
1377 break;
1378 case CT_Encode:
1379 // See if any further constraints here? TODO
1380 ++MIPD_stats[N_POSTs_eq_encode];
1381 break;
1382 default:
1383 MZN_MIPD_assert_hard_msg(0, "Unknown constraint type");
1384 }
1385 }
1386 DBGOUT_MIPD(" Clique domain after proj of " << A << " * " << vd->id()->str() << " + " << B
1387 << ": " << sDomain);
1388 }
1389
1390 static double rndIfInt(VarDecl* vdTarget, double v) {
1391 return vdTarget->type().isint() ? std::round(v) : v;
1392 }
1393 static double rndIfBothInt(VarDecl* vdTarget, double v) {
1394 if (!vdTarget->type().isint()) {
1395 return v;
1396 }
1397 const double vRnd = std::round(v);
1398 return (fabs(v - vRnd) < INT_EPS) ? vRnd : v;
1399 }
1400 static double rndUpIfInt(VarDecl* vdTarget, double v) {
1401 return vdTarget->type().isint() ? std::ceil(v - INT_EPS) : v;
1402 }
1403 static double rndDownIfInt(VarDecl* vdTarget, double v) {
1404 return vdTarget->type().isint() ? std::floor(v + INT_EPS) : v;
1405 }
1406
1407 void makeRangeDomains() {
1408 auto bnds = sDomain.getBounds();
1409 for (auto& iRef1 : cls.mRef1) {
1410 VarDecl* vd = iRef1.first;
1411 // projecting the bounds back:
1412 double lb0 = (bnds.left - iRef1.second.second) / iRef1.second.first;
1413 double ub0 = (bnds.right - iRef1.second.second) / iRef1.second.first;
1414 if (lb0 > ub0) {
1415 MZN_MIPD_assert_hard(iRef1.second.first < 0.0);
1416 std::swap(lb0, ub0);
1417 }
1418 if (vd->type().isint()) {
1419 lb0 = rndUpIfInt(vd, lb0);
1420 ub0 = rndDownIfInt(vd, ub0);
1421 }
1422 setVarDomain(vd, lb0, ub0);
1423 }
1424 }
1425
1426 /// tightens element bounds in the existing eq_encoding of varRef1
1427 /// necessary because if one exists, int_ne is not translated into it
1428 /// Can also back-check from there? TODO
1429 /// And further checks TODO
1430 void syncWithEqEncoding() {
1431 std::vector<Expression*> pp;
1432 auto bnds = sDomain.getBounds();
1433 const long long iMin = mipd.expr2ExprArray(
1434 mipd._vVarDescr[cls.varRef1->payload()].pEqEncoding->e()->dynamicCast<Call>()->arg(1),
1435 pp);
1436 MZN_MIPD_assert_hard(pp.size() >= bnds.right - bnds.left + 1);
1437 MZN_MIPD_assert_hard(iMin <= bnds.left);
1438 long long vEE = iMin;
1439 DBGOUT_MIPD_FLUSH(
1440 " SYNC EQ_ENCODE( "
1441 << (*cls.varRef1) << ", bitflags: "
1442 << *(mipd._vVarDescr[cls.varRef1->payload()].pEqEncoding->e()->dynamicCast<Call>()->arg(
1443 1))
1444 << " ): SETTING 0 FLAGS FOR VALUES: ");
1445 for (const auto& intv : sDomain) {
1446 for (; static_cast<double>(vEE) < intv.left; ++vEE) {
1447 if (vEE >= static_cast<long long>(iMin + pp.size())) {
1448 return;
1449 }
1450 if (pp[vEE - iMin]->isa<Id>()) {
1451 if (pp[vEE - iMin]->dynamicCast<Id>()->decl()->type().isvar()) {
1452 DBGOUT_MIPD_FLUSH(vEE << ", ");
1453 setVarDomain(pp[vEE - iMin]->dynamicCast<Id>()->decl(), 0.0, 0.0);
1454 }
1455 }
1456 }
1457 vEE = static_cast<long long>(intv.right + 1);
1458 }
1459 for (; vEE < static_cast<long long>(iMin + pp.size()); ++vEE) {
1460 if (pp[vEE - iMin]->isa<Id>()) {
1461 if (pp[vEE - iMin]->dynamicCast<Id>()->decl()->type().isvar()) {
1462 DBGOUT_MIPD_FLUSH(vEE << ", ");
1463 setVarDomain(pp[vEE - iMin]->dynamicCast<Id>()->decl(), 0.0, 0.0);
1464 }
1465 }
1466 }
1467 DBGOUT_MIPD("");
1468 }
1469
1470 /// sync varRef1's eq_encoding with those of other variables
1471 void syncOtherEqEncodings() {
1472 // TODO This could be in the var projection? No, need the final domain
1473 }
1474
1475 /// Depending on params,
1476 /// create an equality encoding for an integer variable
1477 /// TODO What if a float's domain is discrete?
1478 void considerDenseEncoding() {
1479 if (cls.varRef1->id()->type().isint()) {
1480 if (sDomain.maxInterval() <= mipd.nMaxIntv2Bits ||
1481 sDomain.cardInt() <= mipd.dMaxNValueDensity * sDomain.size()) {
1482 sDomain.split2Bits();
1483 ++MIPD_stats[N_POSTs_clEEEnforced];
1484 }
1485 }
1486 }
1487
1488 /// if ! eq_encoding, creates a flag for each subinterval in the domain
1489 /// && constrains sum(flags)==1
1490 void createDomainFlags() {
1491 std::vector<Expression*> vVars(sDomain.size()); // flags for each subinterval
1492 std::vector<double> vIntvLB(sDomain.size() + 1);
1493 std::vector<double> vIntvUB_(sDomain.size() + 1);
1494 int i = 0;
1495 double dMaxIntv = -1.0;
1496 for (const auto& intv : sDomain) {
1497 intv.varFlag = addIntVar(0.0, 1.0);
1498 vVars[i] = intv.varFlag->id();
1499 vIntvLB[i] = intv.left;
1500 vIntvUB_[i] = -intv.right;
1501 dMaxIntv = std::max(dMaxIntv, intv.right - intv.left);
1502 ++i;
1503 }
1504 // Sum of flags == 1
1505 std::vector<double> ones(sDomain.size(), 1.0);
1506 addLinConstr(ones, vVars, CMPT_EQ, 1.0);
1507 // Domain decomp
1508 vVars.push_back(cls.varRef1->id());
1509 vIntvLB[i] = -1.0; // var1 >= sum(LBi*flagi)
1510 /// STRICT equality encoding if small intervals
1511 if (dMaxIntv > 1e-6) { // EPS = param? TODO
1512 vIntvUB_[i] = 1.0; // var1 <= sum(UBi*flagi)
1513 addLinConstr(vIntvLB, vVars, CMPT_LE, 0.0);
1514 addLinConstr(vIntvUB_, vVars, CMPT_LE, 0.0);
1515 } else {
1516 ++MIPD_stats[N_POSTs_clEEFound];
1517 addLinConstr(vIntvLB, vVars, CMPT_EQ, 0.0);
1518 }
1519 }
1520
1521 /// deletes them as well
1522 void implementPOSTs() {
1523 auto bnds = sDomain.getBounds();
1524 for (auto& iRef1 : cls.mRef1) {
1525 // DBGOUT_MIPD_FLUSH( " MIPD: implementing constraints of variable " );
1526 // DBGOUT_MIPD_SELF( debugprint(vd) );
1527 VarDecl* vd = iRef1.first;
1528 auto eq1 = iRef1.second;
1529 const double A = eq1.first; // vd = A*arg + B. conversion
1530 const double B = eq1.second;
1531 // process calls. Can use the constr type info.
1532 auto& aCalls = mipd._vVarDescr[vd->payload()].aCalls;
1533 for (Item* pItem : aCalls) {
1534 auto* pCI = pItem->dynamicCast<ConstraintI>();
1535 MZN_MIPD_assert_hard(pCI);
1536 Call* pCall = pCI->e()->dynamicCast<Call>();
1537 MZN_MIPD_assert_hard(pCall);
1538 DBGOUT_MIPD_FLUSH("IMPL CALL ");
1539 DBGOUT_MIPD_SELF(debugprint(pCall));
1540 // check the bounds for bool in reifs? TODO
1541 auto ipct = mipd._mCallTypes.find(pCall->decl());
1542 MZN_MIPD_assert_hard(mipd._mCallTypes.end() != ipct);
1543 const DCT& dct = *ipct->second;
1544 int nCmpType_ADAPTED = dct.nCmpType;
1545 if (A < 0.0) { // negative factor
1546 if (std::abs(nCmpType_ADAPTED) >= 4) { // inequality
1547 nCmpType_ADAPTED = -nCmpType_ADAPTED;
1548 }
1549 }
1550 switch (dct.nConstrType) {
1551 case CT_SetIn:
1552 if (RIT_Reif == dct.nReifType) {
1553 SetOfIntvReal SS;
1554 convertIntSet(pCall->arg(1), SS, cls.varRef1, A, B);
1555 relateReifFlag(pCall->arg(2), SS);
1556 }
1557 break;
1558 case CT_Comparison:
1559 if (RIT_Reif == dct.nReifType) {
1560 const double rhs = (mipd.aux_float_lt_zero_iff_1_POST == pCall->decl())
1561 ? B /* + A*0.0, relating to 0 */
1562 // The 2nd argument is constant:
1563 : A * MIPD::expr2Const(pCall->arg(1)) + B;
1564 const double rhsUp = rndUpIfInt(cls.varRef1, rhs);
1565 const double rhsDown = rndDownIfInt(cls.varRef1, rhs);
1566 const double rhsRnd = rndIfBothInt(
1567 cls.varRef1, rhs); // if the ref var is int, need to round almost-int values
1568 const double delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 3);
1569 switch (nCmpType_ADAPTED) {
1570 case CMPT_LE:
1571 relateReifFlag(pCall->arg(2), {{IntvReal::infMinus(), rhsDown}});
1572 break;
1573 case CMPT_GE:
1574 relateReifFlag(pCall->arg(2), {{rhsUp, IntvReal::infPlus()}});
1575 break;
1576 case CMPT_LT_0:
1577 relateReifFlag(pCall->arg(1), {{IntvReal::infMinus(), rhsDown - delta}});
1578 break;
1579 case CMPT_GT_0:
1580 relateReifFlag(pCall->arg(1), {{rhsUp + delta, IntvReal::infPlus()}});
1581 break;
1582 case CMPT_EQ:
1583 relateReifFlag(pCall->arg(2), {{rhsRnd, rhsRnd}});
1584 break; // ... but if the value is sign. fractional for an int var, the flag is
1585 // set=0
1586 default:
1587 break;
1588 }
1589 } else if (RIT_Static == dct.nReifType) {
1590 // !hing here for NE
1591 MZN_MIPD_assert_hard(CMPT_NE == nCmpType_ADAPTED);
1592 } else { // aux_ relate to 0.0
1593 // But we don't modify domain splitting for them currently
1594 MZN_MIPD_assert_hard(RIT_Halfreif == dct.nReifType);
1595 double rhs = B; // + A*0
1596 const double rhsUp = rndUpIfInt(cls.varRef1, rhs);
1597 const double rhsDown = rndDownIfInt(cls.varRef1, rhs);
1598 const double rhsRnd = rndIfInt(cls.varRef1, rhs);
1599 double delta = 0.0;
1600 if (mipd.aux_float_lt_zero_if_1_POST == pCall->decl()) { // only float && lt
1601 delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 3);
1602 }
1603 if (nCmpType_ADAPTED < 0) {
1604 delta = -delta;
1605 }
1606 if (cls.varRef1->type().isint() && CMPT_EQ_0 != nCmpType_ADAPTED) {
1607 if (nCmpType_ADAPTED < 0) {
1608 rhs = rhsDown;
1609 } else {
1610 rhs = rhsUp;
1611 }
1612 } else {
1613 rhs += delta;
1614 }
1615 // Now we need rhs ! to be in the inner of the domain
1616 bool fUseDD = true;
1617 if (!cls.fRef1HasEqEncode) {
1618 switch (nCmpType_ADAPTED) {
1619 case CMPT_EQ_0: {
1620 auto itLB = sDomain.lower_bound(rhsRnd);
1621 fUseDD = (itLB->left == rhsRnd && itLB->right == rhsRnd); // exactly
1622 } break;
1623 case CMPT_LT_0:
1624 case CMPT_LE_0: {
1625 auto itUB = sDomain.upper_bound(rhsUp);
1626 bool fInner = false;
1627 if (sDomain.begin() != itUB) {
1628 --itUB;
1629 if (itUB->right > rhs) {
1630 fInner = true;
1631 }
1632 }
1633 fUseDD = !fInner;
1634 } break;
1635 case CMPT_GT_0:
1636 case CMPT_GE_0: {
1637 auto itLB = sDomain.lower_bound(rhsDown);
1638 bool fInner = false;
1639 if (sDomain.begin() != itLB) {
1640 --itLB;
1641 if (itLB->right >= rhs) {
1642 fInner = true;
1643 }
1644 }
1645 fUseDD = !fInner;
1646 } break;
1647 default:
1648 MZN_MIPD_assert_hard_msg(0, "Unknown halfreif cmp type");
1649 }
1650 }
1651 if (fUseDD) { // use sDomain
1652 if (CMPT_EQ_0 == nCmpType_ADAPTED) {
1653 relateReifFlag(pCall->arg(1), {{rhsRnd, rhsRnd}}, RIT_Halfreif);
1654 } else if (nCmpType_ADAPTED < 0) {
1655 relateReifFlag(pCall->arg(1), {{IntvReal::infMinus(), rhsDown}}, RIT_Halfreif);
1656 } else {
1657 relateReifFlag(pCall->arg(1), {{rhsUp, IntvReal::infPlus()}}, RIT_Halfreif);
1658 }
1659 } else { // use big-M
1660 DBGOUT_MIPD(" AUX BY BIG-Ms: ");
1661 const bool fLE = (CMPT_EQ_0 == nCmpType_ADAPTED || 0 > nCmpType_ADAPTED);
1662 const bool fGE = (CMPT_EQ_0 == nCmpType_ADAPTED || 0 < nCmpType_ADAPTED);
1663 // Take integer || float indicator version, depending on the constrained var:
1664 const int nIdxInd = // (VT_Int==dct.nVarType) ?
1665 // No: vd->ti()->type().isint() ? 1 : 2;
1666 cls.varRef1->ti()->type().isint()
1667 ? 1
1668 : 2; // need the type of the variable to be constr
1669 MZN_MIPD_assert_hard(static_cast<unsigned int>(nIdxInd) < pCall->argCount());
1670 Expression* pInd = pCall->arg(nIdxInd);
1671 if (fLE && rhs < bnds.right) {
1672 if (rhs >= bnds.left) {
1673 std::vector<double> coefs = {1.0, bnds.right - rhs};
1674 // Use the float version of indicator:
1675 std::vector<Expression*> vars = {cls.varRef1->id(), pInd};
1676 addLinConstr(coefs, vars, CMPT_LE, bnds.right);
1677 } else {
1678 setVarDomain(MIPD::expr2VarDecl(pInd), 0.0, 0.0);
1679 }
1680 }
1681 if (fGE && rhs > bnds.left) {
1682 if (rhs <= bnds.right) {
1683 std::vector<double> coefs = {-1.0, rhs - bnds.left};
1684 std::vector<Expression*> vars = {cls.varRef1->id(), pInd};
1685 addLinConstr(coefs, vars, CMPT_LE, -bnds.left);
1686 } else {
1687 setVarDomain(MIPD::expr2VarDecl(pInd), 0.0, 0.0);
1688 }
1689 }
1690 }
1691 }
1692 break;
1693 case CT_Encode:
1694 // See if any further constraints here? TODO
1695 break;
1696 default:
1697 MZN_MIPD_assert_hard_msg(0, "Unknown constraint type");
1698 }
1699 pItem->remove(); // removing the call
1700 }
1701 // removing the eq_encoding call
1702 if (mipd._vVarDescr[vd->payload()].pEqEncoding != nullptr) {
1703 mipd._vVarDescr[vd->payload()].pEqEncoding->remove();
1704 }
1705 }
1706 }
1707
1708 /// sets varFlag = || <= sum( intv.varFlag : SS )
1709 void relateReifFlag(Expression* expFlag, const SetOfIntvReal& SS, EnumReifType nRT = RIT_Reif) {
1710 MZN_MIPD_assert_hard(RIT_Reif == nRT || RIT_Halfreif == nRT);
1711 // MZN_MIPD_assert_hard( sDomain.size()>=2 );
1712 VarDecl* varFlag = MIPD::expr2VarDecl(expFlag);
1713 std::vector<Expression*> vIntvFlags;
1714 if (cls.fRef1HasEqEncode) { // use eq_encoding
1715 MZN_MIPD_assert_hard(varFlag->type().isint());
1716 std::vector<Expression*> pp;
1717 auto bnds = sDomain.getBounds();
1718 const long long iMin = mipd.expr2ExprArray(
1719 mipd._vVarDescr[cls.varRef1->payload()].pEqEncoding->e()->dynamicCast<Call>()->arg(1),
1720 pp);
1721 MZN_MIPD_assert_hard(pp.size() >= bnds.right - bnds.left + 1);
1722 MZN_MIPD_assert_hard(iMin <= bnds.left);
1723 for (const auto& intv : SS) {
1724 for (long long vv = (long long)std::max(double(iMin), ceil(intv.left));
1725 vv <= (long long)std::min(double(iMin) + pp.size() - 1, floor(intv.right)); ++vv) {
1726 vIntvFlags.push_back(pp[vv - iMin]);
1727 }
1728 }
1729 } else {
1730 MZN_MIPD_assert_hard(varFlag->type().isint());
1731 for (const auto& intv : SS) {
1732 auto it1 = sDomain.lower_bound(intv.left);
1733 auto it2 = sDomain.upper_bound(intv.right);
1734 auto it11 = it1;
1735 // Check that we are looking ! into a subinterval:
1736 if (sDomain.begin() != it11) {
1737 --it11;
1738 MZN_MIPD_assert_hard(it11->right < intv.left);
1739 }
1740 auto it12 = it2;
1741 if (sDomain.begin() != it12) {
1742 --it12;
1743 MZN_MIPD_assert_hard_msg(it12->right <= intv.right,
1744 " relateReifFlag for " << intv << " in " << sDomain);
1745 }
1746 for (it12 = it1; it12 != it2; ++it12) {
1747 if (it12->varFlag != nullptr) {
1748 vIntvFlags.push_back(it12->varFlag->id());
1749 } else {
1750 MZN_MIPD_assert_hard(1 == sDomain.size());
1751 vIntvFlags.push_back(IntLit::a(1)); // just a constant then
1752 }
1753 }
1754 }
1755 }
1756 if (!vIntvFlags.empty()) {
1757 // Could find out if reif is true -- TODO && see above for 1 subinterval
1758 std::vector<double> onesm(vIntvFlags.size(), -1.0);
1759 onesm.push_back(1.0);
1760 vIntvFlags.push_back(varFlag->id());
1761 EnumCmpType nCmpType = (RIT_Reif == nRT) ? CMPT_EQ : CMPT_LE;
1762 addLinConstr(onesm, vIntvFlags, nCmpType, 0.0);
1763 } else { // the reif is false
1764 setVarDomain(varFlag, 0.0, 0.0);
1765 }
1766 }
1767
1768 static void setVarDomain(VarDecl* vd, double lb, double ub) {
1769 // need to check if the new range is in the previous bounds... TODO
1770 if (vd->type().isfloat()) {
1771 // if ( 0.0==lb && 0.0==ub ) {
1772 auto* newDom =
1773 new BinOp(Location().introduce(), FloatLit::a(lb), BOT_DOTDOT, FloatLit::a(ub));
1774 vd->ti()->domain(newDom);
1775 DBGOUT_MIPD(" NULL OUT: " << vd->id()->str());
1776 // }
1777 } else if (vd->type().isint() || vd->type().isbool()) {
1778 auto* newDom = new SetLit(
1779 Location().introduce(),
1780 IntSetVal::a(static_cast<long long int>(lb), static_cast<long long int>(ub)));
1781 // TypeInst* nti = copy(mipd.getEnv()->envi(),varFlag->ti())->cast<TypeInst>();
1782 // nti->domain(newDom);
1783 vd->ti()->domain(newDom);
1784 } else {
1785 MZN_MIPD_assert_hard_msg(0, "Unknown var type ");
1786 }
1787 }
1788
1789 VarDecl* addIntVar(double LB, double UB) {
1790 // GCLock lock;
1791 // Cache them? Only location can be different TODO
1792 auto* newDom =
1793 new SetLit(Location().introduce(),
1794 IntSetVal::a(static_cast<long long int>(LB), static_cast<long long int>(UB)));
1795 auto* ti = new TypeInst(Location().introduce(), Type::varint(), newDom);
1796 auto* newVar = new VarDecl(Location().introduce(), ti, mipd.getEnv()->envi().genId());
1797 newVar->flat(newVar);
1798 mipd.getEnv()->envi().flatAddItem(new VarDeclI(Location().introduce(), newVar));
1799 return newVar;
1800 }
1801
1802 void addLinConstr(std::vector<double>& coefs, std::vector<Expression*>& vars,
1803 EnumCmpType nCmpType, double rhs) {
1804 std::vector<Expression*> args(3);
1805 MZN_MIPD_assert_hard(vars.size() >= 2);
1806 for (auto* v : vars) {
1807 MZN_MIPD_assert_hard(&v);
1808 // throw std::string("addLinConstr: &var=NULL");
1809 MZN_MIPD_assert_hard_msg(v->isa<Id>() || v->isa<IntLit>() || v->isa<FloatLit>(),
1810 " expression at " << (&v) << " eid = " << v->eid()
1811 << " while E_INTLIT=" << Expression::E_INTLIT);
1812 // throw std::string("addLinConstr: only id's as variables allowed");
1813 }
1814 MZN_MIPD_assert_hard(coefs.size() == vars.size());
1815 MZN_MIPD_assert_hard(CMPT_EQ == nCmpType || CMPT_LE == nCmpType);
1816 DBGOUT_MIPD_SELF( // LinEq leq; leq.coefs=coefs; leq.vd=vars; leq.rhs=rhs;
1817 DBGOUT_MIPD_FLUSH(" ADDING " << (CMPT_EQ == nCmpType ? "LIN_EQ" : "LIN_LE") << ": [ ");
1818 for (auto c
1819 : coefs) DBGOUT_MIPD_FLUSH(c << ',');
1820 DBGOUT_MIPD_FLUSH(" ] * [ "); for (auto v
1821 : vars) {
1822 MZN_MIPD_assert_hard(!v->isa<VarDecl>());
1823 if (v->isa<Id>()) DBGOUT_MIPD_FLUSH(v->dynamicCast<Id>()->str() << ',');
1824 // else if ( v->isa<VarDecl>() )
1825 // MZN_MIPD_assert_hard ("addLinConstr: only id's as variables allowed");
1826 else
1827 DBGOUT_MIPD_FLUSH(mipd.expr2Const(v) << ',');
1828 } DBGOUT_MIPD(" ] " << (CMPT_EQ == nCmpType ? "== " : "<= ") << rhs););
1829 std::vector<Expression*> nc_c;
1830 std::vector<Expression*> nx;
1831 bool fFloat = false;
1832 for (auto* v : vars) {
1833 if (!v->type().isint()) {
1834 fFloat = true;
1835 break;
1836 }
1837 }
1838 auto sName = constants().ids.float_.lin_eq; // "int_lin_eq";
1839 FunctionI* fDecl = mipd.float_lin_eq;
1840 if (fFloat) { // MZN_MIPD_assert_hard all vars of same type TODO
1841 for (int i = 0; i < vars.size(); ++i) {
1842 if (fabs(coefs[i]) > 1e-8) /// Only add terms with non-0 coefs. TODO Eps=param
1843 {
1844 nc_c.push_back(FloatLit::a(coefs[i]));
1845 if (vars[i]->type().isint()) {
1846 std::vector<Expression*> i2f_args(1);
1847 i2f_args[0] = vars[i];
1848 Call* i2f = new Call(Location().introduce(), constants().ids.int2float, i2f_args);
1849 i2f->type(Type::varfloat());
1850 i2f->decl(mipd.getEnv()->model()->matchFn(mipd.getEnv()->envi(), i2f, false));
1851 EE ret = flat_exp(mipd.getEnv()->envi(), Ctx(), i2f, nullptr, constants().varTrue);
1852 nx.push_back(ret.r());
1853 } else {
1854 nx.push_back(vars[i]); // ->id(); once passing a general expression
1855 }
1856 }
1857 }
1858 args[2] = FloatLit::a(rhs);
1859 args[2]->type(Type::parfloat(0));
1860 args[0] = new ArrayLit(Location().introduce(), nc_c);
1861 args[0]->type(Type::parfloat(1));
1862 args[1] = new ArrayLit(Location().introduce(), nx);
1863 args[1]->type(Type::varfloat(1));
1864 if (CMPT_LE == nCmpType) {
1865 sName = constants().ids.float_.lin_le; // "float_lin_le";
1866 fDecl = mipd.float_lin_le;
1867 }
1868 } else {
1869 for (int i = 0; i < vars.size(); ++i) {
1870 if (fabs(coefs[i]) > 1e-8) /// Only add terms with non-0 coefs. TODO Eps=param
1871 {
1872 nc_c.push_back(IntLit::a(static_cast<long long int>(coefs[i])));
1873 nx.push_back(vars[i]); //->id();
1874 }
1875 }
1876 args[2] = IntLit::a(static_cast<long long int>(rhs));
1877 args[2]->type(Type::parint(0));
1878 args[0] = new ArrayLit(Location().introduce(), nc_c);
1879 args[0]->type(Type::parint(1));
1880 args[1] = new ArrayLit(Location().introduce(), nx);
1881 args[1]->type(Type::varint(1));
1882 if (CMPT_LE == nCmpType) {
1883 sName = constants().ids.int_.lin_le; // "int_lin_le";
1884 fDecl = mipd.int_lin_le;
1885 } else {
1886 sName = constants().ids.int_.lin_eq; // "int_lin_eq";
1887 fDecl = mipd.int_lin_eq;
1888 }
1889 }
1890 if (mipd.getEnv()->envi().cseMapEnd() != mipd.getEnv()->envi().cseMapFind(args[0])) {
1891 DBGOUT_MIPD_FLUSH(" Found expr ");
1892 DBGOUT_MIPD_SELF(debugprint(args[0]));
1893 }
1894 auto* nc = new Call(Location().introduce(), ASTString(sName), args);
1895 nc->type(Type::varbool());
1896 nc->decl(fDecl);
1897 mipd.getEnv()->envi().flatAddItem(new ConstraintI(Location().introduce(), nc));
1898 }
1899
1900 /// domain / reif set of one variable into that for a!her
1901 void convertIntSet(Expression* e, SetOfIntvReal& s, VarDecl* varTarget, double A, double B) {
1902 MZN_MIPD_assert_hard(A != 0.0);
1903 if (e->type().isIntSet()) {
1904 IntSetVal* S = eval_intset(mipd.getEnv()->envi(), e);
1905 IntSetRanges domr(S);
1906 for (; domr(); ++domr) { // * A + B
1907 IntVal mmin = domr.min();
1908 IntVal mmax = domr.max();
1909 if (A < 0.0) {
1910 std::swap(mmin, mmax);
1911 }
1912 s.insert(IntvReal( // * A + B
1913 mmin.isFinite() ? rndUpIfInt(varTarget, (static_cast<double>(mmin.toInt()) * A + B))
1914 : IntvReal::infMinus(),
1915 mmax.isFinite() ? rndDownIfInt(varTarget, (static_cast<double>(mmax.toInt()) * A + B))
1916 : IntvReal::infPlus()));
1917 }
1918 } else {
1919 assert(e->type().isFloatSet());
1920 FloatSetVal* S = eval_floatset(mipd.getEnv()->envi(), e);
1921 FloatSetRanges domr(S);
1922 for (; domr(); ++domr) { // * A + B
1923 FloatVal mmin = domr.min();
1924 FloatVal mmax = domr.max();
1925 if (A < 0.0) {
1926 std::swap(mmin, mmax);
1927 }
1928 s.insert(IntvReal( // * A + B
1929 mmin.isFinite() ? rndUpIfInt(varTarget, (mmin.toDouble() * A + B))
1930 : IntvReal::infMinus(),
1931 mmax.isFinite() ? rndDownIfInt(varTarget, (mmax.toDouble() * A + B))
1932 : IntvReal::infPlus()));
1933 }
1934 }
1935 }
1936
1937 /// compute the delta for float strict ineq
1938 static double computeDelta(VarDecl* var, VarDecl* varOrig, IntvReal bnds, double A, Call* pCall,
1939 int nArg) {
1940 double delta = varOrig->type().isfloat()
1941 ? MIPD::expr2Const(pCall->arg(nArg))
1942 // * ( bnds.right-bnds.left ) ABANDONED 12.4.18 due to #207
1943 : std::fabs(A); // delta should be scaled as well
1944 if (var->type().isint()) { // the projected-onto variable
1945 delta = std::max(1.0, delta);
1946 }
1947 return delta;
1948 }
1949 }; // class DomainDecomp
1950
1951 /// Vars without explicit clique still need a decomposition.
1952 /// Have !iced all _POSTs, set_in's && eq_encode's to it BEFORE
1953 /// In each clique, relate all vars to one chosen
1954 /// Find all "smallest rel. factor" variables, integer && with eq_encode if avail
1955 /// Re-relate all vars to it
1956 /// Refer all _POSTs && dom() to it
1957 /// build domain decomposition
1958 /// Implement all domain constraints, incl. possible corresp, of eq_encode's
1959 ///
1960 /// REMARKS.
1961 /// ! impose effects of integrality scaling (e.g., int v = int k/3)
1962 /// BUT when using k's eq_encode?
1963 /// And when subdividing into intervals
1964 bool decomposeDomains() {
1965 // for (int iClq=0; iClq<_aCliques.size(); ++iClq ) {
1966 // TClique& clq = _aCliques[iClq];
1967 // }
1968 bool fRetTrue = true;
1969 for (int iVar = 0; iVar < _vVarDescr.size(); ++iVar) {
1970 // VarDescr& var = _vVarDescr[iVar];
1971 if (_vVarDescr[iVar].fDomainConstrProcessed == 0U) {
1972 GCLock lock;
1973 DomainDecomp dd(this, iVar);
1974 dd.doProcess();
1975 _vVarDescr[iVar].fDomainConstrProcessed = 1U;
1976 }
1977 }
1978 // Clean up _POSTs:
1979 for (auto& vVar : _vVarDescr) {
1980 for (auto* pCallI : vVar.aCalls) {
1981 pCallI->remove();
1982 }
1983 if (vVar.pEqEncoding != nullptr) {
1984 vVar.pEqEncoding->remove();
1985 }
1986 }
1987 return fRetTrue;
1988 }
1989
1990 static VarDecl* expr2VarDecl(Expression* arg) {
1991 // The requirement to have actual variable objects
1992 // might be a limitation if more optimizations are done before...
1993 // Might need to flexibilize this TODO
1994 // MZN_MIPD_assert_hard_msg( ! arg->dynamicCast<IntLit>(),
1995 // "Expression " << *arg << " is an IntLit!" );
1996 // MZN_MIPD_assert_hard( ! arg->dynamicCast<FloatLit>() );
1997 // MZN_MIPD_assert_hard( ! arg->dynamicCast<BoolLit>() );
1998 Id* id = arg->dynamicCast<Id>();
1999 // MZN_MIPD_assert_hard(id);
2000 if (nullptr == id) {
2001 return nullptr; // the call using this should be ignored?
2002 }
2003 VarDecl* vd = id->decl();
2004 MZN_MIPD_assert_hard(vd);
2005 return vd;
2006 }
2007
2008 /// Fills the vector of vardecls && returns the least index of the array
2009 template <class Array>
2010 long long expr2DeclArray(Expression* arg, Array& aVD) {
2011 ArrayLit* al = eval_array_lit(getEnv()->envi(), arg);
2012 checkOrResize(aVD, al->size());
2013 for (unsigned int i = 0; i < al->size(); i++) {
2014 aVD[i] = expr2VarDecl((*al)[i]);
2015 }
2016 return al->min(0);
2017 }
2018
2019 /// Fills the vector of expressions && returns the least index of the array
2020 template <class Array>
2021 long long expr2ExprArray(Expression* arg, Array& aVD) {
2022 ArrayLit* al = eval_array_lit(getEnv()->envi(), arg);
2023 checkOrResize(aVD, al->size());
2024 for (unsigned int i = 0; i < al->size(); i++) {
2025 aVD[i] = ((*al)[i]);
2026 }
2027 return al->min(0);
2028 }
2029
2030 static double expr2Const(Expression* arg) {
2031 if (auto* il = arg->dynamicCast<IntLit>()) {
2032 return (static_cast<double>(il->v().toInt()));
2033 }
2034 if (auto* fl = arg->dynamicCast<FloatLit>()) {
2035 return (fl->v().toDouble());
2036 }
2037 if (auto* bl = arg->dynamicCast<BoolLit>()) {
2038 return static_cast<double>(bl->v());
2039 }
2040 MZN_MIPD_assert_hard_msg(0, "unexpected expression instead of an int/float/bool literal: eid="
2041 << arg->eid() << " while E_INTLIT=" << Expression::E_INTLIT);
2042
2043 return 0.0;
2044 }
2045
2046 template <class Container, class Elem = int, size_t = 0>
2047 void checkOrResize(Container& cnt, size_t sz) {
2048 cnt.resize(sz);
2049 }
2050
2051 template <class Elem, size_t N>
2052 void checkOrResize(std::array<Elem, N>& cnt, size_t sz) {
2053 MZN_MIPD_assert_hard(cnt.size() == sz);
2054 }
2055
2056 template <class Array>
2057 void expr2Array(Expression* arg, Array& vals) {
2058 ArrayLit* al = eval_array_lit(getEnv()->envi(), arg);
2059 // if ( typeid(typename Array::pointer) == typeid(typename Array::iterator) ) // fixed
2060 // array
2061 // MZN_MIPD_assert_hard( vals.size() == al->v().size() );
2062 // else
2063 // vals.resize( al->v().size() );
2064 checkOrResize(vals, al->size());
2065 for (unsigned int i = 0; i < al->size(); i++) {
2066 vals[i] = expr2Const((*al)[i]);
2067 }
2068 }
2069
2070 void printStats(std::ostream& os) {
2071 // if ( _aCliques.empty() )
2072 // return;
2073 if (_vVarDescr.empty()) {
2074 return;
2075 }
2076 int nc = 0;
2077 for (auto& cl : _aCliques) {
2078 if (!cl.empty()) {
2079 ++nc;
2080 }
2081 }
2082 for (auto& var : _vVarDescr) {
2083 if (0 > var.nClique) {
2084 ++nc; // 1-var cliques
2085 }
2086 }
2087 // os << "N cliques " << _aCliques.size() << " total, "
2088 // << nc << " final" << std::endl;
2089 MZN_MIPD_assert_hard(nc);
2090 MIPD_stats[N_POSTs_eqNmapsize] = static_cast<double>(_mNViews.size());
2091 double nSubintvAve = MIPD_stats[N_POSTs_NSubintvSum] / nc;
2092 MZN_MIPD_assert_hard(MIPD_stats[N_POSTs_NSubintvSum]);
2093 double dSubSizeAve = MIPD_stats[N_POSTs_SubSizeSum] / MIPD_stats[N_POSTs_NSubintvSum];
2094 os << MIPD_stats[N_POSTs_all]
2095 << " POSTs"
2096#ifdef MZN_MIPDOMAINS_PRINTMORESTATS
2097 " [ ";
2098 MZN_MIPDOMAINS_PRINTMORESTATS
2099 for (int i = N_POSTs_intCmpReif; i <= N_POSTs_floatAux; ++i) {
2100 os << MIPD_stats[i] << ',';
2101 }
2102 os << " ], LINEQ [ ";
2103 for (int i = N_POSTs_eq2intlineq; i <= N_POSTs_eqNmapsize; ++i) {
2104 os << MIPD_stats[i] << ',';
2105 }
2106 os << " ]"
2107#endif
2108 ", "
2109 << MIPD_stats[N_POSTs_varsDirect] << " / " << MIPD_stats[N_POSTs_varsInvolved] << " vars, "
2110 << nc << " cliques, " << MIPD_stats[N_POSTs_NSubintvMin] << " / " << nSubintvAve << " / "
2111 << MIPD_stats[N_POSTs_NSubintvMax] << " NSubIntv m/a/m, " << MIPD_stats[N_POSTs_SubSizeMin]
2112 << " / " << dSubSizeAve << " / " << MIPD_stats[N_POSTs_SubSizeMax] << " SubIntvSize m/a/m, "
2113 << MIPD_stats[N_POSTs_cliquesWithEqEncode] << "+" << MIPD_stats[N_POSTs_clEEEnforced] << "("
2114 << MIPD_stats[N_POSTs_clEEFound] << ")"
2115 << " clq eq_encoded ";
2116 // << std::flush
2117 if (TCliqueSorter::LinEqGraph::dCoefMax > 1.0) {
2118 os << TCliqueSorter::LinEqGraph::dCoefMin << "--" << TCliqueSorter::LinEqGraph::dCoefMax
2119 << " abs coefs";
2120 }
2121 os << " ... ";
2122 }
2123
2124}; // namespace MiniZinc
2125
2126template <class N>
2127template <class N1>
2128void SetOfIntervals<N>::intersect(const SetOfIntervals<N1>& s2) {
2129 if (s2.empty()) {
2130 this->clear();
2131 return;
2132 }
2133 this->cutOut(Interval<N>(Interval<N>::infMinus(), (N)s2.begin()->left));
2134 for (auto is2 = s2.begin(); is2 != s2.end(); ++is2) {
2135 auto is2next = is2;
2136 ++is2next;
2137 this->cutOut(
2138 Interval<N>(is2->right, s2.end() == is2next ? Interval<N>::infPlus() : (N)is2next->left));
2139 }
2140}
2141template <class N>
2142template <class N1>
2143void SetOfIntervals<N>::cutDeltas(const SetOfIntervals<N1>& s2, N1 delta) {
2144 if (this->empty()) {
2145 return;
2146 }
2147 // What if distance < delta? TODO
2148 for (auto is2 : s2) {
2149 if (is2.left > Interval<N1>::infMinus()) {
2150 this->cutOut(Interval<N>(is2.left - delta, is2.left));
2151 }
2152 if (is2.right < Interval<N1>::infPlus()) {
2153 this->cutOut(Interval<N>(is2.right, is2.right + delta));
2154 }
2155 }
2156}
2157template <class N>
2158void SetOfIntervals<N>::cutOut(const Interval<N>& intv) {
2159 DBGOUT_MIPD_FLUSH("Cutting " << intv << " from " << (*this));
2160 if (this->empty()) {
2161 return;
2162 }
2163 auto it1 = (Interval<N>::infMinus() == intv.left)
2164 ? this->lower_bound(Interval<N>(intv.left, intv.right))
2165 : this->upper_bound(Interval<N>(intv.left, intv.right));
2166 auto it2Del1 = it1; // from which to delete
2167 if (this->begin() != it1) {
2168 --it1;
2169 const N it1l = it1->left;
2170 MZN_MIPD_assert_hard(it1l <= intv.left);
2171 if (it1->right > intv.left) { // split it
2172 it2Del1 = split(it1, intv.left).second;
2173 // it1->right = intv.left; READ-ONLY
2174 // this->erase(it1);
2175 // it1 = this->end();
2176 // auto iR = this->insert( Interval<N>( it1l, intv.left ) );
2177 // MZN_MIPD_assert_hard( iR.second );
2178 }
2179 }
2180 DBGOUT_MIPD_FLUSH("; after split 1: " << (*this));
2181 // Processing the right end:
2182 auto it2 = this->lower_bound(Interval<N>(intv.right, intv.right + 1));
2183 auto it2Del2 = it2;
2184 if (this->begin() != it2) {
2185 --it2;
2186 MZN_MIPD_assert_hard(it2->left < intv.right);
2187 const N it2r = it2->right;
2188 if ((Interval<N>::infPlus() == intv.right) ? (it2r > intv.right)
2189 : (it2r >= intv.right)) { // >=: split it
2190 // it2Del2 = split( it2, intv.right ).second;
2191 const bool fEEE = (it2Del1 == it2);
2192 this->erase(it2);
2193 it2 = this->end();
2194 it2Del2 = this->insert(Interval<N>(intv.right, it2r));
2195 if (fEEE) {
2196 it2Del1 = it2Del2;
2197 }
2198 }
2199 }
2200 DBGOUT_MIPD_FLUSH("; after split 2: " << (*this));
2201 DBGOUT_MIPD_FLUSH("; cutting out: " << SetOfIntervals(it2Del1, it2Del2));
2202#ifdef MZN_DBG_CHECK_ITER_CUTOUT
2203 {
2204 auto it = this->begin();
2205 int nO = 0;
2206 do {
2207 if (it == it2Del1) {
2208 MZN_MIPD_assert_hard(!nO);
2209 ++nO;
2210 }
2211 if (it == it2Del2) {
2212 MZN_MIPD_assert_hard(1 == nO);
2213 ++nO;
2214 }
2215 if (this->end() == it) {
2216 break;
2217 }
2218 ++it;
2219 } while (true);
2220 MZN_MIPD_assert_hard(2 == nO);
2221 }
2222#endif
2223 this->erase(it2Del1, it2Del2);
2224 DBGOUT_MIPD(" ... gives " << (*this));
2225}
2226template <class N>
2227typename SetOfIntervals<N>::SplitResult SetOfIntervals<N>::split(iterator& it, N pos) {
2228 MZN_MIPD_assert_hard(pos >= it->left);
2229 MZN_MIPD_assert_hard(pos <= it->right);
2230 Interval<N> intvOld = *it;
2231 this->erase(it);
2232 auto it_01 = this->insert(Interval<N>(intvOld.left, pos));
2233 auto it_02 = this->insert(Interval<N>(pos, intvOld.right));
2234 it = this->end();
2235 return std::make_pair(it_01, it_02);
2236}
2237template <class N>
2238Interval<N> SetOfIntervals<N>::getBounds() const {
2239 if (this->empty()) {
2240 return Interval<N>(Interval<N>::infPlus(), Interval<N>::infMinus());
2241 }
2242 auto it2 = this->end();
2243 --it2;
2244 return Interval<N>(this->begin()->left, it2->right);
2245}
2246template <class N>
2247bool SetOfIntervals<N>::checkFiniteBounds() {
2248 if (this->empty()) {
2249 return false;
2250 }
2251 auto bnds = getBounds();
2252 return bnds.left > Interval<N>::infMinus() && bnds.right < Interval<N>::infPlus();
2253}
2254template <class N>
2255bool SetOfIntervals<N>::checkDisjunctStrict() {
2256 for (auto it = this->begin(); it != this->end(); ++it) {
2257 if (it->left > it->right) {
2258 return false;
2259 }
2260 if (this->begin() != it) {
2261 auto it_1 = it;
2262 --it_1;
2263 if (it_1->right >= it->left) {
2264 return false;
2265 }
2266 }
2267 }
2268 return true;
2269}
2270/// Assumes integer interval bounds
2271template <class N>
2272int SetOfIntervals<N>::cardInt() const {
2273 int nn = 0;
2274 for (auto it = this->begin(); it != this->end(); ++it) {
2275 ++nn;
2276 nn += int(round(it->right - it->left));
2277 }
2278 return nn;
2279}
2280template <class N>
2281N SetOfIntervals<N>::maxInterval() const {
2282 N ll = -1;
2283 for (auto it = this->begin(); it != this->end(); ++it) {
2284 ll = std::max(ll, it->right - it->left);
2285 }
2286 return ll;
2287}
2288/// Assumes integer interval bounds
2289template <class N>
2290void SetOfIntervals<N>::split2Bits() {
2291 Base bsNew;
2292 for (auto it = this->begin(); it != this->end(); ++it) {
2293 for (int v = static_cast<int>(round(it->left)); v <= round(it->right); ++v) {
2294 bsNew.insert(Intv(v, v));
2295 }
2296 }
2297 *(Base*)this = std::move(bsNew);
2298}
2299
2300bool MIPD::fVerbose = false;
2301
2302void mip_domains(Env& env, bool fVerbose, int nmi, double dmd) {
2303 MIPD mipd(&env, fVerbose, nmi, dmd);
2304 if (!mipd.doMIPdomains()) {
2305 GCLock lock;
2306 env.envi().fail();
2307 }
2308}
2309
2310double MIPD::TCliqueSorter::LinEqGraph::dCoefMin = +1e100;
2311double MIPD::TCliqueSorter::LinEqGraph::dCoefMax = -1e100;
2312
2313} // namespace MiniZinc