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