Branch data Line data Source code
1 : : /******************************************************************************
2 : : * Top contributors (to current version):
3 : : * Tim King, Andrew V. Teylu, Gereon Kremer
4 : : *
5 : : * This file is part of the cvc5 project.
6 : : *
7 : : * Copyright (c) 2009-2025 by the authors listed in the file AUTHORS
8 : : * in the top-level source directory and their institutional affiliations.
9 : : * All rights reserved. See the file COPYING in the top-level source
10 : : * directory for licensing information.
11 : : * ****************************************************************************
12 : : *
13 : : * [[ Add one-line brief description here ]]
14 : : *
15 : : * [[ Add lengthier description here ]]
16 : : * \todo document this file
17 : : */
18 : : #include "theory/arith/linear/approx_simplex.h"
19 : :
20 : : #include <math.h>
21 : :
22 : : #include <cfloat>
23 : : #include <cmath>
24 : : #include <unordered_set>
25 : :
26 : : #include "base/cvc5config.h"
27 : : #include "base/output.h"
28 : : #include "proof/eager_proof_generator.h"
29 : : #include "theory/arith/linear/constraint.h"
30 : : #include "theory/arith/linear/cut_log.h"
31 : : #include "theory/arith/linear/matrix.h"
32 : : #include "theory/arith/linear/normal_form.h"
33 : : #include "util/statistics_registry.h"
34 : :
35 : : #ifdef CVC5_USE_GLPK
36 : : #include "theory/arith/linear/partial_model.h"
37 : :
38 : : extern "C" {
39 : : #include <glpk.h>
40 : : } /* extern "C" */
41 : :
42 : : namespace cvc5::internal {
43 : : namespace theory {
44 : : namespace arith::linear {
45 : :
46 : : struct AuxInfo {
47 : : TreeLog* tl;
48 : : int pivotLimit;
49 : : int branchLimit;
50 : : int branchDepth;
51 : : MipResult term; /* terminatation */
52 : : };
53 : :
54 : : enum SlackReplace { SlackUndef=0, SlackLB, SlackUB, SlackVLB, SlackVUB };
55 : :
56 : : struct VirtualBound {
57 : : // Either x <= d * y or x >= d * y
58 : : ArithVar x; // variable being bounded
59 : : Kind k; // either LEQ or GEQ
60 : : Rational d; // the multiple on y
61 : : ArithVar y; // the variable that is the upper bound
62 : : ConstraintP c; // the original constraint relating x and y
63 : :
64 : 0 : VirtualBound()
65 : 0 : : x(ARITHVAR_SENTINEL)
66 : : , k(Kind::UNDEFINED_KIND)
67 : : , d()
68 : : , y(ARITHVAR_SENTINEL)
69 : 0 : , c(NullConstraint)
70 : 0 : {}
71 : 0 : VirtualBound(ArithVar toBound, Kind rel, const Rational& coeff, ArithVar bounding, ConstraintP orig)
72 : 0 : : x(toBound)
73 : : , k(rel)
74 : : , d(coeff)
75 : : , y(bounding)
76 : 0 : , c(orig)
77 : : {
78 : 0 : Assert(k == Kind::LEQ || k == Kind::GEQ);
79 : 0 : }
80 : : };
81 : :
82 : : struct CutScratchPad {
83 : : bool d_failure; // if the construction was unsuccessful
84 : :
85 : : /* GOMORY CUTS Datastructures */
86 : : ArithVar d_basic; // a variable that is basic in the approximate solver
87 : : DenseVector d_tabRow; // a row in the tableau not including d_basic, equal to 0
88 : : DenseMap<ConstraintP> d_toBound; // each variable in toBound maps each variable in tabRow to either an upper/lower bound
89 : :
90 : : /* MIR CUTS Datastructures */
91 : : DenseMap<SlackReplace> d_slacks;// The x'[i] selected for x[i]
92 : : DenseMap<VirtualBound> d_vub; // Virtual upper bounds.
93 : : DenseMap<VirtualBound> d_vlb; // Virtual lower bounds.
94 : : DenseMap<Rational> d_compRanges;
95 : :
96 : : // a sum of rows in the tableau, with possible replacements for fixed
97 : : // sum aggLhs[i] x[i] = aggRhs;
98 : : DenseVector d_agg;
99 : : // Takes agg and replaces x[i] with a slack variable x'[i]
100 : : // Takes agg and replaces x[i] with a slack variable x'[i]
101 : : // sum modLhs[i] x'[i] = modRhs;
102 : : DenseVector d_mod;
103 : :
104 : : // Takes mod, and performs c-Mir on it
105 : : // sum alpha[i] x'[i] <= beta
106 : : DenseVector d_alpha;
107 : :
108 : : /* The constructed cut */
109 : : // sum cut[i] x[i] <= cutRhs
110 : : DenseVector d_cut;
111 : : Kind d_cutKind;
112 : :
113 : : /* The constraints used throughout construction. */
114 : : std::set<ConstraintP> d_explanation; // use pointer equality
115 : 6 : CutScratchPad(){
116 : 6 : clear();
117 : 6 : }
118 : 6 : void clear(){
119 : 6 : d_failure = false;
120 : 6 : d_basic = ARITHVAR_SENTINEL;
121 : 6 : d_tabRow.purge();
122 : 6 : d_toBound.purge();
123 : :
124 : 6 : d_slacks.purge();
125 : 6 : d_vub.purge();
126 : 6 : d_vlb.purge();
127 : 6 : d_compRanges.purge();
128 : :
129 : 6 : d_agg.purge();
130 : 6 : d_mod.purge();
131 : 6 : d_alpha.purge();
132 : :
133 : 6 : d_cut.purge();
134 : 6 : d_cutKind = Kind::UNDEFINED_KIND;
135 : 6 : d_explanation.clear();
136 : 6 : }
137 : : };
138 : :
139 : : struct GmiInfo;
140 : : struct MirInfo;
141 : : class BranchCutInfo;
142 : :
143 : : class ApproxGLPK : public ApproximateSimplex
144 : : {
145 : : public:
146 : : ApproxGLPK(const ArithVariables& v, TreeLog& l, ApproximateStatistics& s);
147 : : ~ApproxGLPK();
148 : :
149 : : LinResult solveRelaxation() override;
150 : 0 : Solution extractRelaxation() const override { return extractSolution(false); }
151 : :
152 : : ArithRatPairVec heuristicOptCoeffs() const override;
153 : :
154 : : MipResult solveMIP(bool al) override;
155 : 6 : Solution extractMIP() const override { return extractSolution(true); }
156 : : void setOptCoeffs(const ArithRatPairVec& ref) override;
157 : : std::vector<const CutInfo*> getValidCuts(const NodeLog& nodes) override;
158 : : ArithVar getBranchVar(const NodeLog& con) const override;
159 : :
160 : : static void printGLPKStatus(int status, std::ostream& out);
161 : :
162 : : virtual void setPivotLimit(int pl) override;
163 : :
164 : : virtual void setBranchingDepth(int bd) override;
165 : :
166 : : virtual void setBranchOnVariableLimit(int bl) override;
167 : :
168 : : virtual std::optional<Rational> estimateWithCFE(double d) const override;
169 : : virtual std::optional<Rational> estimateWithCFE(
170 : : double d, const Integer& D) const override;
171 : :
172 : : private:
173 : : Solution extractSolution(bool mip) const;
174 : : int guessDir(ArithVar v) const;
175 : :
176 : : // get this stuff out of here
177 : : void tryCut(int nid, CutInfo& cut) override;
178 : :
179 : : ArithVar _getArithVar(int nid, int M, int ind) const;
180 : 0 : ArithVar getArithVarFromRow(int nid, int ind) const
181 : : {
182 [ - - ]: 0 : if (ind >= 0)
183 : : {
184 : 0 : const NodeLog& nl = d_log.getNode(nid);
185 : 0 : return nl.lookupRowId(ind);
186 : : }
187 : 0 : return ARITHVAR_SENTINEL;
188 : : }
189 : :
190 : : // virtual void mapRowId(int nid, int ind, ArithVar v){
191 : : // NodeLog& nl = d_log.getNode(nid);
192 : : // nl.mapRowId(ind, v);
193 : : // }
194 : : // virtual void applyRowsDeleted(int nid, const RowsDeleted& rd){
195 : : // NodeLog& nl = d_log.getNode(nid);
196 : : // nl.applyRowsDeleted(rd);
197 : : // }
198 : :
199 : 0 : ArithVar getArithVarFromStructural(int ind) const{
200 [ - - ]: 0 : if(ind >= 0){
201 : 0 : unsigned u = (unsigned) ind;
202 [ - - ]: 0 : if(d_colToArithVar.isKey(u)){
203 : 0 : return d_colToArithVar[u];
204 : : }
205 : : }
206 : 0 : return ARITHVAR_SENTINEL;
207 : : }
208 : :
209 : : /**
210 : : * Attempts to make the row vector vec on the pad.
211 : : * If this is not in the row span of the original tableau this
212 : : * raises the failure flag.
213 : : */
214 : : bool attemptConstructTableRow(int node, int M, const PrimitiveVec& vec);
215 : : bool guessCoefficientsConstructTableRow(int node, int M, const PrimitiveVec& vec);
216 : : bool guessCoefficientsConstructTableRow(int node, int M, const PrimitiveVec& vec, const Integer& D);
217 : : bool gaussianElimConstructTableRow(int node, int M, const PrimitiveVec& vec);
218 : :
219 : : /* This is a guess of a vector in the row span of the tableau.
220 : : * Attempt to cancel out all of the variables.
221 : : * returns true if this is constructable.
222 : : */
223 : : bool guessIsConstructable(const DenseMap<Rational>& guess) const;
224 : :
225 : : /**
226 : : * Loads a vector of statuses into a dense map over bounds.
227 : : * returns true on failure.
228 : : */
229 : : bool loadToBound(int node, int M, int len, int* inds, int* statuses,
230 : : DenseMap<ConstraintP>& toBound) const;
231 : :
232 : : /** checks the cut on the pad for whether it is sufficiently similar to cut. */
233 : : bool checkCutOnPad(int nid, const CutInfo& cut) const;
234 : :
235 : :
236 : : /** turns the pad into a node and creates an explanation. */
237 : : //std::pair<Node, Node> makeCutNodes(int nid, const CutInfo& cut) const;
238 : :
239 : : // true means failure!
240 : : // BRANCH CUTS
241 : : bool attemptBranchCut(int nid, const BranchCutInfo& br);
242 : :
243 : : // GOMORY CUTS
244 : : bool attemptGmi(int nid, const GmiInfo& gmi);
245 : : /** tries to turn the information on the pad into a cut. */
246 : : bool constructGmiCut();
247 : :
248 : : // MIR CUTS
249 : : bool attemptMir(int nid, const MirInfo& mir);
250 : : bool applyCMIRRule(int nid, const MirInfo& mir);
251 : : bool makeRangeForComplemented(int nid, const MirInfo& mir);
252 : : bool loadSlacksIntoPad(int nid, const MirInfo& mir);
253 : : bool loadVirtualBoundsIntoPad(int nid, const MirInfo& mir);
254 : : bool loadRowSumIntoAgg(int nid, int M, const PrimitiveVec& mir);
255 : : bool buildModifiedRow(int nid, const MirInfo& mir);
256 : : bool constructMixedKnapsack();
257 : : bool replaceSlacksOnCuts();
258 : : bool loadVB(int nid, int M, int j, int ri, bool wantUb, VirtualBound& tmp);
259 : :
260 : : double sumInfeasibilities(glp_prob* prob, bool mip) const;
261 : :
262 : : /** UTILITIES FOR DEALING WITH ESTIMATES */
263 : :
264 : : static constexpr double SMALL_FIXED_DELTA = .000000001;
265 : : static constexpr double TOLERENCE = 1 + .000000001;
266 : :
267 : : /** Returns true if two doubles are roughly equal based on TOLERENCE and
268 : : * SMALL_FIXED_DELTA.*/
269 : : static bool roughlyEqual(double a, double b);
270 : :
271 : : /**
272 : : * Converts a rational to a continued fraction expansion representation
273 : : * using a maximum number of expansions equal to depth as long as the
274 : : * expression is not roughlyEqual() to 0.
275 : : */
276 : : static std::vector<Integer> rationalToCfe(const Rational& q, int depth);
277 : :
278 : : /** Converts a continued fraction expansion representation to a rational. */
279 : : static Rational cfeToRational(const std::vector<Integer>& exp);
280 : :
281 : : /** Estimates a rational as a continued fraction expansion.*/
282 : : static Rational estimateWithCFE(const Rational& q, const Integer& K);
283 : :
284 : : private:
285 : : const ArithVariables& d_vars;
286 : : TreeLog& d_log;
287 : : ApproximateStatistics& d_stats;
288 : :
289 : : /* the maximum pivots allowed in a query. */
290 : : int d_pivotLimit;
291 : :
292 : : /* maximum branches allowed on a variable */
293 : : int d_branchLimit;
294 : :
295 : : /* maxmimum branching depth allowed.*/
296 : : int d_maxDepth;
297 : :
298 : : /* Default denominator for diophatine approximation, 2^{26} .*/
299 : : static constexpr uint64_t s_defaultMaxDenom = (1 << 26);
300 : :
301 : : glp_prob* d_inputProb; /* a copy of the input prob */
302 : : glp_prob* d_realProb; /* a copy of the real relaxation output */
303 : : glp_prob* d_mipProb; /* a copy of the integer prob */
304 : :
305 : : DenseMap<int> d_colIndices;
306 : : DenseMap<int> d_rowIndices;
307 : :
308 : : NodeLog::RowIdMap d_rootRowIds;
309 : : // DenseMap<ArithVar> d_rowToArithVar;
310 : : DenseMap<ArithVar> d_colToArithVar;
311 : :
312 : : bool d_solvedRelaxation;
313 : : bool d_solvedMIP;
314 : :
315 : : CutScratchPad d_pad;
316 : :
317 : : std::vector<Integer> d_denomGuesses;
318 : : };
319 : :
320 : 6 : void ApproxGLPK::setPivotLimit(int pl)
321 : : {
322 [ - + ][ - + ]: 6 : Assert(pl >= 0);
[ - - ]
323 : 6 : d_pivotLimit = pl;
324 : 6 : }
325 : :
326 : 6 : void ApproxGLPK::setBranchingDepth(int bd)
327 : : {
328 [ - + ][ - + ]: 6 : Assert(bd >= 0);
[ - - ]
329 : 6 : d_maxDepth = bd;
330 : 6 : }
331 : :
332 : 6 : void ApproxGLPK::setBranchOnVariableLimit(int bl)
333 : : {
334 [ - + ][ - + ]: 6 : Assert(bl >= 0);
[ - - ]
335 : 6 : d_branchLimit = bl;
336 : 6 : }
337 : :
338 : 42 : bool ApproxGLPK::roughlyEqual(double a, double b)
339 : : {
340 [ + + ]: 42 : if (a == 0)
341 : : {
342 [ + - ][ + - ]: 6 : return -SMALL_FIXED_DELTA <= b && b <= SMALL_FIXED_DELTA;
343 : : }
344 [ + + ]: 36 : else if (b == 0)
345 : : {
346 [ + - ][ - + ]: 4 : return -SMALL_FIXED_DELTA <= a && a <= SMALL_FIXED_DELTA;
347 : : }
348 : : else
349 : : {
350 [ + - ][ + + ]: 32 : return std::abs(b / a) <= TOLERENCE && std::abs(a / b) <= TOLERENCE;
351 : : }
352 : : }
353 : :
354 : 0 : Rational ApproxGLPK::cfeToRational(const std::vector<Integer>& exp)
355 : : {
356 [ - - ]: 0 : if (exp.empty())
357 : : {
358 : 0 : return Rational(0);
359 : : }
360 : : else
361 : : {
362 : 0 : Rational result = exp.back();
363 : 0 : std::vector<Integer>::const_reverse_iterator exp_iter = exp.rbegin();
364 : 0 : std::vector<Integer>::const_reverse_iterator exp_end = exp.rend();
365 : 0 : ++exp_iter;
366 [ - - ]: 0 : while (exp_iter != exp_end)
367 : : {
368 : 0 : result = result.inverse();
369 : 0 : const Integer& i = *exp_iter;
370 : 0 : result += i;
371 : 0 : ++exp_iter;
372 : : }
373 : 0 : return result;
374 : : }
375 : : }
376 : :
377 : 0 : std::vector<Integer> ApproxGLPK::rationalToCfe(const Rational& q, int depth)
378 : : {
379 : 0 : std::vector<Integer> mods;
380 [ - - ]: 0 : if (!q.isZero())
381 : : {
382 : 0 : Rational carry = q;
383 [ - - ]: 0 : for (int i = 0; i <= depth; ++i)
384 : : {
385 : 0 : Assert(!carry.isZero());
386 : 0 : mods.push_back(Integer());
387 : 0 : Integer& back = mods.back();
388 : 0 : back = carry.floor();
389 [ - - ]: 0 : Trace("rationalToCfe") << " cfe[" << i << "]: " << back << std::endl;
390 : 0 : carry -= back;
391 [ - - ]: 0 : if (carry.isZero())
392 : : {
393 : 0 : break;
394 : : }
395 [ - - ]: 0 : else if (ApproxGLPK::roughlyEqual(carry.getDouble(), 0.0))
396 : : {
397 : 0 : break;
398 : : }
399 : : else
400 : : {
401 : 0 : carry = carry.inverse();
402 : : }
403 : : }
404 : : }
405 : 0 : return mods;
406 : : }
407 : :
408 : 16 : Rational ApproxGLPK::estimateWithCFE(const Rational& r, const Integer& K)
409 : : {
410 [ + - ]: 32 : Trace("estimateWithCFE") << "estimateWithCFE(" << r << ", " << K << ")"
411 : 16 : << std::endl;
412 : : // references
413 : : // page 4:
414 : : // http://carlossicoli.free.fr/C/Cassels_J.W.S.-An_introduction_to_diophantine_approximation-University_Press(1965).pdf
415 : : // http://en.wikipedia.org/wiki/Continued_fraction
416 [ - + ][ - + ]: 16 : Assert(K >= Integer(1));
[ - - ]
417 [ + - ]: 16 : if (r.getDenominator() <= K)
418 : : {
419 : 16 : return r;
420 : : }
421 : :
422 : : // current numerator and denominator that has not been resolved in the cfe
423 : 0 : Integer num = r.getNumerator(), den = r.getDenominator();
424 : 0 : Integer quot, rem;
425 : :
426 : 0 : unsigned t = 0;
427 : : // For a sequence of candidate solutions q_t/p_t
428 : : // we keep only 3 time steps: 0[prev], 1[current], 2[next]
429 : : // timesteps with a fake timestep 0 (p is 0 and q is 1)
430 : : // at timestep 1
431 : 0 : Integer p[3]; // h
432 : 0 : Integer q[3]; // k
433 : : // load the first 3 time steps manually
434 : 0 : p[0] = 0;
435 : 0 : q[0] = 1; // timestep -2
436 : 0 : p[1] = 1;
437 : 0 : q[1] = 0; // timestep -1
438 : :
439 : 0 : Integer::floorQR(quot, rem, num, den);
440 : 0 : num = den;
441 : 0 : den = rem;
442 : :
443 : 0 : q[2] = q[0] + quot * q[1];
444 : 0 : p[2] = p[0] + quot * p[1];
445 [ - - ]: 0 : Trace("estimateWithCFE") << " cfe[" << t << "]: " << p[2] << "/" << q[2]
446 : 0 : << std::endl;
447 [ - - ]: 0 : while (q[2] <= K)
448 : : {
449 : 0 : p[0] = p[1];
450 : 0 : p[1] = p[2];
451 : 0 : q[0] = q[1];
452 : 0 : q[1] = q[2];
453 : :
454 : 0 : Integer::floorQR(quot, rem, num, den);
455 : 0 : num = den;
456 : 0 : den = rem;
457 : :
458 : 0 : p[2] = p[0] + quot * p[1];
459 : 0 : q[2] = q[0] + quot * q[1];
460 : 0 : ++t;
461 [ - - ]: 0 : Trace("estimateWithCFE")
462 : 0 : << " cfe[" << t << "]: " << p[2] << "/" << q[2] << std::endl;
463 : : }
464 : :
465 : 0 : Integer k = (K - q[0]).floorDivideQuotient(q[1]);
466 : 0 : Rational cand_prev(p[0] + k * p[1], q[0] + k * q[1]);
467 : 0 : Rational cand_curr(p[1], q[1]);
468 : 0 : Rational dist_prev = (cand_prev - r).abs();
469 : 0 : Rational dist_curr = (cand_curr - r).abs();
470 [ - - ]: 0 : if (dist_prev <= dist_curr)
471 : : {
472 [ - - ]: 0 : Trace("estimateWithCFE")
473 : 0 : << cand_prev << " is closer than " << cand_curr << std::endl;
474 : 0 : return cand_prev;
475 : : }
476 : : else
477 : : {
478 [ - - ]: 0 : Trace("estimateWithCFE")
479 : 0 : << cand_curr << " is closer than " << cand_prev << std::endl;
480 : 0 : return cand_curr;
481 : : }
482 : : }
483 : :
484 : 16 : std::optional<Rational> ApproxGLPK::estimateWithCFE(double d,
485 : : const Integer& D) const
486 : : {
487 [ + - ]: 16 : if (std::optional<Rational> from_double = Rational::fromDouble(d))
488 : : {
489 : 32 : return estimateWithCFE(*from_double, D);
490 : : }
491 : 0 : return std::optional<Rational>();
492 : : }
493 : :
494 : 16 : std::optional<Rational> ApproxGLPK::estimateWithCFE(double d) const
495 : : {
496 : 32 : return estimateWithCFE(d, Integer(s_defaultMaxDenom));
497 : : }
498 : :
499 : 0 : Kind glpk_type_to_kind(int glpk_cut_type)
500 : : {
501 [ - - ][ - - ]: 0 : switch (glpk_cut_type)
502 : : {
503 : 0 : case GLP_LO: return Kind::GEQ;
504 : 0 : case GLP_UP: return Kind::LEQ;
505 : 0 : case GLP_FX: return Kind::EQUAL;
506 : 0 : case GLP_DB:
507 : : case GLP_FR:
508 : 0 : default: return Kind::UNDEFINED_KIND;
509 : : }
510 : : }
511 : :
512 : : #ifdef CVC5_ASSERTIONS
513 : 0 : static CutInfoKlass fromGlpkClass(int klass){
514 [ - - ][ - ]: 0 : switch(klass){
515 : 0 : case GLP_RF_GMI: return GmiCutKlass;
516 : 0 : case GLP_RF_MIR: return MirCutKlass;
517 : 0 : case GLP_RF_COV:
518 : : case GLP_RF_CLQ:
519 : 0 : default: return UnknownKlass;
520 : : }
521 : : }
522 : : #endif
523 : :
524 : 6 : ApproxGLPK::ApproxGLPK(const ArithVariables& var,
525 : : TreeLog& l,
526 : 6 : ApproximateStatistics& s)
527 : : : d_vars(var),
528 : : d_log(l),
529 : : d_stats(s),
530 : 6 : d_pivotLimit(std::numeric_limits<int>::max()),
531 : 6 : d_branchLimit(std::numeric_limits<int>::max()),
532 : 6 : d_maxDepth(std::numeric_limits<int>::max()),
533 : : d_inputProb(nullptr),
534 : : d_realProb(nullptr),
535 : : d_mipProb(nullptr),
536 : : d_solvedRelaxation(false),
537 : 6 : d_solvedMIP(false)
538 : : {
539 : 6 : d_denomGuesses.push_back(Integer(1<<22));
540 : 6 : d_denomGuesses.push_back(Integer(ApproxGLPK::s_defaultMaxDenom));
541 : 6 : d_denomGuesses.push_back(Integer(1ul<<29));
542 : 6 : d_denomGuesses.push_back(Integer(1ul<<31));
543 : :
544 : 6 : d_inputProb = glp_create_prob();
545 : 6 : d_realProb = glp_create_prob();
546 : 6 : d_mipProb = glp_create_prob();
547 : 6 : glp_set_obj_dir(d_inputProb, GLP_MAX);
548 : 6 : glp_set_prob_name(d_inputProb, "ApproxGLPK::approximateFindModel");
549 : :
550 : 6 : int numRows = 0;
551 : 6 : int numCols = 0;
552 : :
553 : : // Assign each variable to a row and column variable as it appears in the input
554 [ + + ]: 26 : for(ArithVariables::var_iterator vi = d_vars.var_begin(), vi_end = d_vars.var_end(); vi != vi_end; ++vi){
555 : 20 : ArithVar v = *vi;
556 : :
557 [ + + ]: 20 : if(d_vars.isAuxiliary(v)){
558 : 8 : ++numRows;
559 : 8 : d_rowIndices.set(v, numRows);
560 : : //mapRowId(d_log.getRootId(), numRows, v);
561 : 8 : d_rootRowIds.insert(std::make_pair(numRows, v));
562 : : //d_rowToArithVar.set(numRows, v);
563 [ + - ]: 8 : Trace("approx") << "Row vars: " << v << "<->" << numRows << std::endl;
564 : : }else{
565 : 12 : ++numCols;
566 : 12 : d_colIndices.set(v, numCols);
567 : 12 : d_colToArithVar.set(numCols, v);
568 [ + - ]: 12 : Trace("approx") << "Col vars: " << v << "<->" << numCols << std::endl;
569 : : }
570 : : }
571 [ - + ][ - + ]: 6 : Assert(numRows > 0);
[ - - ]
572 [ - + ][ - + ]: 6 : Assert(numCols > 0);
[ - - ]
573 : :
574 : 6 : glp_add_rows(d_inputProb, numRows);
575 : 6 : glp_add_cols(d_inputProb, numCols);
576 : :
577 : : // Assign the upper/lower bounds and types to each variable
578 [ + + ]: 26 : for(ArithVariables::var_iterator vi = d_vars.var_begin(), vi_end = d_vars.var_end(); vi != vi_end; ++vi){
579 : 20 : ArithVar v = *vi;
580 : :
581 [ - + ]: 20 : if (TraceIsOn("approx-debug"))
582 : : {
583 [ - - ]: 0 : Trace("approx-debug") << v << " ";
584 [ - - ]: 0 : d_vars.printModel(v, Trace("approx-debug"));
585 : : }
586 : :
587 : : int type;
588 : 20 : double lb = 0.0;
589 : 20 : double ub = 0.0;
590 [ + + ][ + + ]: 20 : if(d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
[ + + ]
591 [ + - ]: 2 : if(d_vars.boundsAreEqual(v)){
592 : 2 : type = GLP_FX;
593 : : }else{
594 : 0 : type = GLP_DB;
595 : : }
596 : 2 : lb = d_vars.getLowerBound(v).approx(SMALL_FIXED_DELTA);
597 : 2 : ub = d_vars.getUpperBound(v).approx(SMALL_FIXED_DELTA);
598 [ + + ][ + - ]: 18 : }else if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
[ + + ]
599 : 2 : type = GLP_UP;
600 : 2 : ub = d_vars.getUpperBound(v).approx(SMALL_FIXED_DELTA);
601 [ + - ][ + + ]: 16 : }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
[ + + ]
602 : 6 : type = GLP_LO;
603 : 6 : lb = d_vars.getLowerBound(v).approx(SMALL_FIXED_DELTA);
604 : : }else{
605 : 10 : type = GLP_FR;
606 : : }
607 : :
608 [ + + ]: 20 : if(d_vars.isAuxiliary(v)){
609 : 8 : int rowIndex = d_rowIndices[v];
610 : 8 : glp_set_row_bnds(d_inputProb, rowIndex, type, lb, ub);
611 : : }else{
612 : 12 : int colIndex = d_colIndices[v];
613 : : // is input is correct here
614 [ + - ]: 12 : int kind = d_vars.isInteger(v) ? GLP_IV : GLP_CV;
615 : 12 : glp_set_col_kind(d_inputProb, colIndex, kind);
616 : 12 : glp_set_col_bnds(d_inputProb, colIndex, type, lb, ub);
617 : : }
618 : : }
619 : :
620 : : // Count the number of entries
621 : 6 : int numEntries = 0;
622 [ + + ]: 14 : for(DenseMap<int>::const_iterator i = d_rowIndices.begin(), i_end = d_rowIndices.end(); i != i_end; ++i){
623 : 8 : ArithVar v = *i;
624 : 16 : Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
625 [ + + ]: 24 : for (Polynomial::iterator j = p.begin(), end = p.end(); j != end; ++j)
626 : : {
627 : 16 : ++numEntries;
628 : : }
629 : : }
630 : :
631 [ + - ]: 6 : int* ia = new int[numEntries+1];
632 [ + - ]: 6 : int* ja = new int[numEntries+1];
633 [ + - ]: 6 : double* ar = new double[numEntries+1];
634 : :
635 : 6 : int entryCounter = 0;
636 [ + + ]: 14 : for(DenseMap<int>::const_iterator i = d_rowIndices.begin(), i_end = d_rowIndices.end(); i != i_end; ++i){
637 : 8 : ArithVar v = *i;
638 : 8 : int rowIndex = d_rowIndices[v];
639 : :
640 : 16 : Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
641 : :
642 [ + + ]: 24 : for (Polynomial::iterator j = p.begin(), end = p.end(); j != end; ++j)
643 : : {
644 : 32 : const Monomial& mono = *j;
645 : 16 : const Constant& constant = mono.getConstant();
646 : 16 : const VarList& variable = mono.getVarList();
647 : :
648 : 16 : Node n = variable.getNode();
649 : :
650 [ - + ][ - + ]: 16 : Assert(d_vars.hasArithVar(n));
[ - - ]
651 : 16 : ArithVar av = d_vars.asArithVar(n);
652 : 16 : int colIndex = d_colIndices[av];
653 : 16 : double coeff = constant.getValue().getDouble();
654 : :
655 : 16 : ++entryCounter;
656 : 16 : ia[entryCounter] = rowIndex;
657 : 16 : ja[entryCounter] = colIndex;
658 : 16 : ar[entryCounter] = coeff;
659 : : }
660 : : }
661 : 6 : glp_load_matrix(d_inputProb, numEntries, ia, ja, ar);
662 : :
663 [ + - ]: 6 : delete[] ia;
664 [ + - ]: 6 : delete[] ja;
665 [ + - ]: 6 : delete[] ar;
666 : 6 : }
667 : :
668 : 0 : int ApproxGLPK::guessDir(ArithVar v) const{
669 [ - - ][ - - ]: 0 : if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
[ - - ]
670 : 0 : return -1;
671 [ - - ][ - - ]: 0 : }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
[ - - ]
672 : 0 : return 1;
673 [ - - ][ - - ]: 0 : }else if(!d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
[ - - ]
674 : 0 : return 0;
675 : : }else{
676 : 0 : int ubSgn = d_vars.getUpperBound(v).sgn();
677 : 0 : int lbSgn = d_vars.getLowerBound(v).sgn();
678 : :
679 [ - - ][ - - ]: 0 : if(ubSgn != 0 && lbSgn == 0){
680 : 0 : return -1;
681 [ - - ][ - - ]: 0 : }else if(ubSgn == 0 && lbSgn != 0){
682 : 0 : return 1;
683 : : }
684 : :
685 : 0 : return 1;
686 : : }
687 : : }
688 : :
689 : 6 : ArithRatPairVec ApproxGLPK::heuristicOptCoeffs() const{
690 : 6 : ArithRatPairVec ret;
691 : :
692 : : // Strategies are guess:
693 : : // 1 simple shared "ceiling" variable: danoint, pk1
694 : : // x1 >= c, x1 >= tmp1, x1 >= tmp2, ...
695 : : // 1 large row: fixnet, vpm2, pp08a
696 : : // (+ .......... ) <= c
697 : : // Not yet supported:
698 : : // 1 complex shared "ceiling" variable: opt1217
699 : : // x1 >= c, x1 >= (+ ..... ), x1 >= (+ ..... )
700 : : // and all of the ... are the same sign
701 : :
702 : :
703 : : // Candidates:
704 : : // 1) Upper and lower bounds are not equal.
705 : : // 2) The variable is not integer
706 : : // 3a) For columns look for a ceiling variable
707 : : // 3B) For rows look for a large row with
708 : :
709 : 12 : DenseMap<BoundCounts> d_colCandidates;
710 : 12 : DenseMap<uint32_t> d_rowCandidates;
711 : :
712 : 6 : double sumRowLength = 0.0;
713 : 6 : uint32_t maxRowLength = 0;
714 [ + + ]: 26 : for(ArithVariables::var_iterator vi = d_vars.var_begin(), vi_end = d_vars.var_end(); vi != vi_end; ++vi){
715 : 20 : ArithVar v = *vi;
716 : :
717 [ - + ]: 20 : if (TraceIsOn("approx-debug"))
718 : : {
719 [ - - ]: 0 : Trace("approx-debug") << v << " ";
720 [ - - ]: 0 : d_vars.printModel(v, Trace("approx-debug"));
721 : : }
722 : :
723 : : int type;
724 [ + + ][ + + ]: 20 : if(d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
[ + + ]
725 [ + - ]: 2 : if(d_vars.boundsAreEqual(v)){
726 : 2 : type = GLP_FX;
727 : : }else{
728 : 0 : type = GLP_DB;
729 : : }
730 [ + + ][ + - ]: 18 : }else if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
[ + + ]
731 : 2 : type = GLP_UP;
732 [ + - ][ + + ]: 16 : }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
[ + + ]
733 : 6 : type = GLP_LO;
734 : : }else{
735 : 10 : type = GLP_FR;
736 : : }
737 : :
738 [ + + ][ + + ]: 20 : if(type != GLP_FX && type != GLP_FR){
739 : :
740 [ + + ]: 8 : if(d_vars.isAuxiliary(v)){
741 : 6 : Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
742 : 6 : uint32_t len = p.size();
743 : 6 : d_rowCandidates.set(v, len);
744 : 6 : sumRowLength += len;
745 : 6 : maxRowLength = std::max(maxRowLength, len);
746 [ - + ]: 2 : }else if(!d_vars.isInteger(v)){
747 : 0 : d_colCandidates.set(v, BoundCounts());
748 : : }
749 : : }
750 : : }
751 : :
752 : 6 : uint32_t maxCount = 0;
753 [ + + ]: 14 : for(DenseMap<int>::const_iterator i = d_rowIndices.begin(), i_end = d_rowIndices.end(); i != i_end; ++i){
754 : 8 : ArithVar v = *i;
755 : :
756 [ + + ][ + + ]: 8 : bool lbCap = d_vars.hasLowerBound(v) && !d_vars.hasUpperBound(v);
757 [ + + ][ + - ]: 8 : bool ubCap = !d_vars.hasLowerBound(v) && d_vars.hasUpperBound(v);
758 : :
759 [ + + ][ + + ]: 8 : if(lbCap || ubCap){
760 [ + + ]: 6 : ConstraintP b = lbCap ? d_vars.getLowerBoundConstraint(v)
761 : 2 : : d_vars.getUpperBoundConstraint(v);
762 : :
763 [ + - ]: 6 : if(!(b->getValue()).noninfinitesimalIsZero()){ continue; }
764 : :
765 : 0 : Polynomial poly = Polynomial::parsePolynomial(d_vars.asNode(v));
766 [ - - ]: 0 : if(poly.size() != 2) { continue; }
767 : :
768 : 0 : Polynomial::iterator j = poly.begin();
769 : 0 : Monomial first = *j;
770 : 0 : ++j;
771 : 0 : Monomial second = *j;
772 : :
773 : 0 : bool firstIsPos = first.constantIsPositive();
774 : 0 : bool secondIsPos = second.constantIsPositive();
775 : :
776 [ - - ]: 0 : if(firstIsPos == secondIsPos){ continue; }
777 : :
778 [ - - ]: 0 : Monomial pos = firstIsPos == lbCap ? first : second;
779 [ - - ]: 0 : Monomial neg = firstIsPos != lbCap ? first : second;
780 : : // pos >= neg
781 : 0 : VarList p = pos.getVarList();
782 : 0 : VarList n = neg.getVarList();
783 [ - - ]: 0 : if(d_vars.hasArithVar(p.getNode())){
784 : 0 : ArithVar ap = d_vars.asArithVar(p.getNode());
785 [ - - ]: 0 : if( d_colCandidates.isKey(ap)){
786 : 0 : BoundCounts bc = d_colCandidates.get(ap);
787 : 0 : bc = BoundCounts(bc.lowerBoundCount(), bc.upperBoundCount()+1);
788 : 0 : maxCount = std::max(maxCount, bc.upperBoundCount());
789 : 0 : d_colCandidates.set(ap, bc);
790 : : }
791 : : }
792 [ - - ]: 0 : if(d_vars.hasArithVar(n.getNode())){
793 : 0 : ArithVar an = d_vars.asArithVar(n.getNode());
794 [ - - ]: 0 : if( d_colCandidates.isKey(an)){
795 : 0 : BoundCounts bc = d_colCandidates.get(an);
796 : 0 : bc = BoundCounts(bc.lowerBoundCount()+1, bc.upperBoundCount());
797 : 0 : maxCount = std::max(maxCount, bc.lowerBoundCount());
798 : 0 : d_colCandidates.set(an, bc);
799 : : }
800 : : }
801 : : }
802 : : }
803 : :
804 : : // Attempt row
805 [ + - ]: 6 : double avgRowLength = d_rowCandidates.size() >= 1 ?
806 : 6 : ( sumRowLength / d_rowCandidates.size() ) : 0.0;
807 : :
808 : : // There is a large row among the candidates
809 : 6 : bool guessARowCandidate = maxRowLength >= (10.0 * avgRowLength);
810 : :
811 : 6 : double rowLengthReq = (maxRowLength * .9);
812 : :
813 [ - + ]: 6 : if (guessARowCandidate)
814 : : {
815 [ - - ]: 0 : for (ArithVar r : d_rowCandidates)
816 : : {
817 : 0 : uint32_t len = d_rowCandidates[r];
818 : :
819 : 0 : int dir = guessDir(r);
820 [ - - ]: 0 : if(len >= rowLengthReq){
821 [ - - ]: 0 : if (TraceIsOn("approx-debug"))
822 : : {
823 [ - - ]: 0 : Trace("approx-debug") << "high row " << r << " " << len << " "
824 : 0 : << avgRowLength << " " << dir << std::endl;
825 [ - - ]: 0 : d_vars.printModel(r, Trace("approx-debug"));
826 : : }
827 : 0 : ret.push_back(ArithRatPair(r, Rational(dir)));
828 : : }
829 : : }
830 : : }
831 : :
832 : : // Attempt columns
833 : 6 : bool guessAColCandidate = maxCount >= 4;
834 [ - + ]: 6 : if (guessAColCandidate)
835 : : {
836 [ - - ]: 0 : for (ArithVar c : d_colCandidates)
837 : : {
838 : 0 : BoundCounts bc = d_colCandidates[c];
839 : :
840 : 0 : int dir = guessDir(c);
841 : 0 : double ubScore = double(bc.upperBoundCount()) / maxCount;
842 : 0 : double lbScore = double(bc.lowerBoundCount()) / maxCount;
843 [ - - ][ - - ]: 0 : if(ubScore >= .9 || lbScore >= .9){
844 [ - - ]: 0 : if (TraceIsOn("approx-debug"))
845 : : {
846 [ - - ]: 0 : Trace("approx-debug")
847 : 0 : << "high col " << c << " " << bc << " " << ubScore << " "
848 : 0 : << lbScore << " " << dir << std::endl;
849 [ - - ]: 0 : d_vars.printModel(c, Trace("approx-debug"));
850 : : }
851 : 0 : ret.push_back(ArithRatPair(c, Rational(c)));
852 : : }
853 : : }
854 : : }
855 : :
856 : 12 : return ret;
857 : : }
858 : :
859 : 0 : void ApproxGLPK::setOptCoeffs(const ArithRatPairVec& ref){
860 : 0 : DenseMap<double> nbCoeffs;
861 : :
862 [ - - ]: 0 : for(ArithRatPairVec::const_iterator i = ref.begin(), iend = ref.end(); i != iend; ++i){
863 : 0 : ArithVar v = (*i).first;
864 : 0 : const Rational& q = (*i).second;
865 : :
866 [ - - ]: 0 : if(d_vars.isAuxiliary(v)){
867 : : // replace the variable by its definition and multiply by q
868 : 0 : Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
869 : 0 : Polynomial pq = p * q;
870 : :
871 [ - - ]: 0 : for(Polynomial::iterator j = pq.begin(), jend = pq.end(); j != jend; ++j){
872 : 0 : const Monomial& mono = *j;
873 : 0 : const Constant& constant = mono.getConstant();
874 : 0 : const VarList& variable = mono.getVarList();
875 : :
876 : 0 : Node n = variable.getNode();
877 : :
878 : 0 : Assert(d_vars.hasArithVar(n));
879 : 0 : ArithVar av = d_vars.asArithVar(n);
880 : 0 : int colIndex = d_colIndices[av];
881 : 0 : double coeff = constant.getValue().getDouble();
882 [ - - ]: 0 : if(!nbCoeffs.isKey(colIndex)){
883 : 0 : nbCoeffs.set(colIndex, 0.0);
884 : : }
885 : 0 : nbCoeffs.set(colIndex, nbCoeffs[colIndex]+coeff);
886 : : }
887 : : }else{
888 : 0 : int colIndex = d_colIndices[v];
889 : 0 : double coeff = q.getDouble();
890 [ - - ]: 0 : if(!nbCoeffs.isKey(colIndex)){
891 : 0 : nbCoeffs.set(colIndex, 0.0);
892 : : }
893 : 0 : nbCoeffs.set(colIndex, nbCoeffs[colIndex]+coeff);
894 : : }
895 : : }
896 [ - - ]: 0 : for(DenseMap<double>::const_iterator ci =nbCoeffs.begin(), ciend = nbCoeffs.end(); ci != ciend; ++ci){
897 : 0 : Index colIndex = *ci;
898 : 0 : double coeff = nbCoeffs[colIndex];
899 : 0 : glp_set_obj_coef(d_inputProb, colIndex, coeff);
900 : : }
901 : 0 : }
902 : :
903 : :
904 : : /*
905 : : * rough strategy:
906 : : * real relaxation
907 : : * try approximate real optimization of error function
908 : : * pivot in its basis
909 : : * update to its assignment
910 : : * check with FCSimplex
911 : : * check integer solution
912 : : * try approximate mixed integer problem
913 : : * stop at the first feasible point
914 : : * pivot in its basis
915 : : * update to its assignment
916 : : * check with FCSimplex
917 : : */
918 : :
919 : 0 : void ApproxGLPK::printGLPKStatus(int status, std::ostream& out){
920 [ - - ][ - - ]: 0 : switch(status){
[ - - ][ - ]
921 : 0 : case GLP_OPT: out << "GLP_OPT" << std::endl; break;
922 : 0 : case GLP_FEAS: out << "GLP_FEAS" << std::endl; break;
923 : 0 : case GLP_INFEAS: out << "GLP_INFEAS" << std::endl; break;
924 : 0 : case GLP_NOFEAS: out << "GLP_NOFEAS" << std::endl; break;
925 : 0 : case GLP_UNBND: out << "GLP_UNBND" << std::endl; break;
926 : 0 : case GLP_UNDEF: out << "GLP_UNDEF" << std::endl; break;
927 : 0 : default: out << "Status unknown" << std::endl; break;
928 : : }
929 : 0 : }
930 : :
931 : 12 : ApproxGLPK::~ApproxGLPK(){
932 : 6 : glp_delete_prob(d_inputProb);
933 : 6 : glp_delete_prob(d_realProb);
934 : 6 : glp_delete_prob(d_mipProb);
935 : :
936 : 12 : }
937 : :
938 : 6 : ApproxGLPK::Solution ApproxGLPK::extractSolution(bool mip) const
939 : : {
940 [ - + ][ - + ]: 6 : Assert(d_solvedRelaxation);
[ - - ]
941 [ + - ][ - + ]: 6 : Assert(!mip || d_solvedMIP);
[ - + ][ - - ]
942 : :
943 : 6 : ApproxGLPK::Solution sol;
944 : 6 : DenseSet& newBasis = sol.newBasis;
945 : 6 : DenseMap<DeltaRational>& newValues = sol.newValues;
946 : :
947 [ + - ]: 6 : glp_prob* prob = mip ? d_mipProb : d_realProb;
948 : :
949 [ + + ]: 26 : for(ArithVariables::var_iterator i = d_vars.var_begin(), i_end = d_vars.var_end(); i != i_end; ++i){
950 : 20 : ArithVar vi = *i;
951 : 20 : bool isAux = d_vars.isAuxiliary(vi);
952 [ + + ]: 20 : int glpk_index = isAux ? d_rowIndices[vi] : d_colIndices[vi];
953 : :
954 [ + + ]: 20 : int status = isAux ? glp_get_row_stat(prob, glpk_index)
955 : 12 : : glp_get_col_stat(prob, glpk_index);
956 [ + - ]: 20 : Trace("approx-debug") << "assignment " << vi << std::endl;
957 : :
958 : 20 : bool useDefaultAssignment = false;
959 : :
960 [ + + ][ + + ]: 20 : switch(status){
961 : 8 : case GLP_BS:
962 [ + - ]: 8 : Trace("approx") << "basic" << std::endl;
963 : 8 : newBasis.add(vi);
964 : 8 : useDefaultAssignment = true;
965 : 8 : break;
966 : 8 : case GLP_NL:
967 : : case GLP_NS:
968 [ - + ]: 8 : if(!mip){
969 [ - - ]: 0 : Trace("approx-debug") << "non-basic lb" << std::endl;
970 : 0 : newValues.set(vi, d_vars.getLowerBound(vi));
971 : : }else{// intentionally fall through otherwise
972 : 8 : useDefaultAssignment = true;
973 : : }
974 : 8 : break;
975 : 2 : case GLP_NU:
976 [ - + ]: 2 : if(!mip){
977 [ - - ]: 0 : Trace("approx-debug") << "non-basic ub" << std::endl;
978 : 0 : newValues.set(vi, d_vars.getUpperBound(vi));
979 : : }else {// intentionally fall through otherwise
980 : 2 : useDefaultAssignment = true;
981 : : }
982 : 2 : break;
983 : 2 : default:
984 : : {
985 : 2 : useDefaultAssignment = true;
986 : : }
987 : 2 : break;
988 : : }
989 : :
990 [ + - ]: 20 : if(useDefaultAssignment){
991 [ + - ]: 20 : Trace("approx-debug") << "non-basic other" << std::endl;
992 : :
993 : : double newAssign;
994 [ + - ]: 20 : if(mip){
995 [ + + ]: 32 : newAssign = (isAux ? glp_mip_row_val(prob, glpk_index)
996 : 12 : : glp_mip_col_val(prob, glpk_index));
997 : : }else{
998 [ - - ]: 0 : newAssign = (isAux ? glp_get_row_prim(prob, glpk_index)
999 : 0 : : glp_get_col_prim(prob, glpk_index));
1000 : : }
1001 : 20 : const DeltaRational& oldAssign = d_vars.getAssignment(vi);
1002 : :
1003 : 20 : if (d_vars.hasLowerBound(vi)
1004 [ + + ][ + + ]: 28 : && roughlyEqual(newAssign,
[ + + ]
1005 : 8 : d_vars.getLowerBound(vi).approx(SMALL_FIXED_DELTA)))
1006 : : {
1007 [ + - ]: 2 : Trace("approx") << " to lb" << std::endl;
1008 : :
1009 : 2 : newValues.set(vi, d_vars.getLowerBound(vi));
1010 : : }
1011 : 18 : else if (d_vars.hasUpperBound(vi)
1012 [ + + ][ + - ]: 20 : && roughlyEqual(
[ + + ]
1013 : : newAssign,
1014 : 2 : d_vars.getUpperBound(vi).approx(SMALL_FIXED_DELTA)))
1015 : : {
1016 : 2 : newValues.set(vi, d_vars.getUpperBound(vi));
1017 [ + - ]: 2 : Trace("approx") << " to ub" << std::endl;
1018 : : }
1019 : : else
1020 : : {
1021 : 16 : double rounded = round(newAssign);
1022 [ + - ]: 16 : if (roughlyEqual(newAssign, rounded))
1023 : : {
1024 [ + - ]: 32 : Trace("approx") << "roughly equal " << rounded << " " << newAssign
1025 : 16 : << " " << oldAssign << std::endl;
1026 : 16 : newAssign = rounded;
1027 : : }
1028 : : else
1029 : : {
1030 [ - - ]: 0 : Trace("approx") << "not roughly equal " << rounded << " " << newAssign
1031 : 0 : << " " << oldAssign << std::endl;
1032 : : }
1033 : :
1034 : 32 : DeltaRational proposal;
1035 [ + - ]: 32 : if (std::optional<Rational> maybe_new = estimateWithCFE(newAssign))
1036 : : {
1037 : 16 : proposal = *maybe_new;
1038 : : }
1039 : : else
1040 : : {
1041 : : // failed to estimate the old value. defaulting to the current.
1042 : 0 : proposal = d_vars.getAssignment(vi);
1043 : : }
1044 : :
1045 [ + + ]: 16 : if (roughlyEqual(newAssign, oldAssign.approx(SMALL_FIXED_DELTA)))
1046 : : {
1047 [ + - ]: 4 : Trace("approx") << " to prev value" << newAssign << " " << oldAssign
1048 : 2 : << std::endl;
1049 : 2 : proposal = d_vars.getAssignment(vi);
1050 : : }
1051 : :
1052 [ - + ]: 16 : if (d_vars.strictlyLessThanLowerBound(vi, proposal))
1053 : : {
1054 [ - - ]: 0 : Trace("approx") << " round to lb " << d_vars.getLowerBound(vi)
1055 : 0 : << std::endl;
1056 : 0 : proposal = d_vars.getLowerBound(vi);
1057 : : }
1058 [ - + ]: 16 : else if (d_vars.strictlyGreaterThanUpperBound(vi, proposal))
1059 : : {
1060 [ - - ]: 0 : Trace("approx") << " round to ub " << d_vars.getUpperBound(vi)
1061 : 0 : << std::endl;
1062 : 0 : proposal = d_vars.getUpperBound(vi);
1063 : : }
1064 : : else
1065 : : {
1066 [ + - ]: 32 : Trace("approx") << " use proposal" << proposal << " " << oldAssign
1067 : 16 : << std::endl;
1068 : : }
1069 : 16 : newValues.set(vi, proposal);
1070 : : }
1071 : : }
1072 : : }
1073 : 6 : return sol;
1074 : : }
1075 : :
1076 : 0 : double ApproxGLPK::sumInfeasibilities(glp_prob* prob, bool mip) const{
1077 : : /* compute the sum of dual infeasibilities */
1078 : 0 : double infeas = 0.0;
1079 : :
1080 [ - - ]: 0 : for(ArithVariables::var_iterator i = d_vars.var_begin(), i_end = d_vars.var_end(); i != i_end; ++i){
1081 : 0 : ArithVar vi = *i;
1082 : 0 : bool isAux = d_vars.isAuxiliary(vi);
1083 [ - - ]: 0 : int glpk_index = isAux ? d_rowIndices[vi] : d_colIndices[vi];
1084 : :
1085 : : double newAssign;
1086 [ - - ]: 0 : if(mip){
1087 [ - - ]: 0 : newAssign = (isAux ? glp_mip_row_val(prob, glpk_index)
1088 : 0 : : glp_mip_col_val(prob, glpk_index));
1089 : : }else{
1090 [ - - ]: 0 : newAssign = (isAux ? glp_get_row_prim(prob, glpk_index)
1091 : 0 : : glp_get_col_prim(prob, glpk_index));
1092 : : }
1093 : :
1094 : :
1095 [ - - ]: 0 : double ub = isAux ?
1096 : 0 : glp_get_row_ub(prob, glpk_index) : glp_get_col_ub(prob, glpk_index);
1097 : :
1098 [ - - ]: 0 : double lb = isAux ?
1099 : 0 : glp_get_row_lb(prob, glpk_index) : glp_get_col_lb(prob, glpk_index);
1100 : :
1101 [ - - ]: 0 : if(ub != +DBL_MAX){
1102 [ - - ]: 0 : if(newAssign > ub){
1103 : 0 : double ubinf = newAssign - ub;
1104 : 0 : infeas += ubinf;
1105 [ - - ]: 0 : Trace("approx::soi")
1106 : 0 : << "ub inf" << vi << " " << ubinf << " " << infeas << std::endl;
1107 : : }
1108 : : }
1109 [ - - ]: 0 : if(lb != -DBL_MAX){
1110 [ - - ]: 0 : if(newAssign < lb){
1111 : 0 : double lbinf = lb - newAssign;
1112 : 0 : infeas += lbinf;
1113 : :
1114 [ - - ]: 0 : Trace("approx::soi")
1115 : 0 : << "lb inf" << vi << " " << lbinf << " " << infeas << std::endl;
1116 : : }
1117 : : }
1118 : : }
1119 : 0 : return infeas;
1120 : : }
1121 : :
1122 : 6 : LinResult ApproxGLPK::solveRelaxation(){
1123 [ - + ][ - + ]: 6 : Assert(!d_solvedRelaxation);
[ - - ]
1124 : :
1125 : : glp_smcp parm;
1126 : 6 : glp_init_smcp(&parm);
1127 : 6 : parm.presolve = GLP_OFF;
1128 : 6 : parm.meth = GLP_PRIMAL;
1129 : 6 : parm.pricing = GLP_PT_PSE;
1130 : 6 : parm.it_lim = d_pivotLimit;
1131 : 6 : parm.msg_lev = GLP_MSG_OFF;
1132 [ - + ]: 6 : if (TraceIsOn("approx-debug"))
1133 : : {
1134 : 0 : parm.msg_lev = GLP_MSG_ALL;
1135 : : }
1136 : :
1137 : 6 : glp_erase_prob(d_realProb);
1138 : 6 : glp_copy_prob(d_realProb, d_inputProb, GLP_OFF);
1139 : :
1140 : 6 : int res = glp_simplex(d_realProb, &parm);
1141 [ + - ]: 6 : switch(res){
1142 : 6 : case 0:
1143 : : {
1144 : 6 : int status = glp_get_status(d_realProb);
1145 : 6 : int iterationcount = glp_get_it_cnt(d_realProb);
1146 [ + - ][ - ]: 6 : switch(status){
1147 : 6 : case GLP_OPT:
1148 : : case GLP_FEAS:
1149 : : case GLP_UNBND:
1150 : 6 : d_solvedRelaxation = true;
1151 : 6 : return LinFeasible;
1152 : 0 : case GLP_INFEAS:
1153 : : case GLP_NOFEAS:
1154 : 0 : d_solvedRelaxation = true;
1155 : 0 : return LinInfeasible;
1156 : 0 : default:
1157 : : {
1158 [ - - ]: 0 : if(iterationcount >= d_pivotLimit){
1159 : 0 : return LinExhausted;
1160 : : }
1161 : 0 : return LinUnknown;
1162 : : }
1163 : : }
1164 : : }
1165 : 0 : default:
1166 : 0 : return LinUnknown;
1167 : : }
1168 : : }
1169 : :
1170 : :
1171 : : struct MirInfo : public CutInfo {
1172 : :
1173 : : /** a sum of input rows. */
1174 : : PrimitiveVec row_sum;
1175 : :
1176 : : /* the delta used */
1177 : : double delta;
1178 : :
1179 : : /* all of these are length vars == N+M*/
1180 : : int nvars;
1181 : : char* cset;
1182 : : char* subst;
1183 : : int* vlbRows;
1184 : : int* vubRows;
1185 : 0 : MirInfo(int execOrd, int ord)
1186 : 0 : : CutInfo(MirCutKlass, execOrd, ord)
1187 : : , nvars(0)
1188 : : , cset(NULL)
1189 : : , subst(NULL)
1190 : : , vlbRows(NULL)
1191 : 0 : , vubRows(NULL)
1192 : 0 : {}
1193 : :
1194 : 0 : ~MirInfo(){
1195 : 0 : clearSets();
1196 : 0 : }
1197 : 0 : void clearSets(){
1198 [ - - ]: 0 : if(cset != NULL){
1199 [ - - ]: 0 : delete[] cset;
1200 [ - - ]: 0 : delete[] subst;
1201 [ - - ]: 0 : delete[] vlbRows;
1202 [ - - ]: 0 : delete[] vubRows;
1203 : 0 : cset = NULL;
1204 : 0 : nvars = 0;
1205 : : }
1206 : 0 : }
1207 : 0 : void initSet(){
1208 : 0 : Assert(d_N >= 0);
1209 : 0 : Assert(d_mAtCreation >= 0);
1210 : 0 : clearSets();
1211 : :
1212 : 0 : int vars = 1 + d_N + d_mAtCreation;
1213 : :
1214 : 0 : cset = new char[1+vars];
1215 : 0 : subst = new char[1+vars];
1216 [ - - ]: 0 : vlbRows = new int[1+vars];
1217 [ - - ]: 0 : vubRows = new int[1+vars];
1218 : 0 : }
1219 : : };
1220 : :
1221 : : struct GmiInfo : public CutInfo {
1222 : : int basic;
1223 : : PrimitiveVec tab_row;
1224 : : int* tab_statuses;
1225 : : /* has the length tab_row.length */
1226 : :
1227 : 0 : GmiInfo(int execOrd, int ord)
1228 : 0 : : CutInfo(GmiCutKlass, execOrd, ord)
1229 : : , basic(-1)
1230 : : , tab_row()
1231 : 0 : , tab_statuses(NULL)
1232 : : {
1233 : 0 : Assert(!initialized_tab());
1234 : 0 : }
1235 : :
1236 : 0 : ~GmiInfo(){
1237 [ - - ]: 0 : if(initialized_tab()){
1238 : 0 : clear_tab();
1239 : : }
1240 : 0 : }
1241 : :
1242 : 0 : bool initialized_tab() const {
1243 : 0 : return tab_statuses != NULL;
1244 : : }
1245 : :
1246 : 0 : void init_tab(int N){
1247 [ - - ]: 0 : if(initialized_tab()){
1248 : 0 : clear_tab();
1249 : : }
1250 : 0 : tab_row.setup(N);
1251 [ - - ]: 0 : tab_statuses = new int[1+N];
1252 : 0 : }
1253 : :
1254 : 0 : void clear_tab() {
1255 [ - - ]: 0 : delete[] tab_statuses;
1256 : 0 : tab_statuses = NULL;
1257 : 0 : tab_row.clear();
1258 : 0 : basic = -1;
1259 : 0 : }
1260 : : };
1261 : :
1262 : :
1263 : :
1264 : 0 : static void loadCut(glp_tree *tree, CutInfo* cut){
1265 : : int ord, cut_len, cut_klass;
1266 : : int N, M;
1267 : : int* cut_inds;
1268 : : double* cut_coeffs;
1269 : : int glpk_cut_type;
1270 : : double cut_rhs;
1271 : : glp_prob* lp;
1272 : :
1273 : 0 : lp = glp_ios_get_prob(tree);
1274 : 0 : ord = cut->poolOrdinal();
1275 : :
1276 : 0 : N = glp_get_num_cols(lp);
1277 : 0 : M = glp_get_num_rows(lp);
1278 : :
1279 : 0 : cut->setDimensions(N, M);
1280 : :
1281 : :
1282 : :
1283 : : // Get the cut
1284 : 0 : cut_len = glp_ios_get_cut(tree, ord, NULL, NULL, &cut_klass, NULL, NULL);
1285 : 0 : Assert(fromGlpkClass(cut_klass) == cut->getKlass());
1286 : :
1287 : 0 : PrimitiveVec& cut_vec = cut->getCutVector();
1288 : 0 : cut_vec.setup(cut_len);
1289 : 0 : cut_inds = cut_vec.inds;
1290 : 0 : cut_coeffs = cut_vec.coeffs;
1291 : :
1292 : 0 : cut_vec.len = glp_ios_get_cut(tree, ord, cut_inds, cut_coeffs, &cut_klass, &glpk_cut_type, &cut_rhs);
1293 : 0 : Assert(fromGlpkClass(cut_klass) == cut->getKlass());
1294 : 0 : Assert(cut_vec.len == cut_len);
1295 : :
1296 : 0 : cut->setRhs(cut_rhs);
1297 : :
1298 : 0 : cut->setKind( glpk_type_to_kind(glpk_cut_type) );
1299 : 0 : }
1300 : :
1301 : :
1302 : 0 : static MirInfo* mirCut(glp_tree *tree, int exec_ord, int cut_ord){
1303 [ - - ]: 0 : Trace("approx::mirCut") << "mirCut()" << exec_ord << std::endl;
1304 : :
1305 : : MirInfo* mir;
1306 : 0 : mir = new MirInfo(exec_ord, cut_ord);
1307 : 0 : loadCut(tree, mir);
1308 : 0 : mir->initSet();
1309 : :
1310 : :
1311 : 0 : int nrows = glp_ios_cut_get_aux_nrows(tree, cut_ord);
1312 : :
1313 : 0 : PrimitiveVec& row_sum = mir->row_sum;
1314 : 0 : row_sum.setup(nrows);
1315 : 0 : glp_ios_cut_get_aux_rows(tree, cut_ord, row_sum.inds, row_sum.coeffs);
1316 : :
1317 : 0 : glp_ios_cut_get_mir_cset(tree, cut_ord, mir->cset);
1318 : 0 : mir->delta = glp_ios_cut_get_mir_delta(tree, cut_ord);
1319 : 0 : glp_ios_cut_get_mir_subst(tree, cut_ord, mir->subst);
1320 : 0 : glp_ios_cut_get_mir_virtual_rows(tree, cut_ord, mir->vlbRows, mir->vubRows);
1321 : :
1322 [ - - ]: 0 : if (TraceIsOn("approx::mirCut"))
1323 : : {
1324 [ - - ]: 0 : Trace("approx::mirCut") << "mir_id: " << exec_ord << std::endl;
1325 [ - - ]: 0 : row_sum.print(Trace("approx::mirCut"));
1326 : : }
1327 : :
1328 : 0 : return mir;
1329 : : }
1330 : :
1331 : 0 : static GmiInfo* gmiCut(glp_tree *tree, int exec_ord, int cut_ord){
1332 [ - - ]: 0 : Trace("approx::gmiCut") << "gmiCut()" << exec_ord << std::endl;
1333 : :
1334 : : int gmi_var;
1335 : : int write_pos;
1336 : : int read_pos;
1337 : : int stat;
1338 : : int ind;
1339 : : int i;
1340 : :
1341 : : GmiInfo* gmi;
1342 : : glp_prob* lp;
1343 : :
1344 : 0 : gmi = new GmiInfo(exec_ord, cut_ord);
1345 : 0 : loadCut(tree, gmi);
1346 : :
1347 : 0 : lp = glp_ios_get_prob(tree);
1348 : :
1349 : 0 : int N = gmi->getN();
1350 : 0 : int M = gmi->getMAtCreation();
1351 : :
1352 : : // Get the tableau row
1353 : 0 : int nrows CVC5_UNUSED = glp_ios_cut_get_aux_nrows(tree, gmi->poolOrdinal());
1354 : 0 : Assert(nrows == 1);
1355 : : int rows[1+1];
1356 : 0 : glp_ios_cut_get_aux_rows(tree, gmi->poolOrdinal(), rows, NULL);
1357 : 0 : gmi_var = rows[1];
1358 : :
1359 : 0 : gmi->init_tab(N);
1360 : 0 : gmi->basic = M+gmi_var;
1361 : :
1362 [ - - ]: 0 : Trace("approx::gmiCut") << gmi << " " << gmi->basic << " " << cut_ord << " "
1363 : 0 : << M << " " << gmi_var << std::endl;
1364 : :
1365 : 0 : PrimitiveVec& tab_row = gmi->tab_row;
1366 [ - - ]: 0 : Trace("approx::gmiCut") << "Is N sufficient here?" << std::endl;
1367 : 0 : tab_row.len = glp_eval_tab_row(lp, gmi->basic, tab_row.inds, tab_row.coeffs);
1368 : :
1369 [ - - ]: 0 : Trace("approx::gmiCut") << "gmi_var " << gmi_var << std::endl;
1370 : :
1371 [ - - ]: 0 : Trace("approx::gmiCut") << "tab_pos " << tab_row.len << std::endl;
1372 : 0 : write_pos = 1;
1373 [ - - ]: 0 : for(read_pos = 1; read_pos <= tab_row.len; ++read_pos){
1374 [ - - ]: 0 : if (fabs(tab_row.coeffs[read_pos]) < 1e-10){
1375 : : }else{
1376 : 0 : tab_row.coeffs[write_pos] = tab_row.coeffs[read_pos];
1377 : 0 : tab_row.inds[write_pos] = tab_row.inds[read_pos];
1378 : 0 : ++write_pos;
1379 : : }
1380 : : }
1381 : 0 : tab_row.len = write_pos-1;
1382 [ - - ]: 0 : Trace("approx::gmiCut") << "write_pos " << write_pos << std::endl;
1383 : 0 : Assert(tab_row.len > 0);
1384 : :
1385 [ - - ]: 0 : for(i = 1; i <= tab_row.len; ++i){
1386 : 0 : ind = tab_row.inds[i];
1387 [ - - ]: 0 : Trace("approx::gmiCut") << "ind " << i << " " << ind << std::endl;
1388 [ - - ]: 0 : stat = (ind <= M) ?
1389 : 0 : glp_get_row_stat(lp, ind) : glp_get_col_stat(lp, ind - M);
1390 : :
1391 [ - - ]: 0 : Trace("approx::gmiCut")
1392 : 0 : << "ind " << i << " " << ind << " stat " << stat << std::endl;
1393 [ - - ]: 0 : switch (stat){
1394 : 0 : case GLP_NL:
1395 : : case GLP_NU:
1396 : : case GLP_NS:
1397 : 0 : gmi->tab_statuses[i] = stat;
1398 : 0 : break;
1399 : 0 : case GLP_NF:
1400 : : default:
1401 : 0 : Unreachable();
1402 : : }
1403 : : }
1404 : :
1405 [ - - ]: 0 : if (TraceIsOn("approx::gmiCut"))
1406 : : {
1407 [ - - ]: 0 : gmi->print(Trace("approx::gmiCut"));
1408 : : }
1409 : 0 : return gmi;
1410 : : }
1411 : :
1412 : 0 : static BranchCutInfo* branchCut(glp_tree *tree, int exec_ord, int br_var, double br_val, bool down_bad){
1413 : : //(tree, br_var, br_val, dn < 0);
1414 : : double rhs;
1415 : : Kind k;
1416 [ - - ]: 0 : if(down_bad){
1417 : : // down branch is infeasible
1418 : : // x <= floor(v) is infeasible
1419 : : // - so x >= ceiling(v) is implied
1420 : 0 : k = Kind::GEQ;
1421 : 0 : rhs = std::ceil(br_val);
1422 : : }else{
1423 : : // up branch is infeasible
1424 : : // x >= ceiling(v) is infeasible
1425 : : // - so x <= floor(v) is implied
1426 : 0 : k = Kind::LEQ;
1427 : 0 : rhs = std::floor(br_val);
1428 : : }
1429 : 0 : BranchCutInfo* br_cut = new BranchCutInfo(exec_ord, br_var, k, rhs);
1430 : 0 : return br_cut;
1431 : : }
1432 : :
1433 : 50 : static void glpkCallback(glp_tree *tree, void *info){
1434 : 50 : AuxInfo* aux = (AuxInfo*)(info);
1435 : 50 : TreeLog& tl = *(aux->tl);
1436 : :
1437 : 50 : int exec = tl.getExecutionOrd();
1438 : 50 : int glpk_node_p = -1;
1439 : 50 : int node_ord = -1;
1440 : :
1441 [ - + ]: 50 : if(tl.isActivelyLogging()){
1442 : 0 : switch(glp_ios_reason(tree)){
1443 : 0 : case GLP_LI_DELROW:
1444 : : {
1445 : 0 : glpk_node_p = glp_ios_curr_node(tree);
1446 : 0 : node_ord = glp_ios_node_ord(tree, glpk_node_p);
1447 : :
1448 : 0 : int nrows = glp_ios_rows_deleted(tree, NULL);
1449 [ - - ]: 0 : int* num = new int[1+nrows];
1450 : 0 : glp_ios_rows_deleted(tree, num);
1451 : :
1452 : 0 : NodeLog& node = tl.getNode(node_ord);
1453 : :
1454 : 0 : RowsDeleted* rd = new RowsDeleted(exec, nrows, num);
1455 : :
1456 : 0 : node.addCut(rd);
1457 [ - - ]: 0 : delete[] num;
1458 : : }
1459 : 0 : break;
1460 : 0 : case GLP_ICUTADDED:
1461 : : {
1462 : 0 : int cut_ord = glp_ios_pool_size(tree);
1463 : 0 : glpk_node_p = glp_ios_curr_node(tree);
1464 : 0 : node_ord = glp_ios_node_ord(tree, glpk_node_p);
1465 : 0 : Assert(cut_ord > 0);
1466 [ - - ]: 0 : Trace("approx") << "curr node " << glpk_node_p << " cut ordinal "
1467 : 0 : << cut_ord << " node depth "
1468 : 0 : << glp_ios_node_level(tree, glpk_node_p) << std::endl;
1469 : : int klass;
1470 : 0 : glp_ios_get_cut(tree, cut_ord, NULL, NULL, &klass, NULL, NULL);
1471 : :
1472 : 0 : NodeLog& node = tl.getNode(node_ord);
1473 [ - - ][ - - ]: 0 : switch(klass){
[ - ]
1474 : 0 : case GLP_RF_GMI:
1475 : : {
1476 : 0 : GmiInfo* gmi = gmiCut(tree, exec, cut_ord);
1477 : 0 : node.addCut(gmi);
1478 : : }
1479 : 0 : break;
1480 : 0 : case GLP_RF_MIR:
1481 : : {
1482 : 0 : MirInfo* mir = mirCut(tree, exec, cut_ord);
1483 : 0 : node.addCut(mir);
1484 : : }
1485 : 0 : break;
1486 [ - - ]: 0 : case GLP_RF_COV: Trace("approx") << "GLP_RF_COV" << std::endl; break;
1487 [ - - ]: 0 : case GLP_RF_CLQ: Trace("approx") << "GLP_RF_CLQ" << std::endl; break;
1488 : 0 : default: break;
1489 : : }
1490 : : }
1491 : 0 : break;
1492 : 0 : case GLP_ICUTSELECT:
1493 : : {
1494 : 0 : glpk_node_p = glp_ios_curr_node(tree);
1495 : 0 : node_ord = glp_ios_node_ord(tree, glpk_node_p);
1496 : 0 : int cuts = glp_ios_pool_size(tree);
1497 [ - - ]: 0 : int* ords = new int[1+cuts];
1498 [ - - ]: 0 : int* rows = new int[1+cuts];
1499 : 0 : int N = glp_ios_selected_cuts(tree, ords, rows);
1500 : :
1501 : 0 : NodeLog& nl = tl.getNode(node_ord);
1502 [ - - ]: 0 : Trace("approx") << glpk_node_p << " " << node_ord << " " << cuts << " "
1503 : 0 : << N << std::endl;
1504 [ - - ]: 0 : for(int i = 1; i <= N; ++i){
1505 [ - - ]: 0 : Trace("approx") << "adding to " << node_ord << " @ i= " << i
1506 : 0 : << " ords[i] = " << ords[i]
1507 : 0 : << " rows[i] = " << rows[i] << std::endl;
1508 : 0 : nl.addSelected(ords[i], rows[i]);
1509 : : }
1510 [ - - ]: 0 : delete[] ords;
1511 [ - - ]: 0 : delete[] rows;
1512 : 0 : nl.applySelected();
1513 : : }
1514 : 0 : break;
1515 : 0 : case GLP_LI_BRANCH:
1516 : : {
1517 : : // a branch was just made
1518 : : int br_var;
1519 : : int p, dn, up;
1520 : : int p_ord, dn_ord, up_ord;
1521 : : double br_val;
1522 : 0 : br_var = glp_ios_branch_log(tree, &br_val, &p, &dn, &up);
1523 : 0 : p_ord = glp_ios_node_ord(tree, p);
1524 : :
1525 [ - - ]: 0 : dn_ord = (dn >= 0) ? glp_ios_node_ord(tree, dn) : -1;
1526 [ - - ]: 0 : up_ord = (up >= 0) ? glp_ios_node_ord(tree, up) : -1;
1527 : :
1528 [ - - ]: 0 : Trace("approx::") << "branch: " << br_var << " " << br_val << " tree "
1529 : 0 : << p << " " << dn << " " << up << std::endl;
1530 [ - - ]: 0 : Trace("approx::") << "\t " << p_ord << " " << dn_ord << " " << up_ord
1531 : 0 : << std::endl;
1532 [ - - ][ - - ]: 0 : if(dn < 0 && up < 0){
1533 [ - - ]: 0 : Trace("approx::") << "branch close " << exec << std::endl;
1534 : 0 : NodeLog& node = tl.getNode(p_ord);
1535 : 0 : BranchCutInfo* cut_br = branchCut(tree, exec, br_var, br_val, dn < 0);
1536 : 0 : node.addCut(cut_br);
1537 : 0 : tl.close(p_ord);
1538 [ - - ][ - - ]: 0 : }else if(dn < 0 || up < 0){
1539 [ - - ]: 0 : Trace("approx::") << "branch cut" << exec << std::endl;
1540 : 0 : NodeLog& node = tl.getNode(p_ord);
1541 : 0 : BranchCutInfo* cut_br = branchCut(tree, exec, br_var, br_val, dn < 0);
1542 : 0 : node.addCut(cut_br);
1543 : : }else{
1544 [ - - ]: 0 : Trace("approx::") << "normal branch" << std::endl;
1545 : 0 : tl.branch(p_ord, br_var, br_val, dn_ord, up_ord);
1546 : : }
1547 : : }
1548 : 0 : break;
1549 : 0 : case GLP_LI_CLOSE:
1550 : : {
1551 : 0 : glpk_node_p = glp_ios_curr_node(tree);
1552 : 0 : node_ord = glp_ios_node_ord(tree, glpk_node_p);
1553 [ - - ]: 0 : Trace("approx::") << "close " << glpk_node_p << std::endl;
1554 : 0 : tl.close(node_ord);
1555 : : }
1556 : 0 : break;
1557 : 0 : default:
1558 : 0 : break;
1559 : : }
1560 : : }
1561 : :
1562 [ + + ][ - - ]: 50 : switch(glp_ios_reason(tree)){
[ + ]
1563 : 4 : case GLP_IBINGO:
1564 [ + - ]: 4 : Trace("approx::") << "bingo" << std::endl;
1565 : 4 : aux->term = MipBingo;
1566 : 4 : glp_ios_terminate(tree);
1567 : 4 : break;
1568 : 6 : case GLP_ICUTADDED:
1569 : : {
1570 : 6 : tl.addCut();
1571 : : }
1572 : 6 : break;
1573 : 0 : case GLP_LI_BRANCH:
1574 : : {
1575 : : int p, dn, up;
1576 : 0 : int br_var = glp_ios_branch_log(tree, NULL, &p, &dn, &up);
1577 : :
1578 [ - - ]: 0 : if(br_var >= 0){
1579 : 0 : unsigned v = br_var;
1580 : 0 : tl.logBranch(v);
1581 : 0 : int depth = glp_ios_node_level(tree, p);
1582 : 0 : unsigned ubl = (aux->branchLimit) >= 0 ? ((unsigned)(aux->branchLimit)) : 0u;
1583 : 0 : if(tl.numBranches(v) >= ubl || depth >= (aux->branchDepth)){
1584 : 0 : aux->term = BranchesExhausted;
1585 : 0 : glp_ios_terminate(tree);
1586 : : }
1587 : : }
1588 : : }
1589 : 0 : break;
1590 : 0 : case GLP_LI_CLOSE:
1591 : 0 : break;
1592 : 40 : default:
1593 : : {
1594 : 40 : glp_prob* prob = glp_ios_get_prob(tree);
1595 : 40 : int iterationcount = glp_get_it_cnt(prob);
1596 [ - + ]: 40 : if(exec > (aux->pivotLimit)){
1597 : 0 : aux->term = ExecExhausted;
1598 : 0 : glp_ios_terminate(tree);
1599 [ - + ]: 40 : }else if(iterationcount > (aux->pivotLimit)){
1600 : 0 : aux->term = PivotsExhauasted;
1601 : 0 : glp_ios_terminate(tree);
1602 : : }
1603 : : }
1604 : 40 : break;
1605 : : }
1606 : 50 : }
1607 : :
1608 : 0 : std::vector<const CutInfo*> ApproxGLPK::getValidCuts(const NodeLog& con)
1609 : : {
1610 : 0 : std::vector<const CutInfo*> proven;
1611 : 0 : int nid = con.getNodeId();
1612 [ - - ]: 0 : for(NodeLog::const_iterator j = con.begin(), jend=con.end(); j!=jend; ++j){
1613 : 0 : CutInfo* cut = *j;
1614 : :
1615 [ - - ]: 0 : if(cut->getKlass() != RowsDeletedKlass){
1616 [ - - ]: 0 : if(!cut->reconstructed()){
1617 : 0 : Assert(!cut->reconstructed());
1618 : 0 : tryCut(nid, *cut);
1619 : : }
1620 : : }
1621 : :
1622 [ - - ]: 0 : if(cut->proven()){
1623 : 0 : proven.push_back(cut);
1624 : : }
1625 : : }
1626 : 0 : return proven;
1627 : : }
1628 : :
1629 : 0 : ArithVar ApproxGLPK::getBranchVar(const NodeLog& con) const{
1630 : 0 : int br_var = con.branchVariable();
1631 : 0 : return getArithVarFromStructural(br_var);
1632 : : }
1633 : :
1634 : :
1635 : 6 : MipResult ApproxGLPK::solveMIP(bool activelyLog){
1636 [ - + ][ - + ]: 6 : Assert(d_solvedRelaxation);
[ - - ]
1637 : : // Explicitly disable presolving
1638 : : // We need the basis thus the presolver must be off!
1639 : : // This is default, but this is just being cautious.
1640 : : AuxInfo aux;
1641 : 6 : aux.pivotLimit = d_pivotLimit;
1642 : 6 : aux.branchLimit = d_branchLimit;
1643 : 6 : aux.branchDepth = d_maxDepth;
1644 : 6 : aux.tl = &d_log;
1645 : 6 : aux.term = MipUnknown;
1646 : :
1647 : 6 : d_log.reset(d_rootRowIds);
1648 [ - + ]: 6 : if(activelyLog){
1649 : 0 : d_log.makeActive();
1650 : : }else{
1651 : 6 : d_log.makeInactive();
1652 : : }
1653 : :
1654 : : glp_iocp parm;
1655 : 6 : glp_init_iocp(&parm);
1656 : 6 : parm.presolve = GLP_OFF;
1657 : 6 : parm.pp_tech = GLP_PP_NONE;
1658 : 6 : parm.fp_heur = GLP_ON;
1659 : 6 : parm.gmi_cuts = GLP_ON;
1660 : 6 : parm.mir_cuts = GLP_ON;
1661 : 6 : parm.cov_cuts = GLP_ON;
1662 : 6 : parm.cb_func = glpkCallback;
1663 : 6 : parm.cb_info = &aux;
1664 : 6 : parm.msg_lev = GLP_MSG_OFF;
1665 [ - + ]: 6 : if (TraceIsOn("approx-debug"))
1666 : : {
1667 : 0 : parm.msg_lev = GLP_MSG_ALL;
1668 : : }
1669 : :
1670 : 6 : glp_erase_prob(d_mipProb);
1671 : 6 : glp_copy_prob(d_mipProb, d_realProb, GLP_OFF);
1672 : :
1673 : 6 : int res = glp_intopt(d_mipProb, &parm);
1674 : :
1675 [ + - ]: 12 : Trace("approx::solveMIP")
1676 : 6 : << "res " << res << " aux.term " << aux.term << std::endl;
1677 : :
1678 [ + - ]: 6 : switch(res){
1679 : 6 : case 0:
1680 : : case GLP_ESTOP:
1681 : : {
1682 : 6 : int status = glp_mip_status(d_mipProb);
1683 [ + - ]: 6 : Trace("approx::") << "status " << status << std::endl;
1684 [ + - ][ - ]: 6 : switch(status){
1685 : 6 : case GLP_OPT:
1686 : : case GLP_FEAS:
1687 : 6 : d_solvedMIP = true;
1688 [ + - ]: 6 : Trace("approx::") << "bingo here!" << std::endl;
1689 : 6 : return MipBingo;
1690 : 0 : case GLP_NOFEAS:
1691 : 0 : d_solvedMIP = true;
1692 : 0 : return MipClosed;
1693 : 0 : default:
1694 [ - - ]: 0 : if(aux.term == MipBingo){
1695 : 0 : d_solvedMIP = true;
1696 [ - - ]: 0 : Trace("approx::") << "bingo here?" << std::endl;
1697 : : }
1698 : 0 : return aux.term;
1699 : : }
1700 : : }
1701 : 0 : default:
1702 : 0 : return MipUnknown;
1703 : : }
1704 : : }
1705 : :
1706 : 0 : DeltaRational sumConstraints(const DenseMap<Rational>& xs, const DenseMap<ConstraintP>& cs, bool* anyinf){
1707 [ - - ]: 0 : if(anyinf != NULL){
1708 : 0 : *anyinf = false;
1709 : : }
1710 : :
1711 : 0 : DeltaRational beta(0);
1712 : 0 : DenseMap<Rational>::const_iterator iter, end;
1713 : 0 : iter = xs.begin();
1714 : 0 : end = xs.end();
1715 : :
1716 [ - - ]: 0 : Trace("approx::sumConstraints") << "sumConstraints";
1717 [ - - ]: 0 : for(; iter != end; ++iter){
1718 : 0 : ArithVar x = *iter;
1719 : 0 : const Rational& psi = xs[x];
1720 : 0 : ConstraintP c = cs[x];
1721 : 0 : Assert(c != NullConstraint);
1722 : :
1723 : 0 : const DeltaRational& bound = c->getValue();
1724 : 0 : beta += bound * psi;
1725 [ - - ]: 0 : Trace("approx::sumConstraints") << " +(" << bound << "*" << psi << ")";
1726 [ - - ]: 0 : if(anyinf != NULL ){
1727 : 0 : *anyinf = *anyinf || !bound.infinitesimalIsZero();
1728 : : }
1729 : : }
1730 [ - - ]: 0 : Trace("approx::sumConstraints") << "= " << beta << std::endl;
1731 : :
1732 : 0 : return beta;
1733 : : }
1734 : :
1735 : : // remove fixed variables from the vector
1736 : 0 : void removeFixed(const ArithVariables& vars,
1737 : : DenseVector& dv,
1738 : : std::set<ConstraintP>& exp)
1739 : : {
1740 : 0 : DenseMap<Rational>& vec = dv.lhs;
1741 : 0 : Rational& removed = dv.rhs;
1742 : 0 : std::vector<ArithVar> equal;
1743 : 0 : DenseMap<Rational>::const_iterator vec_iter, vec_end;
1744 : 0 : vec_iter = vec.begin(), vec_end = vec.end();
1745 [ - - ]: 0 : for(; vec_iter != vec_end; ++vec_iter){
1746 : 0 : ArithVar x = *vec_iter;
1747 [ - - ]: 0 : if(vars.boundsAreEqual(x)){
1748 : 0 : equal.push_back(x);
1749 : : }
1750 : : }
1751 : 0 : std::vector<ArithVar>::const_iterator equal_iter, equal_end;
1752 : 0 : equal_iter = equal.begin(), equal_end = equal.end();
1753 [ - - ]: 0 : for(; equal_iter != equal_end; ++equal_iter){
1754 : 0 : ArithVar x = *equal_iter;
1755 : 0 : Assert(vars.boundsAreEqual(x));
1756 : 0 : const DeltaRational& lb = vars.getLowerBound(x);
1757 : 0 : Assert(lb.infinitesimalIsZero());
1758 : 0 : removed -= (vec[x]) * lb.getNoninfinitesimalPart();
1759 : :
1760 : 0 : vec.remove(x);
1761 : :
1762 : 0 : std::pair<ConstraintP, ConstraintP> p = vars.explainEqualBounds(x);
1763 : 0 : exp.insert(p.first);
1764 [ - - ]: 0 : Trace("removeFixed") << "remove fixed " << p.first << std::endl;
1765 [ - - ]: 0 : if(p.second != NullConstraint){
1766 : 0 : exp.insert(p.second);
1767 [ - - ]: 0 : Trace("removeFixed") << "remove fixed " << p.second << std::endl;
1768 : : }
1769 : : }
1770 : 0 : }
1771 : 0 : void removeZeroes(DenseMap<Rational>& v){
1772 : : // Remove Slack variables
1773 : 0 : std::vector<ArithVar> zeroes;
1774 : 0 : DenseMap<Rational>::const_iterator i, iend;
1775 [ - - ]: 0 : for(i = v.begin(), iend = v.end(); i != iend; ++i){
1776 : 0 : ArithVar x = *i;
1777 [ - - ]: 0 : if(v[x].isZero()){
1778 : 0 : zeroes.push_back(x);
1779 : : }
1780 : : }
1781 : :
1782 : 0 : std::vector<ArithVar>::const_iterator j, jend;
1783 [ - - ]: 0 : for(j = zeroes.begin(), jend = zeroes.end(); j != jend; ++j){
1784 : 0 : ArithVar x = *j;
1785 : 0 : v.remove(x);
1786 : : }
1787 : 0 : }
1788 : :
1789 : 0 : void removeZeroes(DenseVector& v){
1790 : 0 : removeZeroes(v.lhs);
1791 : 0 : }
1792 : :
1793 : 0 : void removeAuxillaryVariables(const ArithVariables& vars, DenseMap<Rational>& vec){
1794 : : // Remove auxillary variables
1795 : 0 : std::vector<ArithVar> aux;
1796 : 0 : DenseMap<Rational>::const_iterator vec_iter, vec_end;
1797 : 0 : vec_iter = vec.begin(), vec_end = vec.end();
1798 [ - - ]: 0 : for(; vec_iter != vec_end; ++vec_iter){
1799 : 0 : ArithVar x = *vec_iter;
1800 [ - - ]: 0 : if(vars.isAuxiliary(x)){
1801 : 0 : aux.push_back(x);
1802 : : }
1803 : : }
1804 : :
1805 : 0 : std::vector<ArithVar>::const_iterator aux_iter, aux_end;
1806 : 0 : aux_iter = aux.begin(), aux_end = aux.end();
1807 [ - - ]: 0 : for(; aux_iter != aux_end; ++aux_iter){
1808 : 0 : ArithVar s = *aux_iter;
1809 : 0 : Rational& s_coeff = vec.get(s);
1810 : 0 : Assert(vars.isAuxiliary(s));
1811 : 0 : Assert(vars.hasNode(s));
1812 : 0 : Node sAsNode = vars.asNode(s);
1813 : 0 : Polynomial p = Polynomial::parsePolynomial(sAsNode);
1814 [ - - ]: 0 : for(Polynomial::iterator j = p.begin(), p_end=p.end(); j != p_end; ++j){
1815 : 0 : Monomial m = *j;
1816 : 0 : const Rational& ns_coeff = m.getConstant().getValue();
1817 : 0 : Node vl = m.getVarList().getNode();
1818 : 0 : ArithVar ns = vars.asArithVar(vl);
1819 : 0 : Rational prod = s_coeff * ns_coeff;
1820 [ - - ]: 0 : if(vec.isKey(ns)){
1821 : 0 : vec.get(ns) += prod;
1822 : : }else{
1823 : 0 : vec.set(ns, prod);
1824 : : }
1825 : : }
1826 : 0 : s_coeff = Rational(0); // subtract s_coeff * s from vec
1827 : : }
1828 : 0 : removeZeroes(vec);
1829 : 0 : }
1830 : :
1831 : 0 : ArithVar ApproxGLPK::_getArithVar(int nid, int M, int ind) const{
1832 [ - - ]: 0 : if(ind <= 0){
1833 : 0 : return ARITHVAR_SENTINEL;
1834 [ - - ]: 0 : }else if(ind <= M){
1835 : 0 : return getArithVarFromRow(nid, ind);
1836 : : }else{
1837 : 0 : return getArithVarFromStructural(ind - M);
1838 : : }
1839 : : }
1840 : :
1841 : :
1842 : 0 : bool ApproxGLPK::guessIsConstructable(const DenseMap<Rational>& guess) const {
1843 : : // basic variable
1844 : : // sum g[i] * x_i
1845 : 0 : DenseMap<Rational> g = guess;
1846 : 0 : removeAuxillaryVariables(d_vars, g);
1847 : :
1848 [ - - ]: 0 : if(TraceIsOn("guessIsConstructable")){
1849 [ - - ]: 0 : if(!g.empty()){
1850 [ - - ]: 0 : Trace("approx::guessIsConstructable")
1851 : 0 : << "guessIsConstructable failed " << g.size() << std::endl;
1852 [ - - ]: 0 : DenseVector::print(Trace("approx::guessIsConstructable"), g);
1853 [ - - ]: 0 : Trace("approx::guessIsConstructable") << std::endl;
1854 : : }
1855 : : }
1856 : 0 : return g.empty();
1857 : : }
1858 : :
1859 : 0 : bool ApproxGLPK::loadToBound(int nid, int M, int len, int* inds, int* statuses, DenseMap<ConstraintP>& toBound) const{
1860 [ - - ]: 0 : for(int i = 1; i <= len; ++i){
1861 : 0 : int status = statuses[i];
1862 : 0 : int ind = inds[i];
1863 : 0 : ArithVar v = _getArithVar(nid, M, ind);
1864 : 0 : ConstraintP c = NullConstraint;
1865 [ - - ]: 0 : if(v == ARITHVAR_SENTINEL){ return true; }
1866 : :
1867 [ - - ][ - ]: 0 : switch(status){
1868 : 0 : case GLP_NL:
1869 : 0 : c = d_vars.getLowerBoundConstraint(v);
1870 : 0 : break;
1871 : 0 : case GLP_NU:
1872 : : case GLP_NS: // upper bound sufficies for fixed variables
1873 : 0 : c = d_vars.getUpperBoundConstraint(v);
1874 : 0 : break;
1875 : 0 : case GLP_NF:
1876 : : default:
1877 : 0 : return true;
1878 : : }
1879 [ - - ]: 0 : if(c == NullConstraint){
1880 [ - - ]: 0 : Trace("approx::") << "couldn't find " << v << " @ " << nid << std::endl;
1881 : 0 : return true;
1882 : : }
1883 : 0 : Assert(c != NullConstraint);
1884 : 0 : toBound.set(v, c);
1885 : : }
1886 : 0 : return false;
1887 : : }
1888 : :
1889 : 0 : bool ApproxGLPK::checkCutOnPad(int nid, const CutInfo& cut) const{
1890 [ - - ]: 0 : Trace("approx::checkCutOnPad")
1891 : 0 : << "checkCutOnPad(" << nid << ", " << cut.getId() << ")" << std::endl;
1892 : :
1893 : 0 : const DenseMap<Rational>& constructedLhs = d_pad.d_cut.lhs;
1894 : 0 : const Rational& constructedRhs = d_pad.d_cut.rhs;
1895 : 0 : std::unordered_set<ArithVar> visited;
1896 : :
1897 [ - - ]: 0 : if(constructedLhs.empty()){
1898 [ - - ]: 0 : Trace("approx::checkCutOnPad") << "its empty?" << std::endl;
1899 : 0 : return true;
1900 : : }
1901 [ - - ]: 0 : if(cut.getKind() != d_pad.d_cutKind) {
1902 [ - - ]: 0 : Trace("approx::checkCutOnPad") << "rel doesn't match" << std::endl;
1903 : 0 : return true;
1904 : : }
1905 : :
1906 : 0 : const PrimitiveVec& cv = cut.getCutVector();
1907 [ - - ]: 0 : for(int i = 1; i <= cv.len; ++i){
1908 : 0 : int ind = cv.inds[i]; // this is always a structural variable
1909 : 0 : double coeff = cv.coeffs[i];
1910 : :
1911 : :
1912 : :
1913 [ - - ]: 0 : if(!d_colToArithVar.isKey(ind)){ return true; }
1914 : 0 : ArithVar x = d_colToArithVar[ind];
1915 : : //if(x == ARITHVAR_SENTINEL){ return true; }
1916 : 0 : visited.insert(x);
1917 : :
1918 : :
1919 [ - - ]: 0 : if(!constructedLhs.isKey(x)){
1920 [ - - ]: 0 : if (TraceIsOn("approx::checkCutOnPad"))
1921 : : {
1922 [ - - ]: 0 : Trace("approx::checkCutOnPad")
1923 : 0 : << " didn't find key for " << x << std::endl;
1924 [ - - ]: 0 : cut.print(Trace("approx::checkCutOnPad"));
1925 [ - - ]: 0 : Trace("approx::checkCutOnPad") << std::endl;
1926 [ - - ]: 0 : d_pad.d_cut.print(Trace("approx::checkCutOnPad"));
1927 [ - - ]: 0 : Trace("approx::checkCutOnPad") << std::endl;
1928 : : }
1929 : 0 : return true;
1930 : : }
1931 : :
1932 : 0 : const Rational& onConstructed = constructedLhs[x];
1933 : :
1934 [ - - ]: 0 : Trace("approx::checkCutOnPad") << ind << " " << coeff << " " << std::endl;
1935 [ - - ]: 0 : Trace("approx::checkCutOnPad")
1936 : 0 : << " " << x << " " << onConstructed << std::endl;
1937 : :
1938 [ - - ]: 0 : if(!roughlyEqual(coeff, onConstructed.getDouble())){
1939 [ - - ]: 0 : Trace("approx::checkCutOnPad") << "coeff failure" << std::endl;
1940 : 0 : return true;
1941 : : }
1942 : : }
1943 [ - - ]: 0 : if(visited.size() != constructedLhs.size()){
1944 [ - - ]: 0 : Trace("approx::checkCutOnPad") << "size mismatch" << std::endl;
1945 : 0 : return true;
1946 : : }
1947 : :
1948 : :
1949 [ - - ]: 0 : if(!roughlyEqual(cut.getRhs(), constructedRhs.getDouble())){
1950 [ - - ]: 0 : Trace("approx::checkCutOnPad") << "norm rhs is off " << cut.getRhs() << " "
1951 : 0 : << constructedRhs << std::endl;
1952 : 0 : return true;
1953 : : }
1954 : 0 : return false;
1955 : : }
1956 : :
1957 : :
1958 : :
1959 : 0 : bool ApproxGLPK::attemptBranchCut(int nid, const BranchCutInfo& br_cut){
1960 : 0 : d_pad.clear();
1961 : :
1962 : 0 : const PrimitiveVec& cut_vec = br_cut.getCutVector();
1963 : 0 : int structural = cut_vec.inds[1];
1964 : 0 : Assert(roughlyEqual(cut_vec.coeffs[1], +1.0));
1965 : :
1966 : 0 : ArithVar x = getArithVarFromStructural(structural);
1967 : 0 : d_pad.d_failure = (x == ARITHVAR_SENTINEL);
1968 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
1969 : :
1970 : 0 : Kind brKind = br_cut.getKind();
1971 : :
1972 [ - - ][ - - ]: 0 : d_pad.d_failure = (brKind != Kind::LEQ && brKind != Kind::GEQ);
1973 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
1974 : :
1975 : 0 : d_pad.d_cutKind = brKind;
1976 : 0 : d_pad.d_cut.lhs.set(x, Rational(1));
1977 : :
1978 : 0 : Rational& rhs = d_pad.d_cut.rhs;
1979 : 0 : std::optional<Rational> br_cut_rhs = Rational::fromDouble(br_cut.getRhs());
1980 [ - - ]: 0 : if (!br_cut_rhs)
1981 : : {
1982 : 0 : return true;
1983 : : }
1984 : :
1985 : 0 : rhs = estimateWithCFE(*br_cut_rhs, Integer(1));
1986 : 0 : d_pad.d_failure = !rhs.isIntegral();
1987 [ - - ]: 0 : if(d_pad.d_failure) { return true; }
1988 : :
1989 : 0 : d_pad.d_failure = checkCutOnPad(nid, br_cut);
1990 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
1991 : :
1992 : 0 : return false;
1993 : : }
1994 : :
1995 : 0 : bool ApproxGLPK::attemptGmi(int nid, const GmiInfo& gmi){
1996 : 0 : d_pad.clear();
1997 : :
1998 : 0 : d_pad.d_cutKind = Kind::GEQ;
1999 : :
2000 : 0 : int M = gmi.getMAtCreation();
2001 : 0 : ArithVar b = _getArithVar(nid, M, gmi.basic);
2002 : 0 : d_pad.d_failure = (b == ARITHVAR_SENTINEL);
2003 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2004 : :
2005 : 0 : d_pad.d_failure = !d_vars.isIntegerInput(b);
2006 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2007 : :
2008 : 0 : d_pad.d_basic = b;
2009 : :
2010 : :
2011 : 0 : const PrimitiveVec& tab = gmi.tab_row;
2012 : 0 : d_pad.d_failure = attemptConstructTableRow(nid, M, tab);
2013 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2014 : :
2015 : 0 : int* statuses = gmi.tab_statuses;
2016 : 0 : DenseMap<ConstraintP>& toBound = d_pad.d_toBound;
2017 : 0 : d_pad.d_failure = loadToBound(nid, M, tab.len, tab.inds, statuses, toBound);
2018 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2019 : :
2020 : 0 : d_pad.d_failure = constructGmiCut();
2021 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2022 : :
2023 : 0 : d_pad.d_failure = checkCutOnPad(nid, gmi);
2024 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2025 : :
2026 : 0 : return false;
2027 : : }
2028 : :
2029 : 0 : bool ApproxGLPK::applyCMIRRule(int nid, const MirInfo& mir){
2030 : :
2031 : 0 : const DenseMap<Rational>& compRanges = d_pad.d_compRanges;
2032 : :
2033 : 0 : DenseMap<Rational>& alpha = d_pad.d_alpha.lhs;
2034 : 0 : Rational& b = d_pad.d_alpha.rhs;
2035 : :
2036 : 0 : std::optional<Rational> delta = estimateWithCFE(mir.delta);
2037 [ - - ]: 0 : if (!delta)
2038 : : {
2039 : 0 : return true;
2040 : : }
2041 : :
2042 : 0 : d_pad.d_failure = (delta->sgn() <= 0);
2043 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2044 : :
2045 [ - - ]: 0 : Trace("approx::mir") << "applyCMIRRule() " << delta << " " << mir.delta
2046 : 0 : << std::endl;
2047 : :
2048 : 0 : DenseMap<Rational>::const_iterator iter, iend;
2049 : 0 : iter = alpha.begin(), iend = alpha.end();
2050 [ - - ]: 0 : for(; iter != iend; ++iter){
2051 : 0 : ArithVar v = *iter;
2052 : 0 : const Rational& curr = alpha[v];
2053 : 0 : Rational next = curr / *delta;
2054 [ - - ]: 0 : if(compRanges.isKey(v)){
2055 : 0 : b -= curr * compRanges[v];
2056 : 0 : alpha.set(v, - next);
2057 : : }else{
2058 : 0 : alpha.set(v, next);
2059 : : }
2060 : : }
2061 : 0 : b = b / *delta;
2062 : :
2063 : 0 : Rational roundB = (b + Rational(1,2)).floor();
2064 : 0 : d_pad.d_failure = (b - roundB).abs() < Rational(1,90);
2065 : : // intensionally more generous than glpk here
2066 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2067 : :
2068 : 0 : Rational one(1);
2069 : 0 : Rational fb = b.floor_frac();
2070 : 0 : Rational one_sub_fb = one - fb;
2071 : 0 : Rational gamma = (one / one_sub_fb);
2072 : :
2073 : 0 : DenseMap<Rational>& cut = d_pad.d_cut.lhs;
2074 : 0 : Rational& beta = d_pad.d_cut.rhs;
2075 : :
2076 : 0 : iter = alpha.begin(), iend = alpha.end();
2077 [ - - ]: 0 : for(; iter != iend; ++iter){
2078 : 0 : ArithVar v = *iter;
2079 : 0 : const Rational& a_j = alpha[v];
2080 [ - - ]: 0 : if(d_vars.isIntegerInput(v)){
2081 : 0 : Rational floor_aj = a_j.floor();
2082 : 0 : Rational frac_aj = a_j.floor_frac();
2083 [ - - ]: 0 : if(frac_aj <= fb){
2084 : 0 : cut.set(v, floor_aj);
2085 : : }else{
2086 : 0 : Rational tmp = ((frac_aj - fb) / one_sub_fb);
2087 : 0 : cut.set(v, floor_aj + tmp);
2088 : : }
2089 : : }else{
2090 : 0 : cut.set(v, a_j * gamma);
2091 : : }
2092 : : }
2093 : 0 : beta = b.floor();
2094 : :
2095 : 0 : iter = cut.begin(), iend = cut.end();
2096 [ - - ]: 0 : for(; iter != iend; ++iter){
2097 : 0 : ArithVar v = *iter;
2098 [ - - ]: 0 : if(compRanges.isKey(v)){
2099 : 0 : Rational neg = - cut[v];
2100 : 0 : beta += neg * compRanges[v];
2101 : 0 : cut.set(v, neg);
2102 : : }
2103 : : }
2104 : :
2105 : 0 : return false;
2106 : : }
2107 : :
2108 : 0 : bool ApproxGLPK::attemptMir(int nid, const MirInfo& mir){
2109 : 0 : d_pad.clear();
2110 : :
2111 : 0 : d_pad.d_cutKind = Kind::LEQ;
2112 : :
2113 : : // virtual bounds must be done before slacks
2114 : 0 : d_pad.d_failure = loadVirtualBoundsIntoPad(nid, mir);
2115 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2116 : :
2117 : 0 : d_pad.d_failure = loadSlacksIntoPad(nid, mir);
2118 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2119 : :
2120 : :
2121 : 0 : d_pad.d_failure = loadRowSumIntoAgg(nid, mir.getMAtCreation(), mir.row_sum);
2122 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2123 : :
2124 : 0 : removeFixed(d_vars, d_pad.d_agg, d_pad.d_explanation);
2125 : :
2126 : 0 : d_pad.d_failure = buildModifiedRow(nid, mir);
2127 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2128 : :
2129 : 0 : d_pad.d_failure = constructMixedKnapsack();
2130 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2131 : :
2132 : 0 : d_pad.d_failure = makeRangeForComplemented(nid, mir);
2133 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2134 : :
2135 : 0 : d_pad.d_failure = applyCMIRRule(nid, mir);
2136 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2137 : :
2138 : 0 : d_pad.d_failure = replaceSlacksOnCuts();
2139 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2140 : :
2141 : 0 : removeAuxillaryVariables(d_vars, d_pad.d_cut.lhs);
2142 : :
2143 : 0 : d_pad.d_failure = checkCutOnPad(nid, mir);
2144 [ - - ]: 0 : if(d_pad.d_failure){ return true; }
2145 : :
2146 : 0 : return false;
2147 : : //return makeCutNodes(nid, mir);
2148 : : }
2149 : :
2150 : : /** Returns true on failure. */
2151 : 0 : bool ApproxGLPK::loadVB(int nid, int M, int j, int ri, bool wantUb, VirtualBound& tmp){
2152 [ - - ]: 0 : if(ri <= 0) { return true; }
2153 : :
2154 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() " << std::endl;
2155 : :
2156 : 0 : ArithVar rowVar = _getArithVar(nid, M, ri);
2157 : 0 : ArithVar contVar = _getArithVar(nid, M, j);
2158 [ - - ]: 0 : if(rowVar == ARITHVAR_SENTINEL){
2159 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() "
2160 : 0 : << " rowVar is ARITHVAR_SENTINEL " << rowVar
2161 : 0 : << std::endl;
2162 : 0 : return true;
2163 : : }
2164 [ - - ]: 0 : if(contVar == ARITHVAR_SENTINEL){
2165 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() "
2166 : 0 : << " contVar is ARITHVAR_SENTINEL " << contVar
2167 : 0 : << std::endl;
2168 : 0 : return true; }
2169 : :
2170 [ - - ]: 0 : if(!d_vars.isAuxiliary(rowVar)){
2171 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() "
2172 : 0 : << " rowVar is not auxilliary " << rowVar
2173 : 0 : << std::endl;
2174 : 0 : return true;
2175 : : }
2176 : : // is integer is correct here
2177 [ - - ]: 0 : if(d_vars.isInteger(contVar)){
2178 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() contVar is integer " << contVar
2179 : 0 : << std::endl;
2180 : 0 : return true;
2181 : : }
2182 : :
2183 : 0 : ConstraintP lb = d_vars.getLowerBoundConstraint(rowVar);
2184 : 0 : ConstraintP ub = d_vars.getUpperBoundConstraint(rowVar);
2185 : :
2186 [ - - ][ - - ]: 0 : if(lb != NullConstraint && ub != NullConstraint){
2187 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() lb and ub are both NULL " << lb << " "
2188 : 0 : << ub << std::endl;
2189 : 0 : return true;
2190 : : }
2191 : :
2192 [ - - ]: 0 : ConstraintP rcon = lb == NullConstraint ? ub : lb;
2193 [ - - ]: 0 : if(rcon == NullConstraint) {
2194 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() rcon is NULL " << rcon << std::endl;
2195 : 0 : return true;
2196 : : }
2197 : :
2198 [ - - ]: 0 : if(!rcon->getValue().isZero()){
2199 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() rcon value is not 0 " << rcon->getValue()
2200 : 0 : << std::endl;
2201 : 0 : return true;
2202 : : }
2203 : :
2204 [ - - ]: 0 : if(!d_vars.hasNode(rowVar)){
2205 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() does not have node " << rowVar
2206 : 0 : << std::endl;
2207 : 0 : return true;
2208 : : }
2209 : :
2210 : 0 : Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(rowVar));
2211 [ - - ]: 0 : if (p.size() != 2)
2212 : : {
2213 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() polynomial is not binary: "
2214 : 0 : << p.getNode() << std::endl;
2215 : 0 : return true;
2216 : : }
2217 : :
2218 : 0 : Monomial first = p.getHead(), second = p.getTail().getHead();
2219 : 0 : Rational c1 = first.getConstant().getValue();
2220 : 0 : Rational c2 = second.getConstant().getValue();
2221 : 0 : Node nx1 = first.getVarList().getNode();
2222 : 0 : Node nx2 = second.getVarList().getNode();
2223 : :
2224 [ - - ]: 0 : if(!d_vars.hasArithVar(nx1)) {
2225 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() does not have a variable for nx1: "
2226 : 0 : << nx1 << std::endl;
2227 : 0 : return true;
2228 : : }
2229 [ - - ]: 0 : if(!d_vars.hasArithVar(nx2)) {
2230 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() does not have a variable for nx2 " << nx2
2231 : 0 : << std::endl;
2232 : 0 : return true;
2233 : : }
2234 : 0 : ArithVar x1 = d_vars.asArithVar(nx1), x2 = d_vars.asArithVar(nx2);
2235 : :
2236 : 0 : Assert(x1 != x2);
2237 : 0 : Assert(!c1.isZero());
2238 : 0 : Assert(!c2.isZero());
2239 : :
2240 [ - - ]: 0 : Trace("glpk::loadVB") << " lb " << lb << " ub " << ub << " rcon " << rcon
2241 : 0 : << " x1 " << x1 << " x2 " << x2 << " c1 " << c1
2242 : 0 : << " c2 " << c2 << std::endl;
2243 : :
2244 [ - - ]: 0 : ArithVar iv = (x1 == contVar) ? x2 : x1;
2245 [ - - ]: 0 : Rational& cc = (x1 == contVar) ? c1 : c2;
2246 [ - - ]: 0 : Rational& ic = (x1 == contVar) ? c2 : c1;
2247 : :
2248 [ - - ]: 0 : Trace("glpk::loadVB") << " cv " << contVar << " cc " << cc << " iv " << iv
2249 : 0 : << " c2 " << ic << std::endl;
2250 : :
2251 [ - - ]: 0 : if(!d_vars.isIntegerInput(iv)){
2252 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() iv is not an integer input variable "
2253 : 0 : << iv << std::endl;
2254 : 0 : return true;
2255 : : }
2256 : : // cc * cv + ic * iv <= 0 or
2257 : : // cc * cv + ic * iv <= 0
2258 : :
2259 [ - - ]: 0 : if(rcon == ub){ // multiply by -1
2260 : 0 : cc = -cc; ic = - ic;
2261 : : }
2262 [ - - ]: 0 : Trace("glpk::loadVB") << " cv " << contVar << " cc " << cc << " iv " << iv
2263 : 0 : << " c2 " << ic << std::endl;
2264 : :
2265 : : // cc * cv + ic * iv >= 0
2266 : : // cc * cv >= -ic * iv
2267 : : // if cc < 0:
2268 : : // cv <= -ic/cc * iv
2269 : : // elif cc > 0:
2270 : : // cv >= -ic/cc * iv
2271 : 0 : Assert(!cc.isZero());
2272 : 0 : Rational d = -ic/cc;
2273 [ - - ]: 0 : Trace("glpk::loadVB") << d << " " << cc.sgn() << std::endl;
2274 : 0 : bool nowUb = cc.sgn() < 0;
2275 [ - - ]: 0 : if(wantUb != nowUb) {
2276 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() wantUb is not nowUb " << wantUb << " "
2277 : 0 : << nowUb << std::endl;
2278 : :
2279 : 0 : return true;
2280 : : }
2281 : :
2282 [ - - ]: 0 : Kind rel = wantUb ? Kind::LEQ : Kind::GEQ;
2283 : :
2284 : 0 : tmp = VirtualBound(contVar, rel, d, iv, rcon);
2285 [ - - ]: 0 : Trace("glpk::loadVB") << "loadVB() was successful" << std::endl;
2286 : 0 : return false;
2287 : : }
2288 : :
2289 : 0 : bool ApproxGLPK::loadVirtualBoundsIntoPad(int nid, const MirInfo& mir){
2290 : 0 : Assert(mir.vlbRows != NULL);
2291 : 0 : Assert(mir.vubRows != NULL);
2292 : :
2293 : 0 : int N = mir.getN();
2294 : 0 : int M = mir.getMAtCreation();
2295 : :
2296 : : // Load the virtual bounds first
2297 : 0 : VirtualBound tmp;
2298 [ - - ]: 0 : for(int j=1; j <= N+M; ++j){
2299 [ - - ]: 0 : if(!loadVB(nid, M, j, mir.vlbRows[j], false, tmp)){
2300 [ - - ]: 0 : if(d_pad.d_vlb.isKey(tmp.x)){ return true; }
2301 : 0 : d_pad.d_vlb.set(tmp.x, tmp);
2302 [ - - ]: 0 : }else if(mir.vlbRows[j] > 0){
2303 [ - - ]: 0 : Trace("approx::mir") << "expected vlb to work" << std::endl;
2304 : : }
2305 [ - - ]: 0 : if(!loadVB(nid, M, j, mir.vubRows[j], true, tmp)){
2306 [ - - ]: 0 : if(d_pad.d_vub.isKey(tmp.x)){ return true; }
2307 : 0 : d_pad.d_vub.set(tmp.x, tmp);
2308 [ - - ]: 0 : }else if(mir.vubRows[j] > 0){
2309 [ - - ]: 0 : Trace("approx::mir") << "expected vub to work" << std::endl;
2310 : : }
2311 : : }
2312 : 0 : return false;
2313 : : }
2314 : :
2315 : 0 : bool ApproxGLPK::loadSlacksIntoPad(int nid, const MirInfo& mir){
2316 : 0 : Assert(mir.vlbRows != NULL);
2317 : 0 : Assert(mir.vubRows != NULL);
2318 : :
2319 : 0 : int N = mir.getN();
2320 : 0 : int M = mir.getMAtCreation();
2321 : :
2322 : : bool useVB;
2323 : : // Load the virtual bounds first
2324 : : SlackReplace rep;
2325 : : bool lb;
2326 : : ConstraintP b;
2327 [ - - ]: 0 : Trace("approx::mir") << "loadSlacksIntoPad(): N=" << N << ", M=" << M
2328 : 0 : << std::endl;
2329 [ - - ]: 0 : for(int j=1; j <= N+M; ++j){
2330 : 0 : ArithVar v = _getArithVar(nid, M, j);
2331 [ - - ]: 0 : if(v == ARITHVAR_SENTINEL){
2332 [ - - ]: 0 : Trace("approx::mir") << " for: " << j << " no variable" << std::endl;
2333 : 0 : continue;
2334 : : }
2335 : 0 : rep = SlackUndef;
2336 : 0 : char sub = mir.subst[j];
2337 [ - - ][ - ]: 0 : switch(sub){
2338 : 0 : case 'L':
2339 : : case 'U':
2340 : 0 : lb = (sub == 'L');
2341 [ - - ]: 0 : useVB = lb ? (mir.vlbRows[j] > 0) : (mir.vubRows[j] > 0);
2342 [ - - ]: 0 : if(useVB){
2343 : 0 : if(lb ? d_pad.d_vlb.isKey(v) : d_pad.d_vub.isKey(v)){
2344 [ - - ]: 0 : rep = lb ? SlackVLB : SlackVUB;
2345 : : }
2346 : : }else{
2347 [ - - ]: 0 : b = lb ? d_vars.getLowerBoundConstraint(v)
2348 : 0 : : d_vars.getUpperBoundConstraint(v);
2349 [ - - ]: 0 : if(b != NullConstraint){
2350 [ - - ]: 0 : if(b->getValue().infinitesimalIsZero()){
2351 [ - - ]: 0 : rep = lb ? SlackLB : SlackUB;
2352 : : }
2353 : : }
2354 : : }
2355 : :
2356 [ - - ]: 0 : Trace("approx::mir") << " for: " << j << ", " << v;
2357 : 0 : Trace("approx::mir") << " " << ((rep != SlackUndef) ? "succ" : "fail")
2358 : 0 : << " ";
2359 [ - - ]: 0 : Trace("approx::mir") << sub << " " << rep << " " << mir.vlbRows[j] << " "
2360 : 0 : << mir.vubRows[j] << std::endl;
2361 [ - - ]: 0 : if(rep != SlackUndef){
2362 : 0 : d_pad.d_slacks.set(v,rep);
2363 : : }
2364 : 0 : break;
2365 : 0 : case '?':
2366 : 0 : continue;
2367 : 0 : default:
2368 [ - - ]: 0 : Trace("approx::mir") << " for: " << j << " got subst " << (int)sub
2369 : 0 : << std::endl;
2370 : 0 : continue;
2371 : : }
2372 : : }
2373 : 0 : return false;
2374 : : }
2375 : :
2376 : 0 : bool ApproxGLPK::replaceSlacksOnCuts(){
2377 : 0 : std::vector<ArithVar> virtualVars;
2378 : :
2379 : 0 : DenseMap<Rational>& cut = d_pad.d_cut.lhs;
2380 : 0 : Rational& cutRhs = d_pad.d_cut.rhs;
2381 : :
2382 : 0 : DenseMap<Rational>::const_iterator iter, iend;
2383 : 0 : iter = cut.begin(), iend = cut.end();
2384 [ - - ]: 0 : for(; iter != iend; ++iter){
2385 : 0 : ArithVar x = *iter;
2386 : 0 : SlackReplace rep = d_pad.d_slacks[x];
2387 [ - - ]: 0 : if(d_vars.isIntegerInput(x)){
2388 : 0 : Assert(rep == SlackLB || rep == SlackUB);
2389 : 0 : Rational& a = cut.get(x);
2390 : :
2391 [ - - ]: 0 : const DeltaRational& bound = (rep == SlackLB) ?
2392 : 0 : d_vars.getLowerBound(x) : d_vars.getUpperBound(x);
2393 : 0 : Assert(bound.infinitesimalIsZero());
2394 : 0 : Rational prod = a * bound.getNoninfinitesimalPart();
2395 [ - - ]: 0 : if(rep == SlackLB){
2396 : 0 : cutRhs += prod;
2397 : : }else{
2398 : 0 : cutRhs -= prod;
2399 : 0 : a = -a;
2400 : : }
2401 [ - - ]: 0 : }else if(rep == SlackVLB){
2402 : 0 : virtualVars.push_back(d_pad.d_vlb[x].y);
2403 [ - - ]: 0 : }else if(rep == SlackVUB){
2404 : 0 : virtualVars.push_back(d_pad.d_vub[x].y);
2405 : : }
2406 : : }
2407 : :
2408 [ - - ]: 0 : for(size_t i = 0; i < virtualVars.size(); ++i){
2409 : 0 : ArithVar x = virtualVars[i];
2410 [ - - ]: 0 : if(!cut.isKey(x)){
2411 : 0 : cut.set(x, Rational(0));
2412 : : }
2413 : : }
2414 : :
2415 : 0 : iter = cut.begin(), iend = cut.end();
2416 [ - - ]: 0 : for(; iter != iend; ++iter){
2417 : 0 : ArithVar x = *iter;
2418 [ - - ]: 0 : if(!d_vars.isIntegerInput(x)){
2419 : 0 : SlackReplace rep = d_pad.d_slacks[x];
2420 : 0 : Rational& a = cut.get(x);
2421 [ - - ][ - - ]: 0 : switch(rep){
2422 : 0 : case SlackLB:
2423 : : {
2424 : 0 : const DeltaRational& bound = d_vars.getLowerBound(x);
2425 : 0 : Assert(bound.infinitesimalIsZero());
2426 : 0 : cutRhs += a * bound.getNoninfinitesimalPart();
2427 : : }
2428 : 0 : break;
2429 : 0 : case SlackUB:
2430 : : {
2431 : 0 : const DeltaRational& bound = d_vars.getUpperBound(x);
2432 : 0 : Assert(bound.infinitesimalIsZero());
2433 : 0 : cutRhs -= a * bound.getNoninfinitesimalPart();
2434 : 0 : a = -a;
2435 : : }
2436 : 0 : break;
2437 : 0 : case SlackVLB:
2438 : : case SlackVUB:
2439 : : {
2440 : 0 : bool lb = (rep == SlackVLB);
2441 [ - - ]: 0 : const VirtualBound& vb = lb ?
2442 : 0 : d_pad.d_vlb[x] : d_pad.d_vub[x];
2443 : 0 : ArithVar y = vb.y;
2444 : 0 : Assert(vb.x == x);
2445 : 0 : Assert(cut.isKey(y));
2446 : 0 : Rational& ycoeff = cut.get(y);
2447 [ - - ]: 0 : if(lb){
2448 : 0 : ycoeff -= a * vb.d;
2449 : : }else{
2450 : 0 : ycoeff += a * vb.d;
2451 : 0 : a = -a;
2452 : : }
2453 : : }
2454 : 0 : break;
2455 : 0 : default:
2456 : 0 : return true;
2457 : : }
2458 : : }
2459 : : }
2460 : 0 : removeZeroes(cut);
2461 : 0 : return false;
2462 : : }
2463 : :
2464 : 0 : bool ApproxGLPK::loadRowSumIntoAgg(int nid, int M, const PrimitiveVec& row_sum){
2465 : 0 : DenseMap<Rational>& lhs = d_pad.d_agg.lhs;
2466 : 0 : d_pad.d_agg.rhs = Rational(0);
2467 : :
2468 : 0 : int len = row_sum.len;
2469 [ - - ]: 0 : for(int i = 1; i <= len; ++i){
2470 : 0 : int aux_ind = row_sum.inds[i]; // auxillary index
2471 : 0 : double coeff = row_sum.coeffs[i];
2472 : 0 : ArithVar x = _getArithVar(nid, M, aux_ind);
2473 [ - - ]: 0 : if(x == ARITHVAR_SENTINEL){ return true; }
2474 : 0 : std::optional<Rational> c = estimateWithCFE(coeff);
2475 [ - - ]: 0 : if (!c)
2476 : : {
2477 : 0 : return true;
2478 : : }
2479 : :
2480 [ - - ]: 0 : if (lhs.isKey(x))
2481 : : {
2482 : 0 : lhs.get(x) -= *c;
2483 : : }
2484 : : else
2485 : : {
2486 : 0 : lhs.set(x, -(*c));
2487 : : }
2488 : : }
2489 : :
2490 [ - - ]: 0 : Trace("approx::mir") << "beg loadRowSumIntoAgg() 1" << std::endl;
2491 [ - - ]: 0 : if (TraceIsOn("approx::mir"))
2492 : : {
2493 [ - - ]: 0 : DenseVector::print(Trace("approx::mir"), lhs);
2494 : : }
2495 : 0 : removeAuxillaryVariables(d_vars, lhs);
2496 [ - - ]: 0 : Trace("approx::mir") << "end loadRowSumIntoAgg() 1" << std::endl;
2497 : :
2498 [ - - ]: 0 : if (TraceIsOn("approx::mir"))
2499 : : {
2500 [ - - ]: 0 : Trace("approx::mir") << "loadRowSumIntoAgg() 2" << std::endl;
2501 [ - - ]: 0 : DenseVector::print(Trace("approx::mir"), lhs);
2502 [ - - ]: 0 : Trace("approx::mir") << "end loadRowSumIntoAgg() 2" << std::endl;
2503 : : }
2504 : :
2505 [ - - ]: 0 : for(int i = 1; i <= len; ++i){
2506 : 0 : int aux_ind = row_sum.inds[i]; // auxillary index
2507 : 0 : double coeff = row_sum.coeffs[i];
2508 : 0 : ArithVar x = _getArithVar(nid, M, aux_ind);
2509 : 0 : Assert(x != ARITHVAR_SENTINEL);
2510 : 0 : std::optional<Rational> c = estimateWithCFE(coeff);
2511 [ - - ]: 0 : if (!c)
2512 : : {
2513 : 0 : return true;
2514 : : }
2515 : 0 : Assert(!lhs.isKey(x));
2516 : 0 : lhs.set(x, *c);
2517 : : }
2518 : :
2519 [ - - ]: 0 : if (TraceIsOn("approx::mir"))
2520 : : {
2521 [ - - ]: 0 : Trace("approx::mir") << "loadRowSumIntoAgg() 2" << std::endl;
2522 [ - - ]: 0 : DenseVector::print(Trace("approx::mir"), lhs);
2523 [ - - ]: 0 : Trace("approx::mir") << "end loadRowSumIntoAgg() 3" << std::endl;
2524 : : }
2525 : 0 : return false;
2526 : : }
2527 : :
2528 : 0 : bool ApproxGLPK::buildModifiedRow(int nid, const MirInfo& mir){
2529 : 0 : const DenseMap<Rational>& agg = d_pad.d_agg.lhs;
2530 : 0 : const Rational& aggRhs = d_pad.d_agg.rhs;
2531 : 0 : DenseMap<Rational>& mod = d_pad.d_mod.lhs;
2532 : 0 : Rational& modRhs = d_pad.d_mod.rhs;
2533 : :
2534 [ - - ]: 0 : Trace("approx::mir") << "buildModifiedRow()"
2535 : 0 : << " |agg|=" << d_pad.d_agg.lhs.size()
2536 : 0 : << " |mod|=" << d_pad.d_mod.lhs.size()
2537 : 0 : << " |slacks|=" << d_pad.d_slacks.size()
2538 : 0 : << " |vlb|=" << d_pad.d_vub.size()
2539 : 0 : << " |vub|=" << d_pad.d_vlb.size() << std::endl;
2540 : :
2541 : 0 : mod.addAll(agg);
2542 : 0 : modRhs = aggRhs;
2543 : 0 : DenseMap<Rational>::const_iterator iter, iend;
2544 [ - - ]: 0 : for(iter = agg.begin(), iend = agg.end(); iter != iend; ++iter){
2545 : 0 : ArithVar x = *iter;
2546 : 0 : const Rational& c = mod[x];
2547 [ - - ]: 0 : if(!d_pad.d_slacks.isKey(x)){
2548 [ - - ]: 0 : Trace("approx::mir") << "missed x: " << x << std::endl;
2549 : 0 : return true;
2550 : : }
2551 : 0 : SlackReplace rep = d_pad.d_slacks[x];
2552 [ - - ][ - ]: 0 : switch(rep){
2553 : 0 : case SlackLB: // skip for now
2554 : : case SlackUB:
2555 : 0 : break;
2556 : 0 : case SlackVLB: /* x[k] = lb[k] * x[kk] + x'[k] */
2557 : : case SlackVUB: /* x[k] = ub[k] * x[kk] - x'[k] */
2558 : : {
2559 : 0 : Assert(!d_vars.isIntegerInput(x));
2560 : 0 : bool ub = (rep == SlackVUB);
2561 : : const VirtualBound& vb =
2562 [ - - ]: 0 : ub ? d_pad.d_vub[x] : d_pad.d_vlb[x];
2563 : 0 : Assert(vb.x == x);
2564 : 0 : ArithVar y = vb.y;
2565 : 0 : Rational prod = c * vb.d;
2566 [ - - ]: 0 : if(mod.isKey(y)){
2567 : 0 : mod.get(x) += prod;
2568 : : }else{
2569 : 0 : mod.set(y, prod);
2570 : : }
2571 [ - - ]: 0 : if(ub){
2572 : 0 : mod.set(x, -c);
2573 : : }
2574 : 0 : Assert(vb.c != NullConstraint);
2575 : 0 : d_pad.d_explanation.insert(vb.c);
2576 : : }
2577 : 0 : break;
2578 : 0 : default:
2579 : 0 : return true;
2580 : : }
2581 : : }
2582 : 0 : removeZeroes(mod); /* if something cancelled we don't want it in the explanation */
2583 [ - - ]: 0 : for(iter = mod.begin(), iend = mod.end(); iter != iend; ++iter){
2584 : 0 : ArithVar x = *iter;
2585 [ - - ]: 0 : if(!d_pad.d_slacks.isKey(x)){ return true; }
2586 : :
2587 : 0 : SlackReplace rep = d_pad.d_slacks[x];
2588 [ - - ][ - ]: 0 : switch(rep){
2589 : 0 : case SlackLB: /* x = lb + x' */
2590 : : case SlackUB: /* x = ub - x' */
2591 : : {
2592 : 0 : bool ub = (rep == SlackUB);
2593 [ - - ]: 0 : ConstraintP b = ub ? d_vars.getUpperBoundConstraint(x):
2594 : 0 : d_vars.getLowerBoundConstraint(x);
2595 : :
2596 : 0 : Assert(b != NullConstraint);
2597 : 0 : Assert(b->getValue().infinitesimalIsZero());
2598 : 0 : const Rational& c = mod.get(x);
2599 : 0 : modRhs -= c * b->getValue().getNoninfinitesimalPart();
2600 [ - - ]: 0 : if(ub){
2601 : 0 : mod.set(x, -c);
2602 : : }
2603 : 0 : d_pad.d_explanation.insert(b);
2604 : : }
2605 : 0 : break;
2606 : 0 : case SlackVLB: /* handled earlier */
2607 : : case SlackVUB:
2608 : 0 : break;
2609 : 0 : default:
2610 : 0 : return true;
2611 : : }
2612 : : }
2613 : 0 : removeZeroes(mod);
2614 : 0 : return false;
2615 : : }
2616 : :
2617 : 0 : bool ApproxGLPK::makeRangeForComplemented(int nid, const MirInfo& mir){
2618 : 0 : DenseMap<Rational>& alpha = d_pad.d_alpha.lhs;
2619 : 0 : int M = mir.getMAtCreation();
2620 : 0 : int N = mir.getN();
2621 : 0 : DenseMap<Rational>& compRanges = d_pad.d_compRanges;
2622 : :
2623 : 0 : int complemented = 0;
2624 : :
2625 [ - - ]: 0 : for(int j = 1; j <= M + N; ++j){
2626 [ - - ]: 0 : if(mir.cset[j] != 0){
2627 : 0 : complemented++;
2628 : 0 : ArithVar x = _getArithVar(nid, M, j);
2629 [ - - ]: 0 : if(!alpha.isKey(x)){ return true; }
2630 [ - - ]: 0 : if(!d_vars.isIntegerInput(x)){ return true; }
2631 : 0 : Assert(d_pad.d_slacks.isKey(x));
2632 : 0 : Assert(d_pad.d_slacks[x] == SlackLB || d_pad.d_slacks[x] == SlackUB);
2633 : :
2634 : 0 : ConstraintP lb = d_vars.getLowerBoundConstraint(x);
2635 : 0 : ConstraintP ub = d_vars.getUpperBoundConstraint(x);
2636 : :
2637 [ - - ]: 0 : if(lb == NullConstraint) { return true; }
2638 [ - - ]: 0 : if(ub == NullConstraint) { return true; }
2639 : :
2640 [ - - ]: 0 : if(!lb->getValue().infinitesimalIsZero()){
2641 : 0 : return true;
2642 : : }
2643 [ - - ]: 0 : if(!ub->getValue().infinitesimalIsZero()){
2644 : 0 : return true;
2645 : : }
2646 : :
2647 : 0 : const Rational& uval = ub->getValue().getNoninfinitesimalPart();
2648 : 0 : const Rational& lval = lb->getValue().getNoninfinitesimalPart();
2649 : :
2650 : 0 : d_pad.d_explanation.insert(lb);
2651 : 0 : d_pad.d_explanation.insert(ub);
2652 : :
2653 : 0 : Rational u = uval - lval;
2654 : : // u is the same for both rep == LP and rep == UB
2655 [ - - ]: 0 : if(compRanges.isKey(x)) { return true; }
2656 : 0 : compRanges.set(x,u);
2657 : : }
2658 : : }
2659 : :
2660 [ - - ]: 0 : Trace("approx::mir") << "makeRangeForComplemented()" << complemented
2661 : 0 : << std::endl;
2662 : 0 : return false;
2663 : : }
2664 : :
2665 : :
2666 : 0 : bool ApproxGLPK::constructMixedKnapsack(){
2667 : 0 : const DenseMap<Rational>& mod = d_pad.d_mod.lhs;
2668 : 0 : const Rational& modRhs = d_pad.d_mod.rhs;
2669 : 0 : DenseMap<Rational>& alpha = d_pad.d_alpha.lhs;
2670 : 0 : Rational& beta = d_pad.d_alpha.rhs;
2671 : :
2672 : 0 : Assert(alpha.empty());
2673 : 0 : beta = modRhs;
2674 : :
2675 : 0 : unsigned intVars = 0;
2676 : 0 : unsigned remain = 0;
2677 : 0 : unsigned dropped = 0;
2678 : 0 : DenseMap<Rational>::const_iterator iter, iend;
2679 [ - - ]: 0 : for(iter = mod.begin(), iend = mod.end(); iter != iend; ++iter){
2680 : 0 : ArithVar v = *iter;
2681 : 0 : const Rational& c = mod[v];
2682 : 0 : Assert(!c.isZero());
2683 [ - - ]: 0 : if(d_vars.isIntegerInput(v)){
2684 : 0 : intVars++;
2685 : 0 : alpha.set(v, c);
2686 [ - - ]: 0 : }else if(c.sgn() < 0){
2687 : 0 : remain++;
2688 : 0 : alpha.set(v, c);
2689 : : }else{
2690 : 0 : dropped++;
2691 : : }
2692 : : }
2693 : :
2694 [ - - ]: 0 : Trace("approx::mir") << "constructMixedKnapsack() "
2695 : 0 : << " dropped " << dropped << " remain " << remain
2696 : 0 : << " intVars " << intVars << std::endl;
2697 : 0 : return intVars == 0; // if this is 0 we have failed
2698 : : }
2699 : :
2700 : 0 : bool ApproxGLPK::attemptConstructTableRow(int nid, int M, const PrimitiveVec& vec){
2701 : 0 : bool failed = guessCoefficientsConstructTableRow(nid, M, vec);
2702 [ - - ]: 0 : if(failed){
2703 : 0 : failed = gaussianElimConstructTableRow(nid, M, vec);
2704 : : }
2705 : :
2706 : 0 : return failed;
2707 : : }
2708 : :
2709 : 0 : bool ApproxGLPK::gaussianElimConstructTableRow(int nid, int M, const PrimitiveVec& vec){
2710 : 0 : TimerStat::CodeTimer codeTimer(d_stats.d_gaussianElimConstructTime);
2711 : 0 : ++d_stats.d_gaussianElimConstruct;
2712 : :
2713 : 0 : ArithVar basic = d_pad.d_basic;
2714 : 0 : DenseMap<Rational>& tab = d_pad.d_tabRow.lhs;
2715 : 0 : tab.purge();
2716 : 0 : d_pad.d_tabRow.rhs = Rational(0);
2717 : 0 : Assert(basic != ARITHVAR_SENTINEL);
2718 : 0 : Assert(tab.empty());
2719 : 0 : Assert(d_pad.d_tabRow.rhs.isZero());
2720 : :
2721 [ - - ]: 0 : if(d_vars.isAuxiliary(basic)) { return true; }
2722 : :
2723 [ - - ]: 0 : if(TraceIsOn("gaussianElimConstructTableRow")){
2724 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2725 : 0 : << "1 gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2726 : 0 : << std::endl;
2727 [ - - ]: 0 : vec.print(Trace("gaussianElimConstructTableRow"));
2728 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2729 : 0 : << "match " << basic << "(" << d_vars.asNode(basic) << ")" << std::endl;
2730 : : }
2731 : :
2732 : 0 : std::set<ArithVar> onrow;
2733 [ - - ]: 0 : for(int i = 1; i <= vec.len; ++i){
2734 : 0 : int ind = vec.inds[i];
2735 : 0 : ArithVar var = _getArithVar(nid, M, ind);
2736 [ - - ]: 0 : if(var == ARITHVAR_SENTINEL){
2737 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2738 : 0 : << "couldn't find" << ind << " " << M << " " << nid << std::endl;
2739 : 0 : return true;
2740 : : }
2741 : 0 : onrow.insert(var);
2742 : : }
2743 : :
2744 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2745 : 0 : << "2 gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2746 : 0 : << std::endl;
2747 : :
2748 : 0 : Matrix<Rational> A;
2749 : 0 : A.increaseSizeTo(d_vars.getNumberOfVariables());
2750 : 0 : std::vector< std::pair<RowIndex, ArithVar> > rows;
2751 : : // load the rows for auxiliary variables into A
2752 [ - - ]: 0 : for (ArithVar v : onrow)
2753 : : {
2754 [ - - ]: 0 : if(d_vars.isAuxiliary(v)){
2755 : 0 : Assert(d_vars.hasNode(v));
2756 : :
2757 : 0 : std::vector<Rational> coeffs;
2758 : 0 : std::vector<ArithVar> vars;
2759 : :
2760 : 0 : coeffs.push_back(Rational(-1));
2761 : 0 : vars.push_back(v);
2762 : :
2763 : 0 : Node n = d_vars.asNode(v);
2764 : 0 : Polynomial p = Polynomial::parsePolynomial(n);
2765 : 0 : Polynomial::iterator j = p.begin(), jend=p.end();
2766 [ - - ]: 0 : for(j=p.begin(), jend=p.end(); j!=jend; ++j){
2767 : 0 : Monomial m = *j;
2768 [ - - ]: 0 : if(m.isConstant()) { return true; }
2769 : 0 : VarList vl = m.getVarList();
2770 [ - - ]: 0 : if(!d_vars.hasArithVar(vl.getNode())){ return true; }
2771 : 0 : ArithVar x = d_vars.asArithVar(vl.getNode());
2772 : 0 : const Rational& q = m.getConstant().getValue();
2773 : 0 : coeffs.push_back(q); vars.push_back(x);
2774 : : }
2775 : 0 : RowIndex rid = A.addRow(coeffs, vars);
2776 : 0 : rows.push_back(std::make_pair(rid, ARITHVAR_SENTINEL));
2777 : : }
2778 : : }
2779 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2780 : 0 : << "3 gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2781 : 0 : << std::endl;
2782 : :
2783 [ - - ]: 0 : for(size_t i=0; i < rows.size(); ++i){
2784 : 0 : RowIndex rid = rows[i].first;
2785 : 0 : Assert(rows[i].second == ARITHVAR_SENTINEL);
2786 : :
2787 : : // substitute previous rows
2788 [ - - ]: 0 : for(size_t j=0; j < i; j++){
2789 : 0 : RowIndex prevRow = rows[j].first;
2790 : 0 : ArithVar other = rows[j].second;
2791 : 0 : Assert(other != ARITHVAR_SENTINEL);
2792 : 0 : const Matrix<Rational>::Entry& e = A.findEntry(rid, other);
2793 [ - - ]: 0 : if(!e.blank()){
2794 : : // r_p : 0 = -1 * other + sum a_i x_i
2795 : : // rid : 0 = e * other + sum b_i x_i
2796 : : // rid += e * r_p
2797 : : // : 0 = 0 * other + ...
2798 : 0 : Assert(!e.getCoefficient().isZero());
2799 : :
2800 : 0 : Rational cp = e.getCoefficient();
2801 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2802 : 0 : << "on " << rid << " subst " << cp << "*" << prevRow << " " << other
2803 : 0 : << std::endl;
2804 : 0 : A.rowPlusRowTimesConstant(rid, prevRow, cp);
2805 : : }
2806 : : }
2807 [ - - ]: 0 : if(TraceIsOn("gaussianElimConstructTableRow")){
2808 [ - - ]: 0 : A.printMatrix(Trace("gaussianElimConstructTableRow"));
2809 : : }
2810 : :
2811 : : // solve the row for anything other than non-basics
2812 : 0 : bool solveForBasic = (i + 1 == rows.size());
2813 : 0 : Rational q;
2814 : 0 : ArithVar s = ARITHVAR_SENTINEL;
2815 : 0 : Matrix<Rational>::RowIterator k = A.getRow(rid).begin();
2816 : 0 : Matrix<Rational>::RowIterator k_end = A.getRow(rid).end();
2817 [ - - ]: 0 : for(; k != k_end; ++k){
2818 : 0 : const Matrix<Rational>::Entry& e = *k;
2819 : 0 : ArithVar colVar = e.getColVar();
2820 : 0 : bool selectColVar = false;
2821 [ - - ]: 0 : if(colVar == basic){
2822 : 0 : selectColVar = solveForBasic;
2823 [ - - ]: 0 : }else if(onrow.find(colVar) == onrow.end()) {
2824 : 0 : selectColVar = true;
2825 : : }
2826 [ - - ]: 0 : if(selectColVar){
2827 : 0 : s = colVar;
2828 : 0 : q = e.getCoefficient();
2829 : : }
2830 : : }
2831 : 0 : if(s == ARITHVAR_SENTINEL || q.isZero()){
2832 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2833 : 0 : << "3 fail gaussianElimConstructTableRow(" << nid << ", " << basic
2834 : 0 : << ")" << std::endl;
2835 : 0 : return true;
2836 : : }else{
2837 : : // 0 = q * s + sum c_i * x_i
2838 : 0 : Rational mult = -(q.inverse());
2839 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2840 : 0 : << "selecting " << s << " : complexity " << mult.complexity() << " "
2841 : 0 : << mult << std::endl;
2842 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2843 : 0 : << "selecting " << rid << " " << s << std::endl;
2844 : 0 : A.multiplyRowByConstant(rid, mult);
2845 : 0 : rows[i].second = s;
2846 : : }
2847 : : }
2848 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2849 : 0 : << "4 gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2850 : 0 : << std::endl;
2851 : :
2852 [ - - ]: 0 : if(rows.empty()) {
2853 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2854 : 0 : << "4 fail 1 gaussianElimConstructTableRow(" << nid << ", " << basic
2855 : 0 : << ")" << std::endl;
2856 : 0 : return true;
2857 : : }
2858 : 0 : RowIndex rid_last = rows.back().first;
2859 : 0 : ArithVar rid_var = rows.back().second;
2860 [ - - ]: 0 : if(rid_var != basic){
2861 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2862 : 0 : << "4 fail 2 gaussianElimConstructTableRow(" << nid << ", " << basic
2863 : 0 : << ")" << std::endl;
2864 : 0 : return true;
2865 : : }
2866 : :
2867 : 0 : Assert(tab.empty());
2868 : :
2869 : 0 : Matrix<Rational>::RowIterator k = A.getRow(rid_last).begin();
2870 : 0 : Matrix<Rational>::RowIterator k_end = A.getRow(rid_last).end();
2871 [ - - ]: 0 : for(; k != k_end; ++k){
2872 : 0 : const Matrix<Rational>::Entry& e = *k;
2873 : 0 : tab.set(e.getColVar(), e.getCoefficient());
2874 : : }
2875 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2876 : 0 : << "5 gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2877 : 0 : << std::endl;
2878 [ - - ]: 0 : if(!tab.isKey(basic)){
2879 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2880 : 0 : << "5 fail 1 gaussianElimConstructTableRow(" << nid << ", " << basic
2881 : 0 : << ")" << std::endl;
2882 : 0 : return true;
2883 : : }
2884 [ - - ]: 0 : if(tab[basic] != Rational(-1)){
2885 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2886 : 0 : << "5 fail 2 gaussianElimConstructTableRow(" << nid << ", " << basic
2887 : 0 : << ")" << std::endl;
2888 : 0 : return true;
2889 : : }
2890 : :
2891 : 0 : tab.remove(basic);
2892 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2893 : 0 : << "6 gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2894 : 0 : << std::endl;
2895 : :
2896 [ - - ]: 0 : if(vec.len < 0 ){
2897 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2898 : 0 : << "6 fail 1 gaussianElimConstructTableRow(" << nid << ", " << basic
2899 : 0 : << ")" << std::endl;
2900 : 0 : return true;
2901 : : }
2902 [ - - ]: 0 : if(tab.size() != ((unsigned)vec.len) ) {
2903 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2904 : 0 : << "6 fail 2 gaussianElimConstructTableRow(" << nid << ", " << basic
2905 : 0 : << ")" << tab.size() << " " << vec.len << std::endl;
2906 : 0 : return true;
2907 : : }
2908 : :
2909 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2910 : 0 : << "7 gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2911 : 0 : << std::endl;
2912 : :
2913 [ - - ]: 0 : for(int i = 1; i <= vec.len; ++i){
2914 : 0 : int ind = vec.inds[i];
2915 : 0 : double coeff = vec.coeffs[i];
2916 : 0 : ArithVar var = _getArithVar(nid, M, ind);
2917 : 0 : Assert(var != ARITHVAR_SENTINEL);
2918 [ - - ]: 0 : if(!tab.isKey(var)){
2919 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2920 : 0 : << "7 fail 1 gaussianElimConstructTableRow(" << nid << ", " << basic
2921 : 0 : << ")" << std::endl;
2922 : 0 : return true;
2923 : : }
2924 : :
2925 : 0 : double est = tab[var].getDouble();
2926 : :
2927 [ - - ]: 0 : if (!ApproxGLPK::roughlyEqual(coeff, est))
2928 : : {
2929 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2930 : 0 : << "7 fail 2 gaussianElimConstructTableRow(" << nid << ", " << basic
2931 : 0 : << ")"
2932 : 0 : << " boink on " << ind << " " << var << " " << est << std::endl;
2933 : 0 : return true;
2934 : : }
2935 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2936 : 0 : << var << " cfe " << coeff << std::endl;
2937 : : }
2938 : :
2939 [ - - ]: 0 : Trace("gaussianElimConstructTableRow")
2940 : 0 : << "gaussianElimConstructTableRow(" << nid << ", " << basic << ")"
2941 : 0 : << " superduper" << std::endl;
2942 : :
2943 : 0 : return false;
2944 : : }
2945 : 0 : bool ApproxGLPK::guessCoefficientsConstructTableRow(int nid, int M, const PrimitiveVec& vec){
2946 [ - - ]: 0 : for(size_t i=0; i < d_denomGuesses.size(); ++i){
2947 : 0 : const Integer& D = d_denomGuesses[i];
2948 [ - - ]: 0 : if(!guessCoefficientsConstructTableRow(nid, M, vec, D)){
2949 : 0 : d_stats.d_averageGuesses << i+1;
2950 [ - - ]: 0 : Trace("approx::gmi") << "guesseditat " << i << " D=" << D << std::endl;
2951 : 0 : return false;
2952 : : }
2953 : : }
2954 : 0 : return true;
2955 : : }
2956 : 0 : bool ApproxGLPK::guessCoefficientsConstructTableRow(int nid, int M, const PrimitiveVec& vec, const Integer& D){
2957 : 0 : ArithVar basic = d_pad.d_basic;
2958 : 0 : DenseMap<Rational>& tab = d_pad.d_tabRow.lhs;
2959 : 0 : tab.purge();
2960 : 0 : d_pad.d_tabRow.rhs = Rational(0);
2961 : 0 : Assert(basic != ARITHVAR_SENTINEL);
2962 : 0 : Assert(tab.empty());
2963 : 0 : Assert(d_pad.d_tabRow.rhs.isZero());
2964 : :
2965 [ - - ]: 0 : if(TraceIsOn("guessCoefficientsConstructTableRow")){
2966 [ - - ]: 0 : Trace("guessCoefficientsConstructTableRow")
2967 : 0 : << "attemptConstructTableRow(" << nid << ", " << basic << ",...," << D
2968 : 0 : << ")" << std::endl;
2969 [ - - ]: 0 : vec.print(Trace("guessCoefficientsConstructTableRow"));
2970 [ - - ]: 0 : Trace("guessCoefficientsConstructTableRow")
2971 : 0 : << "match " << basic << "(" << d_vars.asNode(basic) << ")" << std::endl;
2972 : : }
2973 : :
2974 : 0 : tab.set(basic, Rational(-1));
2975 [ - - ]: 0 : for(int i = 1; i <= vec.len; ++i){
2976 : 0 : int ind = vec.inds[i];
2977 : 0 : double coeff = vec.coeffs[i];
2978 : 0 : ArithVar var = _getArithVar(nid, M, ind);
2979 [ - - ]: 0 : if(var == ARITHVAR_SENTINEL){
2980 [ - - ]: 0 : Trace("guessCoefficientsConstructTableRow")
2981 : 0 : << "couldn't find" << ind << " " << M << " " << nid << std::endl;
2982 : 0 : return true;
2983 : : }
2984 [ - - ]: 0 : Trace("guessCoefficientsConstructTableRow")
2985 : 0 : << "match " << ind << "," << var << "(" << d_vars.asNode(var) << ")"
2986 : 0 : << std::endl;
2987 : :
2988 : 0 : std::optional<Rational> cfe = estimateWithCFE(coeff, D);
2989 [ - - ]: 0 : if (!cfe)
2990 : : {
2991 : 0 : return true;
2992 : : }
2993 : 0 : tab.set(var, *cfe);
2994 [ - - ]: 0 : Trace("guessCoefficientsConstructTableRow")
2995 : 0 : << var << " cfe " << cfe << std::endl;
2996 : : }
2997 [ - - ]: 0 : if(!guessIsConstructable(tab)){
2998 [ - - ]: 0 : Trace("guessCoefficientsConstructTableRow")
2999 : 0 : << "failed to construct with " << D << std::endl;
3000 : 0 : return true;
3001 : : }
3002 : 0 : tab.remove(basic);
3003 : 0 : return false;
3004 : : }
3005 : :
3006 : : /* Maps an ArithVar to either an upper/lower bound */
3007 : 0 : bool ApproxGLPK::constructGmiCut(){
3008 : 0 : const DenseMap<Rational>& tabRow = d_pad.d_tabRow.lhs;
3009 : 0 : const DenseMap<ConstraintP>& toBound = d_pad.d_toBound;
3010 : 0 : DenseMap<Rational>& cut = d_pad.d_cut.lhs;
3011 : 0 : std::set<ConstraintP>& explanation = d_pad.d_explanation;
3012 : 0 : Rational& rhs = d_pad.d_cut.rhs;
3013 : :
3014 : 0 : DenseMap<Rational>::const_iterator iter, end;
3015 : 0 : Assert(cut.empty());
3016 : :
3017 : : // compute beta for a "fake" assignment
3018 : : bool anyInf;
3019 : 0 : DeltaRational dbeta = sumConstraints(tabRow, toBound, &anyInf);
3020 : 0 : const Rational& beta = dbeta.getNoninfinitesimalPart();
3021 [ - - ]: 0 : Trace("approx::gmi") << dbeta << std::endl;
3022 : 0 : if(anyInf || beta.isIntegral()){ return true; }
3023 : :
3024 : 0 : Rational one = Rational(1);
3025 : 0 : Rational fbeta = beta.floor_frac();
3026 : 0 : rhs = fbeta;
3027 : 0 : Assert(fbeta.sgn() > 0);
3028 : 0 : Assert(fbeta < one);
3029 : 0 : Rational one_sub_fbeta = one - fbeta;
3030 [ - - ]: 0 : for(iter = tabRow.begin(), end = tabRow.end(); iter != end; ++iter){
3031 : 0 : ArithVar x = *iter;
3032 : 0 : const Rational& psi = tabRow[x];
3033 : 0 : ConstraintP c = toBound[x];
3034 : 0 : const Rational& bound = c->getValue().getNoninfinitesimalPart();
3035 [ - - ]: 0 : if(d_vars.boundsAreEqual(x)){
3036 : : // do not add a coefficient
3037 : : // implictly substitute the variable w/ its constraint
3038 : 0 : std::pair<ConstraintP, ConstraintP> exp = d_vars.explainEqualBounds(x);
3039 : 0 : explanation.insert(exp.first);
3040 [ - - ]: 0 : if(exp.second != NullConstraint){
3041 : 0 : explanation.insert(exp.second);
3042 : : }
3043 : 0 : }else if(d_vars.isIntegerInput(x) && psi.isIntegral()){
3044 : : // do not add a coefficient
3045 : : // nothing to explain
3046 [ - - ]: 0 : Trace("approx::gmi") << "skipping " << x << std::endl;
3047 : : }else{
3048 : 0 : explanation.insert(c);
3049 : 0 : Rational phi;
3050 [ - - ]: 0 : Rational alpha = (c->isUpperBound() ? psi : -psi);
3051 : :
3052 : : // x - ub <= 0 and lb - x <= 0
3053 [ - - ]: 0 : if(d_vars.isIntegerInput(x)){
3054 : 0 : Assert(!psi.isIntegral());
3055 : : // alpha = slack_sgn * psi
3056 : 0 : Rational falpha = alpha.floor_frac();
3057 : 0 : Assert(falpha.sgn() > 0);
3058 : 0 : Assert(falpha < one);
3059 [ - - ]: 0 : phi = (falpha <= fbeta) ?
3060 : 0 : falpha : ((fbeta / one_sub_fbeta) * (one - falpha));
3061 : : }else{
3062 [ - - ]: 0 : phi = (alpha >= 0) ?
3063 : 0 : alpha : ((fbeta / one_sub_fbeta) * (- alpha));
3064 : : }
3065 : 0 : Assert(phi.sgn() != 0);
3066 [ - - ]: 0 : if(c->isUpperBound()){
3067 : 0 : cut.set(x, -phi);
3068 : 0 : rhs -= phi * bound;
3069 : : }else{
3070 : 0 : cut.set(x, phi);
3071 : 0 : rhs += phi * bound;
3072 : : }
3073 : : }
3074 : : }
3075 [ - - ]: 0 : if (TraceIsOn("approx::gmi"))
3076 : : {
3077 [ - - ]: 0 : Trace("approx::gmi") << "pre removeSlackVariables";
3078 [ - - ]: 0 : d_pad.d_cut.print(Trace("approx::gmi"));
3079 [ - - ]: 0 : Trace("approx::gmi") << std::endl;
3080 : : }
3081 : 0 : removeAuxillaryVariables(d_vars, cut);
3082 : :
3083 [ - - ]: 0 : if (TraceIsOn("approx::gmi"))
3084 : : {
3085 [ - - ]: 0 : Trace("approx::gmi") << "post removeAuxillaryVariables";
3086 [ - - ]: 0 : d_pad.d_cut.print(Trace("approx::gmi"));
3087 [ - - ]: 0 : Trace("approx::gmi") << std::endl;
3088 : : }
3089 : 0 : removeFixed(d_vars, d_pad.d_cut, explanation);
3090 : :
3091 [ - - ]: 0 : if (TraceIsOn("approx::gmi"))
3092 : : {
3093 [ - - ]: 0 : Trace("approx::gmi") << "post removeFixed";
3094 [ - - ]: 0 : d_pad.d_cut.print(Trace("approx::gmi"));
3095 [ - - ]: 0 : Trace("approx::gmi") << std::endl;
3096 : : }
3097 : 0 : return false;
3098 : : }
3099 : :
3100 : 0 : void ApproxGLPK::tryCut(int nid, CutInfo& cut)
3101 : : {
3102 : 0 : Assert(!cut.reconstructed());
3103 : 0 : Assert(cut.getKlass() != RowsDeletedKlass);
3104 : 0 : bool failure = false;
3105 [ - - ][ - - ]: 0 : switch(cut.getKlass()){
3106 : 0 : case GmiCutKlass:
3107 : 0 : failure = attemptGmi(nid, static_cast<const GmiInfo&>(cut));
3108 : 0 : break;
3109 : 0 : case MirCutKlass:
3110 : 0 : failure = attemptMir(nid, static_cast<const MirInfo&>(cut));
3111 : 0 : break;
3112 : 0 : case BranchCutKlass:
3113 [ - - ]: 0 : failure = attemptBranchCut(nid, dynamic_cast<const BranchCutInfo&>(cut));
3114 : 0 : break;
3115 : 0 : default:
3116 : 0 : break;
3117 : : }
3118 : 0 : Assert(failure == d_pad.d_failure);
3119 : :
3120 [ - - ]: 0 : if(!failure){
3121 : : // move the pad to the cut
3122 : 0 : cut.setReconstruction(d_pad.d_cut);
3123 : :
3124 [ - - ]: 0 : if(cut.getKlass() != BranchCutKlass){
3125 : 0 : std::set<ConstraintP>& exp = d_pad.d_explanation;
3126 : 0 : ConstraintCPVec asvec(exp.begin(), exp.end());
3127 : 0 : cut.swapExplanation(asvec);
3128 : : }
3129 : : }else{
3130 [ - - ]: 0 : Trace("approx") << "failure " << cut.getKlass() << std::endl;
3131 : : }
3132 : 0 : }
3133 : :
3134 : : } // namespace arith
3135 : : } // namespace theory
3136 : : } // namespace cvc5::internal
3137 : : /* End GPLK implementation. */
3138 : : #endif /*#ifdef CVC5_USE_GLPK */
3139 : :
3140 : : /* Begin GPLK/NOGLPK Glue code. */
3141 : : namespace cvc5::internal {
3142 : : namespace theory {
3143 : : namespace arith::linear {
3144 : :
3145 : 6 : ApproximateSimplex* ApproximateSimplex::mkApproximateSimplexSolver(
3146 : : const ArithVariables& vars, TreeLog& l, ApproximateStatistics& s)
3147 : : {
3148 : : #ifdef CVC5_USE_GLPK
3149 : 6 : return new ApproxGLPK(vars, l, s);
3150 : : #else
3151 : : Unimplemented() << "Approximate simplex solver requires GLPK";
3152 : : #endif
3153 : : }
3154 : :
3155 : 92 : bool ApproximateSimplex::enabled()
3156 : : {
3157 : : #ifdef CVC5_USE_GLPK
3158 : 92 : return true;
3159 : : #else
3160 : : return false;
3161 : : #endif
3162 : : }
3163 : :
3164 : 6 : ApproximateStatistics::ApproximateStatistics(StatisticsRegistry& sr)
3165 : 6 : : d_branchMaxDepth(sr.registerInt("z::approx::branchMaxDepth")),
3166 : 6 : d_branchesMaxOnAVar(sr.registerInt("z::approx::branchesMaxOnAVar")),
3167 : : d_gaussianElimConstructTime(
3168 : 6 : sr.registerTimer("z::approx::gaussianElimConstruct::time")),
3169 : : d_gaussianElimConstruct(
3170 : 6 : sr.registerInt("z::approx::gaussianElimConstruct::calls")),
3171 : 6 : d_averageGuesses(sr.registerAverage("z::approx::averageGuesses"))
3172 : : {
3173 : 6 : }
3174 : :
3175 : 0 : std::ostream& operator<<(std::ostream& out, MipResult res)
3176 : : {
3177 [ - - ][ - - ]: 0 : switch (res)
[ - - ][ - ]
3178 : : {
3179 : 0 : case MipUnknown: out << "MipUnknown"; break;
3180 : 0 : case MipBingo: out << "MipBingo"; break;
3181 : 0 : case MipClosed: out << "MipClosed"; break;
3182 : 0 : case BranchesExhausted: out << "BranchesExhausted"; break;
3183 : 0 : case PivotsExhauasted: out << "PivotsExhauasted"; break;
3184 : 0 : case ExecExhausted: out << "ExecExhausted"; break;
3185 : 0 : default: out << "Unexpected Mip Value!"; break;
3186 : : }
3187 : 0 : return out;
3188 : : }
3189 : :
3190 : : } // namespace arith::linear
3191 : : } // namespace theory
3192 : : } // namespace cvc5::internal
3193 : : /* End GPLK/NOGLPK Glue code. */
|
