GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/approx_simplex.cpp Lines: 4 151 2.6 %
Date: 2021-03-22 Branches: 2 447 0.4 %

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