src/theory/arith/approx_simplex.cpp

   1 #include "cvc4autoconfig.h"
   2
   3 #include "theory/arith/approx_simplex.h"
   4 #include "theory/arith/normal_form.h"
   5 #include "theory/arith/constraint.h"
   6 #include <math.h>
   7 #include <cmath>
   8
   9 using namespace std;
  10
  11 namespace CVC4 {
  12 namespace theory {
  13 namespace arith {
  14
  15 ApproximateSimplex::ApproximateSimplex() :
  16   d_pivotLimit(std::numeric_limits<int>::max())
  17 {}
  18
  19 void ApproximateSimplex::setPivotLimit(int pivotLimit){
  20   Assert(pivotLimit >= 0);
  21   d_pivotLimit = pivotLimit;
  22 }
  23
  24 const double ApproximateSimplex::SMALL_FIXED_DELTA = .000000001;
  25 const double ApproximateSimplex::TOLERENCE = 1 + .000000001;
  26
  27 bool ApproximateSimplex::roughlyEqual(double a, double b){
  28   if (a == 0){
  29     return -SMALL_FIXED_DELTA <= b && b <= SMALL_FIXED_DELTA;
  30   }else if (b == 0){
  31     return -SMALL_FIXED_DELTA <= a && a <= SMALL_FIXED_DELTA;
  32   }else{
  33     return std::abs(b/a) <= TOLERENCE && std::abs(a/b) <= TOLERENCE;
  34   }
  35 }
  36
  37 Rational ApproximateSimplex::cfeToRational(const vector<Integer>& exp){
  38   if(exp.empty()){
  39     return Rational(0);
  40   }else{
  41     Rational result = exp.back();
  42     vector<Integer>::const_reverse_iterator exp_iter = exp.rbegin();
  43     vector<Integer>::const_reverse_iterator exp_end = exp.rend();
  44     ++exp_iter;
  45     while(exp_iter != exp_end){
  46       result = result.inverse();
  47       const Integer& i = *exp_iter;
  48       result += i;
  49       ++exp_iter;
  50     }
  51     return result;
  52   }
  53 }
  54 std::vector<Integer> ApproximateSimplex::rationalToCfe(const Rational& q, int depth){
  55   vector<Integer> mods;
  56   if(!q.isZero()){
  57     Rational carry = q;
  58     for(int i = 0; i <= depth; ++i){
  59       Assert(!carry.isZero());
  60       mods.push_back(Integer());
  61       Integer& back = mods.back();
  62       back = carry.floor();
  63       //cout << "  cfe["<<i<<"]: " << back << endl;
  64       carry -= back;
  65       if(carry.isZero()){
  66         break;
  67       }else if(ApproximateSimplex::roughlyEqual(carry.getDouble(), 0.0)){
  68         break;
  69       }else{
  70         carry = carry.inverse();
  71       }
  72     }
  73   }
  74   return mods;
  75 }
  76
  77 Rational ApproximateSimplex::estimateWithCFE(const Rational& q, int depth){
  78   std::vector<Integer> cfe = rationalToCfe(q,depth);
  79   return cfeToRational(cfe);
  80 }
  81
  82 Rational ApproximateSimplex::estimateWithCFE(double d){
  83   return estimateWithCFE(Rational::fromDouble(d), 10);
  84 }
  85
  86 class ApproxNoOp : public ApproximateSimplex {
  87 public:
  88   ApproxNoOp(const ArithVariables& vars){}
  89   ~ApproxNoOp(){}
  90
  91   virtual ApproxResult solveRelaxation(){
  92     return ApproxError;
  93   }
  94   virtual Solution extractRelaxation() const{
  95     return Solution();
  96   }
  97
  98   virtual ArithRatPairVec heuristicOptCoeffs() const{
  99     return ArithRatPairVec();
 100   }
 101
 102   virtual ApproxResult solveMIP(){
 103     return ApproxError;
 104   }
 105   virtual Solution extractMIP() const{
 106     return Solution();
 107   }
 108
 109   virtual void setOptCoeffs(const ArithRatPairVec& ref){}
 110 };
 111
 112 }/* CVC4::theory::arith namespace */
 113 }/* CVC4::theory namespace */
 114 }/* CVC4 namespace */
 115
 116 /* Begin the declaration of GLPK specific code. */
 117 #ifdef CVC4_USE_GLPK
 118 extern "C" {
 119 /* Sometimes the header is in a subdirectory glpk/, sometimes not.
 120  * The configure script figures it out. */
 121 #ifdef HAVE_GLPK_GLPK_H
 122 #  include <glpk/glpk.h>
 123 #else /* HAVE_GLPK_GLPK_H */
 124 #  include <glpk.h>
 125 #endif /* HAVE_GLPK_GLPK_H */
 126 }/* extern "C" */
 127
 128 namespace CVC4 {
 129 namespace theory {
 130 namespace arith {
 131
 132 class ApproxGLPK : public ApproximateSimplex {
 133 private:
 134   glp_prob* d_prob;
 135   const ArithVariables& d_vars;
 136
 137   DenseMap<int> d_colIndices;
 138   DenseMap<int> d_rowIndices;
 139
 140
 141   int d_instanceID;
 142
 143   bool d_solvedRelaxation;
 144   bool d_solvedMIP;
 145
 146   static int s_verbosity;
 147
 148 public:
 149   ApproxGLPK(const ArithVariables& vars);
 150   ~ApproxGLPK();
 151
 152   virtual ApproxResult solveRelaxation();
 153   virtual Solution extractRelaxation() const{
 154     return extractSolution(false);
 155   }
 156
 157   virtual ArithRatPairVec heuristicOptCoeffs() const;
 158
 159   virtual ApproxResult solveMIP();
 160   virtual Solution extractMIP() const{
 161     return extractSolution(true);
 162   }
 163   virtual void setOptCoeffs(const ArithRatPairVec& ref);
 164
 165   static void printGLPKStatus(int status, std::ostream& out);
 166 private:
 167   Solution extractSolution(bool mip) const;
 168   int guessDir(ArithVar v) const;
 169 };
 170
 171 int ApproxGLPK::s_verbosity = 0;
 172
 173 }/* CVC4::theory::arith namespace */
 174 }/* CVC4::theory namespace */
 175 }/* CVC4 namespace */
 176 #endif /*#ifdef CVC4_USE_GLPK */
 177 /* End the declaration of GLPK specific code. */
 178
 179 /* Begin GPLK/NOGLPK Glue code. */
 180 namespace CVC4 {
 181 namespace theory {
 182 namespace arith {
 183 ApproximateSimplex* ApproximateSimplex::mkApproximateSimplexSolver(const ArithVariables& vars){
 184 #ifdef CVC4_USE_GLPK
 185   return new ApproxGLPK(vars);
 186 #else
 187   return new ApproxNoOp(vars);
 188 #endif
 189 }
 190 bool ApproximateSimplex::enabled() {
 191 #ifdef CVC4_USE_GLPK
 192   return true;
 193 #else
 194   return false;
 195 #endif
 196 }
 197 }/* CVC4::theory::arith namespace */
 198 }/* CVC4::theory namespace */
 199 }/* CVC4 namespace */
 200 /* End GPLK/NOGLPK Glue code. */
 201
 202
 203 /* Begin GPLK implementation. */
 204 #ifdef CVC4_USE_GLPK
 205 namespace CVC4 {
 206 namespace theory {
 207 namespace arith {
 208
 209 ApproxGLPK::ApproxGLPK(const ArithVariables& avars) :
 210   d_vars(avars), d_solvedRelaxation(false), d_solvedMIP(false)
 211 {
 212   static int instance = 0;
 213   ++instance;
 214   d_instanceID = instance;
 215
 216   d_prob = glp_create_prob();
 217   glp_set_obj_dir(d_prob, GLP_MAX);
 218   glp_set_prob_name(d_prob, "ApproximateSimplex::approximateFindModel");
 219
 220   int numRows = 0;
 221   int numCols = 0;
 222
 223   // Assign each variable to a row and column variable as it appears in the input
 224   for(ArithVariables::var_iterator vi = d_vars.var_begin(), vi_end = d_vars.var_end(); vi != vi_end; ++vi){
 225     ArithVar v = *vi;
 226
 227     if(d_vars.isSlack(v)){
 228       ++numRows;
 229       d_rowIndices.set(v, numRows);
 230     }else{
 231       ++numCols;
 232       d_colIndices.set(v, numCols);
 233     }
 234   }
 235   glp_add_rows(d_prob, numRows);
 236   glp_add_cols(d_prob, numCols);
 237
 238   // Assign the upper/lower bounds and types to each variable
 239   for(ArithVariables::var_iterator vi = d_vars.var_begin(), vi_end = d_vars.var_end(); vi != vi_end; ++vi){
 240     ArithVar v = *vi;
 241
 242     if(s_verbosity >= 2){
 243       Message() << v  << " ";
 244       d_vars.printModel(v, Message());
 245     }
 246
 247     int type;
 248     double lb = 0.0;
 249     double ub = 0.0;
 250     if(d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
 251       if(d_vars.boundsAreEqual(v)){
 252         type = GLP_FX;
 253       }else{
 254         type = GLP_DB;
 255       }
 256       lb = d_vars.getLowerBound(v).approx(SMALL_FIXED_DELTA);
 257       ub = d_vars.getUpperBound(v).approx(SMALL_FIXED_DELTA);
 258     }else if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
 259       type = GLP_UP;
 260       ub = d_vars.getUpperBound(v).approx(SMALL_FIXED_DELTA);
 261     }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
 262       type = GLP_LO;
 263       lb = d_vars.getLowerBound(v).approx(SMALL_FIXED_DELTA);
 264     }else{
 265       type = GLP_FR;
 266     }
 267
 268     if(d_vars.isSlack(v)){
 269       int rowIndex = d_rowIndices[v];
 270       glp_set_row_bnds(d_prob, rowIndex, type, lb, ub);
 271     }else{
 272       int colIndex = d_colIndices[v];
 273       int kind = d_vars.isInteger(v) ? GLP_IV : GLP_CV;
 274       glp_set_col_kind(d_prob, colIndex, kind);
 275       glp_set_col_bnds(d_prob, colIndex, type, lb, ub);
 276     }
 277   }
 278
 279   // Count the number of entries
 280   int numEntries = 0;
 281   for(DenseMap<int>::const_iterator i = d_rowIndices.begin(), i_end = d_rowIndices.end(); i != i_end; ++i){
 282     ArithVar v = *i;
 283     Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
 284     for(Polynomial::iterator i = p.begin(), end = p.end(); i != end; ++i){
 285       ++numEntries;
 286     }
 287   }
 288
 289   int* ia = new int[numEntries+1];
 290   int* ja = new int[numEntries+1];
 291   double* ar = new double[numEntries+1];
 292
 293   int entryCounter = 0;
 294   for(DenseMap<int>::const_iterator i = d_rowIndices.begin(), i_end = d_rowIndices.end(); i != i_end; ++i){
 295     ArithVar v = *i;
 296     int rowIndex = d_rowIndices[v];
 297
 298     Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
 299
 300     for(Polynomial::iterator i = p.begin(), end = p.end(); i != end; ++i){
 301
 302       const Monomial& mono = *i;
 303       const Constant& constant = mono.getConstant();
 304       const VarList& variable = mono.getVarList();
 305
 306       Node n = variable.getNode();
 307
 308       Assert(d_vars.hasArithVar(n));
 309       ArithVar av = d_vars.asArithVar(n);
 310       int colIndex = d_colIndices[av];
 311       double coeff = constant.getValue().getDouble();
 312
 313       ++entryCounter;
 314       ia[entryCounter] = rowIndex;
 315       ja[entryCounter] = colIndex;
 316       ar[entryCounter] = coeff;
 317     }
 318   }
 319   glp_load_matrix(d_prob, numEntries, ia, ja, ar);
 320
 321   delete[] ia;
 322   delete[] ja;
 323   delete[] ar;
 324 }
 325 int ApproxGLPK::guessDir(ArithVar v) const{
 326   if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
 327     return -1;
 328   }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
 329     return 1;
 330   }else if(!d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
 331     return 0;
 332   }else{
 333     int ubSgn = d_vars.getUpperBound(v).sgn();
 334     int lbSgn = d_vars.getLowerBound(v).sgn();
 335
 336     if(ubSgn != 0 && lbSgn == 0){
 337       return -1;
 338     }else if(ubSgn == 0 && lbSgn != 0){
 339       return 1;
 340     }
 341
 342     return 1;
 343   }
 344 }
 345
 346 ArithRatPairVec ApproxGLPK::heuristicOptCoeffs() const{
 347   ArithRatPairVec ret;
 348
 349   // Strategies are guess:
 350   // 1 simple shared "ceiling" variable: danoint, pk1
 351   //  x1 >= c, x1 >= tmp1, x1 >= tmp2, ...
 352   // 1 large row: fixnet, vpm2, pp08a
 353   //  (+ .......... ) <= c
 354   // Not yet supported:
 355   // 1 complex shared "ceiling" variable: opt1217
 356   //  x1 >= c, x1 >= (+ ..... ), x1 >= (+ ..... )
 357   //  and all of the ... are the same sign
 358
 359
 360   // Candidates:
 361   // 1) Upper and lower bounds are not equal.
 362   // 2) The variable is not integer
 363   // 3a) For columns look for a ceiling variable
 364   // 3B) For rows look for a large row with
 365
 366   DenseMap<BoundCounts> d_colCandidates;
 367   DenseMap<uint32_t> d_rowCandidates;
 368
 369   double sumRowLength = 0.0;
 370   uint32_t maxRowLength = 0;
 371   for(ArithVariables::var_iterator vi = d_vars.var_begin(), vi_end = d_vars.var_end(); vi != vi_end; ++vi){
 372     ArithVar v = *vi;
 373
 374     if(s_verbosity >= 2){
 375       Message() << v  << " ";
 376       d_vars.printModel(v, Message());
 377     }
 378
 379     int type;
 380     if(d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
 381       if(d_vars.boundsAreEqual(v)){
 382         type = GLP_FX;
 383       }else{
 384         type = GLP_DB;
 385       }
 386     }else if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){
 387       type = GLP_UP;
 388     }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){
 389       type = GLP_LO;
 390     }else{
 391       type = GLP_FR;
 392     }
 393
 394     if(type != GLP_FX && type != GLP_FR){
 395
 396       if(d_vars.isSlack(v)){
 397         Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
 398         uint32_t len = p.size();
 399         d_rowCandidates.set(v, len);
 400         sumRowLength += len;
 401         maxRowLength =std::max(maxRowLength, len);
 402       }else if(!d_vars.isInteger(v)){
 403         d_colCandidates.set(v, BoundCounts());
 404       }
 405     }
 406   }
 407
 408   uint32_t maxCount = 0;
 409   for(DenseMap<int>::const_iterator i = d_rowIndices.begin(), i_end = d_rowIndices.end(); i != i_end; ++i){
 410     ArithVar v = *i;
 411
 412     bool lbCap = d_vars.hasLowerBound(v) && !d_vars.hasUpperBound(v);
 413     bool ubCap = !d_vars.hasLowerBound(v) && d_vars.hasUpperBound(v);
 414
 415     if(lbCap || ubCap){
 416       Constraint b = lbCap ? d_vars.getLowerBoundConstraint(v)
 417         : d_vars.getUpperBoundConstraint(v);
 418
 419       if(!(b->getValue()).noninfinitesimalIsZero()){ continue; }
 420
 421       Polynomial poly = Polynomial::parsePolynomial(d_vars.asNode(v));
 422       if(poly.size() != 2) { continue; }
 423
 424       Polynomial::iterator j = poly.begin();
 425       Monomial first = *j;
 426       ++j;
 427       Monomial second = *j;
 428
 429       bool firstIsPos = first.constantIsPositive();
 430       bool secondIsPos = second.constantIsPositive();
 431
 432       if(firstIsPos == secondIsPos){ continue; }
 433
 434       Monomial pos = firstIsPos == lbCap ? first : second;
 435       Monomial neg = firstIsPos != lbCap ? first : second;
 436       // pos >= neg
 437       VarList p = pos.getVarList();
 438       VarList n = neg.getVarList();
 439       if(d_vars.hasArithVar(p.getNode())){
 440         ArithVar ap = d_vars.asArithVar(p.getNode());
 441         if( d_colCandidates.isKey(ap)){
 442           BoundCounts bc = d_colCandidates.get(ap);
 443           bc = BoundCounts(bc.lowerBoundCount(), bc.upperBoundCount()+1);
 444           maxCount = std::max(maxCount, bc.upperBoundCount());
 445           d_colCandidates.set(ap, bc);
 446         }
 447       }
 448       if(d_vars.hasArithVar(n.getNode())){
 449         ArithVar an = d_vars.asArithVar(n.getNode());
 450         if( d_colCandidates.isKey(an)){
 451           BoundCounts bc = d_colCandidates.get(an);
 452           bc = BoundCounts(bc.lowerBoundCount()+1, bc.upperBoundCount());
 453           maxCount = std::max(maxCount, bc.lowerBoundCount());
 454           d_colCandidates.set(an, bc);
 455         }
 456       }
 457     }
 458   }
 459
 460   // Attempt row
 461   double avgRowLength = d_rowCandidates.size() >= 1 ?
 462     ( sumRowLength / d_rowCandidates.size() ) : 0.0;
 463
 464   // There is a large row among the candidates
 465   bool guessARowCandidate = maxRowLength >= (10.0 * avgRowLength);
 466
 467   double rowLengthReq = (maxRowLength * .9);
 468
 469   if(guessARowCandidate){
 470     for(DenseMap<uint32_t>::const_iterator i = d_rowCandidates.begin(), iend =d_rowCandidates.end(); i != iend; ++i ){
 471       ArithVar r = *i;
 472       uint32_t len = d_rowCandidates[r];
 473
 474       int dir = guessDir(r);
 475       if(len >= rowLengthReq){
 476         if(s_verbosity >= 1){
 477           Message() << "high row " << r << " " << len << " " << avgRowLength << " " << dir<< endl;
 478           d_vars.printModel(r, Message());
 479         }
 480         ret.push_back(ArithRatPair(r, Rational(dir)));
 481       }
 482     }
 483   }
 484
 485   // Attempt columns
 486   bool guessAColCandidate = maxCount >= 4;
 487   if(guessAColCandidate){
 488     for(DenseMap<BoundCounts>::const_iterator i = d_colCandidates.begin(), iend = d_colCandidates.end(); i != iend; ++i ){
 489       ArithVar c = *i;
 490       BoundCounts bc = d_colCandidates[c];
 491
 492       int dir = guessDir(c);
 493       double ubScore = double(bc.upperBoundCount()) / maxCount;
 494       double lbScore = double(bc.lowerBoundCount()) / maxCount;
 495       if(ubScore  >= .9 || lbScore >= .9){
 496         if(s_verbosity >= 1){
 497           Message() << "high col " << c << " " << bc << " " << ubScore << " " << lbScore << " " << dir << endl;
 498           d_vars.printModel(c, Message());
 499         }
 500         ret.push_back(ArithRatPair(c, Rational(c)));
 501       }
 502     }
 503   }
 504
 505
 506   return ret;
 507 }
 508
 509 void ApproxGLPK::setOptCoeffs(const ArithRatPairVec& ref){
 510   DenseMap<double> nbCoeffs;
 511
 512   for(ArithRatPairVec::const_iterator i = ref.begin(), iend = ref.end(); i != iend; ++i){
 513     ArithVar v = (*i).first;
 514     const Rational& q = (*i).second;
 515
 516     if(d_vars.isSlack(v)){
 517       // replace the variable by its definition and multiply by q
 518       Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v));
 519       Polynomial pq = p * q;
 520
 521       for(Polynomial::iterator j = pq.begin(), jend = pq.end(); j != jend; ++j){
 522         const Monomial& mono = *j;
 523         const Constant& constant = mono.getConstant();
 524         const VarList& variable = mono.getVarList();
 525
 526         Node n = variable.getNode();
 527
 528         Assert(d_vars.hasArithVar(n));
 529         ArithVar av = d_vars.asArithVar(n);
 530         int colIndex = d_colIndices[av];
 531         double coeff = constant.getValue().getDouble();
 532         if(!nbCoeffs.isKey(colIndex)){
 533           nbCoeffs.set(colIndex, 0.0);
 534         }
 535         nbCoeffs.set(colIndex, nbCoeffs[colIndex]+coeff);
 536       }
 537     }else{
 538       int colIndex = d_colIndices[v];
 539       double coeff = q.getDouble();
 540       if(!nbCoeffs.isKey(colIndex)){
 541         nbCoeffs.set(colIndex, 0.0);
 542       }
 543       nbCoeffs.set(colIndex, nbCoeffs[colIndex]+coeff);
 544     }
 545   }
 546   for(DenseMap<double>::const_iterator ci =nbCoeffs.begin(), ciend = nbCoeffs.end(); ci != ciend; ++ci){
 547     Index colIndex = *ci;
 548     double coeff = nbCoeffs[colIndex];
 549     glp_set_obj_coef(d_prob, colIndex, coeff);
 550   }
 551 }
 552
 553
 554 /*
 555  * rough strategy:
 556  *  real relaxation
 557  *   try approximate real optimization of error function
 558  *   pivot in its basis
 559  *   update to its assignment
 560  *   check with FCSimplex
 561  *  check integer solution
 562  *   try approximate mixed integer problem
 563  *   stop at the first feasible point
 564  *   pivot in its basis
 565  *   update to its assignment
 566  *   check with FCSimplex
 567  */
 568
 569 void ApproxGLPK::printGLPKStatus(int status, std::ostream& out){
 570   switch(status){
 571   case GLP_OPT:
 572     out << "GLP_OPT" << endl;
 573     break;
 574   case GLP_FEAS:
 575     out << "GLP_FEAS" << endl;
 576     break;
 577   case GLP_INFEAS:
 578     out << "GLP_INFEAS" << endl;
 579     break;
 580   case GLP_NOFEAS:
 581     out << "GLP_NOFEAS" << endl;
 582     break;
 583   case GLP_UNBND:
 584     out << "GLP_UNBND" << endl;
 585     break;
 586   case GLP_UNDEF:
 587     out << "GLP_UNDEF" << endl;
 588     break;
 589   default:
 590     out << "Status unknown" << endl;
 591     break;
 592   }
 593 }
 594
 595 ApproxGLPK::~ApproxGLPK(){
 596   glp_delete_prob(d_prob);
 597 }
 598
 599
 600 ApproximateSimplex::Solution ApproxGLPK::extractSolution(bool mip) const{
 601   Assert(d_solvedRelaxation);
 602   Assert(!mip  || d_solvedMIP);
 603
 604   ApproximateSimplex::Solution sol;
 605   DenseSet& newBasis = sol.newBasis;
 606   DenseMap<DeltaRational>& newValues = sol.newValues;
 607
 608   for(ArithVariables::var_iterator i = d_vars.var_begin(), i_end = d_vars.var_end(); i != i_end; ++i){
 609     ArithVar vi = *i;
 610     bool isSlack = d_vars.isSlack(vi);
 611     int glpk_index = isSlack ? d_rowIndices[vi] : d_colIndices[vi];
 612
 613     int status = isSlack ? glp_get_row_stat(d_prob, glpk_index) : glp_get_col_stat(d_prob, glpk_index);
 614     if(s_verbosity >= 2){
 615       Message() << "assignment " << vi << endl;
 616     }
 617
 618     bool useDefaultAssignment = false;
 619
 620     switch(status){
 621     case GLP_BS:
 622       //Message() << "basic" << endl;
 623       newBasis.add(vi);
 624       useDefaultAssignment = true;
 625       break;
 626     case GLP_NL:
 627     case GLP_NS:
 628       if(!mip){
 629         if(s_verbosity >= 2){ Message() << "non-basic lb" << endl; }
 630         newValues.set(vi, d_vars.getLowerBound(vi));
 631       }else{// intentionally fall through otherwise
 632         useDefaultAssignment = true;
 633       }
 634       break;
 635     case GLP_NU:
 636       if(!mip){
 637         if(s_verbosity >= 2){ Message() << "non-basic ub" << endl; }
 638         newValues.set(vi, d_vars.getUpperBound(vi));
 639       }else {// intentionally fall through otherwise
 640         useDefaultAssignment = true;
 641       }
 642       break;
 643     default:
 644       {
 645         useDefaultAssignment = true;
 646       }
 647       break;
 648     }
 649
 650     if(useDefaultAssignment){
 651       if(s_verbosity >= 2){ Message() << "non-basic other" << endl; }
 652
 653       double newAssign =
 654         mip ?
 655         (isSlack ? glp_mip_row_val(d_prob, glpk_index) :  glp_mip_col_val(d_prob, glpk_index))
 656         : (isSlack ? glp_get_row_prim(d_prob, glpk_index) :  glp_get_col_prim(d_prob, glpk_index));
 657       const DeltaRational& oldAssign = d_vars.getAssignment(vi);
 658
 659
 660       if(d_vars.hasLowerBound(vi) &&
 661          roughlyEqual(newAssign, d_vars.getLowerBound(vi).approx(SMALL_FIXED_DELTA))){
 662         //Message() << "  to lb" << endl;
 663
 664         newValues.set(vi, d_vars.getLowerBound(vi));
 665       }else if(d_vars.hasUpperBound(vi) &&
 666                roughlyEqual(newAssign, d_vars.getUpperBound(vi).approx(SMALL_FIXED_DELTA))){
 667         newValues.set(vi, d_vars.getUpperBound(vi));
 668         // Message() << "  to ub" << endl;
 669       }else{
 670
 671         double rounded = round(newAssign);
 672         if(roughlyEqual(newAssign, rounded)){
 673           // Message() << "roughly equal " << rounded << " " << newAssign << " " << oldAssign << endl;
 674           newAssign = rounded;
 675         }else{
 676           // Message() << "not roughly equal " << rounded << " " << newAssign << " " << oldAssign << endl;
 677         }
 678
 679         DeltaRational proposal = estimateWithCFE(newAssign);
 680
 681
 682         if(roughlyEqual(newAssign, oldAssign.approx(SMALL_FIXED_DELTA))){
 683           // Message() << "  to prev value" << newAssign << " " << oldAssign << endl;
 684           proposal = d_vars.getAssignment(vi);
 685         }
 686
 687
 688         if(d_vars.strictlyLessThanLowerBound(vi, proposal)){
 689           //Message() << "  round to lb " << d_vars.getLowerBound(vi) << endl;
 690           proposal = d_vars.getLowerBound(vi);
 691         }else if(d_vars.strictlyGreaterThanUpperBound(vi, proposal)){
 692           //Message() << "  round to ub " << d_vars.getUpperBound(vi) << endl;
 693           proposal = d_vars.getUpperBound(vi);
 694         }else{
 695           //Message() << "  use proposal" << proposal << " " << oldAssign  << endl;
 696         }
 697         newValues.set(vi, proposal);
 698       }
 699     }
 700   }
 701   return sol;
 702 }
 703
 704 ApproximateSimplex::ApproxResult ApproxGLPK::solveRelaxation(){
 705   Assert(!d_solvedRelaxation);
 706
 707   glp_smcp parm;
 708   glp_init_smcp(&parm);
 709   parm.presolve = GLP_OFF;
 710   parm.meth = GLP_PRIMAL;
 711   parm.pricing = GLP_PT_PSE;
 712   parm.it_lim = d_pivotLimit;
 713   parm.msg_lev = GLP_MSG_OFF;
 714   if(s_verbosity >= 1){
 715     parm.msg_lev = GLP_MSG_ALL;
 716   }
 717
 718   int res = glp_simplex(d_prob, &parm);
 719   switch(res){
 720   case 0:
 721     {
 722       int status = glp_get_status(d_prob);
 723       switch(status){
 724       case GLP_OPT:
 725       case GLP_FEAS:
 726       case GLP_UNBND:
 727         d_solvedRelaxation = true;
 728         return ApproxSat;
 729       case GLP_INFEAS:
 730       case GLP_NOFEAS:
 731         d_solvedRelaxation = true;
 732         return ApproxUnsat;
 733       default:
 734         return ApproxError;
 735       }
 736     }
 737   default:
 738     return ApproxError;
 739   }
 740 }
 741
 742 void stopAtBingoOrPivotLimit(glp_tree *tree, void *info){
 743   int pivotLimit = *((int*)info);
 744   switch(glp_ios_reason(tree)){
 745   case GLP_IBINGO:
 746     glp_ios_terminate(tree);
 747     break;
 748   default:
 749     glp_prob* prob = glp_ios_get_prob(tree);
 750     int iterationcount = lpx_get_int_parm(prob, LPX_K_ITCNT);
 751     if(iterationcount > pivotLimit){
 752       glp_ios_terminate(tree);
 753     }
 754     break;
 755   }
 756 }
 757
 758 ApproximateSimplex::ApproxResult ApproxGLPK::solveMIP(){
 759   Assert(d_solvedRelaxation);
 760   // Explicitly disable presolving
 761   // We need the basis thus the presolver must be off!
 762   // This is default, but this is just being cautious.
 763   glp_iocp parm;
 764   glp_init_iocp(&parm);
 765   parm.presolve = GLP_OFF;
 766   parm.pp_tech = GLP_PP_NONE;
 767   parm.fp_heur = GLP_ON;
 768   parm.gmi_cuts = GLP_ON;
 769   parm.mir_cuts = GLP_ON;
 770   parm.cov_cuts = GLP_ON;
 771   parm.cb_func = stopAtBingoOrPivotLimit;
 772   parm.cb_info = &d_pivotLimit;
 773   parm.msg_lev = GLP_MSG_OFF;
 774   if(s_verbosity >= 1){
 775     parm.msg_lev = GLP_MSG_ALL;
 776   }
 777   int res = glp_intopt(d_prob, &parm);
 778
 779   switch(res){
 780   case 0:
 781   case GLP_ESTOP:
 782     {
 783       int status = glp_mip_status(d_prob);
 784       switch(status){
 785       case GLP_OPT:
 786       case GLP_FEAS:
 787         d_solvedMIP = true;
 788         return ApproxSat;
 789       case GLP_NOFEAS:
 790         d_solvedMIP = true;
 791         return ApproxUnsat;
 792       default:
 793         return ApproxError;
 794       }
 795     }
 796   default:
 797     return ApproxError;
 798   }
 799 }
 800
 801 }/* CVC4::theory::arith namespace */
 802 }/* CVC4::theory namespace */
 803 }/* CVC4 namespace */
 804 #endif /*#ifdef CVC4_USE_GLPK */
 805 /* End GPLK implementation. */