From: Tim King Date: Mon, 6 May 2013 22:38:12 +0000 (-0400) Subject: Adding a heuristic for guessing an optimization function when using glpk. X-Git-Tag: cvc5-1.0.0~7287^2~152 X-Git-Url: https://git.libre-soc.org/?a=commitdiff_plain;h=0ff427f9cb3384f7d85c8f401a14b57b52c87cdd;p=cvc5.git Adding a heuristic for guessing an optimization function when using glpk. --- diff --git a/src/theory/arith/approx_simplex.cpp b/src/theory/arith/approx_simplex.cpp index ef3ea0484..181def552 100644 --- a/src/theory/arith/approx_simplex.cpp +++ b/src/theory/arith/approx_simplex.cpp @@ -2,6 +2,7 @@ #include "theory/arith/approx_simplex.h" #include "theory/arith/normal_form.h" +#include "theory/arith/constraint.h" #include #include @@ -94,6 +95,10 @@ public: return Solution(); } + virtual ArithRatPairVec heuristicOptCoeffs() const{ + return ArithRatPairVec(); + } + virtual ApproxResult solveMIP(){ return ApproxError; } @@ -143,6 +148,8 @@ public: return extractSolution(false); } + virtual ArithRatPairVec heuristicOptCoeffs() const; + virtual ApproxResult solveMIP(); virtual Solution extractMIP() const{ return extractSolution(true); @@ -152,10 +159,11 @@ public: static void printGLPKStatus(int status, std::ostream& out); private: Solution extractSolution(bool mip) const; + int guessDir(ArithVar v) const; }; #warning "set back to 0" -int ApproxGLPK::s_verbosity = 0; +int ApproxGLPK::s_verbosity = 1; }/* CVC4::theory::arith namespace */ }/* CVC4::theory namespace */ @@ -309,6 +317,189 @@ ApproxGLPK::ApproxGLPK(const ArithVariables& avars) : delete[] ja; delete[] ar; } +int ApproxGLPK::guessDir(ArithVar v) const{ + if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){ + return -1; + }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){ + return 1; + }else if(!d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){ + return 0; + }else{ + int ubSgn = d_vars.getUpperBound(v).sgn(); + int lbSgn = d_vars.getLowerBound(v).sgn(); + + if(ubSgn != 0 && lbSgn == 0){ + return -1; + }else if(ubSgn == 0 && lbSgn != 0){ + return 1; + } + + return 1; + } +} + +ArithRatPairVec ApproxGLPK::heuristicOptCoeffs() const{ + ArithRatPairVec ret; + + // Strategies are guess: + // 1 simple shared "ceiling" variable: danoint, pk1 + // x1 >= c, x1 >= tmp1, x1 >= tmp2, ... + // 1 large row: fixnet, vpm2, pp08a + // (+ .......... ) <= c + // Not yet supported: + // 1 complex shared "ceiling" variable: opt1217 + // x1 >= c, x1 >= (+ ..... ), x1 >= (+ ..... ) + // and all of the ... are the same sign + + + // Candidates: + // 1) Upper and lower bounds are not equal. + // 2) The variable is not integer + // 3a) For columns look for a ceiling variable + // 3B) For rows look for a large row with + + DenseMap d_colCandidates; + DenseMap d_rowCandidates; + + double sumRowLength = 0.0; + uint32_t maxRowLength = 0; + for(ArithVariables::var_iterator vi = d_vars.var_begin(), vi_end = d_vars.var_end(); vi != vi_end; ++vi){ + ArithVar v = *vi; + + if(s_verbosity >= 2){ + Message() << v << " "; + d_vars.printModel(v, Message()); + } + + int type; + if(d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){ + if(d_vars.boundsAreEqual(v)){ + type = GLP_FX; + }else{ + type = GLP_DB; + } + }else if(d_vars.hasUpperBound(v) && !d_vars.hasLowerBound(v)){ + type = GLP_UP; + }else if(!d_vars.hasUpperBound(v) && d_vars.hasLowerBound(v)){ + type = GLP_LO; + }else{ + type = GLP_FR; + } + + if(type != GLP_FX && type != GLP_FR){ + + if(d_vars.isSlack(v)){ + Polynomial p = Polynomial::parsePolynomial(d_vars.asNode(v)); + uint32_t len = p.size(); + d_rowCandidates.set(v, len); + sumRowLength += len; + maxRowLength =std::max(maxRowLength, len); + }else if(!d_vars.isInteger(v)){ + d_colCandidates.set(v, BoundCounts()); + } + } + } + + uint32_t maxCount = 0; + for(DenseMap::const_iterator i = d_rowIndices.begin(), i_end = d_rowIndices.end(); i != i_end; ++i){ + ArithVar v = *i; + + bool lbCap = d_vars.hasLowerBound(v) && !d_vars.hasUpperBound(v); + bool ubCap = !d_vars.hasLowerBound(v) && d_vars.hasUpperBound(v); + + if(lbCap || ubCap){ + Constraint b = lbCap ? d_vars.getLowerBoundConstraint(v) + : d_vars.getUpperBoundConstraint(v); + + if(!(b->getValue()).noninfinitesimalIsZero()){ continue; } + + Polynomial poly = Polynomial::parsePolynomial(d_vars.asNode(v)); + if(poly.size() != 2) { continue; } + + Polynomial::iterator j = poly.begin(); + Monomial first = *j; + ++j; + Monomial second = *j; + + bool firstIsPos = first.constantIsPositive(); + bool secondIsPos = second.constantIsPositive(); + + if(firstIsPos == secondIsPos){ continue; } + + Monomial pos = firstIsPos == lbCap ? first : second; + Monomial neg = firstIsPos != lbCap ? first : second; + // pos >= neg + VarList p = pos.getVarList(); + VarList n = neg.getVarList(); + if(d_vars.hasArithVar(p.getNode())){ + ArithVar ap = d_vars.asArithVar(p.getNode()); + if( d_colCandidates.isKey(ap)){ + BoundCounts bc = d_colCandidates.get(ap); + bc = BoundCounts(bc.lowerBoundCount(), bc.upperBoundCount()+1); + maxCount = std::max(maxCount, bc.upperBoundCount()); + d_colCandidates.set(ap, bc); + } + } + if(d_vars.hasArithVar(n.getNode())){ + ArithVar an = d_vars.asArithVar(n.getNode()); + if( d_colCandidates.isKey(an)){ + BoundCounts bc = d_colCandidates.get(an); + bc = BoundCounts(bc.lowerBoundCount()+1, bc.upperBoundCount()); + maxCount = std::max(maxCount, bc.lowerBoundCount()); + d_colCandidates.set(an, bc); + } + } + } + } + + // Attempt row + double avgRowLength = d_rowCandidates.size() >= 1 ? + ( sumRowLength / d_rowCandidates.size() ) : 0.0; + + // There is a large row among the candidates + bool guessARowCandidate = maxRowLength >= (10.0 * avgRowLength); + + double rowLengthReq = (maxRowLength * .9); + + if(guessARowCandidate){ + for(DenseMap::const_iterator i = d_rowCandidates.begin(), iend =d_rowCandidates.end(); i != iend; ++i ){ + ArithVar r = *i; + uint32_t len = d_rowCandidates[r]; + + int dir = guessDir(r); + if(len >= rowLengthReq){ + if(s_verbosity >= 1){ + Message() << "high row " << r << " " << len << " " << avgRowLength << " " << dir<< endl; + d_vars.printModel(r, Message()); + } + ret.push_back(ArithRatPair(r, Rational(dir))); + } + } + } + + // Attempt columns + bool guessAColCandidate = maxCount >= 4; + if(guessAColCandidate){ + for(DenseMap::const_iterator i = d_colCandidates.begin(), iend = d_colCandidates.end(); i != iend; ++i ){ + ArithVar c = *i; + BoundCounts bc = d_colCandidates[c]; + + int dir = guessDir(c); + double ubScore = double(bc.upperBoundCount()) / maxCount; + double lbScore = double(bc.lowerBoundCount()) / maxCount; + if(ubScore >= .9 || lbScore >= .9){ + if(s_verbosity >= 1){ + Message() << "high col " << c << " " << bc << " " << ubScore << " " << lbScore << " " << dir << endl; + d_vars.printModel(c, Message()); + } + ret.push_back(ArithRatPair(c, Rational(c))); + } + } + } + + + return ret; +} void ApproxGLPK::setOptCoeffs(const ArithRatPairVec& ref){ DenseMap nbCoeffs; @@ -354,6 +545,7 @@ void ApproxGLPK::setOptCoeffs(const ArithRatPairVec& ref){ } } + /* * rough strategy: * real relaxation @@ -422,7 +614,7 @@ ApproximateSimplex::Solution ApproxGLPK::extractSolution(bool mip) const{ switch(status){ case GLP_BS: - //cout << "basic" << endl; + //Message() << "basic" << endl; newBasis.add(vi); useDefaultAssignment = true; break; @@ -462,40 +654,40 @@ ApproximateSimplex::Solution ApproxGLPK::extractSolution(bool mip) const{ if(d_vars.hasLowerBound(vi) && roughlyEqual(newAssign, d_vars.getLowerBound(vi).approx(SMALL_FIXED_DELTA))){ - //cout << " to lb" << endl; + //Message() << " to lb" << endl; newValues.set(vi, d_vars.getLowerBound(vi)); }else if(d_vars.hasUpperBound(vi) && roughlyEqual(newAssign, d_vars.getUpperBound(vi).approx(SMALL_FIXED_DELTA))){ newValues.set(vi, d_vars.getUpperBound(vi)); - // cout << " to ub" << endl; + // Message() << " to ub" << endl; }else{ double rounded = round(newAssign); if(roughlyEqual(newAssign, rounded)){ - // cout << "roughly equal " << rounded << " " << newAssign << " " << oldAssign << endl; + // Message() << "roughly equal " << rounded << " " << newAssign << " " << oldAssign << endl; newAssign = rounded; }else{ - // cout << "not roughly equal " << rounded << " " << newAssign << " " << oldAssign << endl; + // Message() << "not roughly equal " << rounded << " " << newAssign << " " << oldAssign << endl; } DeltaRational proposal = estimateWithCFE(newAssign); if(roughlyEqual(newAssign, oldAssign.approx(SMALL_FIXED_DELTA))){ - // cout << " to prev value" << newAssign << " " << oldAssign << endl; + // Message() << " to prev value" << newAssign << " " << oldAssign << endl; proposal = d_vars.getAssignment(vi); } if(d_vars.strictlyLessThanLowerBound(vi, proposal)){ - //cout << " round to lb " << d_vars.getLowerBound(vi) << endl; + //Message() << " round to lb " << d_vars.getLowerBound(vi) << endl; proposal = d_vars.getLowerBound(vi); }else if(d_vars.strictlyGreaterThanUpperBound(vi, proposal)){ - //cout << " round to ub " << d_vars.getUpperBound(vi) << endl; + //Message() << " round to ub " << d_vars.getUpperBound(vi) << endl; proposal = d_vars.getUpperBound(vi); }else{ - //cout << " use proposal" << proposal << " " << oldAssign << endl; + //Message() << " use proposal" << proposal << " " << oldAssign << endl; } newValues.set(vi, proposal); } diff --git a/src/theory/arith/approx_simplex.h b/src/theory/arith/approx_simplex.h index d4680939d..175e56f8e 100644 --- a/src/theory/arith/approx_simplex.h +++ b/src/theory/arith/approx_simplex.h @@ -44,6 +44,8 @@ public: /** Sets a maximization criteria for the approximate solver.*/ virtual void setOptCoeffs(const ArithRatPairVec& ref) = 0; + virtual ArithRatPairVec heuristicOptCoeffs() const = 0; + virtual ApproxResult solveRelaxation() = 0; virtual Solution extractRelaxation() const = 0; diff --git a/src/theory/arith/theory_arith_private.cpp b/src/theory/arith/theory_arith_private.cpp index 08907880e..a49d3776d 100644 --- a/src/theory/arith/theory_arith_private.cpp +++ b/src/theory/arith/theory_arith_private.cpp @@ -115,6 +115,8 @@ TheoryArithPrivate::TheoryArithPrivate(TheoryArith& containing, context::Context d_cutCount(c, 0), d_cutInContext(c), d_likelyIntegerInfeasible(c, false), + d_guessedCoeffSet(c, false), + d_guessedCoeffs(), d_statistics() { srand(79); @@ -1584,6 +1586,14 @@ bool TheoryArithPrivate::solveRealRelaxation(Theory::Effort effortLevel){ ApproximateSimplex* approxSolver = ApproximateSimplex::mkApproximateSimplexSolver(d_partialModel); approxSolver->setPivotLimit(relaxationLimit); + if(!d_guessedCoeffSet){ + d_guessedCoeffs = approxSolver->heuristicOptCoeffs(); + d_guessedCoeffSet = true; + } + if(!d_guessedCoeffs.empty()){ + approxSolver->setOptCoeffs(d_guessedCoeffs); + } + ApproximateSimplex::ApproxResult relaxRes, mipRes; ApproximateSimplex::Solution relaxSolution, mipSolution; relaxRes = approxSolver->solveRelaxation(); diff --git a/src/theory/arith/theory_arith_private.h b/src/theory/arith/theory_arith_private.h index 3da064a80..c8f49ab01 100644 --- a/src/theory/arith/theory_arith_private.h +++ b/src/theory/arith/theory_arith_private.h @@ -572,6 +572,10 @@ private: context::CDO d_likelyIntegerInfeasible; + + context::CDO d_guessedCoeffSet; + ArithRatPairVec d_guessedCoeffs; + /** These fields are designed to be accessible to TheoryArith methods. */ class Statistics { public: