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#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__(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__(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,
530 set_in__POST);
531 _aCT.emplace_back("set_in_reif__POST", t_VISIVI, RIT_Reif, CT_SetIn, CMPT_None, VT_Int,
532 set_in_reif__POST);
533 /// Registering all declared & compatible __POST constraints
534 /// (First, cleanup FunctionIs' payload: -- ! doing now)
535 for (int i = 0; i < _aCT.size(); ++i) {
536 FunctionI* fi = env.model->matchFn(env, ASTString(_aCT[i].sFuncName), _aCT[i].aParams, false);
537 if (fi != nullptr) {
538 _mCallTypes[fi] = _aCT.data() + i;
539 _aCT[i].pfi = fi;
540 // fi->pPayload = (void*)this;
541 // std::cerr << " FOund declaration: " << _aCT[i].sFuncName << std::endl;
542 } else {
543 _aCT[i].pfi = nullptr;
544 DBGOUT_MIPD(" MIssing declaration: " << _aCT[i].sFuncName);
545 return false;
546 }
547 }
548 return true;
549 }
550
551 /// Registering all __POST calls' domain-constrained variables
552 void registerPOSTVariables() {
553 EnvI& env = getEnv()->envi();
554 GCLock lock;
555 Model& mFlat = *getEnv()->flat();
556 // First, cleanup VarDecls' payload which stores index in _vVarDescr
557 for (VarDeclIterator ivd = mFlat.vardecls().begin(); ivd != mFlat.vardecls().end(); ++ivd) {
558 ivd->e()->payload(-1);
559 }
560 // Now add variables with non-contiguous domain
561 for (VarDeclIterator ivd = mFlat.vardecls().begin(); ivd != mFlat.vardecls().end(); ++ivd) {
562 VarDecl* vd0 = ivd->e();
563 bool fNonCtg = false;
564 if (vd0->type().isint()) { // currently only for int vars TODO
565 if (Expression* eDom = vd0->ti()->domain()) {
566 IntSetVal* dom = eval_intset(env, eDom);
567 fNonCtg = (dom->size() > 1);
568 }
569 }
570 if (fNonCtg) {
571 DBGOUT_MIPD(" Variable " << vd0->id()->str() << ": non-contiguous domain "
572 << (*(vd0->ti()->domain())));
573 if (vd0->payload() == -1) { // ! yet visited
574 vd0->payload(static_cast<int>(_vVarDescr.size()));
575 _vVarDescr.emplace_back(vd0, vd0->type().isint()); // can use /prmTypes/ as well
576 if (vd0->e() != nullptr) {
577 checkInitExpr(vd0);
578 }
579 } else {
580 DBGOUT_MIPD__(" (already touched)");
581 }
582 ++MIPD__stats[N_POSTs__domain];
583 ++MIPD__stats[N_POSTs__all];
584 }
585 }
586 // Iterate thru original __POST constraints to mark constrained vars:
587 for (ConstraintIterator ic = mFlat.constraints().begin(); ic != mFlat.constraints().end();
588 ++ic) {
589 if (ic->removed()) {
590 continue;
591 }
592 if (Call* c = ic->e()->dynamicCast<Call>()) {
593 auto ipct = _mCallTypes.find(c->decl());
594 if (ipct != _mCallTypes.end()) {
595 // No ! here because might be deleted immediately in later versions.
596 // ic->remove(); // mark removed at once
597 MZN_MIPD__assert_hard(c->argCount() > 1);
598 ++MIPD__stats[N_POSTs__all];
599 VarDecl* vd0 = expr2VarDecl(c->arg(0));
600 if (nullptr == vd0) {
601 DBGOUT_MIPD__(" Call " << *c << ": 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__(" 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__(" (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__(" 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__(" 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__(" 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()
1040 << " + " << B);
1041 if (fReportRepeat) {
1042 MZN_MIPD__assert_soft(!fVerbose, "LinEqGraph: eqn between "
1043 << (*begV)->id()->str() << " && "
1044 << (*(begV + 1))->id()->str()
1045 << " is repeated. ");
1046 }
1047 return true;
1048 }
1049 }
1050 return false;
1051 }
1052 };
1053
1054 class LinEqGraph : public TMatrixVars {
1055 public:
1056 static double dCoefMin, dCoefMax;
1057 /// Stores the arc (x1, x2) as x1 = a*x2 + b
1058 /// so that a constraint on x2, say x2<=c <-> f,
1059 /// is equivalent to one for x1: x1 <=/>= a*c+b <-> f
1060 //// ( the other way involves division:
1061 //// so that a constraint on x1, say x1<=c <-> f,
1062 //// can easily be converted into one for x2 as a*x2 <= c-b <-> f
1063 //// <=> x2 (care for sign) (c-b)/a <-> f )
1064
1065 template <class ICoef, class IVarDecl>
1066 void addArc(ICoef begC, IVarDecl begV, double rhs) {
1067 MZN_MIPD__assert_soft(!fVerbose || std::fabs(*begC) >= 1e-10,
1068 " Vars " << (*begV)->id()->str() << " to "
1069 << (*(begV + 1))->id()->str() << ": coef=" << (*begC));
1070 // Transform Ax+By=C into x = -B/Ay+C/A
1071 const double negBA = -(*(begC + 1)) / (*begC);
1072 const double CA = rhs / (*begC);
1073 checkExistingArc(begV, negBA, CA);
1074 (*this)[*begV][*(begV + 1)] = std::make_pair(negBA, CA);
1075 const double dCoefAbs = std::fabs(negBA);
1076 if (dCoefAbs < dCoefMin) {
1077 dCoefMin = dCoefAbs;
1078 }
1079 if (dCoefAbs > dCoefMax) {
1080 dCoefMax = dCoefAbs;
1081 }
1082 }
1083 void addEdge(const LinEq2Vars& led) {
1084 addArc(led.coefs.begin(), led.vd.begin(), led.rhs);
1085 addArc(led.coefs.rbegin(), led.vd.rbegin(), led.rhs);
1086 }
1087 /// Propagate linear relations from the given variable
1088 void propagate(iterator itStart, TMapVars& mWhereStore) {
1089 MZN_MIPD__assert_hard(this->end() != itStart);
1090 TMatrixVars mTemp;
1091 mTemp[itStart->first] = itStart->second; // init with existing
1092 DBGOUT_MIPD("Propagation started from " << itStart->first->id()->str() << " having "
1093 << itStart->second.size() << " connections");
1094 propagate2(itStart, itStart, std::make_pair(1.0, 0.0), mTemp);
1095 mWhereStore = mTemp.begin()->second;
1096 MZN_MIPD__assert_hard_msg(
1097 mWhereStore.size() == this->size() - 1,
1098 "Variable " << (*(mTemp.begin()->first))
1099 << " should be connected to all others in the clique, but "
1100 << "|edges|==" << mWhereStore.size() << ", |all nodes|==" << this->size());
1101 }
1102 /// Propagate linear relations from it1 via it2
1103 void propagate2(iterator itSrc, iterator itVia, std::pair<double, double> rel,
1104 TMatrixVars& mWhereStore) {
1105 for (auto itDst = itVia->second.begin(); itDst != itVia->second.end(); ++itDst) {
1106 // Transform x1=A1x2+B1, x2=A2x3+B2 into x1=A1A2x3+A1B2+B1
1107 if (itDst->first == itSrc->first) {
1108 continue;
1109 }
1110 const double A1A2 = rel.first * itDst->second.first;
1111 const double A1B2plusB1 = rel.first * itDst->second.second + rel.second;
1112 bool fDive = true;
1113 if (itSrc != itVia) {
1114 PVarDecl vd[2] = {itSrc->first, itDst->first};
1115 if (!mWhereStore.checkExistingArc(vd, A1A2, A1B2plusB1, false)) {
1116 mWhereStore[vd[0]][vd[1]] = std::make_pair(A1A2, A1B2plusB1);
1117 DBGOUT_MIPD(" PROPAGATING: " << vd[0]->id()->str() << " = " << A1A2 << " * "
1118 << vd[1]->id()->str() << " + " << A1B2plusB1);
1119 } else {
1120 fDive = false;
1121 }
1122 }
1123 if (fDive) {
1124 auto itDST = this->find(itDst->first);
1125 MZN_MIPD__assert_hard(this->end() != itDST);
1126 propagate2(itSrc, itDST, std::make_pair(A1A2, A1B2plusB1), mWhereStore);
1127 }
1128 }
1129 }
1130 };
1131 LinEqGraph leg;
1132
1133 TCliqueSorter(MIPD* pm, int iv) : _mipd(*pm), _iVarStart(iv) {}
1134 void doRelate() {
1135 MZN_MIPD__assert_hard(_mipd._vVarDescr[_iVarStart].nClique >= 0);
1136 const TClique& clq = _mipd._aCliques[_mipd._vVarDescr[_iVarStart].nClique];
1137 for (const auto& eq2 : clq) {
1138 leg.addEdge(eq2);
1139 }
1140 DBGOUT_MIPD(" Clique " << _mipd._vVarDescr[_iVarStart].nClique << ": " << leg.size()
1141 << " variables, " << clq.size() << " connections.");
1142 for (auto& it1 : leg) {
1143 _mipd._vVarDescr[it1.first->payload()].fDomainConstrProcessed = 1U;
1144 }
1145
1146 // Propagate the 1st var's relations:
1147 leg.propagate(leg.begin(), mRef0);
1148
1149 // Find a best main variable according to:
1150 // 1. isInt 2. hasEqEncode 3. abs linFactor to ref0
1151 varRef1 = leg.begin()->first;
1152 std::array<double, 3> aCrit = {
1153 {(double)_mipd._vVarDescr[varRef1->payload()].fInt,
1154 static_cast<double>(_mipd._vVarDescr[varRef1->payload()].pEqEncoding != nullptr), 1.0}};
1155 for (auto& it2 : mRef0) {
1156 VarDescr& vard = _mipd._vVarDescr[it2.first->payload()];
1157 std::array<double, 3> aCrit1 = {{(double)vard.fInt,
1158 static_cast<double>(vard.pEqEncoding != nullptr),
1159 std::fabs(it2.second.first)}};
1160 if (aCrit1 > aCrit) {
1161 varRef1 = it2.first;
1162 aCrit = aCrit1;
1163 }
1164 }
1165 leg.propagate(leg.find(varRef1), mRef1);
1166 }
1167 }; // class TCliqueSorter
1168
1169 /// Build a domain decomposition for a clique
1170 /// a clique can consist of just 1 var without a clique object
1171 class DomainDecomp {
1172 public:
1173 MIPD& mipd;
1174 const int iVarStart;
1175 TCliqueSorter cls;
1176 SetOfIntvReal sDomain;
1177
1178 DomainDecomp(MIPD* pm, int iv) : mipd(*pm), iVarStart(iv), cls(pm, iv) {
1179 sDomain.insert(IntvReal()); // the decomposed domain. Init to +-inf
1180 }
1181 void doProcess() {
1182 // Choose the main variable && relate all others to it
1183 const int nClique = mipd._vVarDescr[iVarStart].nClique;
1184 if (nClique >= 0) {
1185 cls.doRelate();
1186 } else {
1187 cls.varRef1 = mipd._vVarDescr[iVarStart].vd;
1188 }
1189 // Adding itself:
1190 cls.mRef1[cls.varRef1] = std::make_pair(1.0, 0.0);
1191
1192 int iVarRef1 = cls.varRef1->payload();
1193 MZN_MIPD__assert_hard(nClique == mipd._vVarDescr[iVarRef1].nClique);
1194 cls.fRef1HasEqEncode = (mipd._vVarDescr[iVarRef1].pEqEncoding != nullptr);
1195
1196 // First, construct the domain decomposition in any case
1197 // projectVariableConstr( cls.varRef1, std::make_pair(1.0, 0.0) );
1198 // if ( nClique >= 0 ) {
1199 for (auto& iRef1 : cls.mRef1) {
1200 projectVariableConstr(iRef1.first, iRef1.second);
1201 }
1202
1203 DBGOUT_MIPD("Clique " << nClique << ": main ref var " << cls.varRef1->id()->str()
1204 << ", domain dec: " << sDomain);
1205
1206 MZN_MIPD__ASSERT_FOR_SAT(!sDomain.empty(), mipd.getEnv()->envi(), cls.varRef1->loc(),
1207 "clique " << nClique << ": main ref var " << *cls.varRef1->id()
1208 << ", domain decomposition seems empty: " << sDomain);
1209
1210 MZN_MIPD__FLATTENING_ERROR__IF_NOT(sDomain.checkFiniteBounds(), mipd.getEnv()->envi(),
1211 cls.varRef1->loc(),
1212 "variable " << *cls.varRef1->id()
1213 << " needs finite bounds for linearisation."
1214 " Or, use indicator constraints. "
1215 << "Current domain is " << sDomain);
1216
1217 MZN_MIPD__assert_hard(sDomain.checkDisjunctStrict());
1218
1219 makeRangeDomains();
1220
1221 // Then, use equality_encoding if available
1222 if (cls.fRef1HasEqEncode) {
1223 syncWithEqEncoding();
1224 syncOtherEqEncodings();
1225 } else { // ! cls.fRef1HasEqEncode
1226 if (sDomain.size() >= 2) { // need to simplify stuff otherwise
1227 considerDenseEncoding();
1228 createDomainFlags();
1229 }
1230 }
1231 implementPOSTs();
1232
1233 // Statistics
1234 if (sDomain.size() < MIPD__stats[N_POSTs__NSubintvMin]) {
1235 MIPD__stats[N_POSTs__NSubintvMin] = static_cast<double>(sDomain.size());
1236 }
1237 MIPD__stats[N_POSTs__NSubintvSum] += sDomain.size();
1238 if (sDomain.size() > MIPD__stats[N_POSTs__NSubintvMax]) {
1239 MIPD__stats[N_POSTs__NSubintvMax] = static_cast<double>(sDomain.size());
1240 }
1241 for (const auto& intv : sDomain) {
1242 const auto nSubSize = intv.right - intv.left;
1243 if (nSubSize < MIPD__stats[N_POSTs__SubSizeMin]) {
1244 MIPD__stats[N_POSTs__SubSizeMin] = nSubSize;
1245 }
1246 MIPD__stats[N_POSTs__SubSizeSum] += nSubSize;
1247 if (nSubSize > MIPD__stats[N_POSTs__SubSizeMax]) {
1248 MIPD__stats[N_POSTs__SubSizeMax] = nSubSize;
1249 }
1250 }
1251 if (cls.fRef1HasEqEncode) {
1252 ++MIPD__stats[N_POSTs__cliquesWithEqEncode];
1253 }
1254 }
1255
1256 /// Project the domain-related constraints of a variable into the clique
1257 /// Deltas should be scaled but to a minimum of the target's discr
1258 /// COmparison sense changes on negated vars
1259 void projectVariableConstr(VarDecl* vd, std::pair<double, double> eq1) {
1260 DBGOUT_MIPD__(" MIPD: projecting variable ");
1261 DBGOUT_MIPD_SELF(debugprint(vd));
1262 // Always check if domain becomes empty? TODO
1263 const double A = eq1.first; // vd = A*arg + B. conversion
1264 const double B = eq1.second;
1265 // process domain info
1266 double lb = B;
1267 double ub = A + B; // projected bounds for bool
1268 if (vd->ti()->domain() != nullptr) {
1269 if (vd->type().isint() || vd->type().isfloat()) { // INT VAR OR FLOAT VAR
1270 SetOfIntvReal sD1;
1271 convertIntSet(vd->ti()->domain(), sD1, cls.varRef1, A, B);
1272 sDomain.intersect(sD1);
1273 DBGOUT_MIPD(" Clique domain after proj of the init. domain "
1274 << sD1 << " of " << (vd->type().isint() ? "varint" : "varfloat") << A << " * "
1275 << vd->id()->str() << " + " << B << ": " << sDomain);
1276 auto bnds = sD1.getBounds();
1277 lb = bnds.left;
1278 ub = bnds.right;
1279 } else {
1280 MZN_MIPD__FLATTENING_ERROR__IF_NOT(0, mipd.getEnv()->envi(), cls.varRef1->loc(),
1281 "Variable " << vd->id()->str() << " of type "
1282 << vd->type().toString(mipd._env->envi())
1283 << " has a domain.");
1284 }
1285 // /// Deleting var domain:
1286 // vd->ti()->domain( NULL );
1287 } else {
1288 if (nullptr == vd->ti()->domain() && !vd->type().isbool()) {
1289 lb = IntvReal::infMinus();
1290 ub = IntvReal::infPlus();
1291 }
1292 }
1293 auto bnds = sDomain.getBounds(); // can change TODO
1294 // process calls. Can use the constr type info.
1295 auto& aCalls = mipd._vVarDescr[vd->payload()].aCalls;
1296 for (Item* pItem : aCalls) {
1297 auto* pCI = pItem->dynamicCast<ConstraintI>();
1298 MZN_MIPD__assert_hard(pCI != nullptr);
1299 Call* pCall = pCI->e()->dynamicCast<Call>();
1300 MZN_MIPD__assert_hard(pCall != nullptr);
1301 DBGOUT_MIPD__("PROPAG CALL ");
1302 DBGOUT_MIPD_SELF(debugprint(pCall));
1303 // check the bounds for bool in reifs? TODO
1304 auto ipct = mipd._mCallTypes.find(pCall->decl());
1305 MZN_MIPD__assert_hard(mipd._mCallTypes.end() != ipct);
1306 const DCT& dct = *ipct->second;
1307 int nCmpType_ADAPTED = dct.nCmpType;
1308 if (A < 0.0) { // negative factor
1309 if (std::abs(nCmpType_ADAPTED) >= 4) { // inequality
1310 nCmpType_ADAPTED = -nCmpType_ADAPTED;
1311 }
1312 }
1313 switch (dct.nConstrType) {
1314 case CT_SetIn: {
1315 SetOfIntvReal SS;
1316 convertIntSet(pCall->arg(1), SS, cls.varRef1, A, B);
1317 if (RIT_Static == dct.nReifType) {
1318 sDomain.intersect(SS);
1319 ++MIPD__stats[N_POSTs__setIn];
1320 } else {
1321 sDomain.cutDeltas(SS, std::max(1.0, std::fabs(A))); // deltas to scale
1322 ++MIPD__stats[N_POSTs__setInReif];
1323 }
1324 } break;
1325 case CT_Comparison:
1326 if (RIT_Reif == dct.nReifType) {
1327 const double rhs = (mipd.aux_float_lt_zero_iff_1__POST == pCall->decl())
1328 ? B /* + A*0.0, relating to 0 */
1329 // The 2nd argument is constant:
1330 : A * MIPD::expr2Const(pCall->arg(1)) + B;
1331 const double rhsUp = rndUpIfInt(cls.varRef1, rhs);
1332 const double rhsDown = rndDownIfInt(cls.varRef1, rhs);
1333 const double rhsRnd = rndIfInt(cls.varRef1, rhs);
1334 /// Strictly, for delta we should finish domain reductions first... TODO?
1335 const double delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 3);
1336 switch (nCmpType_ADAPTED) {
1337 case CMPT_LE:
1338 sDomain.cutDeltas(IntvReal::infMinus(), rhsDown, delta);
1339 break;
1340 case CMPT_GE:
1341 sDomain.cutDeltas(rhsUp, IntvReal::infPlus(), delta);
1342 break;
1343 case CMPT_LT_0:
1344 sDomain.cutDeltas(IntvReal::infMinus(), rhsDown - delta, delta);
1345 break;
1346 case CMPT_GT_0:
1347 sDomain.cutDeltas(rhsUp + delta, IntvReal::infPlus(), delta);
1348 break;
1349 case CMPT_EQ:
1350 if (!(cls.varRef1->type().isint() && // skip if int target var
1351 std::fabs(rhs - rhsRnd) > INT_EPS)) { // && fract value
1352 sDomain.cutDeltas(rhsRnd, rhsRnd, delta);
1353 }
1354 break;
1355 default:
1356 MZN_MIPD__assert_hard_msg(0, " No other reified cmp type ");
1357 }
1358 ++MIPD__stats[(vd->ti()->type().isint()) ? N_POSTs__intCmpReif
1359 : N_POSTs__floatCmpReif];
1360 } else if (RIT_Static == dct.nReifType) {
1361 // _ne, later maybe static ineq TODO
1362 MZN_MIPD__assert_hard(CMPT_NE == dct.nCmpType);
1363 const double rhs = A * MIPD::expr2Const(pCall->arg(1)) + B;
1364 const double rhsRnd = rndIfInt(cls.varRef1, rhs);
1365 bool fSkipNE = (cls.varRef1->type().isint() && std::fabs(rhs - rhsRnd) > INT_EPS);
1366 if (!fSkipNE) {
1367 const double delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 2);
1368 sDomain.cutOut({rhsRnd - delta, rhsRnd + delta});
1369 }
1370 ++MIPD__stats[(vd->ti()->type().isint()) ? N_POSTs__intNE : N_POSTs__floatNE];
1371 } else { // aux_ relate to 0.0
1372 // But we don't modify domain splitting for them currently
1373 ++MIPD__stats[(vd->ti()->type().isint()) ? N_POSTs__intAux : N_POSTs__floatAux];
1374 MZN_MIPD__assert_hard(RIT_Halfreif == dct.nReifType);
1375 // const double rhs = B; // + A*0
1376 // const double delta = vd->type().isint() ? 1.0 : 1e-5; //
1377 // TODO : eps
1378 }
1379 break;
1380 case CT_Encode:
1381 // See if any further constraints here? TODO
1382 ++MIPD__stats[N_POSTs__eq_encode];
1383 break;
1384 default:
1385 MZN_MIPD__assert_hard_msg(0, "Unknown constraint type");
1386 }
1387 }
1388 DBGOUT_MIPD(" Clique domain after proj of " << A << " * " << vd->id()->str() << " + " << B
1389 << ": " << sDomain);
1390 }
1391
1392 static double rndIfInt(VarDecl* vdTarget, double v) {
1393 return vdTarget->type().isint() ? std::round(v) : v;
1394 }
1395 static double rndIfBothInt(VarDecl* vdTarget, double v) {
1396 if (!vdTarget->type().isint()) {
1397 return v;
1398 }
1399 const double vRnd = std::round(v);
1400 return (fabs(v - vRnd) < INT_EPS) ? vRnd : v;
1401 }
1402 static double rndUpIfInt(VarDecl* vdTarget, double v) {
1403 return vdTarget->type().isint() ? std::ceil(v - INT_EPS) : v;
1404 }
1405 static double rndDownIfInt(VarDecl* vdTarget, double v) {
1406 return vdTarget->type().isint() ? std::floor(v + INT_EPS) : v;
1407 }
1408
1409 void makeRangeDomains() {
1410 auto bnds = sDomain.getBounds();
1411 for (auto& iRef1 : cls.mRef1) {
1412 VarDecl* vd = iRef1.first;
1413 // projecting the bounds back:
1414 double lb0 = (bnds.left - iRef1.second.second) / iRef1.second.first;
1415 double ub0 = (bnds.right - iRef1.second.second) / iRef1.second.first;
1416 if (lb0 > ub0) {
1417 MZN_MIPD__assert_hard(iRef1.second.first < 0.0);
1418 std::swap(lb0, ub0);
1419 }
1420 if (vd->type().isint()) {
1421 lb0 = rndUpIfInt(vd, lb0);
1422 ub0 = rndDownIfInt(vd, ub0);
1423 }
1424 setVarDomain(vd, lb0, ub0);
1425 }
1426 }
1427
1428 /// tightens element bounds in the existing eq_encoding of varRef1
1429 /// necessary because if one exists, int_ne is not translated into it
1430 /// Can also back-check from there? TODO
1431 /// And further checks TODO
1432 void syncWithEqEncoding() {
1433 std::vector<Expression*> pp;
1434 auto bnds = sDomain.getBounds();
1435 const long long iMin = mipd.expr2ExprArray(
1436 mipd._vVarDescr[cls.varRef1->payload()].pEqEncoding->e()->dynamicCast<Call>()->arg(1),
1437 pp);
1438 MZN_MIPD__assert_hard(pp.size() >= bnds.right - bnds.left + 1);
1439 MZN_MIPD__assert_hard(iMin <= bnds.left);
1440 long long vEE = iMin;
1441 DBGOUT_MIPD__(
1442 " SYNC EQ_ENCODE( "
1443 << (*cls.varRef1) << ", bitflags: "
1444 << *(mipd._vVarDescr[cls.varRef1->payload()].pEqEncoding->e()->dynamicCast<Call>()->arg(
1445 1))
1446 << " ): SETTING 0 FLAGS FOR VALUES: ");
1447 for (const auto& intv : sDomain) {
1448 for (; static_cast<double>(vEE) < intv.left; ++vEE) {
1449 if (vEE >= static_cast<long long>(iMin + pp.size())) {
1450 return;
1451 }
1452 if (pp[vEE - iMin]->isa<Id>()) {
1453 if (pp[vEE - iMin]->dynamicCast<Id>()->decl()->type().isvar()) {
1454 DBGOUT_MIPD__(vEE << ", ");
1455 setVarDomain(pp[vEE - iMin]->dynamicCast<Id>()->decl(), 0.0, 0.0);
1456 }
1457 }
1458 }
1459 vEE = static_cast<long long>(intv.right + 1);
1460 }
1461 for (; vEE < static_cast<long long>(iMin + pp.size()); ++vEE) {
1462 if (pp[vEE - iMin]->isa<Id>()) {
1463 if (pp[vEE - iMin]->dynamicCast<Id>()->decl()->type().isvar()) {
1464 DBGOUT_MIPD__(vEE << ", ");
1465 setVarDomain(pp[vEE - iMin]->dynamicCast<Id>()->decl(), 0.0, 0.0);
1466 }
1467 }
1468 }
1469 DBGOUT_MIPD("");
1470 }
1471
1472 /// sync varRef1's eq_encoding with those of other variables
1473 void syncOtherEqEncodings() {
1474 // TODO This could be in the var projection? No, need the final domain
1475 }
1476
1477 /// Depending on params,
1478 /// create an equality encoding for an integer variable
1479 /// TODO What if a float's domain is discrete?
1480 void considerDenseEncoding() {
1481 if (cls.varRef1->id()->type().isint()) {
1482 if (sDomain.maxInterval() <= mipd.nMaxIntv2Bits ||
1483 sDomain.cardInt() <= mipd.dMaxNValueDensity * sDomain.size()) {
1484 sDomain.split2Bits();
1485 ++MIPD__stats[N_POSTs__clEEEnforced];
1486 }
1487 }
1488 }
1489
1490 /// if ! eq_encoding, creates a flag for each subinterval in the domain
1491 /// && constrains sum(flags)==1
1492 void createDomainFlags() {
1493 std::vector<Expression*> vVars(sDomain.size()); // flags for each subinterval
1494 std::vector<double> vIntvLB(sDomain.size() + 1);
1495 std::vector<double> vIntvUB__(sDomain.size() + 1);
1496 int i = 0;
1497 double dMaxIntv = -1.0;
1498 for (const auto& intv : sDomain) {
1499 intv.varFlag = addIntVar(0.0, 1.0);
1500 vVars[i] = intv.varFlag->id();
1501 vIntvLB[i] = intv.left;
1502 vIntvUB__[i] = -intv.right;
1503 dMaxIntv = std::max(dMaxIntv, intv.right - intv.left);
1504 ++i;
1505 }
1506 // Sum of flags == 1
1507 std::vector<double> ones(sDomain.size(), 1.0);
1508 addLinConstr(ones, vVars, CMPT_EQ, 1.0);
1509 // Domain decomp
1510 vVars.push_back(cls.varRef1->id());
1511 vIntvLB[i] = -1.0; // var1 >= sum(LBi*flagi)
1512 /// STRICT equality encoding if small intervals
1513 if (dMaxIntv > 1e-6) { // EPS = param? TODO
1514 vIntvUB__[i] = 1.0; // var1 <= sum(UBi*flagi)
1515 addLinConstr(vIntvLB, vVars, CMPT_LE, 0.0);
1516 addLinConstr(vIntvUB__, vVars, CMPT_LE, 0.0);
1517 } else {
1518 ++MIPD__stats[N_POSTs__clEEFound];
1519 addLinConstr(vIntvLB, vVars, CMPT_EQ, 0.0);
1520 }
1521 }
1522
1523 /// deletes them as well
1524 void implementPOSTs() {
1525 auto bnds = sDomain.getBounds();
1526 for (auto& iRef1 : cls.mRef1) {
1527 // DBGOUT_MIPD__( " MIPD: implementing constraints of variable " );
1528 // DBGOUT_MIPD_SELF( debugprint(vd) );
1529 VarDecl* vd = iRef1.first;
1530 auto eq1 = iRef1.second;
1531 const double A = eq1.first; // vd = A*arg + B. conversion
1532 const double B = eq1.second;
1533 // process calls. Can use the constr type info.
1534 auto& aCalls = mipd._vVarDescr[vd->payload()].aCalls;
1535 for (Item* pItem : aCalls) {
1536 auto* pCI = pItem->dynamicCast<ConstraintI>();
1537 MZN_MIPD__assert_hard(pCI);
1538 Call* pCall = pCI->e()->dynamicCast<Call>();
1539 MZN_MIPD__assert_hard(pCall);
1540 DBGOUT_MIPD__("IMPL CALL ");
1541 DBGOUT_MIPD_SELF(debugprint(pCall));
1542 // check the bounds for bool in reifs? TODO
1543 auto ipct = mipd._mCallTypes.find(pCall->decl());
1544 MZN_MIPD__assert_hard(mipd._mCallTypes.end() != ipct);
1545 const DCT& dct = *ipct->second;
1546 int nCmpType_ADAPTED = dct.nCmpType;
1547 if (A < 0.0) { // negative factor
1548 if (std::abs(nCmpType_ADAPTED) >= 4) { // inequality
1549 nCmpType_ADAPTED = -nCmpType_ADAPTED;
1550 }
1551 }
1552 switch (dct.nConstrType) {
1553 case CT_SetIn:
1554 if (RIT_Reif == dct.nReifType) {
1555 SetOfIntvReal SS;
1556 convertIntSet(pCall->arg(1), SS, cls.varRef1, A, B);
1557 relateReifFlag(pCall->arg(2), SS);
1558 }
1559 break;
1560 case CT_Comparison:
1561 if (RIT_Reif == dct.nReifType) {
1562 const double rhs = (mipd.aux_float_lt_zero_iff_1__POST == pCall->decl())
1563 ? B /* + A*0.0, relating to 0 */
1564 // The 2nd argument is constant:
1565 : A * MIPD::expr2Const(pCall->arg(1)) + B;
1566 const double rhsUp = rndUpIfInt(cls.varRef1, rhs);
1567 const double rhsDown = rndDownIfInt(cls.varRef1, rhs);
1568 const double rhsRnd = rndIfBothInt(
1569 cls.varRef1, rhs); // if the ref var is int, need to round almost-int values
1570 const double delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 3);
1571 switch (nCmpType_ADAPTED) {
1572 case CMPT_LE:
1573 relateReifFlag(pCall->arg(2), {{IntvReal::infMinus(), rhsDown}});
1574 break;
1575 case CMPT_GE:
1576 relateReifFlag(pCall->arg(2), {{rhsUp, IntvReal::infPlus()}});
1577 break;
1578 case CMPT_LT_0:
1579 relateReifFlag(pCall->arg(1), {{IntvReal::infMinus(), rhsDown - delta}});
1580 break;
1581 case CMPT_GT_0:
1582 relateReifFlag(pCall->arg(1), {{rhsUp + delta, IntvReal::infPlus()}});
1583 break;
1584 case CMPT_EQ:
1585 relateReifFlag(pCall->arg(2), {{rhsRnd, rhsRnd}});
1586 break; // ... but if the value is sign. fractional for an int var, the flag is
1587 // set=0
1588 default:
1589 break;
1590 }
1591 } else if (RIT_Static == dct.nReifType) {
1592 // !hing here for NE
1593 MZN_MIPD__assert_hard(CMPT_NE == nCmpType_ADAPTED);
1594 } else { // aux_ relate to 0.0
1595 // But we don't modify domain splitting for them currently
1596 MZN_MIPD__assert_hard(RIT_Halfreif == dct.nReifType);
1597 double rhs = B; // + A*0
1598 const double rhsUp = rndUpIfInt(cls.varRef1, rhs);
1599 const double rhsDown = rndDownIfInt(cls.varRef1, rhs);
1600 const double rhsRnd = rndIfInt(cls.varRef1, rhs);
1601 double delta = 0.0;
1602 if (mipd.aux_float_lt_zero_if_1__POST == pCall->decl()) { // only float && lt
1603 delta = computeDelta(cls.varRef1, vd, bnds, A, pCall, 3);
1604 }
1605 if (nCmpType_ADAPTED < 0) {
1606 delta = -delta;
1607 }
1608 if (cls.varRef1->type().isint() && CMPT_EQ_0 != nCmpType_ADAPTED) {
1609 if (nCmpType_ADAPTED < 0) {
1610 rhs = rhsDown;
1611 } else {
1612 rhs = rhsUp;
1613 }
1614 } else {
1615 rhs += delta;
1616 }
1617 // Now we need rhs ! to be in the inner of the domain
1618 bool fUseDD = true;
1619 if (!cls.fRef1HasEqEncode) {
1620 switch (nCmpType_ADAPTED) {
1621 case CMPT_EQ_0: {
1622 auto itLB = sDomain.lower_bound(rhsRnd);
1623 fUseDD = (itLB->left == rhsRnd && itLB->right == rhsRnd); // exactly
1624 } break;
1625 case CMPT_LT_0:
1626 case CMPT_LE_0: {
1627 auto itUB = sDomain.upper_bound(rhsUp);
1628 bool fInner = false;
1629 if (sDomain.begin() != itUB) {
1630 --itUB;
1631 if (itUB->right > rhs) {
1632 fInner = true;
1633 }
1634 }
1635 fUseDD = !fInner;
1636 } break;
1637 case CMPT_GT_0:
1638 case CMPT_GE_0: {
1639 auto itLB = sDomain.lower_bound(rhsDown);
1640 bool fInner = false;
1641 if (sDomain.begin() != itLB) {
1642 --itLB;
1643 if (itLB->right >= rhs) {
1644 fInner = true;
1645 }
1646 }
1647 fUseDD = !fInner;
1648 } break;
1649 default:
1650 MZN_MIPD__assert_hard_msg(0, "Unknown halfreif cmp type");
1651 }
1652 }
1653 if (fUseDD) { // use sDomain
1654 if (CMPT_EQ_0 == nCmpType_ADAPTED) {
1655 relateReifFlag(pCall->arg(1), {{rhsRnd, rhsRnd}}, RIT_Halfreif);
1656 } else if (nCmpType_ADAPTED < 0) {
1657 relateReifFlag(pCall->arg(1), {{IntvReal::infMinus(), rhsDown}}, RIT_Halfreif);
1658 } else {
1659 relateReifFlag(pCall->arg(1), {{rhsUp, IntvReal::infPlus()}}, RIT_Halfreif);
1660 }
1661 } else { // use big-M
1662 DBGOUT_MIPD(" AUX BY BIG-Ms: ");
1663 const bool fLE = (CMPT_EQ_0 == nCmpType_ADAPTED || 0 > nCmpType_ADAPTED);
1664 const bool fGE = (CMPT_EQ_0 == nCmpType_ADAPTED || 0 < nCmpType_ADAPTED);
1665 // Take integer || float indicator version, depending on the constrained var:
1666 const int nIdxInd = // (VT_Int==dct.nVarType) ?
1667 // No: vd->ti()->type().isint() ? 1 : 2;
1668 cls.varRef1->ti()->type().isint()
1669 ? 1
1670 : 2; // need the type of the variable to be constr
1671 MZN_MIPD__assert_hard(static_cast<unsigned int>(nIdxInd) < pCall->argCount());
1672 Expression* pInd = pCall->arg(nIdxInd);
1673 if (fLE && rhs < bnds.right) {
1674 if (rhs >= bnds.left) {
1675 std::vector<double> coefs = {1.0, bnds.right - rhs};
1676 // Use the float version of indicator:
1677 std::vector<Expression*> vars = {cls.varRef1->id(), pInd};
1678 addLinConstr(coefs, vars, CMPT_LE, bnds.right);
1679 } else {
1680 setVarDomain(MIPD::expr2VarDecl(pInd), 0.0, 0.0);
1681 }
1682 }
1683 if (fGE && rhs > bnds.left) {
1684 if (rhs <= bnds.right) {
1685 std::vector<double> coefs = {-1.0, rhs - bnds.left};
1686 std::vector<Expression*> vars = {cls.varRef1->id(), pInd};
1687 addLinConstr(coefs, vars, CMPT_LE, -bnds.left);
1688 } else {
1689 setVarDomain(MIPD::expr2VarDecl(pInd), 0.0, 0.0);
1690 }
1691 }
1692 }
1693 }
1694 break;
1695 case CT_Encode:
1696 // See if any further constraints here? TODO
1697 break;
1698 default:
1699 MZN_MIPD__assert_hard_msg(0, "Unknown constraint type");
1700 }
1701 pItem->remove(); // removing the call
1702 }
1703 // removing the eq_encoding call
1704 if (mipd._vVarDescr[vd->payload()].pEqEncoding != nullptr) {
1705 mipd._vVarDescr[vd->payload()].pEqEncoding->remove();
1706 }
1707 }
1708 }
1709
1710 /// sets varFlag = || <= sum( intv.varFlag : SS )
1711 void relateReifFlag(Expression* expFlag, const SetOfIntvReal& SS, EnumReifType nRT = RIT_Reif) {
1712 MZN_MIPD__assert_hard(RIT_Reif == nRT || RIT_Halfreif == nRT);
1713 // MZN_MIPD__assert_hard( sDomain.size()>=2 );
1714 VarDecl* varFlag = MIPD::expr2VarDecl(expFlag);
1715 std::vector<Expression*> vIntvFlags;
1716 if (cls.fRef1HasEqEncode) { // use eq_encoding
1717 MZN_MIPD__assert_hard(varFlag->type().isint());
1718 std::vector<Expression*> pp;
1719 auto bnds = sDomain.getBounds();
1720 const long long iMin = mipd.expr2ExprArray(
1721 mipd._vVarDescr[cls.varRef1->payload()].pEqEncoding->e()->dynamicCast<Call>()->arg(1),
1722 pp);
1723 MZN_MIPD__assert_hard(pp.size() >= bnds.right - bnds.left + 1);
1724 MZN_MIPD__assert_hard(iMin <= bnds.left);
1725 for (const auto& intv : SS) {
1726 for (long long vv = (long long)std::max(double(iMin), ceil(intv.left));
1727 vv <= (long long)std::min(double(iMin) + pp.size() - 1, floor(intv.right)); ++vv) {
1728 vIntvFlags.push_back(pp[vv - iMin]);
1729 }
1730 }
1731 } else {
1732 MZN_MIPD__assert_hard(varFlag->type().isint());
1733 for (const auto& intv : SS) {
1734 auto it1 = sDomain.lower_bound(intv.left);
1735 auto it2 = sDomain.upper_bound(intv.right);
1736 auto it11 = it1;
1737 // Check that we are looking ! into a subinterval:
1738 if (sDomain.begin() != it11) {
1739 --it11;
1740 MZN_MIPD__assert_hard(it11->right < intv.left);
1741 }
1742 auto it12 = it2;
1743 if (sDomain.begin() != it12) {
1744 --it12;
1745 MZN_MIPD__assert_hard_msg(it12->right <= intv.right,
1746 " relateReifFlag for " << intv << " in " << sDomain);
1747 }
1748 for (it12 = it1; it12 != it2; ++it12) {
1749 if (it12->varFlag != nullptr) {
1750 vIntvFlags.push_back(it12->varFlag->id());
1751 } else {
1752 MZN_MIPD__assert_hard(1 == sDomain.size());
1753 vIntvFlags.push_back(IntLit::a(1)); // just a constant then
1754 }
1755 }
1756 }
1757 }
1758 if (!vIntvFlags.empty()) {
1759 // Could find out if reif is true -- TODO && see above for 1 subinterval
1760 std::vector<double> onesm(vIntvFlags.size(), -1.0);
1761 onesm.push_back(1.0);
1762 vIntvFlags.push_back(varFlag->id());
1763 EnumCmpType nCmpType = (RIT_Reif == nRT) ? CMPT_EQ : CMPT_LE;
1764 addLinConstr(onesm, vIntvFlags, nCmpType, 0.0);
1765 } else { // the reif is false
1766 setVarDomain(varFlag, 0.0, 0.0);
1767 }
1768 }
1769
1770 static void setVarDomain(VarDecl* vd, double lb, double ub) {
1771 // need to check if the new range is in the previous bounds... TODO
1772 if (vd->type().isfloat()) {
1773 // if ( 0.0==lb && 0.0==ub ) {
1774 auto* newDom =
1775 new BinOp(Location().introduce(), FloatLit::a(lb), BOT_DOTDOT, FloatLit::a(ub));
1776 vd->ti()->domain(newDom);
1777 DBGOUT_MIPD(" NULL OUT: " << vd->id()->str());
1778 // }
1779 } else if (vd->type().isint() || vd->type().isbool()) {
1780 auto* newDom = new SetLit(
1781 Location().introduce(),
1782 IntSetVal::a(static_cast<long long int>(lb), static_cast<long long int>(ub)));
1783 // TypeInst* nti = copy(mipd.getEnv()->envi(),varFlag->ti())->cast<TypeInst>();
1784 // nti->domain(newDom);
1785 vd->ti()->domain(newDom);
1786 } else {
1787 MZN_MIPD__assert_hard_msg(0, "Unknown var type ");
1788 }
1789 }
1790
1791 VarDecl* addIntVar(double LB, double UB) {
1792 // GCLock lock;
1793 // Cache them? Only location can be different TODO
1794 auto* newDom =
1795 new SetLit(Location().introduce(),
1796 IntSetVal::a(static_cast<long long int>(LB), static_cast<long long int>(UB)));
1797 auto* ti = new TypeInst(Location().introduce(), Type::varint(), newDom);
1798 auto* newVar = new VarDecl(Location().introduce(), ti, mipd.getEnv()->envi().genId());
1799 newVar->flat(newVar);
1800 mipd.getEnv()->envi().flatAddItem(new VarDeclI(Location().introduce(), newVar));
1801 return newVar;
1802 }
1803
1804 void addLinConstr(std::vector<double>& coefs, std::vector<Expression*>& vars,
1805 EnumCmpType nCmpType, double rhs) {
1806 std::vector<Expression*> args(3);
1807 MZN_MIPD__assert_hard(vars.size() >= 2);
1808 for (auto* v : vars) {
1809 MZN_MIPD__assert_hard(&v);
1810 // throw std::string("addLinConstr: &var=NULL");
1811 MZN_MIPD__assert_hard_msg(v->isa<Id>() || v->isa<IntLit>() || v->isa<FloatLit>(),
1812 " expression at " << (&v) << " eid = " << v->eid()
1813 << " while E_INTLIT=" << Expression::E_INTLIT);
1814 // throw std::string("addLinConstr: only id's as variables allowed");
1815 }
1816 MZN_MIPD__assert_hard(coefs.size() == vars.size());
1817 MZN_MIPD__assert_hard(CMPT_EQ == nCmpType || CMPT_LE == nCmpType);
1818 DBGOUT_MIPD_SELF( // LinEq leq; leq.coefs=coefs; leq.vd=vars; leq.rhs=rhs;
1819 DBGOUT_MIPD__(" ADDING " << (CMPT_EQ == nCmpType ? "LIN_EQ" : "LIN_LE") << ": [ ");
1820 for (auto c
1821 : coefs) DBGOUT_MIPD__(c << ',');
1822 DBGOUT_MIPD__(" ] * [ "); for (auto v
1823 : vars) {
1824 MZN_MIPD__assert_hard(!v->isa<VarDecl>());
1825 if (v->isa<Id>()) DBGOUT_MIPD__(v->dynamicCast<Id>()->str() << ',');
1826 // else if ( v->isa<VarDecl>() )
1827 // MZN_MIPD__assert_hard ("addLinConstr: only id's as variables allowed");
1828 else
1829 DBGOUT_MIPD__(mipd.expr2Const(v) << ',');
1830 } DBGOUT_MIPD(" ] " << (CMPT_EQ == nCmpType ? "== " : "<= ") << rhs););
1831 std::vector<Expression*> nc_c;
1832 std::vector<Expression*> nx;
1833 bool fFloat = false;
1834 for (auto* v : vars) {
1835 if (!v->type().isint()) {
1836 fFloat = true;
1837 break;
1838 }
1839 }
1840 auto sName = constants().ids.float_.lin_eq; // "int_lin_eq";
1841 FunctionI* fDecl = mipd.float_lin_eq;
1842 if (fFloat) { // MZN_MIPD__assert_hard all vars of same type TODO
1843 for (int i = 0; i < vars.size(); ++i) {
1844 if (fabs(coefs[i]) > 1e-8) /// Only add terms with non-0 coefs. TODO Eps=param
1845 {
1846 nc_c.push_back(FloatLit::a(coefs[i]));
1847 if (vars[i]->type().isint()) {
1848 std::vector<Expression*> i2f_args(1);
1849 i2f_args[0] = vars[i];
1850 Call* i2f = new Call(Location().introduce(), constants().ids.int2float, i2f_args);
1851 i2f->type(Type::varfloat());
1852 i2f->decl(mipd.getEnv()->model()->matchFn(mipd.getEnv()->envi(), i2f, false));
1853 EE ret = flat_exp(mipd.getEnv()->envi(), Ctx(), i2f, nullptr, constants().varTrue);
1854 nx.push_back(ret.r());
1855 } else {
1856 nx.push_back(vars[i]); // ->id(); once passing a general expression
1857 }
1858 }
1859 }
1860 args[2] = FloatLit::a(rhs);
1861 args[2]->type(Type::parfloat(0));
1862 args[0] = new ArrayLit(Location().introduce(), nc_c);
1863 args[0]->type(Type::parfloat(1));
1864 args[1] = new ArrayLit(Location().introduce(), nx);
1865 args[1]->type(Type::varfloat(1));
1866 if (CMPT_LE == nCmpType) {
1867 sName = constants().ids.float_.lin_le; // "float_lin_le";
1868 fDecl = mipd.float_lin_le;
1869 }
1870 } else {
1871 for (int i = 0; i < vars.size(); ++i) {
1872 if (fabs(coefs[i]) > 1e-8) /// Only add terms with non-0 coefs. TODO Eps=param
1873 {
1874 nc_c.push_back(IntLit::a(static_cast<long long int>(coefs[i])));
1875 nx.push_back(vars[i]); //->id();
1876 }
1877 }
1878 args[2] = IntLit::a(static_cast<long long int>(rhs));
1879 args[2]->type(Type::parint(0));
1880 args[0] = new ArrayLit(Location().introduce(), nc_c);
1881 args[0]->type(Type::parint(1));
1882 args[1] = new ArrayLit(Location().introduce(), nx);
1883 args[1]->type(Type::varint(1));
1884 if (CMPT_LE == nCmpType) {
1885 sName = constants().ids.int_.lin_le; // "int_lin_le";
1886 fDecl = mipd.int_lin_le;
1887 } else {
1888 sName = constants().ids.int_.lin_eq; // "int_lin_eq";
1889 fDecl = mipd.int_lin_eq;
1890 }
1891 }
1892 if (mipd.getEnv()->envi().cseMapEnd() != mipd.getEnv()->envi().cseMapFind(args[0])) {
1893 DBGOUT_MIPD__(" Found expr ");
1894 DBGOUT_MIPD_SELF(debugprint(args[0]));
1895 }
1896 auto* nc = new Call(Location().introduce(), ASTString(sName), args);
1897 nc->type(Type::varbool());
1898 nc->decl(fDecl);
1899 mipd.getEnv()->envi().flatAddItem(new ConstraintI(Location().introduce(), nc));
1900 }
1901
1902 /// domain / reif set of one variable into that for a!her
1903 void convertIntSet(Expression* e, SetOfIntvReal& s, VarDecl* varTarget, double A, double B) {
1904 MZN_MIPD__assert_hard(A != 0.0);
1905 if (e->type().isIntSet()) {
1906 IntSetVal* S = eval_intset(mipd.getEnv()->envi(), e);
1907 IntSetRanges domr(S);
1908 for (; domr(); ++domr) { // * A + B
1909 IntVal mmin = domr.min();
1910 IntVal mmax = domr.max();
1911 if (A < 0.0) {
1912 std::swap(mmin, mmax);
1913 }
1914 s.insert(IntvReal( // * A + B
1915 mmin.isFinite() ? rndUpIfInt(varTarget, (static_cast<double>(mmin.toInt()) * A + B))
1916 : IntvReal::infMinus(),
1917 mmax.isFinite() ? rndDownIfInt(varTarget, (static_cast<double>(mmax.toInt()) * A + B))
1918 : IntvReal::infPlus()));
1919 }
1920 } else {
1921 assert(e->type().isFloatSet());
1922 FloatSetVal* S = eval_floatset(mipd.getEnv()->envi(), e);
1923 FloatSetRanges domr(S);
1924 for (; domr(); ++domr) { // * A + B
1925 FloatVal mmin = domr.min();
1926 FloatVal mmax = domr.max();
1927 if (A < 0.0) {
1928 std::swap(mmin, mmax);
1929 }
1930 s.insert(IntvReal( // * A + B
1931 mmin.isFinite() ? rndUpIfInt(varTarget, (mmin.toDouble() * A + B))
1932 : IntvReal::infMinus(),
1933 mmax.isFinite() ? rndDownIfInt(varTarget, (mmax.toDouble() * A + B))
1934 : IntvReal::infPlus()));
1935 }
1936 }
1937 }
1938
1939 /// compute the delta for float strict ineq
1940 static double computeDelta(VarDecl* var, VarDecl* varOrig, IntvReal bnds, double A, Call* pCall,
1941 int nArg) {
1942 double delta = varOrig->type().isfloat()
1943 ? MIPD::expr2Const(pCall->arg(nArg))
1944 // * ( bnds.right-bnds.left ) ABANDONED 12.4.18 due to #207
1945 : std::fabs(A); // delta should be scaled as well
1946 if (var->type().isint()) { // the projected-onto variable
1947 delta = std::max(1.0, delta);
1948 }
1949 return delta;
1950 }
1951 }; // class DomainDecomp
1952
1953 /// Vars without explicit clique still need a decomposition.
1954 /// Have !iced all __POSTs, set_in's && eq_encode's to it BEFORE
1955 /// In each clique, relate all vars to one chosen
1956 /// Find all "smallest rel. factor" variables, integer && with eq_encode if avail
1957 /// Re-relate all vars to it
1958 /// Refer all __POSTs && dom() to it
1959 /// build domain decomposition
1960 /// Implement all domain constraints, incl. possible corresp, of eq_encode's
1961 ///
1962 /// REMARKS.
1963 /// ! impose effects of integrality scaling (e.g., int v = int k/3)
1964 /// BUT when using k's eq_encode?
1965 /// And when subdividing into intervals
1966 bool decomposeDomains() {
1967 // for (int iClq=0; iClq<_aCliques.size(); ++iClq ) {
1968 // TClique& clq = _aCliques[iClq];
1969 // }
1970 bool fRetTrue = true;
1971 for (int iVar = 0; iVar < _vVarDescr.size(); ++iVar) {
1972 // VarDescr& var = _vVarDescr[iVar];
1973 if (_vVarDescr[iVar].fDomainConstrProcessed == 0U) {
1974 GCLock lock;
1975 DomainDecomp dd(this, iVar);
1976 dd.doProcess();
1977 _vVarDescr[iVar].fDomainConstrProcessed = 1U;
1978 }
1979 }
1980 // Clean up __POSTs:
1981 for (auto& vVar : _vVarDescr) {
1982 for (auto* pCallI : vVar.aCalls) {
1983 pCallI->remove();
1984 }
1985 if (vVar.pEqEncoding != nullptr) {
1986 vVar.pEqEncoding->remove();
1987 }
1988 }
1989 return fRetTrue;
1990 }
1991
1992 static VarDecl* expr2VarDecl(Expression* arg) {
1993 // The requirement to have actual variable objects
1994 // might be a limitation if more optimizations are done before...
1995 // Might need to flexibilize this TODO
1996 // MZN_MIPD__assert_hard_msg( ! arg->dynamicCast<IntLit>(),
1997 // "Expression " << *arg << " is an IntLit!" );
1998 // MZN_MIPD__assert_hard( ! arg->dynamicCast<FloatLit>() );
1999 // MZN_MIPD__assert_hard( ! arg->dynamicCast<BoolLit>() );
2000 Id* id = arg->dynamicCast<Id>();
2001 // MZN_MIPD__assert_hard(id);
2002 if (nullptr == id) {
2003 return nullptr; // the call using this should be ignored?
2004 }
2005 VarDecl* vd = id->decl();
2006 MZN_MIPD__assert_hard(vd);
2007 return vd;
2008 }
2009
2010 /// Fills the vector of vardecls && returns the least index of the array
2011 template <class Array>
2012 long long expr2DeclArray(Expression* arg, Array& aVD) {
2013 ArrayLit* al = eval_array_lit(getEnv()->envi(), arg);
2014 checkOrResize(aVD, al->size());
2015 for (unsigned int i = 0; i < al->size(); i++) {
2016 aVD[i] = expr2VarDecl((*al)[i]);
2017 }
2018 return al->min(0);
2019 }
2020
2021 /// Fills the vector of expressions && returns the least index of the array
2022 template <class Array>
2023 long long expr2ExprArray(Expression* arg, Array& aVD) {
2024 ArrayLit* al = eval_array_lit(getEnv()->envi(), arg);
2025 checkOrResize(aVD, al->size());
2026 for (unsigned int i = 0; i < al->size(); i++) {
2027 aVD[i] = ((*al)[i]);
2028 }
2029 return al->min(0);
2030 }
2031
2032 static double expr2Const(Expression* arg) {
2033 if (auto* il = arg->dynamicCast<IntLit>()) {
2034 return (static_cast<double>(il->v().toInt()));
2035 }
2036 if (auto* fl = arg->dynamicCast<FloatLit>()) {
2037 return (fl->v().toDouble());
2038 }
2039 if (auto* bl = arg->dynamicCast<BoolLit>()) {
2040 return static_cast<double>(bl->v());
2041 }
2042 MZN_MIPD__assert_hard_msg(0, "unexpected expression instead of an int/float/bool literal: eid="
2043 << arg->eid() << " while E_INTLIT=" << Expression::E_INTLIT);
2044
2045 return 0.0;
2046 }
2047
2048 template <class Container, class Elem = int, size_t = 0>
2049 void checkOrResize(Container& cnt, size_t sz) {
2050 cnt.resize(sz);
2051 }
2052
2053 template <class Elem, size_t N>
2054 void checkOrResize(std::array<Elem, N>& cnt, size_t sz) {
2055 MZN_MIPD__assert_hard(cnt.size() == sz);
2056 }
2057
2058 template <class Array>
2059 void expr2Array(Expression* arg, Array& vals) {
2060 ArrayLit* al = eval_array_lit(getEnv()->envi(), arg);
2061 // if ( typeid(typename Array::pointer) == typeid(typename Array::iterator) ) // fixed
2062 // array
2063 // MZN_MIPD__assert_hard( vals.size() == al->v().size() );
2064 // else
2065 // vals.resize( al->v().size() );
2066 checkOrResize(vals, al->size());
2067 for (unsigned int i = 0; i < al->size(); i++) {
2068 vals[i] = expr2Const((*al)[i]);
2069 }
2070 }
2071
2072 void printStats(std::ostream& os) {
2073 // if ( _aCliques.empty() )
2074 // return;
2075 if (_vVarDescr.empty()) {
2076 return;
2077 }
2078 int nc = 0;
2079 for (auto& cl : _aCliques) {
2080 if (!cl.empty()) {
2081 ++nc;
2082 }
2083 }
2084 for (auto& var : _vVarDescr) {
2085 if (0 > var.nClique) {
2086 ++nc; // 1-var cliques
2087 }
2088 }
2089 // os << "N cliques " << _aCliques.size() << " total, "
2090 // << nc << " final" << std::endl;
2091 MZN_MIPD__assert_hard(nc);
2092 MIPD__stats[N_POSTs__eqNmapsize] = static_cast<double>(_mNViews.size());
2093 double nSubintvAve = MIPD__stats[N_POSTs__NSubintvSum] / nc;
2094 MZN_MIPD__assert_hard(MIPD__stats[N_POSTs__NSubintvSum]);
2095 double dSubSizeAve = MIPD__stats[N_POSTs__SubSizeSum] / MIPD__stats[N_POSTs__NSubintvSum];
2096 os << MIPD__stats[N_POSTs__all]
2097 << " POSTs"
2098#ifdef __MZN__MIPDOMAINS__PRINTMORESTATS
2099 " [ ";
2100 for (int i = N_POSTs__intCmpReif; i <= N_POSTs__floatAux; ++i) {
2101 os << MIPD__stats[i] << ',';
2102 }
2103 os << " ], LINEQ [ ";
2104 for (int i = N_POSTs__eq2intlineq; i <= N_POSTs__eqNmapsize; ++i) {
2105 os << MIPD__stats[i] << ',';
2106 }
2107 os << " ]"
2108#endif
2109 ", "
2110 << MIPD__stats[N_POSTs__varsDirect] << " / " << MIPD__stats[N_POSTs__varsInvolved]
2111 << " vars, " << nc << " cliques, " << MIPD__stats[N_POSTs__NSubintvMin] << " / "
2112 << nSubintvAve << " / " << MIPD__stats[N_POSTs__NSubintvMax] << " NSubIntv m/a/m, "
2113 << MIPD__stats[N_POSTs__SubSizeMin] << " / " << dSubSizeAve << " / "
2114 << MIPD__stats[N_POSTs__SubSizeMax] << " SubIntvSize m/a/m, "
2115 << MIPD__stats[N_POSTs__cliquesWithEqEncode] << "+" << MIPD__stats[N_POSTs__clEEEnforced]
2116 << "(" << MIPD__stats[N_POSTs__clEEFound] << ")"
2117 << " clq eq_encoded ";
2118 // << std::flush
2119 if (TCliqueSorter::LinEqGraph::dCoefMax > 1.0) {
2120 os << TCliqueSorter::LinEqGraph::dCoefMin << "--" << TCliqueSorter::LinEqGraph::dCoefMax
2121 << " abs coefs";
2122 }
2123 os << " ... ";
2124 }
2125
2126}; // namespace MiniZinc
2127
2128template <class N>
2129template <class N1>
2130void SetOfIntervals<N>::intersect(const SetOfIntervals<N1>& s2) {
2131 if (s2.empty()) {
2132 this->clear();
2133 return;
2134 }
2135 this->cutOut(Interval<N>(Interval<N>::infMinus(), (N)s2.begin()->left));
2136 for (auto is2 = s2.begin(); is2 != s2.end(); ++is2) {
2137 auto is2next = is2;
2138 ++is2next;
2139 this->cutOut(
2140 Interval<N>(is2->right, s2.end() == is2next ? Interval<N>::infPlus() : (N)is2next->left));
2141 }
2142}
2143template <class N>
2144template <class N1>
2145void SetOfIntervals<N>::cutDeltas(const SetOfIntervals<N1>& s2, N1 delta) {
2146 if (this->empty()) {
2147 return;
2148 }
2149 // What if distance < delta? TODO
2150 for (auto is2 : s2) {
2151 if (is2.left > Interval<N1>::infMinus()) {
2152 this->cutOut(Interval<N>(is2.left - delta, is2.left));
2153 }
2154 if (is2.right < Interval<N1>::infPlus()) {
2155 this->cutOut(Interval<N>(is2.right, is2.right + delta));
2156 }
2157 }
2158}
2159template <class N>
2160void SetOfIntervals<N>::cutOut(const Interval<N>& intv) {
2161 DBGOUT_MIPD__("Cutting " << intv << " from " << (*this));
2162 if (this->empty()) {
2163 return;
2164 }
2165 auto it1 = (Interval<N>::infMinus() == intv.left)
2166 ? this->lower_bound(Interval<N>(intv.left, intv.right))
2167 : this->upper_bound(Interval<N>(intv.left, intv.right));
2168 auto it2Del1 = it1; // from which to delete
2169 if (this->begin() != it1) {
2170 --it1;
2171 const N it1l = it1->left;
2172 MZN_MIPD__assert_hard(it1l <= intv.left);
2173 if (it1->right > intv.left) { // split it
2174 it2Del1 = split(it1, intv.left).second;
2175 // it1->right = intv.left; READ-ONLY
2176 // this->erase(it1);
2177 // it1 = this->end();
2178 // auto iR = this->insert( Interval<N>( it1l, intv.left ) );
2179 // MZN_MIPD__assert_hard( iR.second );
2180 }
2181 }
2182 DBGOUT_MIPD__("; after split 1: " << (*this));
2183 // Processing the right end:
2184 auto it2 = this->lower_bound(Interval<N>(intv.right, intv.right + 1));
2185 auto it2Del2 = it2;
2186 if (this->begin() != it2) {
2187 --it2;
2188 MZN_MIPD__assert_hard(it2->left < intv.right);
2189 const N it2r = it2->right;
2190 if ((Interval<N>::infPlus() == intv.right) ? (it2r > intv.right)
2191 : (it2r >= intv.right)) { // >=: split it
2192 // it2Del2 = split( it2, intv.right ).second;
2193 const bool fEEE = (it2Del1 == it2);
2194 this->erase(it2);
2195 it2 = this->end();
2196 it2Del2 = this->insert(Interval<N>(intv.right, it2r));
2197 if (fEEE) {
2198 it2Del1 = it2Del2;
2199 }
2200 }
2201 }
2202 DBGOUT_MIPD__("; after split 2: " << (*this));
2203 DBGOUT_MIPD__("; cutting out: " << SetOfIntervals(it2Del1, it2Del2));
2204#ifdef MZN_DBG_CHECK_ITER_CUTOUT
2205 {
2206 auto it = this->begin();
2207 int nO = 0;
2208 do {
2209 if (it == it2Del1) {
2210 MZN_MIPD__assert_hard(!nO);
2211 ++nO;
2212 }
2213 if (it == it2Del2) {
2214 MZN_MIPD__assert_hard(1 == nO);
2215 ++nO;
2216 }
2217 if (this->end() == it) {
2218 break;
2219 }
2220 ++it;
2221 } while (true);
2222 MZN_MIPD__assert_hard(2 == nO);
2223 }
2224#endif
2225 this->erase(it2Del1, it2Del2);
2226 DBGOUT_MIPD(" ... gives " << (*this));
2227}
2228template <class N>
2229typename SetOfIntervals<N>::SplitResult SetOfIntervals<N>::split(iterator& it, N pos) {
2230 MZN_MIPD__assert_hard(pos >= it->left);
2231 MZN_MIPD__assert_hard(pos <= it->right);
2232 Interval<N> intvOld = *it;
2233 this->erase(it);
2234 auto it_01 = this->insert(Interval<N>(intvOld.left, pos));
2235 auto it_02 = this->insert(Interval<N>(pos, intvOld.right));
2236 it = this->end();
2237 return std::make_pair(it_01, it_02);
2238}
2239template <class N>
2240Interval<N> SetOfIntervals<N>::getBounds() const {
2241 if (this->empty()) {
2242 return Interval<N>(Interval<N>::infPlus(), Interval<N>::infMinus());
2243 }
2244 auto it2 = this->end();
2245 --it2;
2246 return Interval<N>(this->begin()->left, it2->right);
2247}
2248template <class N>
2249bool SetOfIntervals<N>::checkFiniteBounds() {
2250 if (this->empty()) {
2251 return false;
2252 }
2253 auto bnds = getBounds();
2254 return bnds.left > Interval<N>::infMinus() && bnds.right < Interval<N>::infPlus();
2255}
2256template <class N>
2257bool SetOfIntervals<N>::checkDisjunctStrict() {
2258 for (auto it = this->begin(); it != this->end(); ++it) {
2259 if (it->left > it->right) {
2260 return false;
2261 }
2262 if (this->begin() != it) {
2263 auto it_1 = it;
2264 --it_1;
2265 if (it_1->right >= it->left) {
2266 return false;
2267 }
2268 }
2269 }
2270 return true;
2271}
2272/// Assumes integer interval bounds
2273template <class N>
2274int SetOfIntervals<N>::cardInt() const {
2275 int nn = 0;
2276 for (auto it = this->begin(); it != this->end(); ++it) {
2277 ++nn;
2278 nn += int(round(it->right - it->left));
2279 }
2280 return nn;
2281}
2282template <class N>
2283N SetOfIntervals<N>::maxInterval() const {
2284 N ll = -1;
2285 for (auto it = this->begin(); it != this->end(); ++it) {
2286 ll = std::max(ll, it->right - it->left);
2287 }
2288 return ll;
2289}
2290/// Assumes integer interval bounds
2291template <class N>
2292void SetOfIntervals<N>::split2Bits() {
2293 Base bsNew;
2294 for (auto it = this->begin(); it != this->end(); ++it) {
2295 for (int v = static_cast<int>(round(it->left)); v <= round(it->right); ++v) {
2296 bsNew.insert(Intv(v, v));
2297 }
2298 }
2299 *(Base*)this = std::move(bsNew);
2300}
2301
2302bool MIPD::fVerbose = false;
2303
2304void mip_domains(Env& env, bool fVerbose, int nmi, double dmd) {
2305 MIPD mipd(&env, fVerbose, nmi, dmd);
2306 if (!mipd.doMIPdomains()) {
2307 GCLock lock;
2308 env.envi().fail();
2309 }
2310}
2311
2312double MIPD::TCliqueSorter::LinEqGraph::dCoefMin = +1e100;
2313double MIPD::TCliqueSorter::LinEqGraph::dCoefMax = -1e100;
2314
2315} // namespace MiniZinc